1 Introduction↩︎

Reinforcement learning (RL) has emerged as a powerful framework for sequential decision making under uncertainty, bringing together ideas from optimal stochastic control, Markov decision processes (MDPs), and statistical decision theory. In classical MDP theory, the model of the environment is assumed to be fully known and the objective is to optimize a control policy accordingly. In contrast, in reinforcement learning the dynamics of the environment are initially unknown, and the agent must learn about them—either directly or implicitly—through interaction, thereby introducing elements of statistical estimation into the problem. From a statistical perspective, this learning requirement fundamentally changes the performance criteria. Rather than the likelihood-based risk functions common in classical statistical inference, reinforcement learning typically focuses on regret, which quantifies the performance loss incurred by using a learning policy instead of an optimal one from the outset. This regret-based viewpoint, closely connected to online learning, motivates specialized techniques such as optimistic exploration, in which additional exploration is deliberately introduced to gather information that ultimately improves long-term performance.

Within this framework, two broad methodological paradigms have emerged: model-based and model-free approaches. In model-based reinforcement learning, the agent explicitly estimates transition dynamics and rewards and uses these estimates to plan or optimize policies. Model-free methods, by contrast, avoid explicit model estimation and instead learn value functions or policies directly through repeated interaction with the environment. A canonical example is Q-learning, where action-value estimates are updated recursively using observed rewards and successor values without constructing an explicit model of state transitions.

Optimization methods also lie at the heart of reinforcement learning, shaping how agents refine policies and value functions to maximize cumulative rewards. Many algorithms can be interpreted as gradient-based, mirror-descent, or proximal update schemes, thereby connecting reinforcement learning with modern convex and nonconvex optimization theory. When additional regularization terms—such as entropy or Kullback–Leibler (KL) penalties—are incorporated into the objective, the resulting optimization problems often acquire improved stability and favorable geometric properties that facilitate convergence analysis. For example, entropy-regularized policy gradient methods encourage exploration while smoothing the optimization landscape and reducing the risk of convergence to suboptimal deterministic policies. These optimization perspectives provide a principled framework for balancing exploration and exploitation while supporting both theoretical guarantees and practical algorithmic performance.

The mathematical structure of reinforcement learning is further shaped by several foundational analytical tools. Central among these is the Bellman principle, which expresses optimal value functions through recursive fixed-point relations and forms the basis of dynamic programming methods. The performance difference lemma provides a quantitative characterization of policy improvement and plays an important role in the analysis of policy gradient methods. The log-derivative trick (also known as the score-function gradient estimator) enables efficient computation of policy gradients in stochastic environments. In addition, concentration inequalities for dependent data and mixing processes provide the probabilistic control required for finite-sample guarantees in policy evaluation and exploration strategies. Together, these tools form a unified theoretical backbone for analyzing modern reinforcement learning algorithms.

Reinforcement learning problems can be studied under several different time-horizon formulations, each leading to distinct analytical challenges. In finite-horizon episodic settings, interaction with the environment is organized into trajectories of fixed length and performance is typically evaluated over a finite number of steps or episodes. In contrast, infinite-horizon formulations consider ongoing interaction and evaluate performance using discounted or average reward criteria. In the latter setting, analyses frequently rely on ergodic properties of Markov processes. Under suitable assumptions such as irreducibility and aperiodicity, ergodicity ensures convergence to a stationary distribution, which in turn supports theoretical guarantees for long-run performance measures such as average reward or regret. Techniques from Markov chain theory—including drift and minorization conditions and mixing-time estimates—therefore play a central role in establishing convergence and regret bounds in both model-based and model-free reinforcement learning.

Independently of the time-horizon formulation, reinforcement learning problems may be studied in either tabular or function-approximation regimes. When the state–action space is finite and moderate in size, policies and value functions can be represented exactly. In large or continuous spaces, however, function approximation becomes essential. In this setting, tools from approximation theory provide insight into the expressive power and stability of parametric representations such as linear architectures, kernel methods, and deep neural networks. Results such as universal approximation properties and quantitative approximation error bounds help characterize the trade-off between representational flexibility and statistical efficiency, thereby guiding the design and analysis of scalable reinforcement learning algorithms.

This review surveys the mathematical foundations of reinforcement learning with an emphasis on operator-theoretic, probabilistic, and optimization-based perspectives. We highlight how regret-based objectives shape exploration strategies and unify model-based and model-free approaches within a common analytical framework. By tracing the interplay between dynamic programming, stochastic approximation, concentration methods, and function approximation theory, the survey provides a coherent mathematical entry point to the analysis of modern reinforcement learning algorithms.

2 Main setup↩︎

Let \(\mathbb{R}\) denote the set of real numbers, and for any integer \(N\), denote \([N] = \{0,\dots,N-1\}\). We denote by \(\Delta(\mathcal{X})\) the set of probability measures on the measurable space \((\mathcal{X},\mathcal{F}),\) where \(\mathcal{F}\subset 2^{\mathcal{X}}\) denotes a \(\sigma\)-algebra on \(\mathcal{X}\). The central object of the review is the Markov decision process defined by the tuple \(\mathcal{M}=(\mathcal{S},\mathcal{A},P, r)\), where \(\mathcal{S}\) and \(\mathcal{A}\) denote the sets of states and actions, \(P: \mathcal{S}\times\mathcal{A}\rightarrow \Delta(\mathcal{S})\) is the Markov transition kernel, and \(r\) is the reward function. The main part of the review (Chapters 3–6) considers the case where the state and action spaces are finite (tabular setting), with cardinalities \(S\) and \(A\), respectively. Chapter 7 presents results for the case of an infinite state space \(\mathcal{S}\). In the tabular setting, we will often interpret \(P\) as a matrix of dimension \(SA \times S\) and \(R\) as a vector in \(\mathbb{R}^{SA}\). A Markov decision process describes an environment with which an agent interacts sequentially at discrete points in time. At time \(t\), the agent observes a state \(s_t\in \mathcal{S},\) takes an action \(a_t\), receives a reward \(r=r(s,a),\) and transitions to the next state \(s_{t+1} \sim P(\cdot|s_t,a_t)\). The goal of the agent is to maximize the cumulative reward over time. To formalize this objective, we denote by \(\pi : \mathcal{S}\xrightarrow[]{}\Delta(\mathcal{A})\) the agent’s policy that specifies how the agent selects actions. Note that under this definition, the choice of action depends only on the current state. Such policies are called stationary. We denote by \(\Pi\) the set of all stationary policies. In a more general setting, a policy may depend on the entire history of the agent’s observations. However, enlarging the policy class in this way does not improve performance, as there always exists an optimal stationary policy; see, e.g., [1] (Theorem 6.2.10). Several approaches can be used to compare policies; the main ones are outlined below.

2.0.0.1 Finite-horizon setting

In the finite-horizon setting, the interaction between the agent and the environment lasts for \(H > 0\) steps. The measure of policy performance is the expected cumulative reward, given by \[\begin{align} V^{\pi}(s) &= \mathbb{E}_\pi \Bigl[\sum_{h=0}^{H-1}r(s_h,a_h) \, |\, s_0 = s \Bigr], \\ Q^{\pi}(s, a) &= \mathbb{E}_\pi \Bigl[\sum_{h=0}^{H-1}r(s_h,a_h) \, |\, s_0 = s, a_0=a \Bigr]\;. \end{align}\] where the expectation is taken with respect to the policy \(\pi\) and the stochasticity of the transition kernel \(P\). The function \(V^\pi : \mathcal{S}\rightarrow\mathbb{R}\) is referred to as the state-value function, whereas \(Q^\pi : \mathcal{S}\times\mathcal{A}\rightarrow\mathbb{R}\) is known as the state–action value function. Chapter 6 provides a detailed analysis of this formulation.

2.0.0.2 Discounted setting

In many cases, it is important not to restrict attention to a finite planning horizon, but rather to evaluate the performance of a policy over an infinite trajectory. It is often natural to assume that rewards obtained in the distant future have less significance than immediate ones. This leads to the formulation of a discounted Markov Decision process (DMDP), characterized by a discount factor \(\gamma \in (0,1)\). In this setting, the agent aims to maximize the expected discounted cumulative reward \[\label{eq:value95function} V^\pi (s) := \mathbb{E}_\pi \Bigl[ \sum_{t=0}^\infty \gamma^t r(s_t, a_t)\, |\, s_0 = s \Bigr],\tag{1}\] \[\label{eq:action95value} Q^\pi (s, a) := \mathbb{E}_\pi \Bigl[ \sum_{t=0}^\infty \gamma^t r(s_t, a_t)\, |\, s_0 = s, a_0=a \Bigr].\tag{2}\] The main results within this framework are discussed in Chapters 3 - 5.

2.0.0.3 Average reward setting

Alongside the discounted setting, there is the infinite-horizon average-reward setting. This setting is majorly comprised of MDP where the expected return is defined via state value function and state-action value function, averaged along the trajectory rather than computed as a geometrically discounted sum: \[\label{eq:amdp95return} \begin{align} V^{\pi} (s) :=& \lim\limits_{H\rightarrow\infty}\mathbb{E}_{\pi}\Bigl[\frac{1}{H}\sum\limits_{t = 0}^{H - 1}r(s_t, a_t)\mid s_0 = s\Bigr];\\ Q^{\pi} (s, a) :=& \lim\limits_{H\rightarrow\infty}\mathbb{E}_{\pi}\Bigl[\frac{1}{H}\sum\limits_{t = 0}^{H - 1}r(s_t, a_t)\mid s_{0} = s, a_{0} = a\Bigr]. \end{align}\tag{3}\] To guarantee convergence and the finiteness of the expectation in 3 , we assume that the Markov kernel \(P^{\pi}\) (see 5 ) admits a unique invariant distribution \(\mu_\pi\) for all \(\pi\) and is uniformly in \(\pi\) geometrically ergodic. That is, there exists \(t_{\operatorname{mix}}>0\) such that for all \(t \in \mathbb{N}\), \[\bigl\| (P^{\pi})^t - \mathbf{1} \cdot \mu^\top_{\pi} \bigr\|_{\infty} \leq \left(\frac{1}{4}\right)^{\lfloor t / t_{\operatorname{mix}} \rfloor}\;.\] The analysis methods for AMDP are presented in subsections 3.3, 4.1 and 4.2.

Thus far, reinforcement learning has been presented as sequential decision making in which an agent with policy \(\pi\) interacts with an environment modeled as an MDP. In the setting of a partially observable Markov decision process (POMDP), the agent does not have direct access to the underlying states \(s_t\in\mathcal{S}\) but instead receives observations \(o_t\in\mathcal{O}\) conditioned on those states. Decision making must therefore rely on \(o_t\) rather than \(s_t\), producing trajectories with the following joint distribution:

\[\label{eq:pomdp95def} P(s_{0:H}, o_{0:H}, a_{0:H-1}) = P(s_0)\,P(o_0 | s_0) \prod_{t=0}^{H-1}\! \pi(a_t \mid o_{0:t})\, P(s_{t+1} \mid s_t, a_t)\, P(o_{t+1} \mid s_{t+1})\tag{4}\] where \(H\rightarrow\infty\) in the infinite-horizon case. Sometimes the deterministic map \(\hat{O}: \mathcal{S}\to\mathcal{O}\) is introduced to simplify the notation by binding the randomness to the Markov kernel: \(o_t = \hat{O}(s_t)\) [2]. There is also an alternative view on POMDP with explicit conditioning of observations on both states and actions: \(o_{t + 1}\sim P(\cdot|s_{t + 1}, a_t)\) [3]. Some works consider policies in POMDP directly conditioned both on previous observations and actions: \(a_t\sim\pi(\cdot|o_{0:t}, a_{0:t - 1})\), \(t=\overline{1, H}\) [4]. For fully observable MDP, the trajectory likelihood formula simplifies because \(o_t = s_t, t = \overline{0, H}\): \[\begin{align} \tau &:= (s_0, a_0, s_1, a_1, \dots, s_H);\\ P(\tau) &= P(s_0)\cdot\prod_{t = 0}^{H - 1}\pi(a_t|s_t)\cdot P(s_{t + 1}|s_t, a_t). \end{align}\] Here and after we consider exclusively fully observable MDPs unless explicitly stated otherwise. The objective is to choose \(\pi\) to maximize expected return: \[J(\pi) = \mathbb{E}_{P(\tau)}\left[V^{\pi}(s_0)\right].\] The optimization problem as stated is unconstrained. However, one may impose explicit constraints on value functions, leading to a constrained MDP (CMDP): \[\begin{align} &\max\limits_{\pi}J(\pi);\\ &\text{s.t.}\quad J_i(\pi)\geq d_i, \quad i\in[m] \end{align}\] The value functions \(J_i\) can have an arbitrary form. The typical example of these constraints for a DMDP – value functions with different rewards and discounting: \[J_i(\pi) := \mathbb{E}_{\pi}\left[\sum\limits_{t = 0}^{\infty}\lambda_i^t\cdot r_i(s_t, a_t)\right],\quad\lambda_i\in(0, 1].\]

3 Dynamic programming↩︎

In this section, we explore dynamic programming (DP) which is a group of methods aimed at solving an MDP in the case when the dynamics of the environment is known, i.e., transition probabilities \(P(s'|s,a)\) are given for every \(s,s'\in\mathcal{S},\, a\in\mathcal{A}\). DP algorithms are theoretically important and often serve as a foundation for more complex methods that don’t assume availability of a perfect model of the environment. Throughout this section, we assume that \(\mathcal{S}\) and \(\mathcal{A}\) are finite, and rewards are specified by a deterministic function \(r : \mathcal{S}\times\mathcal{A}\to [0,1]\).

3.0.0.1 Stationary policies and induced Markov chains

A stationary policy \(\pi\) naturally induces a Markov chain on the state–action space \(\mathcal{S}\times \mathcal{A}\) with transition kernel \[\label{eq:p94pi95def} P^\pi\bigl((s',a') \mid (s,a)\bigr) = P(s' \mid s,a)\,\pi(a' \mid s').\tag{5}\] Assume that this Markov chain is ergodic (i.e., irreducible and aperiodic). Then it admits a unique stationary distribution \(\mu^\pi\) over \(\mathcal{S}\times \mathcal{A}\), called the state–action visitation distribution, which satisfies \[\mu^\pi(s',a') = \sum_{(s,a)\in \mathcal{S}\times \mathcal{A}} \mu^\pi(s,a)\, P(s' \mid s,a)\,\pi(a' \mid s'),\] By integrating out the action component according to the policy \(\pi\), one obtains the induced transition kernel on the state space \(\mathcal{S}\), \[\label{eq:p95pi95def} P_\pi(s' \mid s) = \sum_{a \in \mathcal{A}} \pi(a \mid s)\, P(s' \mid s,a),\tag{6}\] Assuming that the resulting Markov chain on \(\mathcal{S}\) is ergodic, we refer to its stationary distribution \(\nu_\pi\) as the state visitation distribution. Superscripts are used to indicate quantities induced by the policy \(\pi\), such as the transition kernels \(P^\pi\) and \(P_\pi\), as well as the stationary distribution \(\mu^\pi\) of the state–action Markov chain. In contrast, subscripts are used for distributions over states, such as \(\nu_\pi\), to emphasize their role as marginal (state visitation) distributions derived from \(\mu^\pi\). The two stationary distributions are consistent in the sense that \(\mu^\pi\) extends \(\nu_\pi\) to the product space: \[\mu^\pi(s,a) = \nu_{\pi}(s)\,\pi(a \mid s), \qquad (s,a)\in \mathcal{S}\times\mathcal{A}.\] Thus, \(\nu_\pi\) is the marginal of \(\mu^\pi\) over \(\mathcal{S}\), while \(\mu^\pi\) couples the state visitation distribution with the policy \(\pi\) to form a stationary distribution on \(\mathcal{S}\times\mathcal{A}\).

3.1 Bellman equations↩︎

3.1.1 Consistency equations↩︎

In what follows, we consider stationary policies \(\pi:\mathcal{S}\rightarrow \Delta(\mathcal{A})\). It follows directly from the definitions 1 and 2 , together with the law of total probability, that for all \(s\in\mathcal{S}\) and \(a\in\mathcal{A}\), the value functions \(V^\pi\) and \(Q^\pi\) satisfy the following relations, known as the Bellman consistency equations: \[\begin{align} \tag{7} V^\pi(s) &= \sum_{a\in\mathcal{A}}\pi(a|s)Q^\pi(s, a),\\ \tag{8} Q^\pi(s, a) &= r(s, a) + \gamma \sum_{s'\in\mathcal{S}} P(s'|s,a) V^\pi(s'). \end{align}\] Equation 7 immediately implies a useful identity, commonly employed in reinforcement learning analysis : \(P V^\pi = P^\pi Q^\pi\). Indeed, for any \((s,a)\in \mathcal{S}\times\mathcal{A}\), \[\begin{align} \label{eq:kernel95connection} (P V^\pi)(s,a) &= \sum_{s'\in\mathcal{S}} P(s' | s,a) V^\pi(s') = \sum_{s'} P(s'|s,a) \sum_{a'}\pi(a'|s')Q^\pi(s',a') \\ &= \sum_{(s',a')\in\mathcal{S}\times\mathcal{A}} P^\pi((s',a') | (s,a)) Q^\pi(s',a') = (P^\pi Q^\pi)(s,a). \end{align}\tag{9}\] This identity captures the connection between the state transition kernel \(P\) and the induced kernel \(P^\pi\) on the state–action space. Combining 9 and 8 we obtain the Bellman equations in matrix form: \[\label{eq:bellman95in95matrix} Q^{\pi} = r + \gamma P^\pi Q^\pi\;, \quad V^\pi = r^\pi+ \gamma P_\pi V^\pi \;,\tag{10}\] where \(r^\pi(s) = \sum_{a\in\mathcal{A}}\pi(a|s)r(s,a)\). To establish that 10 admits a unique solution, it suffices to verify that the matrix \(I-\gamma P^\pi\) is invertible. Indeed, for any nonzero vector \(x \in \mathbb{R}^{SA}\) we have by the triangle inequality \[\left\|\left(I-\gamma P^\pi\right) x\right\|_{\infty} =\left\|x-\gamma P^\pi x\right\|_{\infty} \geq\|x\|_{\infty}-\gamma\left\|P^\pi x\right\|_{\infty} \geq\|x\|_{\infty}-\gamma\|x\|_{\infty} > 0\;,\] where the last inequality is due to the fact that each component of \(P^\pi x\) is a convex combination of the components of \(x\). This immediately implies that \(I-\gamma P^\pi\) is invertible. The same argument applies to \(I-\gamma P_\pi\), ensuring that the equation for \(V^\pi\) also admits a unique solution.

3.1.1.1 Bellman Operators

For notational convenience, define the Bellman operators \(\mathcal{B}_\pi: \mathbb{R}^{\mathcal{S}} \to \mathbb{R}^{\mathcal{S}}\) and \(\mathcal{B}^\pi: \mathbb{R}^{SA} \to \mathbb{R}^{SA}\) as \[\label{eq:bellman95operators} \mathcal{B}_\pi V := r^\pi+ \gamma P_\pi V\;,\quad \mathcal{B}^\pi Q =r + \gamma P^\pi Q\;.\tag{11}\] The Bellman consistency equations 10 can now be equivalently rewritten as \[V^\pi = \mathcal{B}_\pi V^\pi, \qquad Q^\pi = \mathcal{B}^\pi Q^\pi \;.\] Hence, \(V^\pi\) and \(Q^\pi\) are fixed points of their respective Bellman operators [1]. Moreover, \(\mathcal{B}^\pi\) and \(\mathcal{B}_\pi\) are \(\gamma\)-contractions in the \(\ell_\infty\)-norm. Let us show it for \(\mathcal{B}_\pi\): \[\begin{align} \mathcal{B}_\pi V_1 - \mathcal{B}_\pi V_2 \stackrel{\eqref{eq:bellman95operators}}{=}& \gamma P_\pi (V_1 - V_2), \\ |(\mathcal{B}_\pi V_1)(s) - (\mathcal{B}_\pi V_2)(s)| =\, & \gamma \left| \sum_{s' \in \mathcal{S}} P_\pi(s' \mid s)\,(V_1(s') - V_2(s')) \right| \\ \le\, & \gamma \sum_{s' \in \mathcal{S}} P_\pi(s' \mid s)\, |V_1(s') - V_2(s')| \\ \le\, & \gamma \|V_1 - V_2\|_\infty, \end{align}\] where the last inequality is due to the fact that \(P_\pi(\cdot \mid s)\) is a probability distribution. Due to Banach’s fixed-point theorem, the contraction property guarantees the existence and uniqueness of the fixed points, as well as the convergence of iterative procedures \(V_{k+1} = \mathcal{B}_\pi V_k\) and \(Q_{k+1} = \mathcal{B}^\pi Q_k\) to \(V^\pi\) and \(Q^\pi\), respectively [1].

3.1.2 Existence of deterministic stationary optimal policy↩︎

Define the optimal value function \(V^\star:\mathcal{S}\to \mathbb{R}\) and the optimal action-value function \(Q^\star:\mathcal{S}\times\mathcal{A}\to \mathbb{R}\) as \[V^{\star}(s) :=\sup_{\pi} V^\pi(s)\;, \quad Q^{\star}(s, a) :=\sup_{\pi} Q^\pi(s, a).\]

It turns out that if an optimal policy exists, there exists a deterministic optimal policy [1]. In particular, it exists when \(\mathcal{S}\) is discrete and \(\mathcal{A}\) is finite [1]. Moreover, such a policy can be obtained as a greedy policy with respect to the optimal action-value function \(Q^\star\) [1]. To characterize optimal value functions and policies, we introduce the following operator.

Similarly to operators \(\mathcal{B}^\pi\) and \(\mathcal{B}_\pi\) defined in 11 , the optimality operators are also \(\gamma\)-contractions [1], i.e., \[\label{eq:B95contraction} \|\mathcal{B}V - \mathcal{B}V'\|_\infty \leq \gamma \| V - V'\|_\infty,\qquad \|\mathcal{B}Q - \mathcal{B}Q'\|_\infty \leq \gamma \| Q - Q'\|_\infty.\tag{14}\] By the Banach’s fixed-point theorem, this implies existence and uniqueness of the solutions to the equations \[\label{eq:bellman95opt} V = \mathcal{B}V,\tag{15}\] referred to as the Bellman optimality equations.

One can similarly derive \(Q\)-function counterparts of 15 and ?? :

3.2 Value iteration and policy iteration↩︎

3.2.0.1 Value iteration

The proof of Proposition 1 naturally suggests a procedure based on iteratively applying the contraction mapping \(\mathcal{B}\) to a current estimate of a value function until it converges to \(V^{\star}\). After that, one uses the equation ?? to compute a near-optimal policy. The resulting algorithm goes by the name of Value Iteration (VI) [5] and runs as follows:

Alternatively, one can use the Bellman equations for \(Q\)-function ?? to similarly construct the \(Q\)-value iteration (QVI) algorithm, which performs updates of the form \[\label{eq:qvi} Q_{j+1}(s,a) \gets r(s,a) + \gamma \sum_{s'\in\mathcal{S}} P(s'|s,a) \max_{a'\in\mathcal{A}} Q_{j}(s',a'), \quad (s,a)\in\mathcal{S}\times\mathcal{A}.\tag{16}\]

Before stating the convergence result, note that if rewards are bounded in \([0,1]\), then for every policy \(\pi\) the corresponding value functions satisfy \[\label{eq:bounded95QV} V^\pi(s) \in \left[0,\frac{1}{1-\gamma}\right],\qquad Q^\pi(s,a) \in \left[0,\frac{1}{1-\gamma}\right].\tag{17}\] The same bound holds for \(V^\star\) and \(Q^\star\). This follows from the definitions 1 , 2 and the geometric progression formula.

3.2.0.2 Policy iteration

We have seen that Value Iteration directly improves the value function estimate until it becomes close to optimal. After that, it outputs a greedy policy w.r.t. that near-optimal value function. Alternatively, given some policy, one could first approximate the respective value or \(Q\)-function using the consistency equations 7 , 8 and after that define an improved policy as a greedy policy w.r.t. the estimated value function. These two steps are called policy evaluation and policy improvement. They are repeated until the policy stabilizes. The following implementation of the described algorithm is called Policy Iteration (PI) [7]:

Policy Iteration also enjoys a linear convergence guarantee; for instance, under the exact policy-evaluation scheme above one has bounds of the form \[\|V^{\pi_k}-V^\star\|_\infty \le C\gamma^k,\] and in particular [8] establishes \[\|V^{\pi_k}-V^\star\|_\infty \le \frac{\gamma^k}{1-\gamma}.\]

3.3 Linear programming↩︎

There exist alternative approaches to finding an optimal policy for a given Markov decision process. Consider the case where the MDP \(\mathcal{M}= (\mathcal{S}, \mathcal{A}, P, R, \gamma)\) is fully known; here we assumed that \(P(s'|s,a),\;r(s,a)\;,\gamma\) are given as rational numbers. Using linear programming (LP) yields a polynomial-time algorithm when the MDP is specified in this way, with the computational complexity depending on the description length of \(\mathcal{M}\), assuming all parameters are rational. Recall that in the AMDP setting, the agent aims to find a policy \(\pi^\star\) that maximizes the expected average reward: \[\begin{gather} \rho^{\star}= \max_{\pi} \lim_{H\rightarrow \infty} \frac{1}{H}\mathbb{E} \Bigl[ \sum_{t=0}^{H-1}r(s_t,a_t) \Bigr ]. \end{gather}\] In the finite-horizon case, the limit is dropped and the maximum finite \(H\) is used. For a policy \(\pi\) one obtains a stationary distribution: \[\begin{gather} \nu_{\pi}(s')=\sum\limits_{(s,a)\in\mathcal{S}\times\mathcal{A}}P(s'|s, a)\pi(a|s)\nu_{\pi}(s),\quad s'\in\mathcal{S}, \end{gather}\] which corresponds to its probability vector \(\nu_{\pi} = \left(\nu_{\pi}(s)\right)_{s\in\mathcal{S}}\). If the Markov chain induced by \(\pi\) is ergodic (i.e., irreducible and aperiodic), then \[\begin{gather} V^\pi=\lim_{H\to\infty} \frac{1}{H} \mathbb{E}_\pi \Bigl[\sum_{t=0}^{H-1} r(s_t,a_t)\Bigr]=\sum\limits_{(s,a)\in\mathcal{S}\times\mathcal{A}}r(s,a)\pi(a|s)\nu_\pi(s). \end{gather}\]

Define the state–action distribution \(\mu(s,a)=\nu_{\pi}(s)\pi(a|s)\). The AMDP optimization problem can then be reformulated as a linear program, expressed in terms of evaluating the policy through a distribution vector \(\mu\) belonging to the probability simplex over the state–action space: \[\label{eq:amdp95lp95formulation} \max_{\mu\in\Delta^{\mathcal{S}\times\mathcal{A}} } \Bigl[ V(\mu)= \langle R,\mu\rangle \;: \;\sum_{a\in\mathcal{A}}\mu(s',a)=\sum_{(s,a)}P(s'|s,a)\mu(s,a)\Bigr]\;.\tag{18}\] Let \(\mu^\star\) denote the optimal solution to the linear program 18 , then the optimal policy \(\pi^\star\) can be recovered according to the rule \[\label{eq:dual95to95policy} \pi^\star(a|s)=\frac{\mu^\star(s,a)}{\sum\limits_{a'\in\mathcal{A}}\mu^\star(s,a')}\;.\tag{19}\] Linear program 18 can be rewritten directly in matrix form: \[\label{eq:amdp95lp95matrix} \max\limits_{\mu\in\Delta^{\mathcal{S}\times\mathcal{A}} }\left \langle R,\mu \right \rangle\quad \text{s.t. } (\widehat I - P)^\top \mu = 0,\tag{20}\] where the matrix \(\widehat{I}\in \mathbb{R}^{SA\times S}\) has a nonstandard structure: in each row indexed by \((s,a)\in\mathcal{S}\times\mathcal{A}\), only the entry corresponding to the state \(s \in \mathcal{S}\) equals one, while all other entries are zero. Equalities and inequalities are understood componentwise for scalars, vectors, and matrices. The dual LP corresponds to evaluating the optimal policy through the \(V\)-function: \[\label{eq:amdp95dual95form} \min\limits_{\overline{V}\in\mathbb{R},V\in\mathbb{R}^{S}}\overline{V} \quad \text{s.t. } \;\overline{V}\cdot \vec{1} + (\widehat{I}-P) V\geq R\;.\tag{21}\] This yields average-reward optimality. The scalar \(\rho^\star\) represents the optimal value of the AMDP \(\mathcal{M}\) if and only if there exists a vector \(h \in \mathbb{R}^S\) such that \[\begin{gather} h(s)=\max\limits_{a \in\mathcal{A}} \Bigl( r(s, a) - \rho^\star + \sum\limits_{s'\in\mathcal{S}} P(s'|s, a) h(s') \Bigr). \end{gather}\] Equations can be derived from inequality constraints of the form \[\widehat{I}\cdot V\geq R - \overline{V}\cdot\vec{1}+ P\, V\;.\] Here we can set \(\overline{V} = \rho^{\star}\), since the LP solution lies on the boundary of the feasible set defined by affine constraints. For discounted MDPs, the LP formulation is \[\label{eq:dual95dmpdp} \min\limits_{V\in\mathbb{R}^{S}}\left \langle q,V \right \rangle \quad \text{s.t. } \; (\widehat{I}-\gamma P)V\geq R\;,\tag{22}\] where \(q\) is the initial state distribution vector. The corresponding dual problem is \[\max\limits_{\mu\in\Delta^{\mathcal{S}\times\mathcal{A}} }\left \langle R,\mu \right \rangle\quad \text{s.t. }\; (\widehat{I}-\gamma P)\mu=q\;.\] Although the LP contains \(\mathcal{O}(|\mathcal{S}||\mathcal{A}|)\) affine constraints, its special structure allows efficient solution methods that exploit sparsity and sampling-based approximations. In the DMDP case, one can equivalently reduce constraints in the primal problem as follows: \[\begin{gather} R-(\widehat{I}-\gamma P)V\leq0\;\Longleftrightarrow\;\max\limits_{(s,a)\in\mathcal{S}\times\mathcal{A}}\Big[r(s,a)- \big\langle \widehat{I}_{(s,a)}-\gamma P_{(s,a)},V \big\rangle\Big]\leq0. \end{gather}\] Here \(\widehat{I}_{(s,a)}\) and \(P_{(s,a)}\) denote the rows corresponding to \((s,a)\in\mathcal{S}\times\mathcal{A}\) of \(\widehat{I}\) and \(P\), respectively. Since the full set of constraints uniquely determines \(V\), it suffices to iterate using stochastic optimization methods, minimizing the left-hand side of the inequality, e.g. via stochastic gradient descent. An unbiased stochastic estimate of the gradient of the constraint function \[r(s,a)- \big\langle \widehat{I}_{(s,a)}-\gamma P_{(s,a)},V\big\rangle\] can be expressed as \(\gamma \,\mathbf{e}_{s'} - \widehat{I}_{(s,a)}\), where \(\mathbf{e}_{s'}\) is a one-hot random vector drawn with distribution \(P(\cdot | s,a)\). This approach applies analogously to the LP formulation for AMDP, but to its dual problem, see subsection 4.2, in the paragraph discussing stochastic mirror descent, for further details. There is also an LP formulation for the constrained discounted MDP (CMDP): \[\begin{gather} \begin{matrix} \max\limits_{\mu\in\Delta^{\mathcal{S}\times\mathcal{A}} }\left \langle R,\mu \right \rangle;\\ s.t.\;(\widehat{I}-\gamma P)\mu=q,\;D\mu\geq c. \end{matrix} \end{gather}\] Compared to previous LPs, an additional affine inequality constraint \(D\mu\geq c\) is introduced. Intuitively, each row of \(D\) corresponds to a constraint cost function, and the inequality \(D\mu\ge c\) enforces lower bounds on their expected values under the policy. This means that the optimized policy \(\pi\) must correlate with each such \(\hat{\pi}\) at least to the level specified by the corresponding entry of \(c\). The constraint \(D\mu\geq c\) can also be used for off-policy policy optimization: we not only maximize the reward under \(\pi\), but also enforce that \(\pi\) remains sufficiently similar to an expert policy \(\hat{\pi}\) in terms of correlation specified by \(c\).

3.4 Accelerated methods for solving MDPs↩︎

In many practical scenarios, it is almost as important to take future rewards into account as it is to consider those in the current period, so the discount factor is close to 1. For instance, when seeking the optimal policy via the policy iteration method, one needs to compute the fixed point of the operator associated with a given policy \(\pi\) during the policy evaluation step. Let us remind the expression for this operator: \[\mathcal{B}_{\pi}(V) := r_{\pi} + \gamma P_{\pi}V\;,\] where \(P_\pi\) and \(r_\pi\) defined in 6 and 10 . Due to the operator being a contraction mapping, the inequality \(\|V^{\pi} - V_t\|_{\infty} \leq \gamma^t \|V^{\pi} - V_0\|_{\infty}\) holds. Thus, there is convergence with a linear rate \(\gamma \in (0, 1)\), and as the discount factor gets closer to 1, the Value Computation algorithm, where on every step \(\mathcal{B}_{\pi}\) is applied to the current approximation of the solution, gets too slow. Therefore, many works have addressed this issue and have tried to develop the idea of accelerated schemes analogous to the accelerated methods applied in convex optimization [9]. The key idea lies in the fact that for any \(x, y \in \mathbb{R}^n\), the following estimations hold: \[\|(I - \mathcal{B}_{\pi})x - (I - \mathcal{B}_{\pi})y\|_{\infty} \leq (1 + \gamma) \|x - y\|_{\infty}\] \[(1 - \gamma) \|x - y\|_{\infty} \leq \|(I - \mathcal{B}_{\pi})x - (I - \mathcal{B}_{\pi})y\|_{\infty},\] where \(I\) is the identity mapping. These estimates suggest an analogy with the gradient of a smooth and strongly convex objective: the residual operator \(I-\mathcal{B}_{\pi}\) is well conditioned, with constants depending on \(1-\gamma\) and \(1+\gamma\). Recalling that the value function of policy \(\pi\) is the solution of the Bellman equation \(V^{\pi} - \mathcal{B}_{\pi}(V^{\pi}) = 0\), we obtain an analogy between finding \(V^{\pi}\) and finding the minimizer of a smooth and strongly convex objective. Note, however, that these estimates alone do not imply that \(I-\mathcal{B}_{\pi}\) is the gradient of such a function, since such an interpretation is valid only under additional structural assumptions. With this idea, many algorithms inspired by convex optimization techniques and that can be applied to solving MDPs were proposed. The most basic schemes are presented in [10]. In particular, these include Accelerated Value Computation (A-VC), inspired by Accelerated Gradient Descent [11], [12], and Momentum Value Computation (M-VC), inspired by Momentum Gradient Descent [13]. In the convergence proofs of these algorithms, it is assumed that the Markov chain defined by the policy \(\pi\) is irreducible and reversible. Unfortunately, the existing works do not prove the same for A-VI and M-VI, which differ from the previous schemes mentioned by replacing \(\mathcal{B}_{\pi}\) with \(\mathcal{B}\). One explanation for this lies in the increased complexity of the Bellman operator, which is neither differentiable nor affine. However, it is piece-wise affine, and in [14], it was examined under what conditions the operator \(I - \mathcal{B}\) can be considered as the gradient of some function. Finally, we briefly note another approach to accelerating Value Computation based on some ideas from control theory. The work [15] presents Value Computation with three controllers: PD, PI, and PID. This approach is more general as it does not require the MDP to be reversible, and a good choice of gains can provide significant acceleration. However, tuning these parameters is a separate task.

3.5 Hidden convexity of MDP↩︎

Although the optimization over policies in a Markov decision process is generally nonconvex, the problem admits an equivalent convex reformulation. To observe it, recall from 22 that the dual LP formulation for discounted MDP takes the form \[\max_{\mu \in \mathcal{D}}\, \sum_{s, a} \mu(s, a) r(s,a),\] where \[\mathcal{D}:=\Bigl\{\mu \in \mathbb{R}^{S \times A} : \mu \geq 0,\, \sum_{a\in\mathcal{A}} \mu(s, a)=q(s)+\gamma \sum_{s^{\prime}, a} P\left(s | s^{\prime}, a\right) \mu\left(s^{\prime}, a\right)\, \forall s \in \mathcal{S}\Bigr\},\] and \(q(s)>0\) are probabilities of initial states, and the variable \(\mu\) is referred to as the discounted state–action occupancy measure. Normalizing \(\mu\) by \((1-\gamma)\) yields a probability measure over the state-action space. It turns out ([1], Corollary 6.9.2), there is a one-to-one correspondence between stationary policies and occupancy measures in \(\mathcal{D}\) on the set of states reachable from the initial distribution \(q\). In particular, a policy \(\pi\) can be associated with an occupancy measure \(\mu_q^\pi\in \mathcal{D}\) defined by \[\mu_q^\pi(s,a):= \sum_{t=0}^\infty \gamma^t \mathbb{P}\left(s_t=s, a_t=a \right),\] where the probability is w.r.t. the random trajectory generated by \(s_0\sim q(\cdot)\), \(a_t\sim \pi(\cdot|s_t)\) and \(s_{t+1} \sim P(\cdot \mid s_t, a_t)\). Thus, \(\mu_q^\pi(s,a)\) is the expected discounted number of visits to state-action pair \((s, a)\) when generating the initial state from \(q\) and following policy \(\pi\). Conversely, occupancy measure \(\mu \in \mathcal{D}\) induces policy \[\pi_\mu(a | s) = \frac{\mu(s, a)}{\sum_{a'\in\mathcal{A}} \mu(s, a')}\] whenever \(\sum_{a'}\mu(s,a')>0.\) In particular, an optimal dual solution \(\mu^\star \in \mathcal{D}\) recovers the optimal policy by \(\pi_{\mu^\star}=\pi^\star\) [1].

When maximizing the expected discounted cumulative reward, it is common to promote safe exploration by adding an entropic regularizer [16], i.e., \[\mathcal{H}(\pi)=-\mathbb{E} \bigg[\sum_{t=0}^{\infty} \gamma^{t} \log \pi\left(a_{t} \mid s_{t}\right) \bigg],\] where the expectation is w.r.t. the random trajectory generated by \(s_0\sim q(\cdot)\), \(a_t\sim \pi(\cdot|s_t)\) and \(s_{t+1} \sim P(\cdot \mid s_t, a_t)\). Remarkably, entropy can be viewed as a strictly concave function of occupancy measures:

Thus, entropy-regularized MDP can be formulated as a convex optimization problem over occupancy measures with an affine feasible set and a strictly concave objective. Another scenario is the imitation learning where the goal is to imitate the expert’s behavior by minimizing the KL-divergence between \(\mu_\pi\) and the occupancy measure learned from the expert’s sampled trajectories, which is also a convex objective [18]. More generally, if an optimization problem admits a convex reformulation via non-linear (but invertible) map \(c(\cdot)\), one talks about hidden convexity. Unfortunately, in many cases the map \(c(\cdot)\) and/or its inverse are unavailable or hard to compute. However, under appropriate structural assumptions, stochastic gradient methods applied to the original policy variables may still enjoy global convergence guarantees; see [19].

4 Generative model setting↩︎

We now transition from the full-knowledge setting to the standard reinforcement learning scenario, where planning must be performed under uncertainty. Specifically, we assume that the reward function is known, while the transition probabilities are unknown, as handling unknown rewards is comparatively straightforward. In this section, we assume access to a generative model:

The generative model setting. A generative model is an oracle that, given a state–action pair \((s,a)\), returns an independent sample \(s'\sim P(\cdot\mid s,a)\) together with the reward \(r(s,a)\) (or a sample of the reward if rewards are stochastic). Throughout this section, unless stated otherwise, we assume that \(\mathcal{S}\) and \(\mathcal{A}\) are finite, and rewards are specified by a deterministic function \(r : \mathcal{S}\times\mathcal{A}\xrightarrow[]{} [0,1]\). The goal in generative model setting is to find an \(\varepsilon\)-optimal policy \(\pi\), i.e., a policy satisfying \[\|Q^\pi - Q^{\pi^{\star}}\|_\infty \le \varepsilon.\] Note that since the rewards \(r(s,a)\) lie in the interval \([0,1]\), the value function \(Q^\pi(s,a)\) takes values in \([0, 1/(1-\gamma)]\), so it is meaningful to consider accuracies \(\varepsilon \in \bigl(0, 1/(1-\gamma)\bigr]\). In this work, we focus on constructing PAC algorithms, which for a given tolerance level \(\delta\) output a policy that is \(\varepsilon\)-optimal with probability at least \(1-\delta\). We primarily focus on how the sample complexity depends on the size of the state-action space \(SA\) and the effective horizon \((1-\gamma)^{-1}\), as well as on the desired accuracy \(\varepsilon\) and the confidence level \(\delta\). The effective horizon reflects the typical time scale over which future rewards influence the value function. Two primary approaches arise in this setting: model-based and model-free:

• In the model-based approach, a sufficient number of transitions is first generated, based on which empirical transition probabilities of the Markov process are estimated. Subsequently, optimization methods are applied to the resulting empirical Markov model to compute an approximate optimal policy.

• In the model-free approach, the algorithm directly constructs an approximation of the optimal policy or value function during its execution, without explicitly estimating the transition probabilities. This approximation is iteratively refined as new samples are obtained from the generative model. In other words, an algorithm is called model-free if it updates value functions or policies directly from sampled transitions without explicitly estimating the transition kernel.

4.1 Model-Based Approach↩︎

The simplest model-based algorithms, such as Empirical QVI or Empirical PI, first collect enough samples using a generative model and construct an empirical MDP \({\mathcal{\widehat{M}}} = (\mathcal{S}, \mathcal{A}, \widehat{P}, r, \gamma)\) based on the estimated transition probabilities, see Algorithm 3. Afterwards, any policy iteration [7] or value iteration (QVI) [5] method can be applied to this empirical model \(\mathcal{M}\).

It was shown in [20] that for \(\varepsilon \in \bigl(0, \tfrac{1}{\sqrt{1-\gamma}}\bigr]\), in order to obtain an \(\varepsilon\)-optimal policy \(\pi\), it is sufficient to construct an empirical model of the process based on \[\label{eq:dmdp95samples} K=\widetilde{\mathcal{O}}\left(\frac{S A}{(1-\gamma)^3\varepsilon^2}\right)\tag{23}\] queries to a generative model.¹ Remarkably, the estimate 23 matches (up to logarithmic factors) the lower bound on the number of generated transitions required to guarantee finding an \(\varepsilon\)-optimal solution with probability at least \(1-\delta\) in the worst case setting. This lower bound was obtained in [21].

Note that the estimate 23 implies that at least on the order of \(\tfrac{S A}{(1-\gamma)^2}\) transitions must be generated, since it holds for accuracies \(\varepsilon \in \bigl(0, \tfrac{1}{\sqrt{1-\gamma}}\bigr]\). In many cases, such a sample size may be infeasible, while lower accuracy requirements may still be acceptable. It turns out that the bound 23 can be extended to the entire range \(\varepsilon \in \bigl(0,\tfrac{1}{1-\gamma}\bigr]\) [22]. To achieve this, one can add a small perturbation to the rewards \(r(s,a)\) and then apply, as before, an optimization method to the empirical model (now with perturbed rewards). Under this approach, the required sample size begins at \(\sim \tfrac{S A}{1-\gamma}\) transitions.

Much less literature has focused on the finite-horizon Markov decision process with horizon \(H\) and a generative model. In this case, the analysis and complexity bounds are analogous to the discounted setting. In particular, the lower bound on the number of queries to the oracle matches the upper bound and is given by \[K=\widetilde{\mathcal{O}}\left(\frac{S A H^3}{\varepsilon^2}\right).\] The lower bound was first established in [23] together with the tight upper bound valid only for the range \(\varepsilon\in (0,1)\). The authors of [24] propose an upper bound for \(\varepsilon\in(0,\sqrt{H}]\) in a manner similar to [21] for the discounted case. Evidently, this bound remains valid over the entire interval \(\varepsilon\in(0,H]\) when the reward perturbation technique of [22] is applied.

4.1.0.1 Optimal Control in AMDP

We conclude the discussion of model-based approaches by considering the problem of optimal control in the average-reward MDP (AMDP). In the discounted setting, the problem is in some sense simplified by introducing an effective planning horizon of order \(1/(1-\gamma)\), which yields optimal policies \(\pi_\gamma^\star\). A natural question is whether there exists a universal optimal policy that does not depend on \(\gamma\). It turns out that such a policy is precisely the optimal policy \(\pi^\star\) in the corresponding AMDP. Moreover, as expected, when the planning horizon is extended (i.e., \(\gamma \to 1\)), the discounted optimal policies converge to the AMDP-optimal policy: \(\pi^\star_\gamma \to \pi^\star\).

This observation motivates reduction-based methods: to solve an AMDP, one may choose \(\gamma\) sufficiently close to \(1\) and instead solve the associated discounted MDP. In particular, [23] established an upper bound of order \(\widetilde{\mathcal{O}}\!\left(\frac{t_{\operatorname{mix}}SA}{{\varepsilon^{3}}}\right),\) and a lower bound of order \(\Omega\!\left(\frac{t_{\operatorname{mix}}SA}{{\varepsilon^{2}}}\right),\) thus leaving a gap in the dependence on the accuracy parameter \(\varepsilon\). This gap was subsequently closed by [25]. Both works rely on the reduction to discounted MDPs by setting \(\gamma = 1 - c(\varepsilon/{t_{\mathrm{mix}}})\) for a suitable constant \(c>0\) and \(\bar{\varepsilon} = {\varepsilon}/({1-\gamma}),\) and solving the corresponding DMDP to accuracy \(\bar{\varepsilon}\).

4.1.0.2 A unified view on sample complexity

The sample complexity bounds for different MDP models (discounted, finite-horizon, and average-reward) admit a common interpretation. After an appropriate normalization of the objective function, the problem can be viewed as estimating an expectation by its empirical counterpart. In particular, the relevant quantities can be interpreted as expectations over state–action pairs (and, in the average-reward setting, over time via an improper uniform measure). From this perspective, the sample complexity takes the generic form \[\widetilde{\mathcal{O}}\!\left( \frac{N_{\mathrm{eff}} \cdot H_{\mathrm{eff}}}{\varepsilon_{\mathrm{eff}}^2} \right),\] where \(N_{\mathrm{eff}}\) is an effective factor proportional to the number of \((s, a)\) pairs required to estimate MDP, \(H_{\mathrm{eff}}\) is the effective horizon and \(\varepsilon_{\mathrm{eff}}\) is the normalized accuracy, accounting for the aggregation of rewards over time. The correspondence between different models is summarized in Table 2.²

In all cases, the complexity can be interpreted as the cost of estimating a collection of \(SA \cdot H_{\mathrm{eff}}\) parameters up to accuracy \(\varepsilon_{\mathrm{eff}}\), corresponding to the deviation between empirical and true expectations. To better understand the lower bounds on complexities presented in Table 2, we introduce hard instances — families of MDPs where different parameter settings lead to nearly identical trajectories yet require distinct optimal policies. We first describe a simple hard instance that underlies many of these lower‑bound proofs. The hard instance used in discounted and average‑reward MDP lower‑bounds is a small, bandit‑like system. It has two controllable states and two absorbing states:

• state 0: taking any action keeps the agent in state 0 with high probability but occasionally moves to state 1. State 0 itself yields no reward.

• state 1: this state hides a multi‑armed bandit. Each action \(a\) leads to a terminal state Good with probability \(p_a\) and to Bad with probability \(1-p_a\). Good returns positive reward in \((0, 1]\) specific for each considered setting, Bad returns reward 0, and both transition back to state 1. The optimal action \(a^*\) is the arm with the largest success probability \(p_{a^*}\). The difference between the best and second‑best arms is \(\varepsilon\) so that distinguishing them requires \(\widetilde{\Omega}(\varepsilon^{-2})\) samples. Once the system reaches state 1, it stays there until the next visit to state 0.

The graph for such a MDP is presented in figure 4.

The intuition is summarized in the discounted‑MDP lower‑bound proof of [28]. They note that once a policy reaches state 1 it should choose the action most likely to lead to Good; the process then returns to state 1 and repeats. Estimating which arm is better thus reduces to identifying the optimal arm of a bandit with arms \(a\) and success probabilities \(p_a\); sample‑complexity results for bandits show that any policy must play the inferior arm at least \(\tilde{\Omega}\left(\tfrac{A}{\varepsilon^2}\right)\) times. The discounted MDP lower bound arises by embedding this bandit in an MDP with an additional state 0 so that the effect of a wrong choice at state 1 propagates back through discounting. Below we analyse this hard instance under three performance criteria and derive the lower‑bound scaling.

Consider an infinite‑horizon MDP with discount factor \(\gamma \in (0,1)\) and rewards scaled to \([0,1]\). In the hard instance, state 0 transitions to state 1 with probability \(\alpha = 1 - \gamma\) or proportional to it and remains in state 0 with probability \(1-\alpha\). Let the bandit in state 1 have two actions: \(a^*\) with success probability \(p_{a^*} = p_a+\varepsilon< 1\) and \(a\) with \(p_a\in (0, 1)\). Both actions lead to state 1 through Good or Bad, so the expected return from state 1 under action \(a\) is \(r(a)=(1-\gamma)\cdot p_a\) because good yields \((1 - \gamma)\) and discounting propagates future rewards. The value difference between \(a^*\) and \(a\) at state 1 is therefore

\[V^{*}(1) - Q(1, a) = (1 - \gamma)(p_{a^*} - p_a) = \varepsilon(1-\gamma).\] State 0 remains in the zero‑reward loop for an average of \(\tfrac{1}{1-\gamma}\) steps. During this period the agent plans its next visit to state 1. If its policy chooses the sub‑optimal action at state 1 in those plans, the error propagates back to state 0 and the expected value at state 0 decreases by roughly \(\gamma^{t}(V^*(1)-Q(1,a))\) each time step. Summing over the geometric time spent in state 0 gives a total value gap on the order of

\[\Delta V(0)\approx \sum\limits_{t = 0}^{\infty}\frac{1}{1-\gamma}\,\gamma^{t}\,(1-\gamma)\cdot \varepsilon= \frac{\varepsilon}{1-\gamma}.\] Distinguishing a difference of size \(\Delta V(0)\) requires \(\Omega\left(\tfrac{1}{\Delta V(0)^2}\right)\) samples. Because the bandit problem at state 1 already requires \(\Omega(\varepsilon^{-2})\) samples to tell whether \(p_{a^*}\) is larger than \(p_a\), the cumulative effect of discounting multiplies the sample requirement by \(\tfrac{1}{(1-\gamma)^2}\). A careful application of Fano’s inequality [29] shows that any algorithm needs at least

\[\widetilde{\Omega}\Bigl(\tfrac{SA}{\varepsilon^{2}\,(1-\gamma)^{3}}\Bigr)\] queries to the generative model to learn an \(\varepsilon\)‑optimal policy [30]. This matches the minimax lower bound for discounted MDPs up to logarithmic factors. Upper bounds based on empirical Q‑value iteration achieve \(\tilde{O}\left(\tfrac{SA}{\varepsilon^2(1-\gamma)^3}\right)\) samples. The key intuition here is discounting attenuates future rewards, but the error incurred when picking the wrong arm at state 1 accumulates while the agent stays in state 0 for roughly \(\tfrac{1}{1-\gamma}\) steps. This yields an additional factor \(\tfrac{1}{1-\gamma}\) in the sample complexity beyond the bandit lower bound.

Average‑reward MDPs optimize the long‑term average reward rather than discounted return. Without discounting, the Markov chain must mix to a stationary distribution. To construct a hard instance, one modifies the previous MDP so that state 0 randomly "restarts" the chain: state 0 transitions to state 1 deterministically \((\alpha = 1)\), but a small probability \(\beta=\Theta\left(\tfrac{1}{t_{\mathrm{mix}}}\right)\) in bandit state causes a reset to a special reset state. This ensures that under any policy the induced Markov chain has mixing time at most \(t_{\mathrm{mix}}\) [30]. The graph of hard AMDP instance is in figure 5. For \(p_a = 1 - \beta\), \(p_{a^*} = p_a + \varepsilon\beta < 1\) one can show the expected average reward difference between choosing \(a^*\) and \(a\) in state 1 according to Lemma 5 in [30]: \[\begin{align} &{\frac{\varepsilon(1 - \beta)}{(1 + \beta)(2 - \beta)(2 - \beta - \varepsilon(1 - \beta))}.}\\ \end{align}\] Basicly, the solution of such AMDP lies on distinguishing between two Bernoulli distributions with means \((1 - \beta)p_a\) and \((1 - \beta)p_{a^*}\). Due to Lemma 7 in [30] it requires the following amount of samples: \[\begin{align} &{\widetilde{\Omega}\left(\frac{1}{(1 - \beta)^2(p_{a^*} - p_a)^2}\right) = \widetilde{\Omega}\left(\frac{1}{\varepsilon^2\beta^2}\right).}\\ \end{align}\] However, in a generative model one can simulate multiple steps of the Markov chain, waiting long enough between samples to let the chain mix so that successive rewards are nearly independent. Because the chain mixes in \(t_{\mathrm{mix}} = \Theta\left(\beta^{-1}\right)\) steps, one can thin the simulation, taking one reward every \(O\left(t_{\mathrm{mix}}\right)\) transitions, and obtain effectively independent observations whose average is \(\mu_{\pi}(1)(1 - \beta)p_{a}\) and \(\mu_{\pi}(1)(1 - \beta)p_{a^*}\) respectively, where \(\mu_{\pi}(1)\) is the stationary probability of being in state 1. To estimate these averages within \(\varepsilon\), , the number of independent observations needed is \(\widetilde{\Omega}\left(\varepsilon^{-2}\right)\). Each independent observation costs \(O(t_{\mathrm{mix}})\) transitions, so the total sample complexity to distinguish between \(a\) and \(a^*\) is \(\widetilde{\Omega}\left(\frac{t_{\mathrm{mix}}}{\varepsilon^2}\right)\). Repeating this analysis for each state-action pair we get the desired estimate: \[\widetilde{\Omega}\bigl(\tfrac{SA\,t_{\mathrm{mix}}}{\varepsilon^{2}}\bigr).\] The work [30] provides matching upper and lower bounds for AMDPs: algorithms based on primal‑dual methods achieve this sample complexity, and their lower bound shows that linear dependence on the mixing time is unavoidable.

In finite‑horizon MDPs each episode lasts exactly \(H\) steps. The agent collects rewards only within an episode; after \(H\) steps the environment resets. To obtain a hard instance, replicate the bandit MDP across the \(H\) stages: each stage \(h\in\{1,\dots,H\}\) has its own state \(h\) where the agent chooses an action, transitions to an absorbing Good or Bad state, and then deterministically moves to stage \(h+1\). The Good state gives immediate reward equal \(\frac{1}{H}\). The optimal arm yields a slightly larger probability of reaching the Good state, producing an extra \(\varepsilon\) expected reward at that stage, for example: \(p_a\in (0, 1)\), \(p_{a^*} = p_a + \varepsilon< 1\). The total reward over an episode under the optimal policy is \(H\) times this stage‑wise reward. A sub‑optimal policy that chooses a wrong arm in stage \(h\) sacrifices \(\frac{\varepsilon}{H}\) immediate reward but does not affect future stages because the horizon is finite. The graph of this hard instance is presented in figure 6. Because the reward gap at a single stage is \(\frac{\varepsilon}{H}\), a bandit lower bound implies that one needs \(\widetilde{\Omega}\left(\frac{H^2}{\varepsilon^2}\right)\) samples per stage to identify the best arm. However, in an episodic setting the value difference accumulates over the \(H\) stages: the total return difference between the optimal and sub‑optimal policies is \(\varepsilon\) per episode. To identify an arm whose expected return is within \(\varepsilon\) of optimal, one must detect a gap of size \(\tfrac{\varepsilon}{H}\) in the per‑stage reward. Fano’s inequality therefore yields a lower bound on the number of episodes of order \(\widetilde{\Omega}\left(\tfrac{H^2}{\varepsilon^2}\right)\). Since each episode contains \(H\) steps, the overall sample complexity (number of state–action pairs sampled) is

\[\widetilde{\Omega}\bigl(\tfrac{SA\,H^{3}}{\varepsilon^{2}}\bigr).\] The work [23] shows the tightness of this bound. For episodic hard instance the intuition means errors do not propagate across episodes but the horizon \(H\) scales the value difference. Finding an \(\varepsilon\)‑optimal policy thus requires \(H^3\) dependence: \(H\) from the accumulated reward gap of size \(\tfrac{\varepsilon}{H}\) at each stage and \(H^2\) from the number of episodes needed to estimate the overall gap.

4.2 Model-free approach↩︎

In this subsection, we discuss the class of model-free algorithms. We first consider Bellman-based methods, most notably \(Q\)-learning and its variants that attain the minimax lower bound in [21]. We then turn to a linear programming perspective, in which the problem admits a saddle-point formulation, and demonstrate that stochastic mirror descent methods provide a natural algorithmic framework in this setting.

4.2.0.1 Q-learning

Assume that a generative model of the MDP is available. The \(Q\)-learning algorithm maintains a sequence of vectors \(Q_t \in \mathbb{R}^{SA}\), initialized arbitrarily with \(\|Q_0\|_\infty \le (1-\gamma)^{-1},\) and for \(t > 0\) updated according to \[\label{eq:q95learning95update} Q_{t} = (1 - \alpha_t) Q_{t-1} + \alpha_t \widehat{\mathcal{B}}_t Q_{t-1}\;,\tag{24}\] where \(\alpha_t \in [0,1)\) denotes the step size, and \(\widehat{\mathcal{B}}_t : \mathbb{R}^{SA} \to \mathbb{R}^{SA}\) is the empirical Bellman operator defined by \[\widehat{\mathcal{B}}_t Q(s,a) = r(s,a) + \gamma \max_{a'\in\mathcal{A}} Q(s_t,a'),\] with \(s_t \sim P(\,\cdot\mid s,a)\). Clearly, \(\widehat{\mathcal{B}}_t\) is an unbiased estimator of the optimal Bellman operator. Equivalently, in matrix form one may write \(\widehat{\mathcal{B}}_t Q_{t-1} = r+\gamma\widehat{P}_t V_{t-1}\), where \(V_{t-1}(s) = \max_{a\in\mathcal{A}} Q_{t-1}(s,a)\), and \(\widehat{P}_t(s'| s,a) = \mathbb{I}\{s' = s_t(s,a)\}\) provides an unbiased estimate of the transition kernel.

We now outline the key ideas underlying the convergence analysis. Letting \(\Delta_t = Q_t - Q^\star_t\) denote the approximation error, one obtains the recursion \[\label{eq:q95learning95martingale95error95recursion} \begin{align} &\Delta_t = A_t\Delta_{t - 1} + \alpha_t\varepsilon_t, \end{align}\tag{25}\] where \(A_t = I - \alpha_t(I - \gamma P^{\pi^\star})\) and \(\pi^\star\) denotes an optimal policy. Iterating this relation with \(\Gamma_{m:n} = \prod_{i = m}^nA_i\) yields \[\label{eq:q95learning95dynamics} \begin{align} &\Delta_{t} = \Gamma_{1:t}\Delta_0 + \sum\limits_{j = 1}^t\alpha_j\Gamma_{j+1:t}\varepsilon_t. \end{align}\tag{26}\] The first term in 26 represents the transient component of the error, which governs the rate at which the initial error \(Q_0 - Q^\star\) is forgotten. The second term corresponds to the fluctuation component, arising from the stochastic variability of the iterates \(Q_t\) around the fixed point \(Q^\star\). The proof of convergence rests on two ingredients: the contraction properties of the matrices \(I - \alpha_t(I - \gamma P^{\pi^\star})\), and concentration inequalities for vector-valued martingales (notably Freedman’s inequality [31]). Here we have an actual martingale: \[\begin{align} &\varepsilon_t = \mathcal{\widehat{B}}_tQ_{t - 1} - \mathcal{B}Q_{t - 1} = \gamma\left(\widehat{P}_t - P\right)V_{t - 1}, \end{align}\] where \(\widehat{P}_t\) is an empirical Markov kernel at iteration \(t\). Below we show it. Let \(\Omega\) be the sample space of all possible sequences of transition samples generated by the generative model, let \(\mathcal{F}\) be the corresponding product \(\sigma\)-algebra, and let \(\mathbb{P}\) be the probability measure induced by the transition kernel \(P\). The natural filtration of \(Q\)-learning algorithm: \[\mathcal{F}_{t - 1} := \sigma\!\left( Q_0,\widehat P_1,\widehat P_2,\ldots,\widehat P_{t - 1} \right), \qquad t\geq1,\] with the convention \[\mathcal{F}_0:=\sigma(Q_0).\] A filtered probability space is a tuple \[(\Omega,\mathcal{F},\{\mathcal{F}_t\}_{t\ge0},\mathbb{P}),\] where \((\Omega,\mathcal{F},\mathbb{P})\) is a probability space and \[\mathcal{F}_0\subseteq \mathcal{F}_1\subseteq \mathcal{F}_2\subseteq \cdots \subseteq \mathcal{F}\] is an increasing sequence of sub-\(\sigma\)-algebras. The \(\sigma\)-algebra \(\mathcal{F}_t\) represents all information available up to time \(t\). A stochastic process \(\{X_t\}_{t\ge0}\), with \(X_t\in\mathbb{R}^d\), is called adapted to the filtration \(\{\mathcal{F}_t\}_{t\ge0}\) if \(X_t\) is \(\mathcal{F}_t\)-measurable for every \(t\ge0\). It is called a martingale with respect to \((\{\mathcal{F}_t\}_{t\ge0},\mathbb{P})\) if \[X_t \text{ is } \mathcal{F}_t\text{-measurable}, \qquad \mathbb{E}[\|X_t\|]<\infty, \qquad \mathbb{E}[X_{t+1}\mid\mathcal{F}_t]=X_t\] for every \(t\ge0\). Equivalently, the conditional expected future value of the process, given the current information, is its current value.

A stochastic process \(\{\xi_{t+1}\}_{t\ge0}\), with \(\xi_{t+1}\in\mathbb{R}^d\), is called a martingale-difference sequence with respect to \(\{\mathcal{F}_t\}_{t\ge0}\) if \[\xi_{t+1} \text{ is } \mathcal{F}_{t+1}\text{-measurable}, \qquad \mathbb{E}[\|\xi_{t+1}\|]<\infty, \qquad \mathbb{E}[\xi_{t+1}\mid\mathcal{F}_t]=0\] for every \(t\ge0\). In this case, \[M_t:=\sum_{k=1}^{t}\xi_k, \qquad M_0:=0,\] is a martingale with respect to \(\{\mathcal{F}_t\}_{t\ge0}\).

In particular, \(Q_t\) is \(\mathcal{F}_t\)-measurable, \(\widehat P_{t}\) is conditionally independent of \(\mathcal{F}_{t - 1}\) given \(P\), and \[\mathbb{E}[\widehat P_{t}\mid\mathcal{F}_{t - 1}]=P.\] Consequently, \[\mathbb{E}[\widehat{\mathcal{B}}_{t}Q_{t - 1}\mid\mathcal{F}_{t - 1}] = r+\gamma P V_{t - 1} = \mathcal{B}Q_{t - 1}.\] Therefore the stochastic Bellman error \[\varepsilon_{t} := \widehat{\mathcal{B}}_{t}Q_{t - 1}-\mathcal{B}Q_{t - 1} = \gamma(\widehat P_{t}-P)V_{t - 1}\] is a bounded martingale-difference sequence with respect to \(\{\mathcal{F}_t\}_{t\ge0}\): \[\mathbb{E}[\varepsilon_{t}\mid\mathcal{F}_{t - 1}]=0,\quad \|\varepsilon_t\|_{\infty}\leq\gamma\left(\|\widehat{P_t}V_{t - 1}\|_{\infty} + \|PV_{t - 1}\|_{\infty}\right)\leq 2\gamma\|V_{t - 1}\|_{\infty}\leq\frac{2\gamma}{1 - \gamma}.\] Equivalently, the cumulative noise process \[M_t:=\sum_{k=1}^{t}\alpha_k\varepsilon_k, \qquad M_0:=0,\] is a martingale whenever the step sizes \(\alpha_k\) are deterministic, or more generally \(\mathcal{F}_{k-1}\)-measurable, and \[\mathbb{E}[\|\alpha_k\varepsilon_k\|]<\infty \qquad \text{for all }k\ge1.\] To establish \(Q\)-learning algorithm convergence almost surely the one thing left is to apply a martingale contraction lemma:

Now, by taking the recursion 25 is of the form ?? with \[u_t=\Delta_t, \qquad g_t=h_t, \qquad \zeta_{t+1}=\xi_{t+1}, \qquad \beta=\gamma, \qquad c_t\equiv0\] all assumptions of Lemma 2 hold. Consequently, for \(t\rightarrow\infty\) \[\Delta_t\to0 \qquad \text{almost surely.}\] Equivalently, when \(t\rightarrow\infty\) \[\|Q_t-Q^\star\|_\infty\to0 \qquad \text{almost surely,}\] which proves \(Q\)-learning convergence. The convergence of \(Q\)-learning algorithm also establishes an almost surely asymptotic optimality of the corresponding greedy policy: \[\pi_t(s)\in\mathop{\mathrm{arg\ max}}_{a\in\mathcal{A}}Q_t(s,a),\qquad\pi_t\rightarrow\pi^\star\text{ as }t\rightarrow\infty.\] or \[\|V^{\pi_t}-V^\star\|_\infty\to0 \qquad \text{almost surely.}\] To prove it one has to bound the norm above. Since \(\pi_t\) is greedy with respect to \(Q_t\): \[\mathcal{B}^{\pi_t}Q_t=\mathcal{B}Q_t.\] Using \(Q^\star=\mathcal{B}Q^\star\) and \(Q^{\pi_t}=\mathcal{B}^{\pi_t}Q^{\pi_t}\), we have \[\begin{align} \|Q^\star-Q^{\pi_t}\|_\infty &\le \|\mathcal{B}Q^\star-\mathcal{B}Q_t\|_\infty + \|\mathcal{B}^{\pi_t}Q_t-\mathcal{B}^{\pi_t}Q^{\pi_t}\|_\infty \\ &\le \gamma\|Q^\star-Q_t\|_\infty + \gamma\|Q_t-Q^{\pi_t}\|_\infty \\ &\le 2\gamma\|Q_t-Q^\star\|_\infty + \gamma\|Q^\star-Q^{\pi_t}\|_\infty. \end{align}\] Rearranging gives \[\|Q^\star-Q^{\pi_t}\|_\infty \le \frac{2\gamma}{1-\gamma}\|Q_t-Q^\star\|_\infty.\] The right-hand side converges to zero almost surely as we recently proved. Since \[\begin{align} \|V^{\pi_t} - V^\star\|_{\infty} &= \max\limits_{s\in\mathcal{S}}\left\{V^{\star}(s) - V^{\pi_t}(s)\right\} \\ &= \max\limits_{s\in\mathcal{S}}\left\{\max\limits_{a\in\mathcal{A}}\left\{Q^\star(s, a)\right\} - \max\limits_{a'\in\mathcal{A}}\left\{Q^{\pi_t}(s, a')\right\}\right\} \\ &\leq\max\limits_{s\in\mathcal{S}}\max\limits_{a\in\mathcal{A}}\left\{Q^\star(s, a) - Q^{\pi_t}(s, a)\right\} = \|Q^{\pi_t} - Q^\star\|_{\infty} \end{align}\] the claim follows.

The work [32] established that \(Q\)-learning achieves a sample complexity of \[\label{eq:q95learning95sample95complexity95generative95model} \begin{align} &\widetilde{\mathcal{O}}\left(\frac{SA}{(1-\gamma)^4 \varepsilon^2}\right)\;, \end{align}\tag{30}\] when using constant or linearly rescaled step sizes. In [33], \(Q\)-learning is analyzed within the more general framework of stochastic approximation with contractive operators, yielding slightly weaker bounds. Finally, we note that classical Q-learning exhibits suboptimal dependence on the effective horizon compared to model-based planning methods, as discussed in [32].

4.2.0.2 Phased Q-learning

In the work [34], a modification of Q-learning algorithm was proposed. The modification is called phased Q-learning and its main distinction from the original algorithm is that it works in phases, as the name suggests, and averages values of the next states observed during a phase. A phase \(\ell\) consists of several calls to the generative model for every state-action pair \((s,a)\) and an update. Let \(m\) be the number of such calls, then the update reads as \[{Q}_{\ell+1}(s, a)=r(s,a)+ \frac{\gamma}{m} \sum_{k=1}^{m} \max_{a'\in\mathcal{A}}Q_{\ell}(s_k^\ell,a'),\] where \(s_1^{\ell}, \ldots, s_{m}^{\ell}\) are the \(m\) next states observed from \((s, a)\) on the calls to the generative model during the \(\ell\)-th phase. Setting the epoch length to \(m = \widetilde{\mathcal{O}}(\varepsilon^{-2}(1-\gamma)^{-4})\) the algorithm computes \(\varepsilon\)-optimal policy, \(\varepsilon \in \bigl(0, \frac{1}{1-\gamma}\bigr]\), after \(\widetilde{\mathcal{O}}\left(\frac{S A}{(1-\gamma)^7\varepsilon^2}\right)\) calls to the generative model, which is one of the first finite-sample convergence results of this type.

4.2.0.3 Polyak-Ruppert Averaged \(Q\)-learning

To attain the theoretical lower bound on sample complexity, several modifications of standard \(Q\)-learning have been proposed. A prominent example is the Polyak-Ruppert averaged \(Q\)-learning algorithm. The averaged iterate is defined as \[\bar{Q}_t = \frac{1}{t}\sum_{i=1}^{t} Q_i,\] where \(\{Q_i\}_{i\ge 1}\) are updated according to the standard \(Q\)-learning recursion 24 . [35] establish that the averaged \(Q\)-learning iterate is minimax optimal in the sense of achieving the lower bound on the expected \(\infty\)-norm error. Furthermore, a central limit theorem is established for the averaged iterates: \[\sqrt{T}\bigl(\bar{Q}_T - Q^\star\bigr) \xrightarrow{d} \mathcal{N}(0, \Sigma),\] where the asymptotic covariance matrix \(\Sigma \in \mathbb{R}^{SA \times SA}\) is given by \[\Sigma = (I - \gamma P^{\pi^\star})^{-1} \operatorname{Var}(Z_1) (I - \gamma P^{\pi^\star})^{-\top}, \quad Z_1 = (P - \widehat{P}_1) V^\star,\] where \(\widehat{P}_1\) is the empirical Markov transition kernel evaluated at the first iteration. The result highlights that averaging stabilizes the stochastic fluctuations of \(Q\)-learning. In particular, the central limit theorem quantifies the asymptotic distribution of the estimation error, providing precise confidence intervals for each component of \(Q^\star\) in the large-sample regime.

4.2.0.4 Variance-Reduced \(Q\)-learning

Another approach to reducing sample complexity is the variance-reduction technique proposed in [36]. We outline the fundamental idea behind variance reduction. Starting from the standard \(Q\)-learning update \[Q_t = (1-\alpha_t) Q_{t-1} + \alpha_t \widehat{\mathcal{B}}_t Q_{t-1},\] suppose that, in principle, we could compute both an empirical \(\widehat{\mathcal{B}}_t(Q^\star)\) and the exact Bellman update \(\mathcal{B}(Q^\star)\). In this case, one could implement the recentered update \[\label{eq:ideal95update} Q_t = (1-\alpha_t) Q_{t-1} + \alpha_t \Big(\widehat{\mathcal{B}}_t Q_{t-1} + \mathcal{B} Q^\star - \widehat{\mathcal{B}}_t Q^\star \Big).\tag{31}\] Defining the error \(\Delta_t = Q_t - Q^\star\), we obtain \[\Delta_t = (1-\alpha_t) \Delta_{t-1} + \alpha_t \widehat{\mathcal{B}}_t( \Delta_{t-1}) .\] By the contractivity of the Bellman operator and the unbiasedness of \(\widehat{\mathcal{B}}_t\), one obtains \(\|\Delta_t\|_{\infty} \leq (1 - (1-\gamma)\alpha_t)\|\Delta_{t-1}\|_{\infty}\). Consequently, if this idealized algorithm could be executed with a constant step size, the error would decay geometrically. This illustrates the key intuition behind variance-reduced \(Q\)-learning: by centering the updates around the optimal \(Q^\star\), the variance of the stochastic updates is significantly diminished, leading to faster convergence. The preceding observations can be formalized as follows. Variance-reduced \(Q\)-learning is executed over \(M\) epochs, each of fixed length \(T\). In the \(m\) -th epoch, the operator \(\widetilde{\mathcal{B}}\) is computed using \(N_m\) samples per state-action pair and serves as an estimate of the Bellman operator \(\mathcal{B}\) in formula 31 . The proxy for \(Q^\star\), denoted \(\bar{Q}\), is set to be the output from the previous epoch, providing a reference for the current epoch’s updates.

Algorithm 7 with appropriate choice of \(N_m\) and \(T\) achieves the minimax-optimal sample complexity 23 .

4.2.0.5 Markovian Sample Trajectory

Unfortunately, the generative model framework is primarily of theoretical interest. In practice, the assumption that one can sample transitions from any state–action pair at will is unrealistic. Typically, an agent interacts with the environment sequentially, collecting data along a Markovian trajectory \(\tau = \{(s_t,a_t,r_t)\}_{t=0}^{\infty}\), which is generated under a fixed stationary behavior policy \(\pi_b\): \[a_t \sim {\pi_b}(s_t), \quad r_t = r(s_t, a_t), \quad s_{t+1} \sim P(\cdot |s_t, a_t).\] The behaviour policy can be arbitrary yet with the positive probability to sample an arbitrary action \(a_t\in\mathcal{A}\) at state \(s_{t}\in\mathcal{S}\) for finite iteration \(t\in\mathbb{Z}_{+}\) to make the whole process convergent. Upon observing the transition tuple \((s_t,a_t,r_t,s_{t+1})\), the update takes the form \[Q_{t+1}(s_{t}, a_{t}) =Q_{t}(s_{t}, a_{t}) + \alpha_t \left(r_t + \gamma \max_{a\in\mathcal{A}} Q_{t}(s_{t+1}, a) - Q_t(s_t,a_t)\right)\;,\] while leaving other entries unchanged. In this online setting, the agent observes correlated transitions along a single trajectory and updates its estimates incrementally. Consequently, algorithms must handle the inherent dependencies between consecutive samples, as well as potential distributional mismatch when learning a target policy \(\pi\) that differs from the behavior policy \(\pi_b\), giving rise to the classical off-policy learning problem.

To ensure that the agent collects sufficient information about every transition \(P(s'|s,a)\), we assume that The Markov chain induced by the behavior policy \(\pi_b\) is uniformly ergodic. Two crucial quantities that dictate the performance of asynchronous Q-learning are the minimum state-action occupancy probability \[\mu_{\text{min}} = \min_{(s,a)\in\mathcal{S}\times\mathcal{A}} \mu_{{\pi_b}}(s,a)\;,\] the mixing time \(t_{\operatorname{mix}}\), which characterizes how quickly the distribution of visited states along the trajectory approaches the stationary distribution and is defined as the smallest integer \(t\) such that \[\sup_{s\in\mathcal{S}} \| P_{\pi_b}^t(\cdot|s) - \nu_{\pi_b} \|_{\mathrm{TV}} \le \tfrac14.\] [32] demonstrates that with high probability, the total sample size needed for asynchronous Q-learning is \[\label{q95upper95bound} \widetilde{\mathcal{O}}\left(\frac{\mu_{\operatorname{min}}^{-1}}{(1-\gamma)^4 \varepsilon^2} + \frac{\mu_{\operatorname{min}}^{-1}t_{\operatorname{mix}}}{(1-\gamma)}\right)\;,\tag{32}\] provided that the learning rates are taken to be some proper constant \(\alpha_{t} = \alpha\). The bound in 32 is comparable to the corresponding result in the generative model setting, with the factor \(\mu_{\operatorname{min}}^{-1}\) playing a role analogous to the size of the state-action space \(SA\).

4.2.0.6 Stochastic Mirror Descent

The first results in this direction were obtained in [37]. Here, however, we focus on the setting studied in [38], with particular attention to the average-reward MDP formulation, while noting that analogous results also hold in the discounted case. We begin by recalling (see 20 , 21 ) that the problem of computing an optimal policy can be cast as the linear program \[\min_{\overline{V},\, V}\overline{V} \quad \text{subject to } (\widehat{I} - P)V + \vec{1}\cdot \overline{V} \;\geq\; R\;,\] whose dual is \[\max_{\mu{\in\Delta^{S\times A}}}\;\langle\mu, R\rangle \quad \text{subject to } {(\widehat{I} - P)^\top\mu} \;=\; 0\;.\] By standard linear programming duality, the problem can be reformulated using Lagrangian multipliers as a bilinear saddle-point problem. For AMDPs, the minimax formulation is given by \[\label{eq:saddle95point95formulation} \min_{V,\overline{V}}\max_{\mu{\in \Delta^{S\times A}}} f(\overline{V}, V, \mu),\tag{33}\] where \[\begin{align} f(\overline{V},V,\mu) &:= \overline{V} + \mu^\top\big( - \overline{V} \cdot\vec{1} + (P - \widehat I)V +R\big).\\ &= \mu^\top\big((P - \widehat I)V +R\big). \end{align}\] The central idea of [38] is to view 33 within the general framework of \(\ell_\infty\)-\(\ell_1\) bilinear games [39]. Such games arise when one player minimizes a bilinear function over a box domain (\(\ell_\infty\)), while the other maximizes over a simplex domain (\(\ell_1\)): \[\label{eq:bilinear} \min_{x \in \mathbb{B}_b^n} \max_{y \in \Delta^m} f(x,y) = y^\top Mx + b^\top x - c^\top y,\tag{34}\] where \(\mathbb{B}_b^n = b \cdot [-1,1]^n\) denotes the \(\ell_\infty\) box constraint, and \(\Delta^m\) denotes the probability simplex. To solve 34 , stochastic mirror descent is applied with divergence functions \[V_x(x') = \tfrac{1}{2}\|x - x'\|^2, \qquad V_y(y') = \operatorname{KL}(y \,\|\, y'),\] corresponding respectively to the Euclidean and Kullback–Leibler divergences. At iteration \(t\), given stochastic gradients \(\mathbb{E}[g_t^x|x_t, y_t] = f_x(x_t,y_t)\) and \(\mathbb{E}[g_t^y|x_t, y_t] = {-f_y(x_t,y_t)}\), the updates take the form \[x_{t+1} = \arg\min_{x \in \mathbb{B}_b^n} \bigl\{\langle \alpha g_t^x, x\rangle + V_{x_t}(x)\bigr\}, \qquad y_{t+1} = \arg\min_{y \in \Delta^m} \bigl\{\langle \beta g_t^y, y\rangle + V_{y_t}(y)\bigr\}.\] It remains to specify how stochastic gradients of 33 can be constructed under access to a generative model. The key observation is that the saddle-point objective admits unbiased gradient estimators obtained via local sampling. For the primal variable \(V\), the gradient is given by \(\nabla_V f(\overline{V}, V, \mu) = {\mu^\top (P - \widehat{I})}\;.\) To obtain a stochastic estimate, we first sample a state-action pair \((s,a) \sim \mu\) and then draw \(s' \sim P(\cdot | s,a)\). An unbiased estimator is \[g^V = \mathbf{e}_{s'} - \mathbf{e}_s\;.\] For the dual variable \(\mu\), the gradient is given by \(\nabla_\mu f(\overline{V}, V, \mu) = {(P - \widehat{I} )V + R\;.}\) Since the dual player maximizes \(f\), while the mirror-descent update above is written in descent form, we use an unbiased estimator of \(-\nabla_\mu f\). In this case, we sample uniformly from the state-action space, \((s,a) \sim \mathrm{Unif}(\mathcal{S}\times\mathcal{A})\), then draw \(s' \sim P(\cdot | s,a)\), and set \[{g^\mu \;=\; SA \bigl(V(s) - V(s') - r(s,a)\bigr)\,\mathbf{e}_{s,a}\;,}\] where the prefactor \(SA\) corrects for uniform sampling and \(\mathbf{e}_{s,a}\) denotes the canonical basis vector for \((s,a)\). This construction illustrates a remarkable feature of the mirror descent framework: despite the apparent complexity of the bilinear formulation, each stochastic gradient step can be implemented using only local one-step samples from the generative model, without the need to reconstruct global transition probabilities.

It is worth remarking on two complementary lines of work. The study [40] examines closely related reductions of planning in MDPs to \(\ell_1\)-regression and linear programming, and develops nearly–linear time algorithms for dense instances; this line of research supplies an algorithmic alternative that is particularly attractive when the underlying LPs are dense and admit efficient dynamic preconditioning. By contrast, [41] investigates a framework that combines linear programming with linear function approximation and employs stochastic primal–dual optimization to produce a compact policy representation via entirely offline computation.

5 Policy Evaluation↩︎

We temporarily shift our focus from the policy control problem to the policy evaluation problem, that is, the estimation of the value function \(V^{\pi}(s)\). Policy evaluation plays a central role in reinforcement learning, as it provides a quantitative assessment of the long-term performance of a given policy and serves as a key component of policy improvement methods. In particular, many algorithms (such as policy iteration and actor-critic methods) rely on repeated evaluation of intermediate policies. It is also essential in off-policy settings, where one seeks to evaluate a target policy using data generated by a different behavior policy, a situation that frequently arises in recommender systems.

5.1 Temporal difference (TD) learning↩︎

TD methods form the basis of many modern reinforcement learning algorithms. They are statistically efficient, update estimates online, and extend naturally to non-tabular approximations. We begin with the classical tabular setting and then move to the more general case of MDPs with linear function approximation.

We begin by recalling that, in the full-knowledge setting, the state–value function \(V^\pi\) is the unique solution of the linear system \[\label{eq:consistency95equation} (I-\gamma P_{\pi})V^\pi = r^\pi\;,\tag{35}\] where \(P_\pi\) and \(r^\pi\) defined in 6 and 10 . Our objective is to recover \(V^{\pi}\) when the kernel \(P_{\pi}\) is unknown and only sample transitions are available. We consider this problem under general Linear Stochastic Approximation (LSA) framework.

5.1.0.1 Linear Stochastic Approximation

LSA aims to solve a linear system \(A\theta = b\) with a unique solution \(\theta^\star\). We do not have access to the coefficients \(A\) and \(b\) but instead we have access to a sequence of observations \(\{\mathbf{A}(X_n), \mathbf{b}(X_n)\}_{n\in\mathbb{N}}\), where \((X_n)_{n\in\mathbb{N}}\) are noise variables taking values in a Polish space \(\mathcal{X}\) and \(\mathbf{A}: \mathcal{X}\rightarrow \mathbb{R}^{d\times d}\), \(\mathbf{b}: \mathcal{X}\rightarrow \mathbb{R}^{d}\) are measurable functions. For a fixed step size \(\alpha>0\) and initialization \(\theta_0\in \mathbb{R}^d\), we consider the sequences of LSA iterates \(\{\theta_n\}_{n\in\mathbb{N}}\) given by \[\theta_{k+1} = \theta_{k} - \alpha (\mathbf{A}(X_{k+1})\theta_{k} - \mathbf{b}(X_{k+1}))\;.\] LSA is a well-studied area of research [42]–[44]. To ensure convergence to the solution, it is standard to assume that the matrix \(-A\) is Hurwitz. Two common types of noise variables are typically considered: independent (i.i.d.) and Markovian.

The analysis of the Markovian case is more challenging; however, it is often worthwhile, as Markov noise sequences arise more frequently in practical applications. The complete TD algorithm is as follows:

To clarify the connection between TD(0) and LSA, we rewrite update rule in matrix form: \[V_{k} = V_{k-1} - \alpha\bigg(\underbrace{\mathbf{e}_{s_k}(\mathbf{e}_{s_k}^\top - \gamma\mathbf{e}_{s_{k+1}}^\top)}_{\mathbf{A}(X_k)}V_{k-1} - \underbrace{\mathbf{e}_{s_k}\mathbf{e}_{s_k,a_k}^\top r}_{\mathbf{b}({X_k})}\bigg)\] It is straightforward to verify that \[\mathbb{E}_{\mu_{\pi}}[\mathbf{A}(X_k)] = \operatorname{diag}(\mu_{\pi})(I - \gamma P_\pi), \quad \mathbb{E}_{\mu_\pi}[\mathbf{b}(X_k)] =\operatorname{diag}(\mu_\pi) r_\pi\;,\] where \(\mu_\pi\) denotes stationary distribution under policy \(\pi\). As long as \(\mu_\pi > 0\), this immediately implies 35 . We note that \(-A\) is Hurwitz due to the stationarity condition \(\mu_\pi^\top P_\pi = \mu_\pi^\top\), which allows the standard linear stochastic approximation theory to be applied directly.

We have considered the on-policy setting. However, similar updates hold in the off-policy regime. Let \(\pi\) be the target policy and \(\pi_b\) the behavior policy. Then, the off-policy TD(0) update is \[\begin{align} \label{eq:importance} V_k(s_k) = {(1-\alpha\rho_k)}V_{k-1}(s_k) + \alpha \rho_k(r_k + \gamma V_{k-1}(s_{k+1})), \quad \rho_k = \frac{\pi(a_k|s_k)}{\pi_b(a_k|s_k)}. \end{align}\tag{36}\] Here, \(\rho_k\) is the importance sampling ratio, which corrects for the distributional shift between the behavior policy \(\pi_b\) and the target policy \(\pi\). It reweights the observed transitions so that the expected update corresponds to the dynamics under \(\pi\), despite sampling from \(\pi_b\). Since the behavior policy distribution \(\mu_{\pi_b}\) does not correspond to the target policy transition kernel \(P_\pi\), the matrix \(A = \mathbb{E}[\mathbf{A}(X_k)] = \operatorname{diag}(\mu_{\pi_b})(I-\gamma P_\pi)\) cannot be guaranteed to be Hurwitz. Consequently, the standard LSA approach cannot be applied directly.

5.1.0.2 Generalizations of TD

Beyond the classical TD(0), there exists a broader family of algorithms denoted by TD(\(k\)), where the update is based on the \(k\)-step return \[V_{t+1}(s_t) = V_t(s_t) + \alpha \big( G_t^{(k)} - V_t(s_t) \big),\] \[G_t^{(k)} = \sum_{\ell=0}^{k-1} \gamma^\ell r_{t+\ell} + \gamma^k V_t(s_{t+k}) .\]

The extreme case \(k \to \infty\) gives the Monte Carlo estimator, which uses the full return until the end of the trajectory. Thus, TD(\(k\)) interpolates between TD(0) (\(k=1\)) and Monte Carlo (\(k=\infty\)), allowing one to balance bias and variance.

A more general approach is TD(\(\lambda\))[45], which combines all \(k\)-step returns using exponentially decaying weights: \[G_t^{(\lambda)} = (1-\lambda) \sum_{k=1}^\infty \lambda^{k-1} G_t^{(k)},\] where \(\lambda \in [0,1]\). This unifies TD(0) (\(\lambda=0\)) and Monte Carlo (\(\lambda=1\)).

5.2 Linear function approximation↩︎

Temporal-difference methods based on gradient descent and linear function approximation [46] form a core part of the modern field of reinforcement learning. We consider the linear function approximation for \(V^{\pi}(s)\), defined for \(s \in \mathcal{S}\), \(\theta \in \mathbb{R}^{d}\), and a feature mapping \(\varphi : \mathcal{S}\to \mathbb{R}^{d}\) as \(V_{\theta}^{}(s) = \varphi^\top(s) \theta\). Our goal is to find a parameter \(\theta^\star\), which defines the best linear approximation to \(V^{\pi}\) (the precise notion of “best” will be specified below). In this context, it is natural to consider \(\mathbb{R}^{S}\) as the vector space of functions \(V:\mathcal{S}\to \mathbb{R}\). The set \(= \{ V_\theta : \theta \in \mathbb{R}^d \}\) is then a \(d\)-dimensional subspace of \(\mathbb{R}^{S}\), spanned by the feature functions \(\{\varphi_i\}_{i=1}^d\). We denote by \(\Phi\in \mathbb{R}^{S\times d}\) the feature matrix with entries \(\Phi(s,i) = \varphi_i(s)\). Using this notation, the value function can be written compactly as \(V_\theta = \Phi\theta\). The TD(0) algorithm with linear function approximation in the on-policy setting updates the parameter \(\theta_k\) as \[\label{eq:td95with95approx} \theta_{k} = \theta_{k-1} - \alpha \varphi(s_k) \delta_k,\quad \delta_k = \varphi(s_k)^\top \theta_{k-1} - r_k -\gamma\varphi(s_{k+1})^\top \theta_{k-1} \;.\tag{37}\] Note that the update rule in Algorithm 8 can be viewed as a special case of 37 when the feature map is chosen to be the one-hot encoding, i.e. \(\varphi(s) = \mathbf{e}_s\). In the general case, consider the linear system \(A \theta^\star = b\) where \[\begin{align} A &= \mathbb{E}_{\mu_\pi} \big[ \varphi(s)\,(\varphi(s) - \gamma \varphi(s'))^\top \big] = \Phi^\top \operatorname{diag}(\mu_{\pi} )(I - \gamma P_\pi) \Phi\;\;,\\ b &= \mathbb{E}_{\mu_\pi} \big[ \varphi(s)\, r(s,\pi(s))\big] = \Phi^\top \operatorname{diag}(\mu_{\pi}) r_\pi\;. \end{align}\] As in the on-policy TD without function approximation, the matrix \(A\) can be shown to be Hurwitz, ensuring that the sequence \(\{\theta_k\}\) converges to \(\theta^\star\), the point at which \(\mathbb{E}[\delta \phi] = 0\), i.e., with zero temporal difference error.

It is well-known that the TD(0) algorithm with linear function approximation may diverge in the off-policy setting. A classical counterexample demonstrating this phenomenon was presented by [47], where TD(0) updates can lead to unbounded growth of the estimated value function, even if the feature space can perfectly represent the true value function. Divergence may be hidden in the matrix \(A = \Phi^\top \operatorname{diag}(\mu_{\pi_b}) (I - \gamma P_\pi) \Phi\), which is not necessarily Hurwitz when the distribution of the behavior policy differs significantly from the target distribution \(\mu_\pi\). To address this, advanced policy evaluation methods such as Generalized Temporal Difference learning (GTD2) and Temporal-Difference learning with Gradient Correction (TDC) were developed in [48], [49]. These algorithms can be interpreted as true gradient methods, designed to optimize the Mean Squared Bellman Projection Error (MSBPE) given by \[\mathbb{E}_{\mu_\pi}[\|V_\theta - \Pi_{\varphi}\mathcal{B}_\pi V_\theta\|^2]\;,\] where \(\Pi_\varphi\) denotes the projection onto the feature subspace. The MSBPE is a natural loss function because the Bellman operator generally maps \(V_\theta\) outside of the feature space. The convergence of these methods can be analyzed within the Two-Time-Scale Stochastic Approximation framework; see [50] for details. Importantly, GTD and TDC remain stable even in off-policy settings, with their iterates converging to the fixed point of the projected Bellman operator \(\Pi_\varphi \mathcal{B}_\pi\).

6 Forward model setting↩︎

6.1 Introduction to Stochastic Multi-Armed Bandits↩︎

The multi-armed bandit (MAB) problem is a fundamental model in sequential decision making, capturing the trade-off between exploration and exploitation. The setting can be informally described as a gambler facing several slot machines (“arms”), each providing stochastic rewards from an unknown distribution, and aiming to maximize the total reward over time.

The concept of the MAB problem was first introduced by [51], who formulated it in the context of the sequential design of experiments. Since then, the problem has served as a canonical benchmark framework for developing and analyzing learning algorithms with provable performance guarantees.

A comprehensive overview of the field is provided in the modern survey by [52], which covers stochastic, adversarial, contextual, and structured variants of the problem. Another influential tutorial is the work of Bubeck and Cesa-Bianchi [53], which presents a unified regret analysis framework for both stochastic and non-stochastic settings. We now provide a formal definition of the problem:

6.1.1 Definitions↩︎

This problem is quite general and have several clear applications.

Let us first summarize the notation. The set of arms (or actions) is denoted by \(\mathcal{A}\), and individual arms are typically represented by \(a\), with or without indices. The underlying reward distribution associated with arm \(a\) is denoted by \(\mathcal{D}_a\). The expected reward (mean) of arm \(a\) is defined as \(\mu(a) \triangleq \mathbb{E}[\mathcal{D}_a]\), and the optimal (best) expected reward is given by \(\mu^\star \triangleq \max_{a \in \mathcal{A}} \mu(a)\). The quantity \(\Delta(a) \triangleq \mu^\star - \mu(a)\) measures how much worse arm \(a\) is compared to the optimal reward \(\mu^\star\); this value is referred to as the gap of arm \(a\). An optimal arm \(a^\star\) is an arm with \(\mu(a^\star) = \mu^\star\) or equivalently, \(\Delta(a^\star) = 0\). Note that the optimal arm need not be unique. Next, we define a performance measure as a cumulative regret (or just regret) at round \(T\)⁴ \[\mathcal{R}(T) \triangleq \sum_{t=1}^T \mu^\star - \mu(a_t) = T \mu^\star - \sum_{t=1}^T \mu(a_t)\;.\] Notice that \(\mathcal{R}(T)\) is a random variable that depends on actions \(a_t\) that may depend on the randomness in obtained rewards and in the algorithm itself. During this chapter, we will focus on the expected regret \(\mathbb{E}[\mathcal{R}(T)]\).

6.1.2 Exploration-Exploitation Tradeoff↩︎

A natural question arises: why do we choose this particular notion of performance, instead of simply considering the gap of the final selected action \(\mu^\star - \mu(a_T)\)? The practical reason is connected to the fact that it is assumed to use the answers of bandit algorithms from the beginning without waiting a lot of time while it converges well.

Another reason is theoretical: considering regret allows us to face the so-called exploration-exploitation dilemma: at each round \(t\), the algorithm faces a trade-off between exploring less frequently chosen arms to gather more information (exploration) and exploiting the arm with the highest estimated mean reward based on current knowledge (exploitation).

In this scenario, it is necessary to ensure that the expected regret grows sublinearly, i.e., \(\mathbb{E}[\mathcal{R}(T)/T] \to 0\). In particular, this implies that \(\mu^\star - \mu(a_T) \to 0\), so by considering regret we can only make our task harder. Next we consider the following naive scheme for regret minimization.

6.1.2.1 Explore-First Algorithm

• Exploration phase: Try each arm \(N\) times;

• Select arm \(\hat{a}\) with the highest average rewards;

• Exploitation phase: play arm \(\hat{a}\) in all remaining \(T - NK\) rounds.

To provide a proof of regret bound, we first recall a fundamental concentration inequality:

Now we are ready establish a regret bound for Explore-First Algorithm.

The presented algorithm has sublinear regret that means that this algorithm actually learn the optimal arm asymptotically but not so efficient as it possible. Next we introduce the main principle that allows us to derive much better regret.

6.1.2.2 Optimism in the Face of Uncertainty (OFU)

The idea of OFU is a key tool for solving the MAB problem. It first appeared in the [54] and was later formalized in the finite-time UCB algorithm by [55]. The main principle is to choose the arm with the highest upper confidence bound on its mean reward. This concept extends naturally to reinforcement learning, where it has been used to design near-optimal algorithms for episodic MDPs.

Let us define \(\overline{\mu}_t(a)\) as an upper confidence bound for arm \(a\), which means that, with high probability we have \(\overline{\mu}_t(a) \geq \mu(a)\). By Hoeffding’s inequality, this upper confidence bound can be chosen as \[\overline{\mu}_t(a) = \widehat{\mu}_t(a) + \beta_t(a) := \frac{1}{n_t(a)} \sum_{t: a_t = a} r_t + \sqrt{\frac{2\log(T)}{n_t(a)}}\;,\] where \(n_t(a)\) is a number of times when the arm \(a\) was picked up to a timestamp \(t\). But importantly, this bound follows not directly from Hoeffding inequality because number of pulls \(n_t(a)\) is random and this randomness should be handled properly. The UCB algorithm is defined as follows: at each round \(t\), select the arm \[a_t = \mathop{\mathrm{arg\ max}}_{a \in \mathcal{A}} \overline{\mu}_t(a)\;,\] observe the reward \(r_t\;\), update the empirical mean \(\widehat{\mu}_t(a_t)\) and the count \(n_t(a_t),\) and then proceed to the next round. Why does this make sense? There are two main reasons to choose arm \(a\) at round \(t\):

• Arm \(a\) has a high mean reward \(\widehat{\mu}_t(a)\), which implies that it is likely to have a high mean reward \(\mu(a)\).

• Arm \(a\) has a large confidence interval \(\beta_t(a)\), which implies that this arms is not explored properly.

In the following we call \(\beta_t(a)\) as an exploration bonus.

6.1.2.3 Thompson Sampling

Thompson Sampling is one of the earliest heuristics developed for multi-armed bandit problems, originally proposed by [56]. For a comprehensive review of modern results, we refer to [57], while a complete Bayesian decision-theoretic foundation can be found in [58]. This randomized algorithm, rooted in Bayesian principles, has gained considerable attention, as empirical studies have demonstrated its strong performance compared to OFU-based methods, see [59].

For simplicity of discussion we consider Thompson Sampling using Beta priors in the Bernoulli bandit problem, i.e., when the rewards are either \(0\) or \(1\), and the probability of reward \(1\) for arm \(i\) (the probability of success) is \(\mu_i\). Let us briefly recall that the Beta distribution is a family of continuous probability distributions supported on the interval \((0, 1)\). The pdf of \(\mathrm{Beta}(\alpha, \beta)\), with parameters \(\alpha > 0\) and \(\beta > 0\), is given by \[f(x; \alpha, \beta) = \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)\,\Gamma(\beta)} \, x^{\alpha - 1} (1 - x)^{\beta - 1}, \quad x \in (0, 1),\] where \(\Gamma(\cdot)\) denotes the gamma function. The mean of \(\mathrm{Beta}(\alpha, \beta)\) is \(\alpha / (\alpha + \beta)\). As can be seen from the pdf, larger values of \(\alpha\) and \(\beta\) result in a tighter concentration of the distribution around its mean. Using Beta priors is useful for Bernoulli rewards because if the prior is a \(\text{Beta}(\alpha, \beta)\) distribution, then after observing a Bernoulli trial, the posterior distribution is simply \(\text{Beta}(\alpha+1,\beta)\) or \(\text{Beta}(\alpha,\beta + 1)\), depending on whether the trial resulted in a success or failure, respectively. TS initially assumes arm \(i\) to have prior \(\text{Beta}(1,1)\) on \(\mu_i\), which is natural because \(\text{Beta}(1,1)\) is the uniform distribution on \((0,1)\). At time \(t\), having observed \(S_i(t)\) successes (reward = 1) and \(F_i(t)\) failures (reward = 0) in \(k_i(t) = S_i(t) + F_i(t)\) plays of arm \(i\), the algorithm updates the distribution on \(\mu_i\) as \(\text{Beta}(S_i(t)+1, F_i(t)+ 1)\). The algorithm then generates independent samples from these posterior distributions of the \(\mu_i\)’s, and plays the arm with the largest sample value.

The first minimax regret bounds for Thompson Sampling were established in [60].

To handle the general case, we slightly modify Algorithm 9. After observing the reward \(r_t\), we perform a Bernoulli trial with success probability \(r_t\) and obtain an outcome \(\tilde{r}_t \in \{0, 1\}\). If \(\tilde{r}_t = 1\), we increment \(S_i \leftarrow S_i + 1\); otherwise, we increment \(F_i \leftarrow F_i + 1\). For further details, see [61]. Later, we will show how to generalize Thompson Sampling to the general case of MDPs.

Besides Thompson Sampling, a closely related idea is Bayes-UCB, which constructs upper confidence bounds using quantiles of the posterior distributions. In [62], the authors analyze the algorithm using \((1 - 1/t)\)-quantiles and show that it achieves asymptotically optimal regret. Empirical results indicate that Bayes-UCB often converges faster than classical UCB methods, benefiting from tighter posterior-based confidence intervals.

6.2 General Problem Settings↩︎

We consider the reinforcement learning problem of an agent interacting with an environment in order to maximize its cumulative rewards through time. We model the environment as a Markov decision process whose transition dynamics are unknown from the agent. For simplicity we assume that the immediate reward function \(r\) is available; however, note that unknown reward functions can be handled without additional difficulties. As the agent interacts with the environment it observes the states, actions and rewards generated by the system dynamics. This leads to a fundamental trade off: should the agent explore poorly-understood states and actions to gain information and improve future performance, or exploit its knowledge to optimize short-run rewards. The most common approach to this learning problem is to separate the process of estimation and optimization. In this paradigm, point estimates of the unknown quantities are used in place of the unknown parameters and a plan is made with respect to these estimates. Naive optimization with respect to these point estimates can lead to premature exploitation and so may never learn the optimal policy. To overcome these difficulties, we will consider two classes of algorithms that simultaneously combine exploration and exploitation and have strong convergence guarantees.

In the following, we focus on tabular finite-horizon MDPs \(\mathcal{M}= (\mathcal{S}, \mathcal{A}, P, R)\). Let us note that both the episodic and non-episodic settings of the problem are considered. In the non-episodic settings, agent starts at state \(s_0\) and selects an action \(a_0\). Next, agent observes a new state \(s_1 \sim P(\cdot | s_0, a_0)\) generated by transition probability kernel and reward \(r_0 \sim R(\cdot | s_0,a_0)\), then agent selects next action \(a_1\) and so on. The interaction with the environment terminates at fixed time horizon \(T\). The agent’s goal is to minimize regret, defined as: \[\mathcal{R}(T) = TV^{\star} - \sum_{t=0}^{T-1}r(s_t,a_t),\] where \(V^{\star}\) is the optimal value of AMDP \(M\), see section 2.

In episodic settings, the agent’s interaction with the environment is divided into episodes, each lasting \(H\) steps. At the end of a round, the interaction with the environment stops, and the agent returns to the initial state \(s_0 = s_0^k\).⁵ . In such settings, the agent minimizes the total regret over \(T = K H\) steps, defined as: \[\mathcal{R}(T) = KV^{\star}(s_0) - \sum_{k=0}^{K-1}\sum_{h=0}^{H-1}r(s_h^k, a_h^k),\] where \((s_h^k, a_h^k)\) is the state and action of the agent in the \(k\)-th episode at step \(h\). Very important feature that here we do not have discounting factor. In the following, we will focus on the episodic case, but note that most algorithms for the episodic model, with minor adjustments, can be transferred to the non-episodic case, and the same regret estimates hold for them, with the only difference being that the horizon length \(H\) is replaced by the diameter of the MDP \(d(\mathcal{M})\):

In applications, it is sometimes useful to consider the case where the transition probabilities of the Markov kernel may change during the episode. So in the following we concentrate on a finite episodic MDPs \(\mathcal{M}= (\mathcal{S}, \mathcal{A}, H, \{P_h\}_{h\in[H]}, \{r_h\}_{h\in[H]})\) , where \(P_h(s' |s, a)\) is the probability transition from state \(s\) to state \(s'\) by taking action a at step \(h\), and \(r_h(s, a)\) is the immediate reward received after taking the action \(a\) in state \(s\) at step \(h.\) As a rule, the analysis of algorithms in the case of a stage-dependent MDP does not differ from the stage-independent MDP case, except for an additional factor \(\sqrt{H}\) in the regret bounds. Note that a stage-dependent MDP \(\mathcal{M}= (\mathcal{S}, \mathcal{A}, P_h, r_h)\) can be equivalently represented as a stage-independent MDP. To this end, define \(\widetilde{\mathcal{M}} = (\widetilde{\mathcal{S}},\mathcal{A}, \widetilde{P}, \widetilde{R} )\), where \(\widetilde{\mathcal{S}} = \bigsqcup_{h=0}^{H-1} \mathcal{S}_h\) denotes the disjoint union of \(H\) copies of \(\mathcal{S}\), which we refer to as levels. The transition kernel \(\widetilde{P}\) is nonzero only between consecutive levels; more precisely, \[\widetilde{P}(s'|s, a) = P_h(s'|s,a), \quad \tilde{r}(s,a) = r(s,a)\quad \text{where}\;s\in\mathcal{S}_h\;,\;s'\in\mathcal{S}_{h+1}\;.\] Note that \(|\widetilde{S}| = HS\) which explains the appearance of an additional factor of \(\sqrt{H}\) in subsequent bounds. The main technical difference in the finite-horizon setting is that the policy depends explicitly on the current time step \(h\in[H]\). This dependence is necessary because, at early stages \(h\), the agent aims to move towards regions with potentially high rewards, while at later stages it prefers to exploit immediate rewards. Consequently, in the finite-horizon setting, we consider a set of policies \((\pi_0, \ldots, \pi_{H-1})\). Slightly abusing notation, we denote this entire sequence simply as \(\pi\). Analogously to the discounted setting, we define the value and action-value functions as follows: \[V^\pi_h(s) = \mathbb{E}\left[ \sum_{t=h}^{H-1} r_t \mid s_{h} = s \right], \qquad Q^\pi_h(s,a) = \mathbb{E}\left[ \sum_{t=h}^{H-1} r_t \mid s_{h} = s, a_h = a \right],\] where \(s_{t+1} \sim P_t(\cdot|s_t,a_t),\;a_t \sim \pi_t(\cdot | s_t), \;r_t =r_t(s_t, a_t)\). Importantly, in this setting value function depends on the current step \(h\). We may define in a similar manner optimal \(V\)- and \(Q\)-functions as a supremum over all possible sets of policies, \[V^\star_h(s) = \sup_{\pi \in \Pi} V^\pi_h(s)\;, \qquad Q^\star_h(s,a) = \sup_{\pi \in \Pi} Q^\pi_h(s,a)\;,\] and, moreover, Bellman and optimal Bellman equations holds for this model \[\begin{align} \label{bel95eq} &Q^\pi_h(s,a) = r_h(s,a) + P_hV^\pi_{h+1}(s,a) , & &Q^\star_h(s,a) = r_h(s,a) + P_hV^\star_{h+1}(s,a)\;, \\ &V^\pi_h(s) = \sum_{a\in\mathcal{A}}\pi_h(a|s)Q^\pi_h(s,a)\;, & &V^\star_h(s) = \max_{a \in \mathcal{A}} Q^\star_h(s, a)\;, \\ &V^\pi_{H}(s) = 0, & &V^\star_{H}(s) = 0\;. \end{align}\tag{40}\] When proving regret bounds, it is often crucial to compare the value function of a fixed policy \(\pi\) when evaluated in two different MDPs (for example, the true environment and its empirical estimate). [The following lemma, commonly referred to as the Simulation Lemma, provides a precise decomposition of this difference and serves as a key analytical tool.

The optimal Bellman equations provide an approcah to computing the optimal \(Q\)-value function via backward dynamic programming. Moreover, the optimal policy \(\pi^\star\) can be represented as the greedy policy with respect to the optimal state-action value function \(Q^\star_h:\) \[\pi_h^\star(s) = \arg \max_{a}Q^\star_h(s,a)\;.\]

6.2.1 Lower Bounds for Episodic Reinforcement Learning↩︎

Understanding the minimax limits of performance in episodic settings is crucial for assessing the optimality of algorithms. Recent results establish sharp lower bounds on the expected regret. In particular, the work [63] shows that for an episodic MDP with horizon \(H\), \(S\) states, and \(A\) actions, any algorithm must incur expected regret at least \[\label{eq:lower95regret95bound} \Omega\!\left(\sqrt{H^2SA T}\right),\tag{41}\] where \(T\) is the total number of interaction steps. The main distinction of this episodic estimate from the finite-horizon estimate in table 2 is that HMDP assumes a generative model assumption while problems in work [63] are solved in an online manner with estimates formed over the trajectories therefore a regret bound here is more natural. The bound 41 matches the upper bounds achieved by several recent algorithms, which we discuss later in this section.

6.2.2 UCRL↩︎

The Upper Confidence Reinforcement Learning (UCRL) algorithm is a foundational optimistic algorithm designed for regret minimization. It was first introduced by [64], who provided near-optimal regret guarantees in the non-episodic average reward setting. Before presenting the algorithm, we slightly abuse notation and define \(V^\pi_h(s\mid P')\) corresponding to policy \(\pi\) evaluated in an MDP with transition kernel \(P',\) i.e. \(\mathcal{M}' = (\mathcal{S}, \mathcal{A}, P', R, H)\). We are now ready to formally define the UCRL algorithm, see Algorithm 10.

Notice that we do not claim that this algorithm is computationally friendly but it requires only solve finite-dimensional convex optimization problem that is doable. Next we are going to define a set of plausible models. First, we define empirical estimate of the model using available data⁶

\[\widehat{P}^k_h(s'|s,a) = \begin{cases} \frac{{N^k_h(s,a, s')}}{N_h^k(s,a)}, & \text{if } N_h^k(s,a) \geq 1 \\ \frac{1}{S}, & \text{otherwise.} \end{cases}\] where \(N_h^k(s, a)\) and \(N_h^k(s,a,s')\) are empirical counts, constructed based on the agent’s observations during the first \(k\) episodes: \[N_h^k(s,a)=\sum_{i=0}^{k-1}\mathbb{I}\{(s_h^i,a_h^i)= (s,a)\}\;,\] \[N_h^k(s,a,s')=\sum_{i=0}^{k-1}\mathbb{I}\{(s_h^i,a_h^i,s_{h+1}^i)= (s,a,s')\}\;.\] where \(a^{k}_h = \pi^{k}_h(s^k_h)\) is an action that was selected at episode \(k\) and step \(h\). In other words, it is a number of visits of state-action pair \((s,a)\) and transitions from state-action pair \((s,a)\) to \(s'\) respectively. To construct the confidence sets, we rely on concentration inequalities for sums of independent random variables. The following lemma will be useful:

Using Lemma 5, we can derive the following result:

Finally, we are ready to state regret bounds for this algorithm.

For a comprehensive survey of UCRL-like algorithms, we refer the reader to [65].

6.2.3 UCBVI↩︎

6.2.3.1 Value Iteration

The UCBVI (Upper Confidence Bound Value Iteration) algorithm updates the policy for the next episode based on information collected from previous episodes. To achieve this, an MDP is constructed with the Markov kernel \(\widehat{P}_h^k\) and rewards defined as \[\widehat{r}_h(s,a) = r_h^k(s,a) + b_h^k(s,a)\;.\] Then the optimal policy \(\pi_h^k\) is found using Bellman equations and dynamic programming. The additional terms added to the true rewards are called bonuses and defined as: \[\label{eq:hoeffding95bonus} b_h^k(s,a) = 2H\sqrt{\frac{\ln{(SAHK/\delta)}}{N_h^k(s,a)}}\;.\tag{43}\] Bonuses allow the algorithm to construct a policy that balances exploration and exploitation, where poorly explored states are assigned higher bonuses. We now present the formal description of the dynamic value iteration algorithm. \[\widehat{V}_H^k(s) = 0, \quad \widehat{Q}_h^k =\min \{ r_h + b_h^k + \widehat{P}_h^k \widehat{V}_{h+1}^{k}, \; H\}\;,\] \[\label{dp95eq952} \widehat{V}_h^k(s) = \max_{a\in\mathcal{A}}\widehat{Q}_h^k(s,a), \; \pi_h^k(s) = \mathop{\mathrm{arg\ max}}_{a\in\mathcal{A}} \widehat{Q}_h^k(s,a)\;.\tag{44}\] We are now ready to formally define the UCBVI algorithm:

UCBVI was first introduced in [66] with two distinct types of exploration bonuses. The bonus specified in equation 43 is referred to as a Hoeffding-type bonus, since the corresponding confidence intervals are derived via the Chernoff–Hoeffding inequality. The following result holds:

The second type of bonus is more intricate and is based on the Bernstein concentration inequality. This more refined bonus design eliminates the extra \(\sqrt{H}\) factor from the regret bound, thereby matching the lower and upper bounds up to logarithmic factors:

We next describe techniques that lead to a regret bound for UCBVI-H that is slightly weaker than the one stated above.

We begin the proof with an important lemma which is proven using a careful application of the Hoeffding inequality:

Note that the bound in lemma 7 could be obtained from lemma 5, but the determinism of the vector \(V^{\star}\) allows us to eliminate the factor \(\sqrt{S}\) from the right-hand side. The above lemmas allow us to define the high-probability event \(\mathcal{E}\), satisfying \(P(\mathcal{E}) \geq 1 - 2\delta\), where inequalities ?? and ?? hold. In the subsequent analysis, we temporarily condition on \(\mathcal{E}\) and prove the key lemma.

We now provide an outline of the proof for the weakened regret bound.

6.2.4 Model-free approaches↩︎

There are two main approaches to RL: model-based and model-free. Model-based algorithms (UCBVI, UCRL) make use of a model for the environment, forming a control policy based on this learned model. Model-free approaches dispense with the model and directly update the value function—the expected reward starting from each state. Advantages of model-free algorithms in reinforcement learning:

• No need to model the environment – these algorithms learn directly from experience without constructing an explicit model of transitions and rewards.

• Flexibility – can be applied in complex and partially observable environments where building a model is difficult.

• Lower computational complexity – do not require storing and updating a transition table like model-based methods.

• Better performance in large state spaces – especially useful for high-dimensional problems.

• Suitable for real-time applications – can adapt to changes in the environment without recomputing the entire model.

Due to these advantages, model-free algorithms became the foundation of deep RL and laid the groundwork for state-of-the-art algorithms such as DQN [67], A2C [68], and TRPO [69]. On the other hand it is believed that model-free algorithms suffer from a higher sample complexity compared to model-based approaches. We draw an analogy with a generative model setting in discounted MDPs. The minimax lower bound on the sample complexity for a generative model is \(\widetilde{{\Omega}}(SA/(1-\gamma)^3\varepsilon^2)\), whereas for Q-learning, a hard example has been constructed with a proven lower bound of \(\widetilde{\Omega}(SA/(1-\gamma)^4\varepsilon^2)\) although a variance reduction technique can close this gap [36]. The simplest and most commonly used algorithm learns the Q-function by updating it at each step, similar to classical Q-learning with a generative model, with the only difference being that the form of the Bellman operator is adapted for a finite-horizon MDP. Exploration of the environment is achieved by allowing the possibility of taking a random action with probability \(\varepsilon\).

Unfortunately, as shown in the example below, an unstructured exploration strategy requires exponentially many episodes to learn.

To design an algorithm that achieves good upper bounds on regret, a model-free algorithm 13 was developed, following the idea of constructing bonuses as in UCBVI. The key difference is that instead of adding bonuses to the rewards, we now add them directly to the Q-function. Using UCB (Upper Confidence Bound) exploration instead of \(\varepsilon\)-greedy exploration in the model-free setting helps manage uncertainties more effectively for different states and actions [71].

The analysis implies that it is important to apply a learning rate of \(\alpha_t = O(H/t)\) rather than \(1/t\) when updating a state-action pair for the \(t\)-th time. The former approach places more emphasis on recent updates, as opposed to treating all past updates equally. This careful adjustment in reweighting is key to the significant difference in sample efficiency, as it leads to a guarantee that is much more efficient compared to earlier methods, which required exponentially many samples in terms of \(H\).

As in optimistic model-based algorithms, the regret analysis follows from the fact that \(Q_h^k\) is an upper bound for \(Q_h^\star\) with high probability. Note that the modification of the algorithm, called Q-Learning-Bernstein, with more accurate bonus modeling leads to an upper bound \(\widetilde{\mathcal{O}}(\sqrt{H^3 SAT})\). The work [72] introduces a novel model-free algorithm, UCB-Advantage, which achieves a regret bound of \(\widetilde{\mathcal{O}}(\sqrt{H^2SAT})\), improving previous results and matching the information-theoretic lower bound up to logarithmic factors.

6.3 Randomized Algorithms↩︎

One of the primary challenges in achieving sample-efficiency in forward-model settings lies in the problem of pessimistic estimation of critical states. Specifically, due to noise in the observed transitions and rewards, an algorithm may assign an unjustifiably low value to the action-value function in certain regions of the state space. As a result, the agent will tend to avoid such states, leading to insufficient visitation and thereby preventing correction of the underestimated values.

To address this issue, early sample-efficient algorithms focused on constructing upper confidence bound on the value function to ensure optimism in the face of uncertainty. These approaches, which we discussed in detail earlier, guarantee that with high probability, the estimated action-values are optimistic, thereby encouraging sufficient exploration.

We now shift focus to an alternative and equally important paradigm: randomized algorithms. The core idea is conceptually simple — rather than explicitly modeling an upper confidence bound on the value function, one perturbs the collected data with appropriately scaled random noise, and then plans using the randomized dataset. This form of data perturbation induces an implicit exploration-exploitation trade-off. Remarkably, this approach not only offers strong theoretical guarantees but also tends to outperform deterministic methods in many practical scenarios.

6.3.1 PSRL and Bayes regret↩︎

Algorithms based on the idea of optimism were deterministic in the following sense: previous observations obtained from the environment uniquely determined the policy for the new round. Now, we will turn to the consideration of another class of so-called randomized algorithms, which trace their roots to the previously discussed Thompson Sampling for MAB. Recall that \(\mathcal{H}_k\) denote the history of observations made prior to episode k. A reinforcement learning algorithm is a deterministic sequence \(\{\mu_k\}\) of functions, each mapping \(\mathcal{H}_k\) to a probability distribution \(\mu_k(\mathcal{H}_k)\) over policies. At the start of the kth episode, the algorithm samples a policy \(\pi_k\) from the distribution \(\mu_k(\mathcal{H}_k)\). The algorithm then execute policy \(\pi^k\) during the kth episode.

We model the agent’s initial uncertainty over the environment through a prior distribution. Posterior sampling for reinforcement learning (PSRL) selects policy through two simple steps. First, a single instance of the environment is sampled from the posterior distribution at the start of an episode. Then, PSRL solves for and executes the policy that is optimal under the sampled environment over the episode. PSRL randomly selects policies according to the probability they are optimal. With this approach exploration is guided by the variance of sampled policies as opposed to optimism.

Let’s list some advantages of posterior sampling relative to optimistic algorithms:

• since PSRL only requires solving for an optimal policy for a single sampled MDP, it is computationally efficient both relative to many optimistic methods, which require simultaneous optimization across a family of plausible environments.

• the presence of an explicit prior allows an agent to incorporate known environment structure in a natural way. This is crucial for most practical applications, as learning without prior knowledge requires exhaustive experimentation in each possible state.

• PSRL is naturally suited to more complex settings where design of an efficiently optimistic algorithm might not be possible.

The first analysis of a randomized algorithm was conducted in the work [73]. In that work, a bound was obtained for the so-called Bayesian regret of the algorithm, which we will define below. Consider an episodic stage-independent MDP model \(M = (S, A, R^M, P^M, H)\), in which the agent only has unknown transition probabilities \(P^M\), but the distribution \(f\) over the set of admissible MDPs \(\Omega\) is assumed to be known. Let the sequence of policies \(\{\pi^0, \ldots, \pi^{K-1}\}\) be obtained during the execution of the algorithm \(\mu\). We define the Bayesian regret of the algorithm \(\mu\) over \(K\) rounds as \[\mathcal{R}(K, \mu) = \sum_{k=0}^{K-1} \Delta_k,\] where \(\Delta_k\) denotes the regret in round \(k\), and is given by \[\Delta_k = V^{\pi^{\star}}_{M^{\star}}(s_0) - {V^{\pi^k}_{M^\star}(s_0)}.\]

Note that the regret is a random variable, which now also depends on the true underlying MDP \(M^{\star}\) and the samples \(\{M_k\}\), which are also random variables (compare with the previous formulations). The probability space consists of the following elements: \[[M^{\star}, M_0, \tau_0, M_1,\tau_1,...,M_{K-1}, \tau_{K-1}]\in \Omega^{K+1}\times\mathcal{T}^K,\] where \(\Omega\) is the set of admissible MDPs, and \(\mathcal{T}\) is the set of trajectories of length \(H\). The probability measure is defined by the density function \[p([M^{\star}, M_0, \tau_0,...,M_{K-1}, \tau_{K-1}]) = f(M^{\star})\prod_{i=0}^{K-1}f(M_i|\mathcal{H}_i)\mathbb{P}(\tau_i|M^{\star}, M_i),\] where \(\mathbb{P}(\tau_i | M^{\star}, M_i)\) is the probability of observing trajectory \(\tau_i\) while using the optimal policy for \(M_i\) within \(M^{\star}\).It is straightforward to show that the sequence of MDPs \(\{M_k\}\) converges almost surely to the true underlying MDP \(M^{\star}\) as the number of observations grows: \[\mathbb{E}[M_k \mid \mathcal{H}_k] \xrightarrow[]{a.s.} M^{\star}.\] However, if we consider the random variables \(M_k\) without conditioning on the history \(\mathcal{H}_k\), then \(M_k\) and \(M^{\star}\) are identically distributed: \[M_k \overset{d}{=} M^{\star}.\] This fact allows us to relate quantities that depend on the true, but unknown, MDP \(M^{\star}\), to those of the sampled MDP \(M_k\), which is fully observed by the agent. The following lemma is an immediate consequence of this observation

Recall that in the definition of the regret, \(\Delta_k = V^{\pi^{\star}}_{M^{\star}}(s_0) - V^{\pi^k}_{M_k}(s_0),\) the optimal policy \(\pi^{\star}\) is not accessible to the agent. This complicates the analysis of the value function \(V^{\pi^{\star}}_{M^{\star}}\) for states encountered during learning. To overcome this difficulty, we introduce an alternative decomposition: \[\widetilde{\Delta}_k = V^{\pi^k}_{M_k}(s_0) - V^{\pi^k}_{M^{\star}}(s_0).\]

The subsequent regret analysis is based on Lemma 10; for further details, we refer the reader to the work [73].

Theorem 8 presents the first regret bound established for the class of randomized algorithms. However, the Bayesian regret analysis only provides guarantees in expectation over a prior distribution on the MDP class, thereby allowing for the existence of hard MDP instances on which the algorithm may exhibit suboptimal learning efficiency. To address this limitation, subsequent works proposed modified algorithms such as SOS-PSRL [74] and OPSRL [75], Bayes-UCBVI [27] which combine optimism-based exploration with posterior sampling. These approaches achieve tight regret bounds even under worst-case MDP instances.

6.3.2 RLSVI↩︎

Randomized Least-Squares Value Iteration (RLSVI) is a model-free reinforcement learning method that promotes deep exploration in episodic MDPs via randomized value function estimates. In the tabular setting, [70] showed that this approach achieves nearly optimal Bayesian regret, outperforming naive exploration strategies such as \(\varepsilon\)-greedy. The idea was extended in [76] to large-scale problems with function approximation, including neural networks. Theoretical analysis in [77] established the first near-optimality bound in the worst-case setting. Subsequent refinements [78] introduced a clipping-based modification that improves worst-case regret bounds.

In this section, we provide a detailed overview of the RLSVI algorithm, highlighting its connection to optimistic approaches and its extension to the non-tabular setting. Conceptually, RLSVI is closely related to UCBVI. At the beginning of each episode \(k\), the agent constructs an empirical model of the environment \(\widehat{\mathcal{M}}_k = ({\mathcal{S}, \mathcal{A}}, \widehat{P}_k, \widehat{r}_k, H)\). Recall that UCBVI adds optimistic bonuses to rewards via \(\bar{r}_k = \widehat{r}_k + b_h^k\). In contrast, RLSVI incorporates randomness by sampling a noise vector \(w^k \in \mathbb{R}^{HSA}\), where each component \(w^k(h, s, a) \sim \mathcal{N}(0, \sigma_k(h,s,a))\). The variance is defined as \(\sigma_k(h,s,a) = (\beta_k / N_h^k(s,a))^{1/2}\), where \(\beta_k\) is a tuning parameter. The random perturbation \(w^k\) should be sufficiently large to dominate the error introduced by performing approximate Bellman updates. A common choice is given by \(\beta_k = \widetilde{\mathcal{O}}(SH^3)\). The agent then computes an optimal policy \(\pi^k\) for the perturbed MDP \(\bar{\mathcal{M}} = ({\mathcal{S}, \mathcal{A}}, \widehat{P}_k, \widehat{r}_k + w^k, H)\).

An alternative view of RLSVI. To clarify the origin of the algorithm’s name, we present an equivalent formulation that generalizes naturally to non-tabular MDPs. Informally, RLSVI constructs a randomized dataset by perturbing the observed transitions and then fits value functions \(Q_h\) via least-squares regression to minimize temporal difference (TD) error. Let the agent’s experience before episode \(k\) be \[D = \sqcup_h D_h, \quad D_h = \{(s_h^i, a_h^i, r_h^i, s_{h+1}^i)\}_{i=0}^{k-1}.\] Perturb each tuple in \(D_h\) by adding independent Gaussian noise \(w_h^i \sim \mathcal{N}(0, \beta_k)\), yielding the randomized dataset \[\widetilde{D} = \sqcup_h \widetilde{D}_h, \quad \widetilde{D}_h = \{(s_h^i, a_h^i, r_h^i + w_h^i, s_{h+1}^i)\}_{i=0}^{k-1}.\] Define the TD-error loss as \[\mathcal{L}(Q \mid Q_{\text{next}}, {D}) = \sum_{(s,a,r,s') \in {D}} \left(Q(s,a) - r - \max_{a'} Q_{\text{next}}(s', a')\right)^2.\] We iteratively solve the regression problem: \[\widehat{Q}_h = \mathop{\mathrm{arg\ min}}_{Q \in \mathbb{R}^{SA}} \mathcal{L}(Q \mid Q_{h+1}, \widetilde{D}_h).\] Each regression step corresponds to one iteration of Value Iteration over \(\widehat{\mathcal{M}}_k\) with injected noise into the Q-values. This view enables generalization to continuous or high-dimensional state spaces.

The key technical contribution in proving this result is the optimism lemma, which ensures that with constant probability, the value of the sampled policy \(\pi^k\) is optimistic relative to the optimal policy \(\pi^{\star}\), conditioned on the agent’s history up to episode \(k\):

Note that the regret bound is suboptimal in its dependence on the horizon \(H\) and state space size \(S\). One challenge identified in [77] is that the range of the Q-function can become unbounded. To address this, [78] propose a simple clipping technique to stabilize the algorithm. After computing \(\widehat{Q}_h\), they clip its values as follows: \[\bar{Q}_{h}(s,a) = \begin{cases} \widehat{Q}_h(s,a) & \text{if } N_h^k(s,a) > \alpha_k \\ H - h + 1 & \text{otherwise} \end{cases}\] where \(\alpha_k = \widetilde{\mathcal{O}}(H^3 S)\). The idea is to retain optimistic estimates until a sufficient amount of data has been collected. This leads to an improved regret bound: \[\mathbb{E}[\mathcal{R}(K)] = \widetilde{\mathcal{O}}( \sqrt{H^4S^2AT}).\]

Finally, the work [79] proposed a “single-seed” version of RLSVI, called the SSR algorithm, which in worst-case settings achieves regret \(\widetilde{\mathcal{O}}(\sqrt{H^2SAT})\), matching the theoretical lower bound. The key idea is to use a single Gaussian sample \(z_k\sim\mathcal{N}(0,1)\) per episode to generate perturbations, instead of sampling independent noises as in standard RLSVI. SSR runs a Value Iteration with the modified rewards \[\widetilde{r}_h(s,a) = r(s,a) + z_k\sigma_k(s,a)\;,\] where \(\sigma_k(s,a)\) corresponds to the Bernstein-type magnitude of the noise. The algorithm also incorporates clipping techniques from [78]. Using different random seeds at each time step, as in RLSVI, can lead to optimism at some steps and pessimism at others, causing perturbation effects to cancel out. SSR avoids this problem by applying the same perturbation throughout an episode, ensuring sufficient exploration without requiring excessively large noise. Numerical simulations show that SSR significantly outperforms RLSVI, consistent with the theoretical regret analysis.

7 Continuous setting↩︎

This chapter moves beyond the tabular setting, focusing on algorithms whose regret guarantees do not scale explicitly with the size of the state space. However, the lower bound in 41 shows that even in finite episodic MDPs, regret must scale at least as \(\sqrt{H^2SA T}\), making such dependence unavoidable in general. This highlights the need for additional structural assumptions to obtain sample-efficient guarantees in large or continuous domains.

Before we describe some of the most influential papers considering the case of continuous state (or state-action) space, let us introduce the following notions.

7.1 Sampling-based fitted value iteration↩︎

Consider the case of an infinite-horizon \(\gamma\)-discounted MDP with a continuous state space represented by a compact subset of some Euclidean space, and a finite action space. Let a generative model of the environment be available (see Section 4). We will assume that rewards are bounded by some positive number \(\hat{R}_{\text{max }}\), w.p. 1 (with probability one).

Recall that Value Iteration (Algorithm 1) is a fundamental algorithm based on applying the backup operator 12 to the current approximation of the optimal value function, \(V_{k+1}:= \mathcal{B}V_k\). In the case of continuous \(\mathcal{S}\), the backup operator has the form \[\label{eq:backup95cont} \mathcal{B}: V \mapsto \mathcal{B}V,\quad (\mathcal{B}V)(s) = \max_{a\in\mathcal{A}} \Bigl( r(s, a) + \gamma \int_{\mathcal{S}} V(s') P(ds'|s,a) \Bigr),\; s\in\mathcal{S}.\tag{45}\] Unfortunately, computation of the integral is often infeasible. One solution is to calculate \(\mathcal{B}V(s)\) approximately using Monte-Carlo integration at a finite number of states \(s\in \mathcal{S}\). After that, one can find a best fit to the computed values in a chosen family of measurable functions \(\mathcal{F}\) on \(\mathcal{S}\). It will be assumed that functions from \(\mathcal{F}\) are bounded by \(V_{\text{max}}< \infty\). This results in fitted value iteration (FVI) [80] presented in Algorithm 16. There, \(M\) represents the number of points in Monte-Carlo integration, \(N\) represents the number of basepoints, i.e., states at which the backup operator is computed approximately, \(\mu \in \mathcal{P}(\mathcal{S})\) is a distribution for sampling basepoints, \(K\) represents the number of iterations, \(p\geq 1\) is a parameter that controls the norm used to compute the best fit in \(\mathcal{F}\).

Before we state the convergence result, let us introduce some notions.

• To measure the approximation power of the function class \(\mathcal{F}\), define its inherent Bellman error as \[d_{p,\mu}(\mathcal{B} \mathcal{F}, \mathcal{F})=\sup _{g \in \mathcal{F}} \inf _{f \in \mathcal{F}}\|f-\mathcal{B} g\|_{p,\mu},\] where \(p\geq 1\), \(\|h\|_{p,\mu}^p:=\int |h(s)|^p \mu(ds)\) for a real-valued measurable function \(h\) defined over \(\mathcal{S}\). This measure captures how closely the function space \(\mathcal{F}\) aligns with the Bellman operator, that is, the dynamics of the MDP.

• Given a deterministic stationary policy \(\pi\), consider a linear operator that maps a distribution \(\nu\) on \(\mathcal{S}\) to the distribution \[\left(\nu P^\pi\right)(d s'):=\int P^\pi(d s' \mid s) \nu(d s),\] where \(P^\pi(d s' \mid s):=P(d s' \mid s, \pi(s) )\). This is the distribution of states if the system is started from \(s_0 \sim \nu\) and policy \(\pi\) was followed for a single time-step. Similarly, \(\nu P^{\pi_1} \ldots P^{\pi_m}\) is the distribution of states if the system is started from \(s_0 \sim \nu\), policy \(\pi_1\) was followed for the first step, policy \(\pi_2\) was followed for the second step, and so on.

• Given distributions \(\nu\) and \(\mu\), and a sequence of stationary policies \(\{\pi_k\}_{k=1}^m\), assume that \(\nu P^{\pi_1} \ldots P^{\pi_m}\) is absolutely continuous w.r.t. \(\mu\), and define the \(m\)-step concentrability of a future-state distribution as \[c(m) := \sup _{\pi_1, \ldots, \pi_m}\left\|\frac{d\left(\nu P^{\pi_1} P^{\pi_2} \ldots P^{\pi_m}\right)}{d \mu}\right\|_{\infty}.\] This number captures how much \(\nu\) can grow in \(m\) steps compared to the reference distribution \(\mu\).

• One way to impose a growth rate condition on \(c(m)\) is by introducing the discounted-average concentrability coefficient of the future-state distributions \[C_{\nu, \mu} := (1-\gamma)^2 \sum_{m \geq 1} m \gamma^{m-1} c(m).\] and assuming that it is finite, \(C_{\nu, \mu} < \infty\). Thanks to discounting, this holds for a reasonably large class of systems, e.g., for uniformly stochastic transitions [81] and in the cases when the top-Lyapunov exponent of the MDP is finite [80].

• If \(s^{1: N} := \left(s_1, \ldots, s_N\right)\) are i.i.d. random variables in \(\mathcal{S}\) with common underlying distribution \(\mu\), then \(\mathbb{E}\left[\mathcal{N}\left(\sigma, \mathcal{F}\left(s^{1: N}\right), \rho_1 \right)\right]\) shall be denoted by \(\tilde{\mathcal{N}}(\sigma, \mathcal{F}, N, \mu)\), where \(\rho_1(x,y)=\|x-y\|_1\), and \(\mathcal{N}\) is the covering number, see Definition 6.

The convergence of Algorithm 16 is characterized as follows. Let \(\delta>0\) be the failure probability, \(\varepsilon>0\) be a positive number. Fix the distributions \(\nu, \mu\), and the function \(V_0 \in \mathcal{F}\). For a sufficient number of iterations \(K\) linear in \(\log (1 / \varepsilon)\), \(\log V_{\text{max}}\) and \(\log (1 /(1-\gamma))\), and numbers of basepoints \(N\) and integration points \(M\) polynomial in \(1 / \varepsilon\), \(\log (1 / \delta)\), \(\log (1 /(1-\gamma))\), \(V_{\text{max}}\), \(\hat{R}_{\text{max }}\), \(\log (|\mathcal{A}|)\), \(\log \tilde{\mathcal{N}}\bigl(\frac{c \varepsilon(1-\gamma)^2}{C_{\nu, \mu}^{1 / p} \gamma}, \mathcal{F}, N, \mu\bigr)\) for some constant \(c>0\), it holds w.p. at least \(1-\delta\), \[\left\|V^*-V^{\pi_K}\right\|_{p, \nu} \leq \frac{2 \gamma}{(1-\gamma)^2} C_{\nu, \mu}^{1 / p}\, d_{p, \mu}(\mathcal{B} \mathcal{F}, \mathcal{F})+\varepsilon,\] where \(\pi_K\) is a policy greedy w.r.t. the \(K\)-th iterate, i.e., \[\pi_K(s) = \arg\max_a \Bigl( r(s, a) + \gamma \int_{\mathcal{S}} V_K(s') P(ds'|s,a) \Bigr)\quad \forall s\in\mathcal{S}.\]

7.2 Kernel UCBVI↩︎

In [82], the authors proposed an approach to continuous RL with a generative model based on smoothing kernels, and provided asymptotic convergence guarantees. When a generative model of the environment is unavailable, agent interacts with the environment and thereby collects data about the system dynamics, like in Section 6. The interactions come in form of \(K\) episodes of length \(H\) (horizon). Consider the case when \(\mathcal{S}\) and \(\mathcal{A}\) are endowed with metrics \(\rho_{\mathcal{S}}\) and \(\rho_{\mathcal{A}}\), respectively, and let the metric \(\rho\) on \(\mathcal{S}\times \mathcal{A}\) be given by their sum. The work [83] proposes Kernel-UCBVI, an algorithm that can be seen as a version of the UCBVI (Algorithm 11) for metric state-action space that uses kernel smoothing to estimate the rewards and transition. To proceed, we make the following assumption:

Let \(u, v \in \mathcal{S}\times \mathcal{A}\). For some function \(g: \mathbb{R}_{+} \rightarrow[0,1]\), define the kernel function as \[\psi_\sigma(u, v) := g(\rho[u, v] / \sigma)\] where \(\sigma\) is the bandwidth parameter (larger bandwidths introduce more smoothing). Moreover, let \((s_h^i, a_h^i, s_{h+1}^i, r_h^i)\) be the state, the action, the next state and the reward at step \(h\) of episode \(i\). The kernel function allows us to define weights \[w_h^i(s, a) \stackrel{\text{ def }}{=} \psi_\sigma\left((s, a),\left(s_h^i, a_h^i\right)\right),\quad (i, h) \in[K] \times[H]\] that determine how much a past transition sample \((s_h^i, a_h^i, s_{h+1}^i, r_h^i)\) should contribute to our estimate of the reward and transitions at a query point \((s, a)\). Specifically, fix an episode \(k\) and define normalized weights as \[\widetilde{w}_h^i(s, a) \stackrel{\text{ def }}{=} \frac{w_h^i(s, a)}{\mathbf{C}_h^k(s, a)}\] where \(\mathbf{C}_h^k(s, a) \stackrel{\text{ def }}{=} \beta+\sum_{i=1}^{k-1} w_h^i(s, a)\) are generalized counts that act like a proxy for the number of visits to \((s, a)\), and \(\beta>0\) is a regularization term. An estimate of the rewards and transitions for each state-action pair can now be defined as \[\begin{align} & \widehat{r}_h^k(s, a) \stackrel{\text{ def }}{=} \sum_{i=1}^{k-1} \widetilde{w}_h^i(s, a) r_h^i, \\ & \widehat{P}_h^k(s' \mid s, a) \stackrel{\text{ def }}{=} \sum_{i=1}^{k-1} \widetilde{w}_h^i(s, a) \delta_{s_{h+1}^i}(s'). \end{align}\] Here, \(\delta_s\) denotes the Dirac measure with mass at \(s\). Like the usual UCBVI, Kernel-UCBVI computes an optimistic Q-function \(\widetilde{Q}_h^k\) through value iteration, a.k.a. backward induction: \[\widetilde{Q}_h^k(s, a)=\widehat{r}_h^k(s, a)+\widehat{P}_h^k V_{h+1}^k(s, a)+\mathrm{B}_h^k(s, a),\] where \(V_{H+1}^k(x)=0\) for all \(x \in \mathcal{S}\) and \(\mathrm{B}_h^k(s, a)\) is an exploration bonus that equals, up to constants and logarithmic terms, \[\frac{H}{\sqrt{\mathbf{C}_h^k(s, a)}}+\frac{\beta H}{\mathbf{C}_h^k(s, a)}+L_1 \sigma.\] Here, \(L_1 \stackrel{\text{ def }}{=} \sum_{h^{\prime}=1}^H \lambda_r \lambda_p^{H-h^{\prime}}\), which is shown to be a Lipschitz constant of the optimal Q-function \(Q^*_1(s, a)\) [83]. More generally, \(L_h \stackrel{\text{ def }}{=} \sum_{h^{\prime}=h}^H \lambda_r \lambda_p^{H-h^{\prime}}\) is a Lipschitz constant of \(Q^*_h(s, a)\).

We denote by \(\mathcal{D}_h=\left\{\left(s_h^i, a_h^i, s_{h+1}^i, r_h^i\right)\right\}_{i \in[k-1]}\) for \(h \in[H]\) the samples collected at step \(h\) before episode \(k\), and give the pseudocode of Kernel-UCBVI and its subroutine in Algorithms 17 and 18.

7.3 Q-learning with adaptive discretization↩︎

When dealing with a continuous state-action space, a straightforward idea is to discretize it and apply known algorithms. However, complexity of naive discretization grows exponentially with the dimension, therefore, a more sophisticated approach is required. The work [84] introduces an approach based on data-driven adaptive discretization, meaning that a partition of regions of the state-action space is refined if they are frequently visited in historical trajectories, and have higher payoff estimates. The approach is based on Q-learning and is thus model-free. Similarly to the previous subsection, it is assumed that \(\mathcal{S}\times \mathcal{A}\) is a subset of a metric space with a metric \(\rho\), and that the optimal Q-function is Lipschitz continuous with respect to the metric. Additionally, \(\mathcal{S}\times \mathcal{A}\) is assumed to be bounded. Moreover, the algorithm requires a covering oracle which takes a finite collection of balls and a set \(X\) and either declares that they cover \(X\) or outputs an uncovered point. An alternative assumption is to assume the covering oracle is able to take a set \(X\) and value \(r\) and output an \(r\)-net of \(X\), i.e., a set of points \(\mathcal{G}\) such that the distance between any two distinct points in \(\mathcal{G}\) is at least \(r\) and \(\cup_{x \in \mathcal{G}} B(x, r) \supseteq X\).

The algorithms begins step \(h\) of episode \(k\) with a collection of balls \(\mathcal{G}^k_h\) covering \(\mathcal{S}\times \mathcal{A}\) and an upper confidence value \(Q^k_h(B)\) for each ball \(B \in \mathcal{G}^k_h\). Domain of a ball \(B\) is defined as a subset of \(B\) which excludes all balls in \(\mathcal{G}^k_h\) having a strictly smaller radius. A ball \(B\) is called relevant for a point \(x \in S\) if \((x,a) \in \operatorname{dom}(B)\) for some \(a \in A\). The algorithm proceeds with the following three steps:

1. select a ball \(B \in \mathcal{G}^k_h\) relevant for the current state \(s^k_h\) which has maximum \(Q^k_h(B)\);

2. update the upper confidence value;

3. if the number of times \(B\) or its ancestors have been selected at step \(h\) in previous episodes is large enough, re-partition the space by splitting \(B\).

The expression for updating \(Q^k_h(B)\) and other detail can be found in the original work [84]. Over \(K\) episodes of length \(H\), the algorithm achieves a regret bound \[\mathcal{R}(K)=\widetilde{\mathcal{O}}\left(H^{5 / 2} K^{(d+1) /(d+2)}\right)\] with probability at least \(1-\delta\), where \(d\) is the covering dimension of \(\mathcal{S}\times\mathcal{A}\) (see Definition 6).

7.4 Linear Function Approximation↩︎

An important approach to handling infinite state (or state-action) spaces is based on function approximation. Perhaps the most studied setting is the case of a linear MDP. One of the ways to define linear MDP [85] is as follows.

Recall that (action) value iteration in tabular case (see subsection 3.2) performs Bellman updates of the form \[Q_h(s, a) \leftarrow r_h(s,a) + \sum_{s'\in\mathcal{S}}P_h(s'|s,a) \max_{a'\in\mathcal{A}} Q_{h+1}(s', a')\quad \forall (s,a) \in \mathcal{S}\times \mathcal{A}.\] These updates are infeasible in the case of infinite state-action space and unknown transitions. A key property of the linear MDP is that, for all policies, the action-value functions are linear in the feature map \(\boldsymbol{\phi}\) [85]. This allows us to parametrize \(Q_h^{\star}(s, a)\) by a linear form \(\mathbf{w}_h^{\top} \boldsymbol{\phi}(s, a)\) and replace the Bellman update by solving for \(\mathbf{w}_h\) in the following (regularized) least-squares problem: \[\mathbf{w}_h \leftarrow \underset{\mathbf{w} \in \mathbb{R}^d}{\operatorname{argmin}} \sum_{i=1}^{k-1}\left[r_h\left(s_h^i, a_h^i\right)+\max _{a \in \mathcal{A}} Q_{h+1}\left(s_{h+1}^i, a\right)-\mathbf{w}^{\top} \boldsymbol{\phi}\left(s_h^i, a_h^i\right)\right]^2+\lambda\|\mathbf{w}\|^2 .\] This results in the algorithm known as Least-Square Value Iteration (LSVI) [70], [86]. Algorithm 19 implements an optimistic modification of LSVI by adding an Upper Confidence Bound (UCB) bonus term of form \(\beta\left(\boldsymbol{\phi}^{\top} \Lambda_h^{-1} \boldsymbol{\phi}\right)^{1 / 2}\) to encourage exploration, where \(\Lambda_h\) is the Gram matrix of the regularized least-squares problem, and \(\beta\) is a scalar.

A regret bound for the algorithm is given by the following theorem, where \(T = KH\) is the total number of steps.

The theorem can also be generalized to the case when the MDP is approximately linear [85].

8 Policy Gradient-Based Reinforcement Learning↩︎

Policy gradient (PG) methods optimize a differentiable, parameterized policy \(\pi_\theta(a\mid s)\) directly by ascending an estimate of the gradient of the expected return. The direct policy-search viewpoint is particularly convenient when policies must be differentiable in parameters (e.g., softmax policies) and when policies are represented by flexible function approximators. In a finite-horizon MDP \((\mathcal{S},\mathcal{A},P,r,\gamma,H)\), define a trajectory \(\tau=(s_0,a_0,\dots,s_H)\) generated by \(s_0\sim \rho\), \(a_t\sim\pi_\theta(\cdot\mid s_t)\), and \(s_{t+1}\sim P(\cdot\mid s_t,a_t)\). The canonical objective is \[\label{eq:pg95objective} J(\theta):=\mathbb{E}_{\tau\sim\pi_\theta}\left[\sum_{t=0}^{H-1}\gamma^t r(s_t,a_t)\right].\tag{46}\] Under mild regularity conditions, the policy gradient theorem yields \[\label{eq:policy95gradient95theorem} \nabla_\theta J(\theta) = \mathbb{E}_{\tau\sim\pi_\theta}\left[\sum_{t=0}^{H-1}\gamma^t\nabla_\theta\log \pi_\theta(a_t\mid s_t)\,R_t^{\pi_\theta}(s_t,a_t)\right],\tag{47}\] where \(R_t^{\pi}(s,a)=\mathbb{E}_{\pi}\!\left[\sum_{k=t}^{H-1}\gamma^{k-t}r(s_k,a_k)\mid s_t=s,a_t=a\right]\) provides a credit assignment signal for each decision [87].

The REINFORCE estimator replaces \(R_t^{\pi_\theta}(s_t,a_t)\) with a Monte Carlo return \(G_t:=\sum_{k=t}^{H-1}\gamma^{k-t}r(s_k,a_k)\) and optionally subtracts a baseline \(b(s_t)\) (which leaves the estimator unbiased) to reduce variance: \[\label{eq:reinforce95estimator} \widehat{\nabla_\theta J(\theta)} = \sum_{t=0}^{H-1}\gamma^t\nabla_\theta\log \pi_\theta(a_t\mid s_t)\,\big(G_t-b(s_t)\big).\tag{48}\] Stochastic gradient ascent updates \(\theta\leftarrow\theta+\alpha\,\widehat{\nabla_\theta J(\theta)}\) [88].

Actor–critic methods maintain an actor (the policy \(\pi_\theta\)) and a critic (a value-function approximation) and use the critic to construct low-variance advantage estimates. Let \(V_w(s)\approx V^{\pi_\theta}(s)\), a neural network based approximation of \(V^{\pi_\theta}\), and define the temporal-difference (TD) residual \[\label{eq:td95error} \delta_t:=r(s_t,a_t)+\gamma V_w(s_{t+1})-V_w(s_t).\tag{49}\] A one-step advantage estimate is \(\hat{A}_t\approx \delta_t\), and the actor update takes the generic form \[\label{eq:actor95update95adv} \widehat{\nabla_\theta J(\theta)}\approx \sum_{t=0}^{H-1}\nabla_\theta\log\pi_\theta(a_t\mid s_t)\,\hat{A}_t.\tag{50}\] Two-time-scale analyses formalize convergence for classes of actor–critic algorithms under linear function approximation and suitable step sizes [89]. Natural-gradient variants (natural actor–critic) further precondition updates using the Fisher information, improving invariance properties and empirical stability [90].

Generalized Advantage Estimation (GAE) defines a family of advantage estimators that interpolates between high-variance Monte Carlo returns and low-variance but biased TD estimates. With TD residuals 49 and parameter \(\lambda\in[0,1]\), the (truncated) GAE estimator is \[\label{eq:gae} \hat{A}_t^{\mathrm{GAE}(\gamma,\lambda)}:=\sum_{l=0}^{H-1-t}(\gamma\lambda)^l\,\delta_{t+l}.\tag{51}\] Hence \(\lambda\to 0\) approaches one-step TD, while \(\lambda\to 1\) approaches Monte Carlo-style returns [91]. The whole framework allows us to consider \(H\rightarrow\infty\) whenever \(\gamma\in(0, 1)\), but in practice we usually work with finite \(H\) in policy gradient based reinforcement learning, sometimes taking \(\gamma=1\).

Advantage Actor–Critic (A2C) and Asynchronous Advantage Actor–Critic (A3C) implementations typically optimize a combined objective consisting of an advantage-weighted policy loss using \(\hat{A}_t\) (often GAE), a squared-error value loss for the critic \(V_w\), and an entropy regularizer \(\mathcal{H}(\pi_\theta(\cdot\mid s_t))\) to encourage exploration. Asynchronous variants (A3C) parallelize data collection and gradient updates over multiple workers, trading off bias due to stale parameters against throughput gains [91].

Large policy updates can be unstable because the sampling distribution induced by \(\pi_\theta\) changes with \(\theta\). Trust-region policy optimization stabilizes learning by solving a surrogate optimization that limits the average change between the previous policy \(\pi_{\theta_{\mathrm{old}}}\) and the updated policy \(\pi_\theta\), commonly measured by KL divergence. Using an advantage estimate \(\hat{A}_t\), define the importance ratio \[\label{eq:is95ratio} \kappa_t(\theta):=\frac{\pi_\theta(a_t\mid s_t)}{\pi_{\theta_{\mathrm{old}}}(a_t\mid s_t)}.\tag{52}\] A trust-region update can be written as \[\label{eq:trpo} \max_{\theta}\;\mathbb{E}\!\left[\kappa_t(\theta)\hat{A}_t\right] \quad\text{s.t.}\quad \mathbb{E}\!\left[D_{\mathrm{KL}}\!\left(\pi_{\theta_{\mathrm{old}}}(\cdot\mid s_t)\,\|\,\pi_{\theta}(\cdot\mid s_t)\right)\right]\le \delta,\tag{53}\] where expectations are over samples collected under \(\pi_{\theta_{\mathrm{old}}}\) [92]. The method with steps of form 53 is also called Natural policy gradient (NPG), it defines steps in the information geometry induced by the Fisher information matrix of \(\pi_\theta\), which is closely related to KL–divergence–based trust regions and mirror descent in policy space. The work [93] constructs a unifying mathematical viewpoint is formalized by the policy mirror descent (PMD) framework while establishing sample complexity for tabular discounted MDPs under generative model setting similar to model–free Q-learning 30 .

Proximal Policy Optimization (PPO) replaces the constrained problem 53 with an unconstrained objective that heuristically prevents excessively large policy updates. The widely used clipped surrogate objective is \[\label{eq:ppo95clip} \begin{align} &L^{\mathrm{CLIP}}(\theta) := \mathbb{E}\!\left[\min\!\Big(\kappa_t(\theta)\hat{A}_t,\;\mathrm{clip}\big(\kappa_t(\theta),1-\varepsilon,1+\varepsilon\big)\hat{A}_t\Big)\right],\\ &\mathrm{clip}\big(\kappa_t(\theta),1-\varepsilon,1+\varepsilon\big) := \min\left\{\max\left\{\kappa_t(\theta), 1 - \varepsilon\right\}, 1 + \varepsilon\right\}, \end{align}\tag{54}\] often augmented with value-function and entropy terms as in actor–critic [94], [95]. PPO is one of the most used RL algorithm in general, also it is the principal RL algorithm in some specific Natural Language Processing tasks, such as reinforcement learning from human feedback (RLHF) for model alignment due to its simplicity and reliability [96]. PPO operates in the actor–critic framework:

• The actor represents the current policy \(\pi_\theta\) and produces an action \(a_t\) based on the current state \(s_t\) (the input).

• The critic (or value model) estimates the long–term expected return of the policy in the current state. This makes it possible to judge whether the actual reward for the actor’s action was better or worse than expected, and to update the policy accordingly.

Without a critic we cannot know whether a reward of, say, \(10\) is good or bad. Ten points may be an excellent reward if in that state we usually obtain only zero or one points, or very poor if we usually obtain \(100\) or \(500\) points. The critic provides a baseline by predicting the expected value of each state; we then compare the received reward with this baseline to determine the advantage and to update the actor effectively. In PPO the training proceeds schematically as follows:

1. All experienced trajectories—tuples of the form \((s_t, a_t, r_t, s_{t+1})\)—are stored as lists.

2. The environment accepts actions from the actor and returns the reward \(r_{t}\) and the next state \(s_{t+1}\).

3. A loss function is defined for updating the critic’s parameters.

4. We update the parameters of both the actor and the critic using standard gradient descent.

5. The actor’s loss employs a clipping coefficient that prevents the new policy \(\pi_{\theta}\) from deviating too far from the old policy \(\pi_{\theta_{\text{old}}}\). This "clip" coefficient is a distinctive feature of PPO.

6. The advantage — which may be positive or negative — is computed by comparing the actual reward with the critic’s predicted value for the state.

Unlike vanilla PG, A2C and A3C, which are purely on–policy methods, TRPO and PPO at each step estimate average of \(\kappa_t(\theta)\) and \(\hat{A}_t\) based on trajectories obtained using \(\pi_{\theta_{\mathrm{old}}}\). This importance–weighting–based estimate makes these mehthods behave similar to off-policy methods such as \(Q\)-learning [97].

Despite being simple and practically reliable, PPO method has its drawbacks: the critic model is expensive to train and resource–intensive, which slows down learning. To mitigate this, the Group Relative Policy Optimization (GRPO) algorithm was proposed, a more efficient variant of PPO that quickly gained popularity [98], [99]. Its key ideas are:

1. The model generates not just a single answer but several answers to the same question.

2. GRPO eliminates the value model entirely; there is no critic. To evaluate rewards correctly we use the average reward of the group of answers to the same question to determine how good each action is.

3. Each answer in the group receives its own reward. Instead of comparing the reward of a particular action to the critic’s expectation, we compare the reward to the group’s average reward.

4. The policy is updated on the basis of the relative advantages within the group.

In language-model post-training and other sequence-level settings, group-relative variants of on-policy PG estimate advantages by comparing multiple rollouts generated for the same input forming a group relative policy optimization via group-relative baselines. Given an input \(x\), sample \(K\) outputs \(\{y^{(k)}\}_{k=1}^K\) and compute rewards \(R(x,y^{(k)})\). A simple group-relative baseline is the within-group mean [100]: \[\label{eq:dr95grpo95adv} \begin{align} \hat{A}^{\mathrm{Dr.GRPO}}(x,y^{(k)}) &:=R(x,y^{(k)})-\frac{1}{K}\sum_{j=1}^K R(x,y^{(j)}), \end{align}\tag{55}\] possibly followed by normalization [98]: \[\label{eq:grpo95adv} \begin{align} \hat{A}^{\mathrm{GRPO}}(x,y^{(k)}) &:=\frac{1}{\sigma(x, \{y^{(j)}\}_{j = 1}^{K}\})}\left(R(x,y^{(k)})-\frac{1}{K}\sum_{j=1}^K R(x,y^{(j)})\right),\\ \sigma(x, \{y^{(j)}\}_{j = 1}^{K}\}) &:=\sqrt{\frac{1}{K - 1}\sum\limits_{k = 1}^{K}\left(R(x, y^{(k)}) - \frac{1}{K}\sum\limits_{j = 1}^{K}R(x, y^{(j)})\right)^{2}}. \end{align}\tag{56}\] Despite simpler form, \(\hat{A}^{\mathrm{Dr.GRPO}}(x,y^{(k)})\) was introduced after \(\hat{A}^{\mathrm{GRPO}}(x,y^{(k)})\), Dr.GRPO stands for GRPO done right as it was introduced to tackle an optimization bias in GRPO training procedure and to combat an artificial increase in response length during training as well. The initial advantage estimate \(\hat{A}^{\mathrm{GRPO}}\) has a bias introduced by a standard normalization term, \(\hat{A}^{\mathrm{Dr.GRPO}}(x,y^{(k)})\) removes this scaling to make a less biased advantage estimate, which can speed up the optimization procedure [100]. Policy updates then mirror 54 with \(\hat{A}_t\) replaced by \(\hat{A}^{\mathrm{GRPO}}\) or \(\hat{A}^{\mathrm{Dr.GRPO}}\), yielding a critic-free on-policy optimization rule based on relative comparisons [101], [102].

Another approach is based on preference optimization, it starts from pairwise comparisons of outputs. Given tuples \((x,y^+,y^-)\) from preference dataset \(\mathcal{P}=\{(x,y^+,y^-)\}\), where \(y^+\) is preferred over \(y^-\) for input \(x\), a common probabilistic model assumes \[\label{eq:bt95preference} \begin{align} \mathbb{P}(y^+\succ y^-\mid x)&=\sigma\!\big(r_\phi(x,y^+)-r_\phi(x,y^-)\big),\\ \sigma(z) &:= \frac{1}{1 + \exp(-z)}, \end{align}\tag{57}\] which links preference learning to reward modeling via a latent reward function \(r_\phi\) [103], [104]. Direct Preference Optimization (DPO) eliminates explicit on-policy rollouts by optimizing the policy parameters directly against preference pairs relative to a fixed reference policy \(\pi_{\mathrm{ref}}\) [105]. DPO was proposed as a simpler and more efficient alternative to PPO. Its key distinction is that DPO does not require a separate reward model or an explicit reinforcement–learning phase. Instead, it directly optimises the policy using pairs of answers ranked by human preference:

1. Unlike classical RLHF approaches, DPO completely avoids training a reward model and the associated RL step. Instead it directly optimises the model (policy) on a small set of human preference data specified via pairs of answers.

2. For each prompt in the dataset there are two candidate answers from the model, one of which humans deemed more preferable.

3. These answer pairs are used directly to optimise the model. The objective is to maximise the probability of the preferred answers and minimise the probability of less preferred answers. In other words, the model learns directly to choose answers that best match human preferences without explicitly constructing a reward model or computing advantages through a critic.

DPO representative objective has the form \[\label{eq:dpo95objective} \max_{\theta}\;\mathbb{E}_{(x,y^+,y^-)\sim\mathcal{P}}\left[\log \sigma\!\Big(\beta\,\Delta_{\theta,\mathrm{ref}}(x,y^+,y^-)\Big)\right],\tag{58}\] with \(\beta>0\) and \[\Delta_{\theta,\mathrm{ref}}(x,y^+,y^-) = \log\frac{\pi_\theta(y^+\mid x)}{\pi_\theta(y^-\mid x)} - \log\frac{\pi_{\mathrm{ref}}(y^+\mid x)}{\pi_{\mathrm{ref}}(y^-\mid x)}.\] Thus, alignment can be cast as a stable, offline (supervised) optimization procedure that retains an implicit reward interpretation, and admits reward-guided extensions [104], [106].

9 Reinforcement Learning applications in Natural Language Processing↩︎

Reinforcement learning (RL) contributes to natural language processing (NLP) both as a framework for sequence-level optimization of generative models and a specific fine-tuning mechanism — post-training for aligning large language models (LLMs) with human preferences, safety constraints, and interactive task success [107]. Classical fine‑tuning does an excellent job of adapting a model to specific scenarios. However, it has an important limitation: fine‑tuning teaches the model to reproduce correct answers from a pre‑prepared dataset, but it does not always give the model an understanding of which answers are genuinely good, useful and appropriate from a human point of view. Therefore, in industry, after the fine‑tuning stage the next player is usually RL. Unlike simple fine‑tuning, RL does not merely show the model examples of good answers but directly provides feedback on how useful or correct its answer was. Thanks to this feedback, the model gradually learns which generations are most appropriate and why.

In the context of generative model training (eg. LLMs) the learning from a reward signal combines two main entities – agent and environment as follows:

• agent is the large language model. It takes actions, for example generating a response to a user question.

• the environment reacts to the agent’s action and produces a reward. For LLMs the environment is any external evaluator; most often this is a special model that imitates human evaluations (commonly called the reward model).

After receiving feedback from the environment the agent adjusts its behaviour — this behavioural function is called the policy — so that it receives greater rewards in the future.

The classical learning framework utilizes sequence generation as an episodic decision process. For an input \(x\sim\mathcal{D}\) (e.g., a document to summarize, a dialogue history, or a user prompt), conditional generation can be modeled as an episodic MDP whose state at time \(t\) is the prefix \(s_t=(x,y_{<t})\), action \(a_t\in\mathcal{A}\) is the next token, and the transition is deterministic \(s_{t+1}=(x,y_{\le t})\). A policy \(\pi_\theta(a_t\mid s_t)\) induces a distribution over complete sequences \(y=y_{1:T}\) and the objective is to maximize expected sequence-level utility \[\label{eq:nlp95seq95rl95objective} J_{\mathrm{NLP}}(\theta):=\mathbb{E}_{x\sim\mathcal{D}}\;\mathbb{E}_{y\sim\pi_\theta(\cdot\mid x)}\!\big[R(x,y)\big],\tag{59}\] where \(R(x,y)\) may encode task-specific metrics, human preference scores, constraint violations, or automatically verifiable rewards. This formulation enables optimization of non-differentiable sequence-level criteria using Monte Carlo policy gradients and actor–critic methods [87], [88].

In more modern LLM post-training, Reinforcement Learning from Human Feedback (RLHF) is used. This name means that in the case of LLMs the reward is mainly defined through human evaluations. Unlike classical RL, where a model typically receives numerical rewards from a predefined reward function, in RLHF the model learns to generate answers directly from human preferences. RLHF is used to make the language model better conform to human expectations and preferences; this process is known as alignment — that is, aligning the model’s behaviour with human values and desires. In practice RLHF is implemented as a multi–stage pipeline:

1. People compare several candidate answers produced by the model for the same question — this may take the form of pairwise comparisons or ranking. It is also possible to collect numerical or binary ratings.

2. The human–annotated data are used to train a reward model. The reward model learns to reproduce human evaluations; in other words it automatically predicts which answers people prefer over others.

3. The reward model then serves as the environment for further training of the primary language model. The LLM generates answers and receives numerical rewards from the reward model, based on which it gradually improves its responses. This stage is usually performed via policy gradient-based methods.

4. The model’s best answers can be sent back to humans for additional evaluation and validation. By iterating this feedback cycle one can collect more data and improve the reward model.

RLHF combines preference learning with RL-style policy optimization. A typical pipeline starts from a preference dataset \(\mathcal{P}=\{(x,y^+,y^-)\}\) and fits a reward model \(r_\phi(x,y)\) by maximum likelihood under a Bradley–Terry–Luce (BTL) or logistic assumption [103]: \[\label{eq:rm95loss} \min_{\phi}\;\mathbb{E}_{(x,y^+,y^-)\sim\mathcal{P}} \left[ -\log \sigma\!\big(r_\phi(x,y^+)-r_\phi(x,y^-)\big) \right].\tag{60}\] The subsequent policy optimization stage can be written as a KL-regularized RL problem \[\label{eq:rlhf95kl} \max_{\theta}\;\mathbb{E}_{x\sim\mathcal{D}}\;\mathbb{E}_{y\sim\pi_\theta(\cdot\mid x)} \left[ r_\phi(x,y)-\beta\,D_{\mathrm{KL}}\!\left(\pi_\theta(\cdot\mid x)\,\|\,\pi_{\mathrm{ref}}(\cdot\mid x)\right) \right],\tag{61}\] where \(\pi_{\mathrm{ref}}\) is a reference (often supervised fine-tuned) policy and \(\beta>0\) controls the strength of regularization. While 61 is frequently introduced in a one-step abstraction by treating an entire response \(y\) as a single action, it is compatible with the token-level episodic view by distributing rewards and KL penalties across timesteps. In practice, on-policy surrogates such as PPO or group-relative variants such as GRPO are used to optimize 61 , while direct preference methods (e.g., DPO) bypass explicit on-policy updates by solving a closely related offline objective on preference pairs [2], [101], [102], [106]. RLHF is intended to make language models more "human–like": models trained with human feedback can conduct polite dialogue or show empathy because they learn not merely to predict text but to account for human evaluations and preferences. RLHF is also indispensable for safety and appropriateness, helping to minimize undesirable or harmful outputs.

For more complex tasks including reasoning an Agentic RL approach entangled with tool-augmented environments is usually employed. Beyond single-turn alignment, agentic RL treats an LLM as a policy that acts over multiple steps in a (partially observable) environment. The agent receives observations \(o_t\) (dialogue context, retrieved documents, tool outputs), chooses structured actions \(u_t\) (token emissions, tool calls, API actions), and the environment transitions according to a kernel \(P(o_{t+1}\mid o_{\le t},u_{\le t})\). With task reward \(r_t\) (e.g., successful completion, verifiable correctness, cost or latency penalties), the objective becomes \[\label{eq:agentic95objective} J_{\mathrm{agent}}(\theta):=\mathbb{E}_{\pi_\theta}\Big[\sum_{t\ge 0}\gamma^t r_t\Big],\tag{62}\] which supports long-horizon credit assignment, planning, memory, and learned tool use beyond hand-crafted heuristics [2], [91].

Many recent multi-task NLP systems instantiate multi-agent formulations in which \(N\) interacting LLM agents coordinate via messages and actions. This setting can be formalized as a Markov game with state \(s_t\), joint actions \(a_t=(a_t^1,\dots,a_t^N)\), transition \(P(s_{t+1}\mid s_t,a_t)\), and agent-specific rewards \(r_t^i\). Each agent has a (possibly shared) policy \(\pi_{\theta_i}(a_t^i\mid o_t^i)\) and optimizes an expected return; for cooperative tasks, a common objective is \[\label{eq:marl95objective} \max_{\theta_1,\dots,\theta_N}\;\mathbb{E}\Big[\sum_{t\ge 0}\gamma^t\,\sum_{i=1}^N r_t^i\Big],\tag{63}\] while mixed cooperative–competitive settings may seek equilibria [108]. Within LLM-based architectures, this formalism underlies RL-driven multi-agent LLM systems along two axes — whether agent configurations adapt per task (dynamic) and whether the LLM backbones are trained—and catalogs RL-based approaches such as GPTSwarm [109], MaAS [110], G-Designer [111], MAPoRL [112], MLPO [113], ReMA [114], FlowReasoner [115], CURE [116], MMedAgent-RL [117], Chain-of-Agents (COA) [118], RLCCF [119], and MAGRPO [120], together with their optimization choices (e.g., PPO/GRPO or bespoke variants). Representative designs separate higher-level meta-planning from execution and share parameters across roles (ReMA), use multi-dimensional rewards tied to execution feedback for query-level meta-agents (FlowReasoner), and combine MARL with LLM-generated hybrid rewards plus evolutionary search to address credit assignment and partial observability (LERO) [121]. For domain-specific cooperation, code-centric systems co-train a generator and tester to enrich reward signals and generalize across benchmarks (CURE), while medical Visual Question Answering frameworks coordinate generalists and specialists with a curriculum, outperforming prior Medical Large Vision–Language Models and yielding more human-like diagnostic behavior (MMedAgent-RL). At the foundation-model end, COA distills trajectories from strong multi-agent systems and then applies agentic RL with carefully designed rewards to produce Agent Foundation Models, and SPIRAL [122] advances fully online, multi-turn self-play in zero-sum games using a shared policy with role-conditioned advantage estimation, showing that gameplay transfers to stronger mathematical and general reasoning.

The current state of the art in reinforcement learning applications shows us RL is used not only to improve answer quality and adherence to human preferences. It is of inestimable value to the industry because of reasoning. Reasoning models differ from ordinary models in that they maintain explicit chains of thought. Such models are believed to be capable of thinking in some sense: they do not merely reproduce information but can reason step by step and justify their answers. RL plays a key role in training reasoning models: during training the model receives a reward not only for the correct answer but also for the logical and consistent chain of reasoning that leads to that answer. The most interesting aspect of reasoning is that it constitutes a new paradigm for scaling LLMs. Until recently it was believed that LLMs could only become smarter by increasing the amount of training data and computational resources. However, with the advent of the first large reasoning models it was discovered that reasoning itself can scale models. The longer the model reasons and the more tokens it devotes to its chains of thought, the higher the quality of the answers. Because this scaling occurs only at inference time, it is known as test‑time scaling.

10 Acknowledgements↩︎

We thank Mark Obozov from Innopolis University for fruitful discussions and valuable suggestions on the text that helped lead to this work.

References↩︎

---

1. Here we hide poly-logarithmic factors in asymptotic bound by using a notation \(\widetilde{\mathcal{O}}(f(x))\) that is an upper bound on \(f(x)\) up to constant and poly-logarithmic factors for sufficiently large \(x\).↩︎

2. Here we also hide poly-logarithmic factors in asymptotic bound by using a notation \(\widetilde{\Omega}(f(x))\) that is a lower bound on \(f(x)\) up to constant and poly-logarithmic factors for sufficiently large \(x\).↩︎

3. When more data is available, classical algorithms for known users (such as collaborative filtering) can be applied. However, these approaches do not work well for newly registered “cold” users. To address this, the problem needs to be formulated in a contextual setting.↩︎

4. In the literature this quantity is often called pseudo-regret.↩︎

5. The setting with a random initial state sampled from initial distribution \(\mu\) could be modelled by new artificial state with prescribed transition from this state.↩︎

6. Moreover, it will be maximum likelihood estimate↩︎

7. Since the set of all models is a product of simplices, it is compact and therefore maximum attains.↩︎