Score-Aware Policy-Gradient and Performance
Guarantees using Local Lyapunov Stability

Céline Comte, Matthieu Jonckheere, Jaron Sanders, Albert Senen-Cerda

► To cite this version:

Céline Comte, Matthieu Jonckheere, Jaron Sanders, Albert Senen-Cerda. Score-Aware Policy-Gradient and Performance Guarantees using Local Lyapunov Stability. Journal of Machine Learning Research, 2025, 26 (132), pp.1-74. hal-04329790v3

HAL Id: hal-04329790

https://hal.science/hal-04329790v3

Submitted on 13 Oct 2025

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Score-Aware Policy-Gradient and Performance Guarantees using Local Lyapunov Stability

Céline Comte

Matthieu Jonckheere

LAAS-CNRS, Université de Toulouse, CNRS, Toulouse, France
7 Avenue du Colonel Roche, 31400 Toulouse, France

Jaron Sanders

Eindhoven University of Technology, Eindhoven, The Netherlands
MetaForum, Groene Loper 5, 5612 AZ Eindhoven, The Netherlands

Albert Senen-Cerda

LAAS-CNRS, IRIT, and Université de Toulouse, Toulouse, France
7 Avenue du Colonel Roche, 31400 Toulouse, France

CELINE.COMTE@CNRS.FR

MATTHIEU.JONCKHEERE@LAAS.FR

JARON.SANDERS@TUE.NL

ALBERT.SENEN-CERDA@IRIT.FR

**Editor:** Nan Jiang

## Abstract

In this paper, we introduce a policy-gradient method for model-based reinforcement learning (RL) that exploits a type of stationary distributions commonly obtained from Markov decision processes (MDPs) in stochastic networks, queueing systems, and statistical mechanics. Specifically, when the stationary distribution of the MDP belongs to an exponential family that is parametrized by policy parameters, we can improve existing policy gradient methods for average-reward RL. Our key identification is a family of gradient estimators, called score-aware gradient estimators (SAGEs), that enable policy gradient estimation without relying on value-function estimation in the aforementioned setting. We show that SAGE-based policy-gradient locally converges, and we obtain its regret. This includes cases when the state space of the MDP is countable and unstable policies can exist. Under appropriate assumptions such as starting sufficiently close to a maximizer and the existence of a local Lyapunov function, the policy under SAGE-based stochastic gradient ascent has an overwhelming probability of converging to the associated optimal policy. Furthermore, we conduct a numerical comparison between a SAGE-based policy-gradient method and an actor-critic method on several examples inspired from stochastic networks, queueing systems, and models derived from statistical physics. Our results demonstrate that a SAGE-based method finds close-to-optimal policies faster than an actor-critic method.

**Keywords:** reinforcement learning, policy-gradient method, exponential families, product-form stationary distribution, stochastic approximation

## 1 Introduction

Reinforcement learning (RL) has become the primary tool for optimizing controls in uncertain environments. Model-free RL, in particular, can be used to solve generic Markov decision processes (MDPs) with unknown dynamics with an agent that learns to maximize a reward incurred upon acting on the environment. In stochastic systems, examples of pos-
sible applications of RL can be found in stochastic networks, queueing systems, and particle
systems, where an optimal policy is desirable. For example, a policy yielding a good routing
policy, an efficient scheduling, or an annealing schedule to reach a desired state.

As stochastic systems expand in size and complexity, however, the RL agent must deal
with large state and action spaces. This leads to several computational concerns, namely,
the combinatorial explosion of action choices, the computationally intensive exploration and
evaluation of policies (Qian et al., 2019), and a more complex optimization landscape.

One way to circumvent issues pertaining to large state spaces and/or nonconvex objective functions is to include features of the underlying MDP in the RL algorithm. If the model class of the environment is known, a model-based RL approach estimates first an approximate model of the environment in the class that can later be used to solve an MDP describing its approximate dynamics. This approach is common in queueing networks (Liu et al., 2022; Anselmi et al., 2024). Nevertheless, solving an approximate MDP adds a computational burden if the number of states is large.

Policy-gradient methods are learning algorithms that instead directly optimize policy
parameters through stochastic gradient ascent (SGA) (Sutton and Barto, 2018). These
methods have gained attention and popularity due to their perceived ability to handle
large state and action spaces in model-free settings (Daneshmand et al., 2018; Khadka and
Tumer, 2018). Policy-gradient methods rely on the estimation of value functions, which
encode reward-weighted representations of the underlying model dynamics. Computing
such functions, however, is challenging in high-dimensional settings and, different from a
model-based approach, key model features are initially unknown.

In this paper, we improve policy-gradient methods for some stochastic systems by incor-
porating model-specific information of the MDP into the gradient estimator. Specifically,
we exploit the fact that long-term average behavior of such systems are described using
exponential families of distributions. In the context of stochastic networks and queueing
systems, this typically means that the Markov chains associated to fixed policies have a
product-form stationary distribution. This structural assumption holds in various relevant
scenarios, including Jackson and Whittle networks (Serfozo, 1999, Chapter 1), BCMP net-
works (Baskett et al., 1975), and more recent models arising in datacenter scheduling and
online matching (Gardner and Righter, 2020). By encoding this key model feature into
policy-gradient methods, we aim to expand the current model-based RL techniques for
control policies of stochastic systems.

Our primary contributions are the following:

* We present a new gradient estimator for policy-gradient methods that incorporates information from the stationary measure of the MDP. Under an average-reward and infinite-horizon learning setting, we namely consider policy parametrizations such that there is a known relationship between the policy on the one hand, and the MDP’s stationary distribution on the other hand. In practice, this translates to assuming that the stationary distribution forms an exponential family explicitly depending on the policy parameters. Using this structure, we define score-aware gradient estimators (SAGEs), a class of estimators that exploit the aforementioned assumption to estimate the policy gradient *without* relying on value or action-value functions.
* We show the local convergence and bound the regret of a SAGE-based policy-gradient under broad assumptions, such as a countable state space, nonconvex objective functions, and unbounded rewards. To do so, we first generalize the approach of Fehrman et al. (2020) to a general RL setting that includes Markovian updates with a countable state space and does not require the stationary distribution to be exponential. We show local convergence by first, using a local Lyapunov function that guarantees stability (i.e., positive recurrence) of the Markov chain as long as the iterates are close to the optimum, and second, by using the nondegeneracy of the Hessian at the optimum, which allows to keep track of the updates locally. These two key elements allow us to show convergence if the bias of the gradient estimator and its variance can be controlled. Remarkably, our local assumptions may be satisfied even when unstable policies exist. For policy gradient with a SAGE in particular, we can then crucially estimate its bias and variance explicitly due to the exponential family assumption, and show convergence with large probability by using the aforementioned approach, whenever the trajectory of the iterates gets close enough to an optimum. The convergence proof approach is of independent interest, and can also be adapted to other policy gradient-based methods as long as the bias and variance of the policy-gradient estimator can be controlled.

* We numerically evaluate the performance of SAGE-based policy-gradient on several models from stochastic networks, queueing systems, and statistical physics. We observe that, compared to an actor-critic algorithm, SAGE-based policy-gradient methods exhibit faster convergence and lower variance.

Our results suggest that exploiting model-specific information is a promising approach to improve RL algorithms, especially for stochastic networks and queueing systems. Sections 1.1 and 1.2 below describe our contributions in more details.

## 1.1 Score-Aware Gradient Estimators (SAGEs)

We introduce SAGEs for MDPs following the exponential-family assumption in Section 4. These estimators leverage the structure of the stationary distribution, with the goal of reducing variance and favoring stable learning. Notably, their usage requires neither knowledge nor explicit estimation of model parameters, ensuring practical applicability. The key step of the derivation exploits information on the form of the score of exponential families—that is, the gradient of the logarithm of the probability mass function.

We can illustrate the working principle using a toy example on a countable state space $S$: given a sufficient statistic $x: S \to \mathbb{R}^n$, the associated exponential family in canonical form is the family of distributions with probability mass functions $p(\cdot|\theta) \propto \exp(\theta^T x(\cdot))$ parametrized by $\theta \in \mathbb{R}^n$. Observe now that these distributions satisfy the relation

$$ \frac{d \log(p(s|\theta))}{d\theta} = x(s) - \mathbb{E}_{S \sim p(s|\theta)}[x(S)], \quad s \in S, \qquad (1) $$

and that (1) gives an exact expression for the gradient of the score.

Now, a more general version of (1) that is also applicable beyond this toy example—see Theorem 1 below—allows us to bypass the commonly used policy-gradient theorem (Sutton
and Barto, 2018, Section 13.2), which ties the estimation of the gradient with that of first estimating value or action-value functions. A key aspect that SAGEs practically exploit is that, in the models from queueing and statistical physics that we will study, we know the sufficient statistic $x$ fully. Furthermore, such models commonly possess an 'effective dimension' that is much lower than the size of the state space and is reflected by the sufficient statistic. For example, in a load-balancing model we consider in the numerical section, an agnostic model-free RL algorithm would learn a value function defined over the complete state space, whose size grows exponentially in the number $n$ of servers. On the contrary, a SAGE only requires learning the expected value of a single $n$-dimensional vector—the sufficient statistics.

## 1.2 Convergence of Policy-Gradient Methods

We examine the convergence properties of the SAGE-based policy-gradient method theoretically in Section 5. Specifically, we consider the setting of policy-gradient RL with average rewards, which consists of finding a parameter $\theta$ such that the parametric policy $\pi(\theta) = \pi(\cdot|\cdot, \theta)$ maximizes

$$J(\theta) = \lim_{T \to \infty} \frac{1}{T} \mathbb{E} \left[ \sum_{t=1}^{T} R_t \right]. \qquad (2)$$

Here, $R_{t+1}$ denotes the reward that is given after choosing action $A_t$ while being in state $S_t$, which happens with probability $\pi(A_t|S_t, \theta)$. As is common in episodic RL, we consider epochs, that is, time intervals where the parameter $\theta$ is fixed and a trajectory of the Markov chain is observed. For each epoch $m$, and under the exponential-family assumption for the stationary distribution, a SAGE yields a gradient estimator $H_m$ from a trajectory of state-action-reward tuples $(S_t, A_t, R_{t+1})$ sampled from a policy with $\Theta_m$ as an epoch-dependent parameter. Convergence analysis of the SAGE-based policy-gradient method aligns with ascent algorithms like SGA by considering updates at the end of epoch $m$ with step-size $\alpha_m > 0$,

$$\Theta_{m+1} = \Theta_m + \alpha_m H_m. \qquad (3)$$

Convergence analyses for policy-gradient RL and SGA are quite standard; see Section 2. Our work specifically aligns with the framework of Fehrman et al. (2020), who study local convergence of unbiased stochastic gradient descent (SGD), that is, when the conditional estimator $H_m$ of $\nabla J(\Theta_m)$ on the past $\mathcal{F}$ is unbiased. This occurs typically in supervised learning. A main contribution in our work consists in expanding the results of Fehrman et al. (2020) to the case of Markovian data, leading to biased estimators (i.e., $\mathbb{E}[H_m|\mathcal{F}] \neq \nabla J(\Theta_m)$). In our RL setting, we handle potentially unbounded rewards and unbounded state spaces as well as the existence of unstable policies. We also assume an online application of the policy-gradient method, where restarts are impractical or costly: the last state of the prior epoch is used as the initial state for the next, distinguishing our work from typical episodic RL setups where an initial state $S_0$ is sampled from a predetermined distribution.

Our main result in Section 5 shows convergence of the SAGE parameter update in (3) to the set $\mathcal{M}$ that attains the maximum $J^*$ of (2), assuming nondegeneracy of $J$ on $\mathcal{M}$ and existence of a local Lyapunov function. If the trajectory of SGA ends up within a
sufficiently small neighborhood $V$ of a maximizer $\theta^* \in \mathcal{M}$, then with appropriate epoch length and step-sizes, convergence to $\mathcal{M}$ occurs with large probability: for any epoch $m > 0$ and $\epsilon > 0$, if $\Theta_0$ is the first iterate in $V$,

$$ \mathbb{P}[J^* - J(\Theta_m) > \epsilon | \Theta_0 \in V] \le O\left(\epsilon^{-2}m^{-\sigma-\kappa} + m^{1-\sigma/2-\kappa/2} + m^{-\kappa/2} + \frac{\alpha^2}{\ell}\right), \quad (4) $$

where the parameters $\sigma \in (2/3, 1)$, $\kappa > 0$, $\alpha \in (0, \alpha_0]$, and $\ell \in [\ell_0, \infty)$ depend on the step and batch sizes and can be tuned to make the bound in (4) arbitrarily small. While our focus is on the global optimum, the bound (4) also holds under the same assumptions in case $J^*$ is a local optimum instead.

Our key assumption relies on the existence of a local Lyapunov function in the neighborhood $V$, which ensures stability. Crucially, this assumption is required only for policies that are close to an optimal policy. This sets our work further apart from others in the RL literature, which typically require the existence of a global Lyapunov function and/or the state space to be finite. In fact, our numerical results in Section 6 show an instance where stability around a global optimizer is sufficient, highlighting the benefits of SAGE. The set $\mathcal{M}$ of global maxima is also not required to be finite or convex, thanks to the local nondegeneracy assumption.

For a large batch index $m$, the bound in (4) can be made arbitrarily small by setting the initial step size $\alpha$ small and batch size $\ell$ large. In (4), the chance that the policy escapes the set $V$, outside which stability cannot be guaranteed, does not vanish when $m \to \infty$; it remains as $\alpha^2/\ell$. We show that this term is inherent to the statistical estimation error of the gradient and cannot be avoided. Specifically, for any $\beta > 0$, there are objective functions $f$ such that $\mathbb{P}[f(\theta^*) - f(\Theta_m) > \epsilon | \Theta_0 \in V] > c\alpha^{2+\beta}/\ell$ for some $c > 0$. Hence, a lower bound shows that the proof method cannot be improved without using additional structure of $H_m$ or $J$. Furthermore, our proof can be adapted to other generic policy-gradients that have similar bounds on the gradient estimator $H_m$ as those of SAGE, thus showing that such phenomenon not just happens with SAGE but also with any other policy-gradient algorithm with similar properties of the gradient estimator.

Denoting by $T \ge 1$ the total number of samples drawn by the algorithm, we obtain from (4) a regret bound of our algorithm when reaching the set $V$. In the case of bounded rewards, we namely show that for any $0 < \epsilon < 1$, if $\Theta_0$ is the first iterate in $V$, then

$$ \mathbb{E}[T J^* - \sum_{t=1}^{T} r(S_t, A_t) | \Theta_0 \in V] = O\left((\mathcal{L}^*)^{\frac{1}{3}} T^{\frac{2}{3}+\epsilon} + \frac{\alpha^2}{\ell} T\right). \quad (5) $$

The linear term in (5) arises from the estimation error of the policy gradient and has been seen in other recent works (Abbasi-Yadkori et al., 2019), and also for countable state spaces (Murthy et al., 2024). The other term is sublinear, and its coefficient $\mathcal{L}^*$ characterizes an ‘effective’ size of the state space, and directly depends on the local Lyapunov function. For a given $T$, we can find $\ell$—a term related to the batchsize and thus estimation error—such that the regret becomes $O(T^{3/4+\epsilon})$. Remarkably, while we start from a stable policy, the expectation in (5) includes trajectories where policies may be unstable. When the reward is unbounded, we similarly obtain a bound without a linear term if we restrict to trajectories in $V$.
For cases where the optimum is reached only as $\Theta_m \to \infty$, as with deterministic policies, we show that adding a small relative entropy regularization term to $J(\theta)$ ensures that maxima are bounded and that $M$ satisfies the nondegeneracy assumption required to show local convergence.

## 1.3 Numerical Experiments

We finally assess the applicability of the SAGE-based policy-gradient algorithm in Section 6 by comparing its performance with that of the actor-critic algorithm on three models from queueing systems, stochastic networks, and statistical physics. Specifically, we consider an admission control problem in the M/M/1 queue, a load balancing system, and the Ising model with Glauber dynamics.

These numerical results suggest that, when applicable, SAGEs can expedite convergence towards an optimal policy (compared to actor-critic) by leveraging the structure of the stationary distribution. Furthermore, the lower variance of SAGE becomes decisive when stability is not guaranteed for all policies. Namely, we observe in an example that the SAGE-based policy-gradient method converges to a close-to-optimal policy even if some policies are unstable, provided that a stable policy is used as initialization. This behavior contrasts with actor-critic, whose output policies are not always stable. SAGE also reproduces a well-known phenomenon in annealing schedules for Ising models. Specifically, the agent momentarily increases the temperature in order to escape stable states that do not correspond to the global optimum.

# 2 Related Works

The work in the present manuscript resides at the intersection of distinct lines of research. We therefore broadly review, relate, and position our work to other research in this section.

## 2.1 Gradient Estimation, Exponential Families, and Product Forms

Operations on high-dimensional probability distributions, such as marginalization and inference, are numerically intractable in general. Exponential families—see Section 4.1 for a definition—are parametric sets of distributions that lead to more tractable operations and approximations while also capturing well-known probability distributions, such as probabilistic graphical models (Wainwright and Jordan, 2018), popular in machine learning. In the context of stochastic networks and queueing systems, the stationary distribution of many product-form systems can be seen as forming an exponential family.

Our first contribution is related to several works on exponential families, product-form distributions, and probabilistic graphical models. Key performance metrics in these distributions are numerically intractable a priori, but can be expressed as expectations of random vectors that can be sampled by simulation. The most basic and well-known result, which appears in Section 1.1 and will be exploited in Section 4.2, rewrites the gradient of the logarithm of the normalizing constant (a.k.a. the log-partition function) as the expectation of the model's sufficient statistics. In probabilistic graphical models, this relation has been mainly used to learn a distribution that best describes a data set via SGD (Wainwright and Jordan, 2018; Koller and Friedman, 2009). In stochastic networks, this relation
has been applied to analyze systems with *known* parameters, for instance to predict their performance (de Souza e Silva and Muntz, 1988; Zachary and Ziedins, 1999; Bonald and Virtamo, 2004; Shah, 2011; Shah and de Veciana, 2015), to characterize their asymptotic behavior in scaling regimes (Shah, 2011; Shah and de Veciana, 2015), for sensitivity analysis (de Souza e Silva and Muntz, 1988; Liu and Nain, 1991), and occasionally to optimize control parameters via gradient ascent (Liu and Nain, 1991; de Souza e Silva and Gerla, 1991; Shah, 2011).

To the best of our knowledge, an approach similar to ours is that introduced by *Sanders et al.* (2016). This work derives a gradient estimator and performs SGA in a class of product-form reversible networks. However, the procedure requires first estimating the stationary distribution, convergence is proven only for convex objective functions, and the focus is more on developing a distributed algorithm than on canonical RL. The algorithm of *Jiang and Walrand* (2009) is similarly noteworthy, although the focus there is on developing a distributed control algorithm specifically for wireless networks and not general product-form networks.

## 2.2 Stochastic Gradient Ascent (SGA) and Policy-Gradient Methods

When a gradient is estimated using samples from a Markov chain, methods from Markov Chain Monte Carlo (MCMC) are commonly used (Mohamed et al., 2020). In our case, we have moreover bias from being unable to restart the chain at each epoch. Convergence of biased SGD to approximate stationary points of smooth nonconvex functions—points $\theta$ such that $|\nabla J(\theta)| < \epsilon$ for some $\epsilon > 0$—has been addressed in the literature (Tadic and Doucet, 2017; Atchadé et al., 2017; Karimi et al., 2019; Doan et al., 2020). The asymptotic conditions for local convergence to a stationary point were first investigated by Tadic and Doucet (2017), who assumed conditions for the asymptotic stochastic variance of the gradient estimator and bias (see Assumptions 2.1–2.3, Tadic and Doucet 2017). Karimi et al. (2019) showed a nonasymptotic analysis of biased SGD. Under Lipschitz assumptions on the transition probabilities and bounded variance of the gradient estimator $H_m$, Karimi et al. (2019) showed that, under appropriate step-sizes, for some $m^* \le M$, we have $E[|\nabla J(\Theta_{m^*})|^2] = O(\log(M)/\sqrt{M})$, where $M$ is a time horizon. Tadic and Doucet (2017) and Karimi et al. (2019) applied these results in an RL context. While these works demonstrate convergence to stationary points, our contribution lies in proving convergence to a maximum, albeit locally. This approach is essential for addressing scenarios with only local assumptions and potentially unstable policies (i.e., policies such that the corresponding Markov chain is not positive recurrent).

Finally, several recent works build on gradient domination for policy-gradient methods, addressing convexity limitations and ensuring global convergence (Fazel et al., 2018; Agarwal et al., 2021; Xiao, 2022; Kumar et al., 2024). Notable differences to our work are their use of a finite state space and that we do not consider an episodic setting with restarts, as well as distinct structural assumptions on policy parametrization like natural gradients. Murthy et al. (2024) considered a convergent natural policy-gradient (NPG) method for countable state spaces where the cost is the norm of the state, i.e., the queue length in the queueing setting. In this paper, a stable max-weight policy is applied with a probability that tends to 1 as the queue size increases, therefore avoiding instability issues. We do not use such a
stabilizing policy, and instead require the algorithm to start sufficiently close to an optimal policy. Another unique aspect of our contribution lies in specialized gradient estimation schemes based on the exponential family assumption on the stationary distribution, which crucially avoids estimating or learning value functions, common in all previous works.

A succinct non-exhaustive table is provided in Table 2 to guide interested readers to key references within this vast body of literature.

<table><caption>Table 2: This table provides a non-exhaustive list of convergence results to help contextualize our contribution. Compared to the (few) results that can be applied to MDPs with a countably-infinite state space, our convergence result only requires local convexity and positive-recurrence assumptions in the parameter space.</caption><thead><tr><th>Method</th><th>Context/State Space</th><th>Main Assumptions</th><th>Convergence</th><th>References</th></tr></thead><tbody><tr><td>SGD</td><td>Markovian data with convex objectives</td><td>Lipschitz transition probabilities, bounded variance estimators</td><td>Local convergence</td><td>Tadic and Doucet (2017); Daneshmand et al. (2018); Karimi et al. (2019)</td></tr><tr><td>SGD</td><td>Id data with non convex objectives</td><td>Local convexity</td><td>Local convergence</td><td>Fehrman et al. (2020)</td></tr><tr><td>Policy Gradient</td><td>RL, Finite state space</td><td>Gradient domination</td><td>Global convergence</td><td>Agarwal et al. (2021), Kumar et al. (2024)</td></tr><tr><td>Policy Gradient</td><td>RL, Finite state space</td><td>ABC condition</td><td>Global convergence</td><td>Yuan et al. (2022)</td></tr><tr><td>Natural Policy Gradient</td><td>RL, Finite state space</td><td>Gradient domination / entropy regularization</td><td>Global convergence</td><td>Liu et al. (2020); Agarwal et al. (2021), Cen et al. (2022)</td></tr><tr><td>Natural Policy Gradient</td><td>RL, Countable state space</td><td>Stabilization via a known policy independently of the model and policy parameters</td><td>Global convergence</td><td>Murthy et al. (2024)</td></tr><tr><td>Policy gradient for product-form Networks</td><td>Countable state space</td><td>Convex objectives, requires knowing the functional form of the gradient in terms of the stationary distribution</td><td>Global convergence for convex objectives</td><td>Sanders et al. (2016)</td></tr><tr><td>SAGE-based policy gradient</td><td>Countable state space</td><td>Lyapunov assumptions, Non-degenerate Hessian, Application to exponential family stationary distributions</td><td>Local convergence</td><td>This paper</td></tr></tbody></table>

# 3 Problem Formulation

After introducing basic notation in Section 3.1, we introduce MDPs in Section 3.2 and the infinite-horizon average-reward optimality criterion in Section 3.3. Section 3.4 gives a brief introduction to policy-gradient algorithms.

## 3.1 Basic Notation

The sets of nonnegative integers, positive integers, reals, and nonnegative reals are denoted by $\mathbb{N}$, $\mathbb{N}_+$, $\mathbb{R}$, and $\mathbb{R}_{\ge 0}$, respectively. For a differentiable function $f: \theta \in \mathbb{R}^n \to f(\theta) \in \mathbb{R}$, $\nabla f(\theta)$ denotes the gradient of $f$ taken at $\theta \in \mathbb{R}^n$, that is, the $n$-dimensional column vector whose $j$-th component is the partial derivative of $f$ with respect to $\theta_j$, for $j \in \{1, 2, \dots, n\}$. If $f$ is twice differentiable, $\text{Hess}_\theta f$ denotes the Hessian of $f$ at $\theta$, that is, the $n \times n$ matrix of second partial derivatives. For a differentiable vector function $f: \theta \in \mathbb{R}^n \to f(\theta) = (f_1(\theta), \dots, f_d(\theta)) \in \mathbb{R}^d$, $Df(\theta)$ is the Jacobian matrix of $f$ taken at $\theta$, that is, the $d \times n$ matrix whose $i$-th row is $\nabla f_i(\theta)^\top$, for $i \in \{1, 2, \dots, d\}$. For a vector $x = (x_1, \dots, x_n) \in \mathbb{R}^n$,
we denote its $l_2$-norm by $|x| = \sqrt{x_1^2 + \cdots + x_n^2}$. We define the operator norm of a matrix $A \in \mathbb{R}^{a \times b}$ as $|A|_{\text{op}} = \sup_{x \in \mathbb{R}^b:|x|=1} |Ax|$. We use uppercase to denote random variables and vectors, and a calligraphic font for their sets of outcomes.

## 3.2 Markov Decision Process (MDP)

We consider a Markov decision process (MDP) with countable state, action, and reward spaces $S$, $\mathcal{A}$, and $\mathcal{R}$, respectively, and transition probability kernel $P: (s, a, r, s') \in S \times \mathcal{A} \times \mathcal{R} \times S \mapsto P(r, s'|s, a) \in [0, 1]$. Thus $P(r, s'|s, a)$ gives the conditional probability that the next reward-state pair is $(r, s')$ given that the current state-action pair is $(s, a)$. With a slight abuse of notation, we introduce

$$P(r|s, a) = \sum_{s' \in S} P(r, s'|s, a), \quad \text{for each } s \in S, a \in \mathcal{A}, \text{ and } r \in \mathcal{R}, \text{ and}$$

$$P(s'|s, a) = \sum_{r \in \mathcal{R}} P(r, s'|s, a), \quad \text{for each } s, s' \in S \text{ and } a \in \mathcal{A}.$$

Our results also generalize to absolutely continuous rewards; an example will appear in Section 6.1.

Following the framework of policy-gradient algorithms (Sutton and Barto, 2018, Chapter 13), we assume that the agent is given a random policy parametrization $\pi: (s, \theta, a) \in S \times \mathbb{R}^n \times \mathcal{A} \to \pi(a|s, \theta) \in (0, 1)$, such that $\pi(a|s, \theta)$ is the conditional probability that the next action is $a \in \mathcal{A}$ given that the current state is $s \in S$ and the parameter vector $\theta \in \mathbb{R}^n$. We assume that the function $\theta \mapsto \pi(a|s, \theta)$ is differentiable for each $(s, a) \in S \times \mathcal{A}$. The goal of the learning algorithm will be to find a parameter (vector) that maximizes the long-run average reward; see Section 3.3.

As a concrete example, we will often consider a class of softmax policies that depend on a feature extraction map $\xi: S \times \mathcal{A} \to \mathbb{R}^n$ as follows:

$$\pi(a|s, \theta) = \frac{e^{\theta^\top \xi(s, a)}}{\sum_{a' \in \mathcal{A}} e^{\theta^\top \xi(s, a')}}, \quad s \in S, a \in \mathcal{A}. \qquad (6)$$

The feature extraction map $\xi$ may leverage prior information on the system dynamics. In queueing systems for instance, we may decide to make similar decisions in large states, as these states are typically visited rarely, and it may be beneficial to aggregate the information collected about them.

## 3.3 Stationary Analysis and Optimality Criterion

Given $\theta \in \mathbb{R}^n$, if the agent applies the policy $\pi(\theta): (s, a) \in S \times \mathcal{A} \mapsto \pi(a|s, \theta)$ at every time step, the random state-action-reward sequence $((S_t, A_t, R_{t+1}), t \in \mathbb{N})$ obtained by running this policy is a Markov chain such that, for each $s, s' \in S$, $a \in \mathcal{A}$, and $r \in \mathcal{R}$, $P[A_t = a|S_t = s] = \pi(a|s, \theta)$ and $P[R_{t+1} = r, S_{t+1} = s'|S_t = s, A_t = a] = P(r, s'|s, a)$. The dependency of the random variables on the parameter is left implicit to avoid cluttering notation. Leaving aside actions and rewards, the state sequence $(S_t, t \in \mathbb{N})$ also defines a Markov chain, with transition probability kernel $P(\theta): (s, s') \in S \times S \mapsto P(s'|s, \theta)$ given by

$$P(s'|s, \theta) = \sum_{a \in \mathcal{A}} \pi(a|s, \theta)P(s'|s, a), \quad s, s' \in S.$$
In the remainder, we assume that Assumptions 1 and 2 below are satisfied.

**Assumption 1** There exists an open set $\Omega \subseteq \mathbb{R}^n$ such that, for each $\theta \in \Omega$, the Markov chain $(S_t, t \in \mathbb{N})$ with transition probability kernel $P(\theta)$ is irreducible and positive recurrent.

In the remainder, we use the words positive recurrent and stable interchangeably. Also, with a slight abuse of language, we say that a policy is stable if the corresponding Markov chain is stable. Thanks to Assumption 1, for each $\theta \in \Omega$, the corresponding Markov chain $(S_t, t \in \mathbb{N})$ has a unique stationary distribution $p(\cdot|\theta)$. We say that a triplet $(S, A, R)$ of random variables is a stationary state-action-reward triplet, and we write $(S, A, R) \sim \text{STAT}(\theta)$, if $(S, A, R)$ follows the stationary distribution of the Markov chain $((S_t, A_t, R_{t+1}), t \in \mathbb{N})$ given by

$$ \mathbb{P}[S = s, A = a, R = r] = p(s|\theta)\pi(a|s, \theta)P(r|s, a), \quad s \in \mathcal{S}, \ a \in \mathcal{A}, \ r \in \mathcal{R}. \quad (\text{STAT}(\theta)) $$

**Assumption 2** For each $\theta \in \Omega$, the stationary state-action-reward triplet $(S, A, R) \sim \text{STAT}(\theta)$ is such that $|R|, |R \nabla \log p(S|\theta)|$, and $|R \nabla \log \pi(A|S, \theta)|$ have a finite expectation.

By ergodicity (Brémaud, 1999, Theorem 4.1), the running average reward $\frac{1}{T} \sum_{t=1}^T R_t$ tends to $J(\theta)$ as defined in (2) almost surely as $T$ tends to infinity. $J(\theta)$ is called the long-run average reward and is also given by

$$ J(\theta) = \mathbb{E}[R] = \sum_{s \in \mathcal{S}} \sum_{a \in \mathcal{A}} \sum_{r \in \mathcal{R}} p(s|\theta)\pi(a|s, \theta)P(r|s, a)r, \quad \theta \in \Omega. \quad (7) $$

Our end goal, further developed in Section 3.4, is to find a learning algorithm that maximizes the objective function $J$. For now, we only observe that the objective function $J: \theta \in \Omega \mapsto J(\theta)$ is differentiable thanks to Assumption 2, and that its gradient is given by

$$ \nabla J(\theta) = \sum_{s \in \mathcal{S}} \sum_{a \in \mathcal{A}} \sum_{r \in \mathcal{R}} p(s|\theta)\pi(a|s, \theta)P(r|s, a)r(\nabla \log p(s|\theta) + \nabla \log \pi(a|s, \theta)), \quad \theta \in \Omega. \quad (8) $$

In general, computing $\nabla J(\theta)$ using (8) is challenging: (i) computing $\nabla \log p(s|\theta)$ is in itself challenging because $p(s|\theta)$ depends in a complex way on the unknown transition kernel $P(r, s'|s, a)$ and the parameter $\theta$ via the policy $\pi(\theta)$, and (ii) enumerating and thus summing over the state space $\mathcal{S}$ is often practically infeasible (for instance, when the state space $\mathcal{S}$ is infinite and/or high-dimensional). Our first contribution, in Section 4, is precisely a new family of estimators for the gradient (8).

## 3.4 Learning Algorithm

In Section 3.3, we defined the objective function $J$ by considering trajectories where the agent applied a policy $\pi(\theta)$ parametrized by a constant vector $\theta$. Going back to a learning setting, we consider a state-action-reward sequence $((S_t, A_t, R_{t+1}), t \in \mathbb{N})$ and a parameter sequence $(\Theta_m, m \in \mathbb{N})$ obtained by updating the parameter periodically according to the gradient-ascent step $\Theta_{m+1} = \Theta_m + \alpha_m H_m$ introduced in (3). Here, $H_m$ is provided by a family of learning algorithms, called policy gradient. The pseudocode of a generic policy-gradient algorithm, shown in Algorithm 1, is parametrized by a sequence $0 \triangleq t_0 < t_1 < t_2 <$
**Algorithm 1** Generic policy-gradient algorithm. Examples of GRADIENT procedure, based on different estimators for the gradient $\nabla J$, are given in Algorithms 2 and 3. All variables of Algorithm 1 are accessible within the GRADIENT procedure.

1: **Input:** • Observation times $0 \triangleq t_0 < t_1 < t_2 < \dots$
• Step size sequence $\alpha_0, \alpha_1, \alpha_2, \dots > 0$
• Positive and differentiable policy parametrization $(s, \theta, a) \mapsto \pi(a|s, \theta)$

2: **Initialization:** Policy parameter $\Theta_0 \in \Omega$ and initial state $S_0 \in \mathcal{S}$

3: **Main loop:**

4: **for** $m = 0, 1, 2, \dots$ **do**

5: &nbsp;&nbsp;&nbsp;&nbsp;**for** $t = t_m, \dots, t_{m+1} - 1$ **do**

6: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Sample $A_t \sim \pi(\cdot|S_t, \Theta_m)$

7: &nbsp;&nbsp;&nbsp;&nbsp;Take action $A_t$ and observe $R_{t+1}, S_{t+1}$

8: &nbsp;&nbsp;&nbsp;&nbsp;**end for**

9: &nbsp;&nbsp;&nbsp;&nbsp;Update $\Theta_{m+1} \leftarrow \Theta_m + \alpha_m \text{GRADIENT}(m)$

10: **end for**

... of observation times and a sequence $\alpha_0, \alpha_1, \alpha_2, \dots > 0$ of step sizes. For each $m \in \mathbb{N}$, $\mathcal{D}_m$ denotes batch $m$, obtained by applying policy $\pi(\Theta_m)$ at epoch $m$, given by

$$ \mathcal{D}_m = ((S_t, A_t, R_{t+1}), t \in \{t_m, \dots, t_{m+1} - 1\}). \quad (9) $$

Given some initialization $\Theta_0$, Algorithm 1 calls the GRADIENT procedure that computes an estimate $H_m$ of $\nabla J(\Theta_m)$ from $\mathcal{D}_m$, and it updates the parameter according to (3).

As discussed at the end of Section 3.3, finding an estimator $H_m$ for $\nabla J(\Theta_m)$ directly from (7) is difficult in general. A common way to obtain $H_m$ follows from the policy-gradient theorem (Sutton and Barto, 2018, Chapter 13), which instead writes the gradient $\nabla J(\theta)$ using the action-value function $q$:

$$ \nabla J(\theta) = \mathbb{E}[q(S, A) \nabla \log \pi(A|S, \theta)], $$

where $(S, A, R) \sim \text{STAT}(\theta)$, for each $\theta \in \Omega$. Consistently, in a model-free setting, policy-gradient methods like the actor-critic algorithm (recalled in Appendix A.1) estimate $\nabla J(\Theta_m)$ by first estimating a value function. However, this approach can suffer from high-variance of the estimator, which slows down convergence, as described in Section 1. Some of these problems can be circumvented by exploiting the problem structure, as we will see now.

## 4 Score-Aware Gradient Estimator (SAGE)

We now define the key structural assumption in our paper. Namely, that we have information on the impact of the policy parameter $\theta$ on the stationary distribution $p$. In Section 4.2, we will use this assumption to build SAGEs, a new family of estimators for the gradient $\nabla J$ that do not involve the state-value function, contrary to actor-critic. In Section 4.3, we will further explain how to use this insight to design a SAGE-based policy-gradient method.
## 4.1 Product-Form and Exponential Family

As announced in the introduction, our end goal is to design a gradient estimator capable of exploiting information on the stationary distribution $p(\cdot|\theta)$ of the MDP when such information is available. Assumption 3 below formalizes this idea by assuming that the stationary distribution forms an exponential family parametrized by the policy parameter $\theta$.

**Assumption 3 (Stationary Distribution)** There exist a scalar function $\Phi: S \to \mathbb{R}_{>0}$, an integer $d \in \mathbb{N}_+$, a differentiable vector function $\rho: \Omega \to \mathbb{R}_{>0}^d$, and a vector function $x: S \to \mathbb{R}^d$ such that the following two equivalent equations are satisfied:

$$p(s|\theta) = \frac{1}{Z(\theta)} \Phi(s) \prod_{i=1}^{d} \rho_i(\theta)^{x_i(s)}, \quad s \in S, \quad \theta \in \Omega, \quad (10-\text{PF})$$

$$\log p(s|\theta) = \log \Phi(s) + \log \rho(\theta)^{\top} x(s) - \log Z(\theta), \quad s \in S, \quad \theta \in \Omega, \quad (10-\text{EF})$$

where the partition function $Z: \Omega \to \mathbb{R}_{>0}$ follows by normalization:

$$Z(\theta) = \sum_{s \in S} \Phi(s) \prod_{i=1}^{d} \rho_i(\theta)^{x_i(s)} = \sum_{s \in S} e^{\log \Phi(s) + \log \rho(\theta)^{\top} x(s)}, \quad \theta \in \Omega. \quad (11)$$

We will call $\Phi$ the balance function, $\rho$ the load function, and $x$ the sufficient statistics.

(10-PF) is the product-form variant of the stationary distribution, classical in queueing theory. (10-EF) is the exponential-family description of the distribution. This latter representation is more classical in machine learning (Wainwright and Jordan, 2018) and will simplify our derivations. Let us briefly discuss the implications of this assumption as well as examples where this assumption is satisfied.

Assumption 3 implies that the stationary distribution $p$ depends on the policy parameter $\theta$ only via the load function $\rho$. Yet, this assumption may not seem very restrictive a priori. Assuming for instance that the state space $S$ is finite, with $S = \{s_1, s_2, \dots, s_N\}$, we can write the stationary distribution in the form (10) with $d = N$, $\rho_i(\theta) = p(s_i|\theta)$, $x_i(s) = 1[s = s_i]$, and $\Phi(s) = Z(\theta) = 1$, for each $\theta \in \mathbb{R}^n$, $s \in S$, and $i \in \{1, 2, \dots, N\}$. However, writing the stationary distribution in this form is not helpful, in the sense that in general the function $\rho$ will be prohibitively intricate. As we will see in Section 4.2, what will prove important in Assumption 3 is that the load function $\rho$ is simple enough so that we can evaluate its Jacobian matrix function $D \log \rho$ numerically.

There is much literature on stochastic networks and queueing systems with a stationary distribution of the form (10-PF). Most works focus on performance evaluation, that is, evaluating $J(\theta)$ for some parameter $\theta \in \Omega$, assuming that the MDP's transition probability kernel is known. In this context, the product-form (10-PF) arises in Jackson and Whittle networks (Serfozo, 1999, Chapter 1), BCMP networks (Baskett et al., 1975), as well as more recent models arising in datacenter scheduling and online matching (Gardner and Righter,
2020)<fnref n="1" />. Building on this literature, in Section 6, we will consider policy parametrizations for control problems that also lead to a stationary distribution of the form (10).

In the next section, we exploit Assumption 3 to construct a gradient estimator that requires knowing the functions D log ρ and x but not the functions ρ, Φ, and Z.

## 4.2 Score-Aware Gradient Estimator (SAGE)

As our first contribution, Theorem 1 below gives simple expressions for $∇ log p(s|θ)$ and $∇J(θ)$ under Assumptions 1 to 3. Gradient estimators that will be formed using (13) will be called score-aware gradient estimators (SAGEs), to emphasize that the estimators rely on the simple expression (12) for the score $∇ log p(s|θ)$. Particular cases of this result have been obtained by de Souza e Silva and Muntz (1988); de Souza e Silva and Gerla (1991); Liu and Nain (1991) for specific stochastic networks; our proof is shorter and more general thanks to the exponential form (10-EF).

**Theorem 1** Suppose that Assumptions 1 to 3 hold. For each $θ ∈ Ω$, we have

$$∇ log p(s|θ) = D log ρ(θ)^T(x(s) - E[x(S)]), \quad (12)$$

$$∇J(θ) = D log ρ(θ)^T \mathrm{Cov}[R, x(S)] + E[R∇ log π(A|S, θ)], \quad (13)$$

where $(S, A, R) \sim \mathrm{STAT}(θ)$, $\mathrm{Cov}[R, x(S)] = (\mathrm{Cov}[R, x_1(S)], \dots, \mathrm{Cov}[R, x_d(S)])^T$, and the gradient and Jacobian operators, $∇$ and $D$ respectively, are taken with respect to $θ$.

**Proof** Applying the gradient operator to the logarithm of (11) and simplifying yields

$$∇ log Z(θ) = D log ρ(θ)^T E[x(S)]. \quad (14)$$

This equation is well-known and was already discussed in Section 2.1. Equation (12) follows by applying the gradient operator to (10-EF) and injecting (14). Equation (13) follows by injecting (12) into (8) and simplifying. ■

Assuming that the functions D log ρ and x are known in closed-form, Theorem 1 allows us to construct an estimator of $∇J(θ)$ from a state-action-reward sequence $((S_t, A_t, R_{t+1}), t \in \{0, 1, \dots, T\})$ obtained by applying policy $π(θ)$ at every time step as follows:

$$H = D log ρ(θ)^T \bar{C} + \bar{E}, \quad (15)$$

where $\bar{C}$ and $\bar{E}$ are estimators of $\mathrm{Cov}[R, x(S)]$ and $E[R∇ log π(A|S, θ)]$, respectively, obtained for instance by taking the sample mean and sample covariance. An estimator of the form (15) will be called a score-aware gradient estimator (SAGE). This idea will form the basis of the SAGE-based policy-gradient method that will be introduced in Section 4.3. Observe that such an estimator will typically be biased since the initial state $S_0$ is not

<footnote n="1">1. Although the distributions recalled in (Gardner and Righter, 2020, Theorems 3.9, 3.10, 3.13) do not seem to fit the framework of (10) a priori because the number of factors in the product can be arbitrarily large, some of these distributions can be rewritten in the form (10) by using an expanded state descriptor, as in (Adan et al., 2017, Equation 4, Corollary 2, and Theorem 6) and (Moyal et al., 2021, Equation 7 and Proposition 3.1).</footnote>
stationary. Nonetheless, we will show in the proof of the convergence result in Section 5 that this bias does not prevent convergence.

The advantage of using a SAGE is twofold. First, the challenging task of estimating $\nabla J(\theta)$ is reduced to the simpler task of estimating the $d$-dimensional covariance $\text{Cov}[R, x(S)]$ and the $n$-dimensional expectation $\mathbb{E}[R \nabla \log \pi(A|S, \theta)]$, for which leveraging estimation techniques in the literature is possible. Also recall that the gradient estimator used in the actor-critic algorithm (Appendix A.1) relies on the state-value function, so that it requires estimating $|S|$ values; we therefore anticipate SAGEs to yield better performance when $\max(n, d) \ll |S|$; see examples from Sections 6.2 and 6.3. Second, as we will also observe in Section 6, SAGEs can exploit information on the structure of the policy and stationary distribution “by design”. Actor-critic exploits this information only indirectly due to its dependency on the state-value function.

## 4.3 SAGE-Based Policy-Gradient Algorithm

Algorithm 2 introduces a SAGE-based policy-gradient method based on Theorem 1. For each $m \in \mathbb{N}$, the GRADIENT($m$) procedure is called in the gradient-update step (Line 9) of Algorithm 1, at the end of epoch $m$, and returns an estimate of $\nabla J(\Theta_m)$ based on batch $\mathcal{D}_m$, defined in (9). Algorithm 2 can be understood as follows. According to Theorem 1, we have $\nabla J(\Theta_m) = \text{D log } \rho(\Theta_m)^{\top} \text{Cov}[R, x(S)] + \mathbb{E}[R \nabla \log \pi(A|S, \Theta_m)]$ with $(S, A, R) \sim \text{STAT}(\Theta_m)$. Lines 4 to 6 estimate $\text{Cov}[R, x(S)]$ using the usual sample covariance estimator. Line 7 estimates $\mathbb{E}[R \nabla \log \pi(A|S, \theta)]$ using the usual sample mean estimator. To simplify the signature of GRADIENT($m$), we assume that all variables from Algorithm 1, in particular batch $\mathcal{D}_m$, are accessible within Algorithm 2. The variable $N_m$ computed on Line 3 is the batch size, i.e., the number of samples used to estimate the gradient $\nabla J(\Theta_m)$, and we assume that it is greater than or equal to 2. An alternate implementation of the SAGE-based policy-gradient method that allows for batch sizes equal to 1 is given in Appendix A.2.

Recall that our initial goal was to exploit information on the stationary distribution, when such information is available. Consistently, compared to actor-critic (Appendix A.1), the SAGE-based method of Algorithm 2 requires as input the Jacobian matrix function

**Algorithm 2** SAGE-based policy-gradient method, to be called on Line 9 of Algorithm 1.

1: **Input:** • Positive and differentiable policy parametrization $(s, \theta, a) \mapsto \pi(a|s, \theta)$
• Jacobian matrix function $\theta \mapsto \text{D log } \rho(\theta)$
• Feature function $s \mapsto x(s)$

2: **procedure** GRADIENT($m$)

3: $N_m \leftarrow t_{m+1} - t_m$

4: $\bar{X}_m \leftarrow \frac{1}{N_m} \sum_{t=t_m}^{t_{m+1}-1} x(S_t)$

5: $\bar{R}_m \leftarrow \frac{1}{N_m} \sum_{t=t_m}^{t_{m+1}-1} R_{t+1}$

6: $\bar{C}_m \leftarrow \frac{1}{N_{m-1}} \sum_{t=t_m}^{t_{m+1}-1} (x(S_t) - \bar{X}_m)(R_{t+1} - \bar{R}_m)$

7: $\bar{E}_m \leftarrow \frac{1}{N_m} \sum_{t=t_m}^{t_{m+1}-1} R_{t+1} \nabla \log \pi(A_t|S_t, \Theta_m)$

8: **return** $\text{D log } \rho(\Theta_m)\bar{C}_m + \bar{E}_m$

9: **end procedure**
$D \log \rho$ and the sufficient statistics $x$. In return, as we will see in Sections 5 and 6, the SAGEs-based method relies on a lower-dimensional estimator whenever $\max(n, d) \ll |S|$, which can lead to improved convergence properties.

## 5 A Local Convergence Result

Our goal in this section is to study the limiting behavior of Algorithm 2. To do so, we will consider this algorithm as an SGA algorithm that uses biased gradient estimates. The gradient estimates are biased because they arise from the MCMC estimations from Lines 4–7 in Algorithm 2. Throughout the proof, we assume that the reward is a deterministic function $r: S \times A \to \mathbb{R}$ for simplicity. Under this assumption, for each $m \in \mathbb{N}$, Algorithm 2 follows the gradient ascent step (3), with

$$ H_m = D \log \rho(\Theta_m)^\top \bar{C}_m + \bar{E}_m, \text{ where } \left\{ \begin{array}{ll} \bar{X}_m = \frac{\sum_{t=t_m}^{t_{m+1}-1} x(S_t)}{t_{m+1} - t_m}, & \bar{R}_m = \frac{\sum_{t=t_m}^{t_{m+1}-1} r(S_t, A_t)}{t_{m+1} - t_m}, \\ \bar{C}_m = \frac{\sum_{t=t_m}^{t_{m+1}-1} (x(S_t) - \bar{X}_m)(r(S_t, A_t) - \bar{R}_m)}{t_{m+1} - t_m - 1}, & \\ \bar{E}_m = \frac{\sum_{t=t_m}^{t_{m+1}-1} r(S_t, A_t) \nabla \log \pi(A_t | S_t, \Theta_m)}{t_{m+1} - t_m}. & \end{array} \right. \quad (16) $$

The estimates $\bar{X}_m$, $\bar{R}_m$, and $\bar{C}_m$ are functions of $D_m$, while $H_m$ and $\bar{E}_m$ are functions of $D_m$ and $\Theta_m$. We will additionally apply decreasing step sizes and increasing batch sizes of the form

$$ \alpha_m = \frac{\alpha}{(m+1)^\sigma} \quad \text{and} \quad t_{m+1} = t_m + \ell m^{\frac{\sigma}{2}+\kappa}, \quad \text{for each } m \in \mathbb{N}, \quad (17) $$

for some parameters $\alpha \in (0, \infty)$, $\ell \in (1, \infty)$, $\sigma \in (2/3, 1)$, and $\kappa \in [0, \infty)$.

Our goal—studying the limiting algorithmic behavior of Algorithm 2—is equivalent to studying the limiting algorithmic behavior of the stochastic recursion (3). In particular, we will focus on the local convergence of the iterates of (3) and (16) to the following set of global maximizers:

$$ \mathcal{M} = \{\theta \in \Omega : J(\theta) = J^*\}, \quad \text{where } J^* = \sup_{\theta \in \Omega} J(\theta). \quad (18) $$

We will assume that $\mathcal{M}$ is nonempty, that is, $\mathcal{M} \neq \emptyset$ and at least one $\theta \in \Omega$ satisfies $J(\theta) = J^*$. Note that Assumption 7 below allows $\mathcal{M}$ to be a manifold just locally. Consequently, $J$ can be nonconvex with noncompact level-subsets, and $J$ is even allowed not to exist outside the local neighborhood, for instance if the policy is unstable. If the policy $\pi(\theta)$ is unstable, then necessarily we have $\theta \notin \Omega$, and we adopt the convention that $J(\theta) = \inf_{s \in S, a \in A} r(s, a) \ge -\infty$. While the previous assumptions allow for general objective functions, the convergence will be guaranteed close to the set of maxima $\mathcal{M}$, or to a set of local maxima that satisfy equivalent assumptions.
## 5.1 Assumptions Pertaining to Algorithmic Convergence

We use the Markov chain of state-action pairs. Specifically, consider the pairs $\{(S_t, A_t)\}_{t \ge 0} \subset S \times A$, where $A_t$ is generated according to policy $\pi(\cdot|S_t, \theta)$. For a given $\theta \in \Omega$, the one-step transition probability and the stationary distribution of this Markov chain are

$$P((s', a')|(s, a), \theta) = \pi(a'|s', \theta)P(s'|s, a), \quad \text{for } (s, a), (s', a') \in S \times A, \quad (19)$$

$$\tilde{p}((s, a)|\theta) = p(s|\theta)\pi(a|s, \theta) \quad \text{for } (s, a) \in S \times A. \quad (20)$$

The following are assumed:

**Assumption 4** There exists a function $\mathcal{L}: S \times A \to [1, \infty)$ such that, for any $\theta^* \in \mathcal{M}$, there exist a neighborhood $U$ of $\theta^*$ in $\Omega$ and four constants $\lambda \in (0, 1)$, $C > 0$, $b \in \mathbb{R}_{>0}$, and $v \ge 16$ such that, for each $\theta \in U$, the policy $\pi(\theta)$ is such that

$$\sum_{(s', a') \in S \times A} P((s', a')|(s, a), \theta)(\mathcal{L}(s', a'))^v \le \lambda(\mathcal{L}(s, a))^v + b, \quad \text{for each } (s, a) \in S \times A,$$

and, for each $\ell \in \mathbb{N}_+$ and $(s, a), (s', a') \in S \times A$,

$$|P^{\ell}((s', a')|(s, a), \theta) - \tilde{p}((s', a')|\theta)| \le C\lambda^{\ell}\mathcal{L}(s, a),$$

where $P^{\ell}(\theta)$ is the $\ell$-step transition probability kernel of the Markov chain with transition probability kernel (19).

**Assumption 5** There exists a constant $C > 0$ such that $|D \log \rho(\theta)|_{op} < C$ for each $\theta \in \Omega$.

**Assumption 6** Let $\mathcal{L}$ be the Lyapunov function from Assumption 4. For any $\theta^* \in \mathcal{M}$, if $U$ is a local neighborhood satisfying the conditions of Assumption 4, then there exists a constant $C > 0$ such that for any $\theta \in U$ and $(s, a) \in S \times A$,

$$|x(s)| < C\mathcal{L}(s, a), \quad |r(s, a)| < C\mathcal{L}(s, a), \quad |r(s, a)\nabla \log \pi(a|s, \theta)| < C\mathcal{L}(s, a). \quad (21)$$

**Assumption 7** There exist an integer $n \in \{0, 1, \dots, n-1\}$ and an open subset $U \subseteq \Omega$ such that (i) $\mathcal{M} \cap U$ is a nonempty $n$-dimensional $\mathbb{R}^2$-submanifold of $\mathbb{R}^n$, and (ii) the Hessian of $J$ at $\theta^*$ has rank $n-n$, for each $\theta^* \in \mathcal{M} \cap U$.

These assumptions have the following interpretation. Assumption 4 formalizes that the Markov chain is geometrically ergodic in a neighborhood around the optimizer, which implies in particular that policies in this neighborhood are stable (i.e., positive recurrent). Remarkably, it does not require that the chain is geometrically ergodic for all policies, only for those close to an optimal policy when $\theta \in \Omega$. This stability is guaranteed by a local Lyapunov function $\mathcal{L}$ uniformly over some neighborhood close to a maximizer. In the notation for $b$, and $\lambda \in (0, 1)$ of Assumption 4, if $S_0$ is the initial state, a term that will later bound the size of an 'effective' state space in the regret of the algorithm is

$$\mathcal{L}^* = \max\left(\frac{b}{(1-\lambda)^2}, \max_{a \in A}\mathcal{L}(S_0, a)^v\right). \quad (22)$$
Assumptions 5 and 6 together guarantee that the estimator $H_m$ concentrates around $\nabla J(\Theta_m)$ at an appropriate rate. Assumption 5 is easy to verify in our examples since $\rho$ is always positive and bounded. Assumption 6 guarantees that the empirical estimators in (16) converge fast enough to their expected values. In many applications from queueing, Assumption 6 holds. Namely, $S$ is usually a normed space and the order of the Lyapunov function $\mathcal{L}(s, a)$ is exponential in the norm of the state $s \in S$, compared to the sufficient statistic $x$ which has an order linear in the norm of $s$. We remark that, in a setting with a bounded reward function $r$ and a bounded map $x$ or with a finite state space, Assumption 6 becomes trivial.

Assumption 7 is a geometric condition. It guarantees that, locally around the set of maxima $\mathcal{M}$ or set of local maxima satisfying the same assumptions, in directions perpendicular to $\mathcal{M}$, $J$ behaves approximately in a convex manner. Concretely, this means that $\text{Hess}_\theta J$ has strictly negative eigenvalues in the directions normal to $\mathcal{M}$—also referred to as the Hessian being *nondegenerate*. Thus, there is a one-to-one correspondence between local directions around $\theta \in \mathcal{M}$ that decrease $J$ and directions that do not belong to the tangent space of $\mathcal{M}$. Strictly concave functions satisfy that $n = 0$ and Assumption 7 is thus automatically satisfied in such cases. If $\mathcal{M} \cap U = \{\theta^*\}$ is a singleton, Assumption 7 reduces to assuming that $\text{Hess}_{\theta^*} J$ is negative definite. Assumption 7 in a general setting can be difficult to verify, but by adding a regularization term, it can be guaranteed to hold in a broad sense (see Section 5.5).

## 5.2 Local Convergence Results

This is our main convergence result for the case that the set of maxima is not necessarily bounded.

**Theorem 2 (Noncompact Case)** Suppose that Assumptions 1 to 7 hold. For every maximizer $\theta^* \in \mathcal{M} \cap U$, there exist constants $c > 0$ and $\alpha_0 > 0$ such that, for each $\alpha \in (0, \alpha_0]$, there exists a nonempty neighborhood $V$ of $\theta^*$ and $\ell_0 \ge 1$ such that, for each $\ell \in [\ell_0, \infty)$, $\sigma \in (2/3, 1)$, $\kappa \in [0, \infty)$ with $\sigma + \kappa > 1$, we have, for each $m \in \mathbb{N}_+$,

$$ \mathbb{P}[J(\Theta_m) < J^* - \epsilon | \Theta_0 \in V] \le c \left( \epsilon^{-2} \mathcal{L}^* m^{-\sigma-\kappa} + \frac{m^{1-\sigma-\kappa}}{\ell} + \frac{\alpha^2}{\ell} + \alpha m^{-\kappa/2} + \frac{\alpha m^{1-(\sigma+\kappa)/2}}{\sqrt{\ell}} \right), \quad (23) $$

where $(\Theta_m, m \in \mathbb{N})$ is a random sequence with $\mathbb{P}[\Theta_0 \in V] > 0$, and built by recursively applying the gradient ascent step (3) with the gradient update (16) and the step and batch sizes (17) parameterized by these values of $\alpha$, $\ell$, $\sigma$, and $\kappa$.

In Theorem 2, by setting the parameters $\alpha$, $\ell$, $\sigma$, and $\kappa$ in (17) appropriately, we can make the probability of $\Theta_m$ being $\epsilon$-suboptimal arbitrarily small. Specifically, the step and batch sizes for each epoch allow us to control the variance of the estimators in (16). This shows that the SAGE-based policy-gradient method converges with large probability. The bound can be understood as follows. The term in (23) on the bound depending on $\epsilon$ characterizes the convergence rate assuming that all iterates up to time $m$ remain in $V$. The remaining terms in (23) estimate the probability that the iterates escape the set $V$,
which can be made small by tuning parameters that diminish the variance of the estimator $H_m$, such as setting $\kappa$ or $\ell$ large—the batch size becomes larger.

Theorem 2 extends the result of (Fehrman et al., 2020, Theorem 25) to a Markovian setting with inability to restart. In our case, the bias can be controlled by using a longer batch size with exponent at least $\sigma/2$. Furthermore, we also use the Lyapunov function to keep track of the state of the MDP as we update the parameter in $V$ and ensure stability. The proof sketch of Theorem 2 can be found in Section 5.6 and the full proof in Appendix C. In Appendix D, we also consider the case that $M \cap U$ is compact, which can be used to improve Theorem 2. Note that the sequence $(\Theta_m, m \in \mathbb{N})$ from Theorem 2 is well defined even if unstable policies occur, since the update $H_m$ from (16) is finite. In this case, recall that we have the convention that if $\theta$ yields an unstable policy then $J(\theta) = \inf_{s \in S, a \in A} r(s, a) \ge -\infty$. In Theorem 2, we can thus assume instead of initializing $\Theta_0 \in V$ that we restrict to trajectories of SGA that end up in the neighborhood $V$—the first iterate satisfying this being $\Theta_0$. In this alternative description, we can assume that the trajectory is $\{\tilde{\Theta}_t\}_{t \in [0, T+t_0]}$ for some $t_0 \in \mathbb{N}$, and $\Theta_0 = \tilde{\Theta}_{t_0}$ reaches $V$.

Theorem 2 also holds for any estimator $\tilde{H}_m$ of the gradient $J(\Theta_m)$ provided that this estimator satisfies appropriate bias and variance bounds typical for estimators using Markov chains (see Lemma 8 and Proposition 9 in Section 5.6 below). Thus, Theorem 9 and its consequences in the following sections hold for a wide range of policy-gradient methods. Similarly, Theorem 2 also holds when $M$ is a manifold of local maxima instead of global maxima. Indeed, the assumptions are all local and the proof is equivalent.

From Theorem 2, we immediately obtain a typical sample complexity bound.

**Corollary 3 (Sample Complexity)** Under the same assumptions and notation as in Theorem 2, there exists a constant $c > 0$ such that for any $1 > \epsilon > 0$ and $\delta > 0$, if we fix $\ell \ge \alpha^2/(5\delta c)$ and $\sigma + \kappa > 2$ then, for any $m \in \mathbb{N}$ satisfying

$$m \ge m(\epsilon, \delta) = c \max\left( (\epsilon^2 \delta)^{-\frac{1}{\sigma+\kappa}}, \delta^{-\frac{1}{\sigma+\kappa-1}}, \delta^{-\frac{2}{\kappa}}, \delta^{-\frac{1}{(\sigma+\kappa)/2-1}} \right), \quad (24)$$

we have

$$\mathbb{P}[J(\Theta_m) < J^* - \epsilon | \Theta_0 \in V] < \delta. \quad (25)$$

## 5.3 Lower Bound

As noted in Theorem 2, the rate in (23) includes the probability that the iterates escape $V$, outside which convergence cannot be guaranteed. Indeed, there is a term $O(\alpha^2/\ell)$ that characterizes the probability that the iterates escape the basin of attraction. For general settings, this term cannot be avoided, even in the unbiased case. In fact, the proposition below shows that for any $\beta > 0$ there are cases where there is a positive lower bound depending on $\alpha^{2+\beta}/\ell$. In Theorem 4 below, we consider an SGA setting with i.i.d. data, where the target is to maximize a function $f$ using estimators $H_m$ for the gradient $\nabla f(\Theta_m)$ at epoch $m$. In a non-RL setting, we usually have $H_m = H_m(\Theta_m, Z_m)$, where $Z_m$ is a collection of i.i.d. random variables and $\mathcal{F}_m$ denotes the sigma algebra of the random variables $\Theta_0, \dots, \Theta_m$ as well as $Z_0, \dots, Z_{m-1}$. For this result, we consider an RL setting where the iterates $\Theta_m$ satisfy (3), and $\eta_m = H_m - \nabla f(\Theta_m)$ satisfies the following unbiased
conditional concentration bounds for some $C > 0$:

$$ \mathbb{E}[\eta_m | \mathcal{F}_m] = 0 \quad \text{and} \quad \mathbb{E}[|\eta_m|^2 | \mathcal{F}_m] \le \frac{C}{t_{m+1} - t_m}. \qquad (26) $$

Proposition 4 below shows that Theorem 2 is almost sharp and characterizes the limitations of using a statistical estimator for the gradient, which can lead to instability. As we will see in Section 6.1, however, there are examples where only local convergence can be expected. The proof of Proposition 4 can be found in Appendix E.

**Proposition 4** For any $\beta > 0$, there are functions $f \in C^\infty(\mathbb{R}^n)$ with a maximum $f^* = f(\theta^*)$ satisfying Assumption 7, such that if the iterates $\Theta_m$ satisfy (3) and the gradient estimator $H_m = \nabla f(\Theta_m) + \eta_m$ satisfies (26), there exists a constant $c > 0$ depending on $f$ and independent of $m$ such that for any $\epsilon \in (0, 1)$, $1 > \alpha > 0$, $\delta > 0$, $\ell \ge 1$ and any $\sigma \ge 0, \kappa \ge 0$, in (17) we have that

$$ \mathbb{P}[f(\Theta_m) < f^* - \epsilon | \Theta_0 \in V] \ge c \frac{\alpha^{2+\beta}}{\ell} \quad \text{for any } m \ge 1. \qquad (27) $$

## 5.4 Performance Gap and Regret

A performance gap bound that we can obtain from Theorem 2 is not fully satisfactory for the epoch number $m$. Indeed, we can set the batch size very large ($\kappa$ large) since the cost of exploration is not factored in. We will therefore obtain a performance gap depending on the number of samples drawn. In particular, for a time-step $T \ge 1$, we will define the parameter at this time-step as $\Theta_{m(T)}$ where

$$ m(T) = \min\{n \in \mathbb{N} : \sum_{i=1}^{n} \ell i^{\sigma/2+\kappa} \ge T\}, \qquad (28) $$

that is, the corresponding epoch of the sample drawn at time $T$. We show the bounds on the performance gap in terms of the total number of samples $T$. The proof of Proposition 5 can be found in Appendix F.

**Proposition 5 (Performance Gap)** Under the same assumptions and notation as in Theorem 2, we fix $\alpha$ and $\delta$. Then for any $1 > \zeta > 0$ there is $\kappa(\zeta) \ge 0$, $c > 0$ and $\ell_0 > 0$ such that for any $\ell \ge \ell_0$ and $T \ge 1$ we have the following.

(i) If $\sup_{(s,a)} |r(s,a)| < \infty$, then

$$ \mathbb{E}[J^* - J(\Theta_{m(T)}) | \Theta_0 \in V] \le c \left( (\mathcal{L}^*)^{\frac{1}{3}} \ell^{\frac{1}{3}+\zeta} T^{-\frac{1}{3}+\zeta} + \ell^{\frac{1}{2}+\zeta} T^{-\frac{1}{2}+\zeta} + \ell^{\frac{2}{3}+\zeta} T^{-1+\zeta} + \frac{\alpha^2}{\ell} \right). \qquad (29) $$

(ii) If $\sup_{(s,a)} |r(s,a)|$ is unbounded, let $\mathcal{B}_{m(T)} = \{\Theta_n \in V, n \in [m(T)]\}$ be the event that all iterates up to the epoch of sample $T$ stay in $V$. Then, we have that

$$ \mathbb{P}[\mathcal{B}_{m(T)}] \ge 1 - c \left( \frac{\alpha^2}{\ell} + \ell^{1/2+\zeta} T^{-1/2-\zeta} + \ell^{2/3+\zeta} T^{-1-\zeta} \right), \qquad (30) $$

and

$$ \mathbb{E}[J^* - J(\Theta_{m(T)}) | \mathcal{B}_{m(T)}] \le c(\mathcal{L}^*)^{\frac{1}{3}} \ell^{\frac{1}{3}+\zeta} T^{-\frac{1}{3}+\zeta}. \qquad (31) $$
From Proposition 5 we obtain the regret of SAGE and, in general, of other policy-gradient algorithms that satisfy typical bounds on bias and variance bounds (see Lemma 8 and Proposition 9 in Section 5.6). The proof of Corollary 6 can be found in Appendix G.

**Corollary 6 (Regret)** Suppose the assumptions and notation as in Theorem 2 hold. Then for any $1 > \zeta > 0$ there exist $\kappa(\zeta) \ge 0$, $c > 0$, $\ell_0$ such that if $\ell \ge \ell_0$, when $\Theta_0$ is the first iterate of (3) in V

(i) If $\sup_{(s,a)} |r(s,a)| < \infty$, then for any $T > 1$,

$$ \mathbb{E}[TJ^* - \sum_{t=1}^{T} r(S_t, A_t) | \Theta_0 \in V] \le c \left( (\mathcal{L}^*)^{\frac{1}{3}} \ell^{\frac{1}{3}+\zeta} T^{\frac{2}{3}+\zeta} + \frac{\alpha^2}{\ell} T \right). \quad (32) $$

(ii) If $\sup_{(s,a)} |r(s,a)|$ is unbounded, then for any $T \ge 1$

$$ \mathbb{E}[TJ^* - \sum_{t=1}^{T} r(S_t, A_t) | \mathcal{B}_{m(T)}] \le c(\mathcal{L}^*)^{\frac{1}{3}} T^{\frac{2}{3}+\zeta}. \quad (33) $$

Besides the term $\alpha^2/\ell$ in Corollary 6 that captures the instability due to the estimation error of the gradient, the term $T^{2/3+\zeta}$ in (32) cannot be easily compared with other common regret bounds that assume global features for $J(\theta)$ or the gradient estimator $H_m$. However, the coefficient $(\mathcal{L}^*)^{\frac{1}{3}}$ plays an analogous role to the size of the state space in regret bounds for finite state space MDPs, and directly depends on the Lyapunov function. In Appendix B.1, we find explicitly the value of $\mathcal{L}^*$ for the single-server queue example of Section 6 and see that it behaves as $\mathcal{L}^* \sim O(\mathrm{Vol}(V))$, that is, it encodes the volume of the parameter space where the iterates are confined.

Remarkably, in the expectation of Corollary 6(i) policies that are unstable are not avoided. Indeed, note that we only condition on initializing in a stable policy in V but afterwards the trajectory may escape the set V and encounter unstable policies.

In (32), by setting $\ell = T^{1/4}$ we would obtain a horizon-dependent regret of $O(T^{3/4+\zeta})$, which is sublinear but far from the optimum $T^{1/2}$. This is most likely due to the decreasing step and increasing batch-sizes. While they allow for asymptotic convergence, they are slower to reach a fixed suboptimality gap compared to, e.g., using a constant step and batchsize algorithm. It may then be possible to use horizon-dependent step, and batch sizes together with a 'doubling trick' (Besson and Kaufmann, 2018) argument to achieve an anytime optimal suboptimality in the sublinear term. In this case, however, it is unclear if we would still obtain an equivalent factor $\alpha^2/\ell$ in (32) that is optimal as shown with Proposition 4.

## 5.5 Local Convergence with Entropy Regularization

A well-known phenomenon that can occur when using the softmax policy (6) is that, if the optimal policy is deterministic, the iterates converge to this optimal policy only when $\Theta_m \to \infty$. Problems where this occurs will thus not satisfy Assumption 7: the set of maxima will be empty. This phenomenon is illustrated in the example of Section 6.1. One prevalent method to mitigate the occurrence of maxima at the boundary involves incorporating a
regularization term, often linked to relative entropy KL[ ˜ π ‖ π ] of the policy π compared to a given ˜ π , defined below in (34).

Let ˜ π be a policy of the same type as those defined in (6) and let ζ be a distribution on S such that ζ ( h − $^{1}$( i )) > 0 for any i ∈ I , where h is the index map defined for the class of policies that we use in (6). We define the regularization term as

$$
\mathcal{R}_{\tilde{\pi}}(\theta) = \mathbb{E}_{S \sim \zeta}[\text{KL}[\tilde{\pi}(\cdot | S) \parallel \pi(\cdot | S, \theta)]] = \sum_{s \in S} \zeta(s) \mathbb{E}_{A \sim \tilde{\pi}(\cdot | s)} \left[ \log \left( \frac{\tilde{\pi}(A|s)}{\pi(A|s, \theta)} \right) \right]. \quad (34)
$$

For some b > 0 we define

$$
J_{\tilde{\pi}}(\theta) = J(\theta) - b \mathcal{R}_{\tilde{\pi}}(\theta). \tag{35}
$$

We can show that adding (34) to J ( θ ) defined in (7) not only prevents maxima from being at the boundary, but also allows us to avoid using Assumption 7 altogether. The next proposition is proved in Appendix H.

Proposition 7 Assume that we use the softmax policy from (6) and let J ( θ ) be defined as in (7). Then for almost every policy ˜ π in the class of (6) with respect to its Lebesgue measure,

1. the function J$_{˜}$ $_{π}$( θ ) in (35) satisfies Assumption 7 and the set of maximizers is bounded, and

2. Theorem 2 for J$_{˜}$ $_{π}$( θ ) holds without Assumption 7.

By using the regularization in (34) we are changing the original objective. Nevertheless, we can explicitly bound the difference between J and J$_{˜}$ π at their respective optima. Let J ∗ be the optimal value of J , and θ ∗ r be an optimum of J$_{˜}$ $_{π}$. For any θ we have the inequalities

$$
J_{\bar{\pi}}(\theta) = J(\theta) - b \mathcal{R}_{\bar{\pi}}(\theta) \leq J_{\bar{\pi}}(\theta_r^*) \leq J(\theta_r^*), \quad (36)
$$

so that rearranging and letting θ tend to θ ∗ with J ( θ $^{∗}$) = J ∗ we have

$$
J^* - J(\theta_r^*) \le b \mathbb{E}_{S \sim \zeta} [\mathrm{KL}[\tilde{\pi}(\cdot | S) \parallel \pi(\cdot | S, \theta^*)]]. \quad (37)
$$

5.6 Proof Outline for Theorem 2

We extend the local approach presented in (Fehrman et al., 2020, Section 5), that deals with convergence of SGD where the samples used to estimate the gradient are i.i.d. We consider instead an RL setting where data is Markovian and thus presents a bias. Fortunately, we can overcome its presence by adding an increasing batch size while tracking the states of the Markov chain via the local Lyapunov function from Assumption 4, which guarantees a stable MDP trajectory as long as the parameter is in a neighborhood close to the maximum. Below we give an outline of the technique employed. For the full proof we refer to Appendix C.
### 5.6.1 STRUCTURE OF THE PROOF

The proof of Theorems 2 consists of several parts. To show a bound on the probability that $\Theta_m$ is $\epsilon$-suboptimal, we consider the event $\mathcal{B}_m$ that all previous iterates $\Theta_0, \dots, \Theta_m$ belong to a local neighborhood $V$, and the complementary event $\bar{\mathcal{B}}_m$. We bound these separately. Firstly, on the event $\mathcal{B}_m$, we show in Lemma 10 that the iterates converge to $M$, and we obtain a bound on the $\epsilon$-suboptimal probability for this case. Secondly, we bound the probability of the complement $\bar{\mathcal{B}}_m$ by using a recursive identity relating the probability of $\bar{\mathcal{B}}_{m-1}$ and the sum of the probabilities of two disjoint events, namely, (i) $\Theta_m \notin V$ and the distance of $\Theta_m$ to $M$ is larger than $\delta$, and (ii) $\Theta_m \notin V$ and the distance of $\Theta_m$ to $M$ is less than $\delta$. Intuitively, these events group the cases when $\Theta_m$ escapes $V$ in ‘normal directions’ to $M$ and in ‘tangent directions’ to $M$, respectively. We can bound the former by using concentration inequalities, but for the latter we need a maximal excursion bound (Lemma 11 below). Combining all bounds results in an upper bound on $\mathbb{P}[\mathcal{B}_m]$ (Lemma 12). The local properties of $J$ are then used to complete the proof. Crucially, we use throughout the proof that the local Lyapunov function guarantees stability of the Markov chain and the gradient estimator within $V$, as well as keeps track of the initial state for each epoch. Geometric ergodicity also allows us to precisely quantify the bias and variance of the estimators. Note that if we assumed only stability instead of geometric ergodicity, only asymptotic control of the empirical estimators would be possible.

### 5.6.2 PRELIMINARY STEP: DEFINITION OF THE LOCAL NEIGHBORHOOD AND BOUND STRATEGY

For $\theta^* \in M \cap U$ and two positive numbers $r > 0$, and $\delta > 0$, we now define a neighborhood $V_{r,\delta}(\theta^*)$ of $\theta^*$ where the algorithm will eventually operate by choosing $\delta$ and $r$ appropriately. Let $\bar{B}_r(\theta^*) := \{\theta \in \Omega : |\theta - \theta^*| \le r\}$ denote a closed ball around $\theta^*$ with radius $r$ and $\mathrm{dist}(\theta, L) = \sup_{\theta' \in L} |\theta - \theta'|$ for an open set $L$. Let $U$ be the neighborhood of $\theta^*$ described in Assumptions 4 and 7. We define a tubular neighborhood of $\theta^*$ as follows:

$$V_{r,\delta}(\theta^*) := \{\theta \in \Omega \cap U : \mathrm{dist}(\theta, M \cap U) = \mathrm{dist}(\theta, \bar{B}_r(\theta^*) \cap M \cap U) < \delta\}. \quad (38)$$

Crucially, Assumption 7 implies that there exist $\delta_0, \mathbf{r}_0 > 0$ such that for any $\delta \in (0, \delta_0]$ and $\mathbf{r} \in (0, \mathbf{r}_0]$ an equivalent definition of the set is then

$$V_{r,\delta}(\theta^*) = \{y+v : y \in (\bar{B}_r(\theta^*) \cap M \cap U) \text{ and } v \in (T_y(M \cap U))^{\perp} \text{ with } |v| < \delta, \mathbf{p}(y+v) = y\}. \quad (39)$$

Here, $\mathbf{p}$ is the unique local projection onto $M \cap U$, and $T_y(M \cap U)^{\perp}$ denotes the cotangent space of $M \cap U$ at $y$. For further details on this geometric statement, we refer to (Fehrman et al., 2020, Proposition 13) or (Lee, 2013, Theorem 6.24).

In the following, we let $U$ denote the intersection of the neighborhoods from Assumptions 4 and 7, and $\mathcal{L}$ the Lyapunov function from Assumption 4. For any $m \in \mathbb{N}_+$ define the event and filtration

$$\mathcal{B}_m := \bigcap_{l=1}^{m} \{\Theta_l \in V_{r,\delta}(\theta^*)\}, \quad (40)$$

$$\mathcal{F}_m := \sigma(\mathcal{D}_1 \cup \dots \cup \mathcal{D}_{m-1} \cup \{\Theta_0, \dots, \Theta_m\}). \quad (41)$$
Due to the local properties of $J$, Theorem 2 can be shown by bounding $\mathbb{P}[\text{dist}(\Theta_m, \mathcal{M} \cap U) \ge \epsilon|\mathcal{B}_0]$. By separating into the event $\mathcal{B}_m$ and its complement, we can show that

$$ \mathbb{P}[\text{dist}(\Theta_m, \mathcal{M} \cap U) \ge \epsilon|\mathcal{B}_0] \le \mathbb{P}[\text{dist}(\Theta_m, \mathcal{M} \cap U)1[\mathcal{B}_{m-1}] \ge \epsilon] + \mathbb{P}[\overline{\mathcal{B}_m}]. \quad (42) $$

The remaining steps of the proof consist of bounding both terms in the right-hand side of (42).

### 5.6.3 STEP 1: THE VARIANCE OF THE GRADIENT ESTIMATOR DECREASES, IN SPITE OF THE BIAS

For each $m \in \mathbb{N}_+$, let

$$ \eta_m := H_m - \nabla J(\Theta_m), \quad (43) $$

denote the difference between the gradient estimator $H_m$ in (16) and the true gradient $\nabla J(\Theta_m)$. Lemma 8 below implies that the difference in (43) is, ultimately, small. From Assumption 4, since the state-action chain $\{(S_t, A_t)\}_{t \ge 0}$ has a Lyapunov function $\mathcal{L}$, so does the chain $\{S_t\}_{t>0}$ with

$$ \mathcal{L}_v(s) = \sum_{a \in \mathcal{A}} \mathcal{L}(s, a)^v \pi(a|s, \theta), \quad (44) $$

where $v \ge 16$ is the exponent from Assumption 4. We can define $\mathcal{L}_4(s)$ similarly. The following lemma bounds the variance of $\eta_m$ on the event $\mathcal{B}_m$, which can be controlled with the local Lyapunov function. The proof of Lemma 8 is deferred to Appendix C.3.

**Lemma 8** Suppose that Assumptions 1–7 hold. There exists a constant $C > 0$ that depends on $\theta^*$, $U$, and $J$ such that for every $m \in \mathbb{N}_+$,

$$ |\mathbb{E}[\eta_m 1[\mathcal{B}_m]|\mathcal{F}_m]| \le \frac{C}{t_{m+1} - t_m} \mathcal{L}_4(S_{t_m})^{1/2}, \quad (45) $$

$$ \mathbb{E}[|\eta_m|^l 1[\mathcal{B}_m]|\mathcal{F}_m] \le \frac{C}{(t_{m+1} - t_m)^{l/2}} \mathcal{L}_4(S_{t_m})^{l/2}, \quad \text{for every } l \in \{1, 2\}. \quad (46) $$

Lemma 8 helps to determine the bias incurred when starting at a different state than that of stationarity, and is used to bound the term $\text{dist}(\Theta_m, \mathcal{M} \cap U)$ from (42) in Lemma 10 below. Note that the definition of SAGE and Assumptions 5 and 6 are used.

As a matter of fact, any other estimator $\tilde{H}_m$ of $\nabla J$ satisfying (45) and (46) from Lemma 8 will yield similar guarantees. In particular, instead of using the aforementioned assumptions for the proof of Theorem 2 that involve the structure of SAGE or the stationary distribution, we may repeat the proof using instead Lemma 8.

**Proposition 9** Suppose Assumptions 1, 2, 4 and 7 hold. Suppose moreover that the estimator of the gradient $H_m$ satisfies Lemma 8 for each $m \ge 1$. Then, the results of Theorem 2, Section 5.4 and 5.5 also hold for the step-size from (17) and policy-gradient update in (3).
### 5.6.4 STEP 2: CONVERGENCE ON THE EVENT $\mathcal{B}_{m-1}$

We turn to the first term on the right-hand side of (42) and examine, on the event $\mathcal{B}_{m-1}$, if the iterates converge. Using a similar proof strategy as that of (Fehrman et al., 2020, Proposition 20) for the unbiased non-Markovian case, we prove in Lemma 10 that the variance of the distance to the set of minima decreases under the appropriate step and batch sizes. Moreover, the rate depends on an effective size of the state space, characterized by the Lyapunov function, and (22).

The proof of Lemma 10 is in Appendix C.4.

**Lemma 10** Suppose that Assumptions 1–7 hold. There then exist $\mathbf{r}_0$, $\alpha_0$, $\ell_0 > 0$ and $c > 0$ such that for any $\mathbf{r} \in (0, \mathbf{r}_0]$, $\alpha \in (0, \alpha_0]$ and $\ell \in [\ell_0, \infty)$ there also exists $\delta_0 > 0$ such that for any $\delta \in (0, \delta_0]$ and $m \in \mathbb{N}_+$,

$$ \mathbb{E}[(\mathrm{dist}(\Theta_m, \mathcal{M} \cap U) \wedge \delta)^2 \mathbf{1}[\mathcal{B}_{m-1}]] \le c\mathcal{L}^*m^{-\sigma-\kappa}. \quad (47) $$

Compared to the unbiased case of Fehrman et al. (2020), Lemma 10 needs to use a larger batch size to deal with the bias of Lemma 8. A key result required is that on the event $\mathcal{B}_{m-1}$, the Lyapunov function is bounded in expectation by $\mathcal{L}^*$, which captures the size of the ‘effective’ state space for the policies around an optimum. With Lemma 10 together with Markov’s inequality, a bound of order $\epsilon^{-2}m^{-\sigma-\kappa}$ for the first term in (42) follows.

### 5.6.5 STEP 3: EXCURSION AND THE PROBABILITY OF STAYING IN $V_{\mathbf{r},\delta}(\theta^*)$

We next focus on $\mathbb{P}[\bar{\mathcal{B}}_m]$. Since

$$ \mathbb{P}[\mathcal{B}_m] \ge \mathbb{P}[\mathcal{B}_{m-1}] - \mathbb{P}[\Theta_m \notin V_{\mathbf{r},\delta}(\theta^*), \mathcal{B}_{m-1}], \quad (48) $$

we can use a recursive argument to obtain a lower bound, if we can bound first the probability

$$ \mathbb{P}[\Theta_m \notin V_{\mathbf{r},\delta}(\theta^*), \mathcal{B}_{m-1}] = \mathbb{P}[\mathrm{dist}(\Theta_m, \mathcal{M} \cap U) > \delta, \mathcal{B}_{m-1}] \\ + \mathbb{P}[\mathrm{dist}(\Theta_m, \mathcal{M} \cap U) \le \delta, \Theta_m \notin V_{\mathbf{r},\delta}(\theta^*), \mathcal{B}_{m-1}]. \quad (49) $$

The first term in (49) represents the event that the iterand $\Theta_m$ escapes the set $V_{\mathbf{r},\delta}(\theta^*)$ in directions ‘normal’ to $\mathcal{M}$, while the second term represents the escape in directions ‘tangent’ to $\mathcal{M}$—intuition derived from the fact that, in that latter event, we still have $\mathrm{dist}(\Theta_m, \mathcal{M} \cap U) \le \delta$.

The first term in (49) can be bounded by using the local geometric properties around minima in the set $U$ and associating the escape probability with the probability that on the event $\mathcal{B}_{m-1}$ escape can only occur if $|\eta_m|$ is large enough. The probability of this last event happening can then be controlled with the variance estimates from Lemma 8.

After a recursive argument, we have to consider the second term in (49) for all $l \le m$. Fortunately, this term can be bounded by first looking at the maximal excursion event for the iterates $\{\Theta_l\}_{l=1}^m$. The proof can be found in Appendix C.5. Here, the Lyapunov function again plays a crucial role to control the variance of the gradient estimator on the events $\mathcal{B}_l$ for $l \le m$, compared to an unbiased and non-Markovian case.
**Lemma 11** Suppose that Assumptions 1-7 hold. Then there exist $\mathbf{r}_0$, $\alpha_0$, $\ell_0 > 0$, and $c > 0$ such that for any $\mathbf{r} \in (0, \mathbf{r}_0]$, $\alpha \in (0, \alpha_0]$ and $\ell \in [\ell_0, \infty)$, there exist $\delta_0 > 0$ such that for any $\delta \in (0, \delta_0]$ and $m \ge 1$,

$$ \mathbb{E}\left[\max_{1 \le l \le m} |\Theta_l - \Theta_0| \mathbf{1}[\mathcal{B}_{l-1}]\right] < c\alpha(m^{1-3\sigma/2-\kappa/2} + \sqrt{\frac{1}{\ell}}m^{1-5\sigma/8-\kappa/2}). \quad (50) $$

Finally, with the previous steps we obtain a bound on $\mathbb{P}[\mathcal{B}_m]$ in Lemma 12 below<fnref n="2" />. The proof of Lemma 12 can be found in Appendix C.6.

**Lemma 12** Suppose that Assumptions 1-7 hold and $\sigma + \kappa > 1$. There exist $\mathbf{r}_0$, $\alpha_0$, such that for any $\mathbf{r} \in (0, \mathbf{r}_0]$, $\alpha \in (0, \alpha_0]$, there also exists a constant $c > 0$, $\delta_0 > 0$ such that for any $\delta \in (0, \delta_0]$, if $\Theta_0 \in V_{\mathbf{r}/2, \delta}(\theta^*)$, there exists $\ell_0 > 0$ such that for any $\ell \in [\ell_0, \infty)$ and $m \in \mathbb{N}_+$,

$$ \mathbb{P}[\mathcal{B}_m] \ge \exp\left(-\frac{c\alpha^2}{\delta^2\ell}\right) - \frac{c}{\delta^4\ell}m^{1-\sigma-\kappa} - c\alpha\frac{\left(m^{1-3\sigma/2-\kappa/2} + \ell^{-1/2}m^{1-5\sigma/8-\kappa/2}\right)}{(\mathbf{r}/2 - 2\delta)_+}. \quad (51) $$

### 5.6.6 STEP 4: COMBINING THE BOUNDS IN (42)

The proof of Theorem 2 follows the same steps as are used to prove (Fehrman et al., 2020, Theorem 25) by substituting the modified bounds that we have obtained from Lemmas 10 and 12 in (42). The details can be found in Appendix C.2.

## 6 Examples and Numerical Results

In Section 5, we have shown convergence of a SAGE-based policy-gradient method under the assumptions in Section 5.1. We now numerically assess<fnref n="3" /> its performance in examples from stochastic networks and statistical physics that go beyond these assumptions. Specifically, we examine a single-server queue with admission control in Section 6.1, a load-balancing system in Section 6.2, and the Ising model with Glauber dynamics in Section 6.3. These examples satisfy the assumptions of Section 4, required to implement the SAGE-based policy-gradient method, but not necessarily those of Section 5.1 used to prove convergence; more discussion about these assumptions is provided in Appendix B. Section 6.1 can be seen as a toy warm-up example whose simple structure allows us to gain insight into the compared behavior of SAGE and actor-critic. The larger-scale examples of Sections 6.2 and 6.3 show the superiority of SAGE, in the sense that the dimension of these examples makes them out of reach for existing approaches without the use of function approximation.

The simulation setup is as follows. Plots are obtained by averaging 10 independent simulation runs, each lasting $T_{max} = 10^6$ time steps. The initial parameter vector $\Theta_0$ is taken to be the zero vector, yielding in each example a uniform policy over the action

<footnote n="2">2. In (50), if the Lyapunov function has only smaller moments than order $\nu$, then condition on $\kappa \ge 0$ will become stricter. In particular, $\kappa$ tunes the batch size required to sample from the tails of the stationary distribution and may be required to be positive depending the moments of the Lyapunov function. The terms $\sigma$ and $\kappa$ can be tuned to control the bias coming from variance and nonstationarity, and finite batch size, respectively.</footnote>

<footnote n="3">3. The code is available at https://gitlab.laas.fr/ccomte/policy-gradient-reinforcement-learning.</footnote>
space. The SAGE-based algorithm (Algorithm 2) is run with batch size 100 and step size $\alpha_m = 10^{-1}$. The actor-critic algorithm (Appendix A.1) is run with batch size 1 and step sizes $\alpha_m = 10^{-3}$ and $\alpha_v = \alpha_R = 10^{-2}$. It uses a tabular value function which we initially treat as containing all-zeros, and which in practice is expanded as states are visited.

## 6.1 Admission Control in a Single-Server Queue

Consider a queueing system where jobs arrive according to a Poisson process with rate $\lambda > 0$, service times are independent and exponentially distributed with rate $\mu > 0$, and the server applies an arbitrary nonidling nonanticipating scheduling policy such as first-come-first-served or processor-sharing. This model is also commonly known as an M/M/1 queue in the literature. When a job arrives, the agent decides to either admit or reject it: in the former case, the job is added to the queue, otherwise it is lost permanently. The agent receives a one-time reward $\gamma > 0$ for each admitted job (incentive to accept jobs) and incurs a holding cost $\eta > 0$ per job per time unit (incentive to reject jobs). The goal is to find an admission-control policy that achieves a trade-off between these two conflicting objectives.

The problem can be related to the framework of Section 3 as follows. For $t \in \mathbb{N}$, let $S_t$ denote the number of jobs in the system right before the arrival of the $(t+1)$th job, and let $A_t$ denote the decision of either admitting or rejecting this job. We have $S = \mathbb{N}$ and $\mathcal{A} = \{\text{admit}, \text{reject}\}$. Also, let $(\Sigma_{\tau}, \tau \in \mathbb{R}_{\ge 0})$ denote the continuous-time process that describes the evolution of the number of jobs over time and $(T_t, t \in \mathbb{N})$ the sequence of job arrival times, so that $S_0 = \Sigma_0$ and $S_t = \lim_{\sigma \uparrow T_t} \Sigma_{\tau}$ for $t \in \mathbb{N}_+$. Rewards are given by

$$R_{t+1} = r_{\text{disc}}(S_t, A_t) + \int_{T_t}^{T_{t+1}} r_{\text{cont}}(\Sigma_{\tau}) d\tau,$$

where $r_{\text{disc}}(s, a) = \gamma \mathbb{1}[a = \text{admit}]$ represents the one-time admission reward and $r_{\text{cont}}(s) = -\eta s$ the holding cost incurred continuously over time. We use this common reward structure in this example, but we remark that arbitrary reward functions $r_{\text{disc}}$ and $r_{\text{cont}}$ are possible.

For each $k \in \mathbb{N}$, we define a random policy parametrization $\pi_k$ with threshold $k$ and parameter vector$^4 \theta = (\theta_0, \theta_1, \dots, \theta_k) \in \mathbb{R}^{k+1}$ as follows. Under policy $\pi_k$, an incoming job finding $s$ jobs in the system is accepted with probability

$$\pi_k(\text{admit}|s, \theta) = \frac{1}{1 + e^{-\theta_{\min(s,k)}}}, \quad s \in \mathbb{N}. \qquad (52)$$

Taking $k=0$ yields a static (i.e., state-independent) random policy, while letting $k$ tend to infinity yields a fully state-dependent random policy. We believe this parametrization makes intuitive sense because, in a stable queueing system, small states tend to be visited more frequently than large states.

Under policy parametrization $\pi_k$, Assumptions 1 to 3 are satisfied with $n = d = k + 1$, $\Omega = \{\theta \in \mathbb{R}^{k+1} : \pi_k(\text{admit}|s, \theta) < \frac{\mu}{\lambda}\}$, $\Phi(s) = (\frac{\lambda}{\mu})^s$ for each $s \in S$, $x_i(s) = \mathbb{1}[s \ge i+1]$ for each $i \in \{0, 1, \dots, k-1\}$ and $x_k(s) = \max(s-k, 0)$, and $\rho_i(\theta) = \pi_k(\text{admit}|i, \theta)$ for each $i \in \{0, 1, \dots, k\}$. It follows that $\nabla \log \rho_i = \nabla \log \pi_k(\text{admit}|i, \cdot)$ for each $i \in \{0, 1, \dots, k\}$. We refer to Appendix B.1 for further details.

<footnote n="4">4. In this example, vectors and matrices are indexed starting at 0 (instead of 1) for notational convenience.</footnote>
6.1.1 NUMERICAL RESULTS IN A STABLE QUEUE

We study the impact of the policy threshold $k \in \mathbb{N}$ on the performance of SAGE and actor-critic. The parameters are $\lambda = 0.7$, $\mu = 1$, $\gamma = 5$, and $\eta = 1$, and we consider random policies $\pi_k$ with various thresholds. We have $\Omega = \mathbb{R}^{k+1}$ because $\lambda < \mu$, i.e., the queue is always stable. As we can verify using Appendix B.1, if $k \le 2$ the best policy is random, while if $k \ge 3$, the best policy (deterministically) admits incoming jobs if and only if there are at most 2 jobs in the system. Thus, if $k \ge 3$, the best policy is approximated when $\theta_i \to +\infty$ if $i \in \{0, 1, 2\}$ and $\theta_i \to -\infty$ if $i \in \{3, 4, \dots, k\}$. This deterministic policy is optimal among all Markovian policies. The initial policy is $\pi_k(\text{admit}|s, \Theta_0) = \frac{1}{2}$ for each $s \in \mathbb{N}$, and the system is initially empty, i.e., $S_0 = 0$ with probability 1.

Figure 1 depicts the impact of the threshold $k$ on the evolution of the long-run average reward $J(\Theta_t)$ (defined in (7)) and computed using the formulas of Appendix B.1) under SAGE and actor-critic. Figure 2 shows the admission probabilities $\pi_3(\text{admit}|i, \Theta_t)$ for each $i \in \{0, 1, 2, 3\}$ (i.e., the admission probabilities under the policy with threshold $k=3$). In both plots, the x-axis has a logarithmic scale starting at time $t = 10^2$, lines are obtained by averaging the results over 10 independent simulations, and transparent areas show the standard deviation. Both SAGE and actor-critic eventually converge to the maximal attainable long-run average reward, and under both algorithms the convergence is initially faster under policy $\pi_0$ than under $\pi_1$, $\pi_3$, $\pi_{100}$, and $\pi_{1000}$. For a particular threshold $k$, the convergence is initially faster under actor-critic than under SAGE. However, the long-run average reward under SAGE increases monotonically from its initial value to its maximal value while, under actor-critic, there is a time period (comprised between $10^3$ and $10^5$ time steps) where the long-run average reward stagnates or even decreases. Similar qualitative remarks can be made when looking at the running average reward $\frac{1}{t} \sum_{t'=1}^t R_{t'}$ instead of the long-run average reward $J(\Theta_t)$. Figure 2b suggests that, under $\pi_3$, this is because actor-critic first “overshoots” by increasing $\pi_3(\text{admit}|3, \Theta_t)$ too much and then decreasing $\pi_3(\text{admit}|2, \Theta_t)$ too much before eventually converging to the best admission probabilities. This overshooting is more pronounced with a small threshold $k$, but it is still visible with $k = 100$ and $k = 1000$.

Figures 1 and 2 suggest actor-critic has more difficulty to correctly estimate the policy update compared to SAGE, especially under parametrizations $\pi_k$ with small thresholds $k$. We conjecture this is due to the combination of two phenomena which reaches a peak when $k$ is small. First, a close examination of the evolution of the value function under $\pi_3$ and $\pi_{10}$ (not shown here) reveals that there is a transitory bias in the estimate of the value function. For instance, right after increasing the admission probability in state 0, the estimate of the value function at states 2 and 3 becomes negative, even if the optimal value function at these states is positive. Second, due to the policy parametrization, parameter $\theta_k$ is updated whenever a state $s \in \{k, k+1, k+2, \dots\}$ is visited (while, for each $i \in \{0, 1, \dots, k-1\}$, parameter $\theta_i$ is updated only when state $i$ is visited). As a result, the correlated biases in the estimates of the value function at states $k, k+1, k+2, \dots$ add up and lead actor-critic to overshoot the update of $\theta_k$, which has a knock-on effect on other states.
<figure><img src="image_2.png" /><figcaption>Figure 1: Long-run average reward $J(\Theta_t)$ in the admission-control problem with $\lambda = 0.7$, $\mu = 1$, $\gamma = 5$, and $\eta = 1$. Using Appendix B.1, we can verify that the long-run average reward under the best policy is approximately 2.183 if $k = 0$, 2.566 if $k = 1$, and 2.795 if $k \ge 3$.</figcaption></figure>

### 6.1.2 NUMERICAL RESULTS IN A POSSIBLY-UNSTABLE QUEUE

Figure 3 is the counterpart of Figure 1 when the arrival rate is $\lambda = 1.4 > 1 = \mu$. Now the set of policy parameters for which the system is stable is $\Omega = \{\theta \in \mathbb{R}^{k+1} : \pi_k(\text{admit}|k, \theta) < \frac{\mu}{\lambda}\} \subset \mathbb{R}^{k+1}$, with $\frac{\mu}{\lambda} \approx 0.714$. For simplicity, we will say that a policy is *stable* if the Markov chain defined by the system state under this policy is positive recurrent (i.e., if $\pi_k(\text{admit}|k, \theta) < \frac{\mu}{\lambda}$), and *unstable* otherwise. This is an example where convergence can only be guaranteed
<figure><img src="image_3.png" /><figcaption>Figure 2: Admission probabilities under policy parametrization $\pi_3$.</figcaption></figure>

locally, as not all policies are stable. Again using Appendix B.1, we can verify that if $k \le 1$, the best policy is random, while if $k \ge 2$, the best policy (deterministically) admits incoming jobs if and only if there are fewer than 2 jobs in the system. This deterministic policy is optimal among all Markovian policies. The initial policy is again the (stable) uniform policy, and the system is initially empty, i.e., $S_0 = 0$ with probability 1.

The first take-away of Figure 3 is that SAGE converges to a close-to-optimal policy despite the fact that some policies are unstable. The convergence of SAGE is actually faster under $\lambda = 1.4$ compared to $\lambda = 0.7$ (Figure 1). By looking at the evolution of the admission probability (not shown here), we conjecture this is due to the fact that the admission probability in states larger than or equal to 2 decreases much faster when $\lambda = 1.4$ compared to $\lambda = 0.7$, and that this probability has a significant impact on the long-run average reward. In none of the simulations does SAGE reach an unstable policy. This suggests that the updates of SAGE have lower chance of reaching unstable regions of the policy space per observed sample.

The second take-away of Figure 3 is that, on the contrary, actor-critic has difficulties coping with instability in this example. In all simulation runs used to plot this figure, the long-run average reward $J(\Theta_t)$ first decreases before possibly increasing again and converging to the best achievable long-run average reward. Under parametrizations $\pi_0, \pi_2$, and $\pi_4$, unstable policies are visited for thousands of steps in all simulation runs, and a stable policy is eventually reached in only 7 out of 10 runs. Under parametrization $\pi_0$, the long-run average reward under the last policy is close to the best only in 2 out of 10 runs. Under $\pi_{100}$ and $\pi_{1000}$, the policy remains stable throughout all runs, but the long-run average reward transitorily decreases before increasing again.

## 6.2 Load-Balancing System

Consider a cluster of $n$ servers. Jobs arrive according to a Poisson process with rate $\lambda > 0$, and a new job is admitted if and only if there are fewer than $c \in \mathbb{N}_+$ jobs in the system. Each server $i \in \{1, 2, \dots, n\}$ processes jobs in its queue according to a nonidling, nonanticipating policy. The service time of each job at server $i$ is exponentially distributed with rate $\mu_i > 0$, independently of all other random variables. The agent aims to maximize the admission probability by adequately distributing load across servers.
<figure><img src="image_4.png" /><figcaption>Figure 3: Long-run average reward in the admission-control problem with parameters $\lambda = 1.4$, $\mu = 1$, $\gamma = 5$, and $\eta = 1$. Using Appendix B.1, we can verify that the maximal value of the long-run average reward is approximately 1.091 if $k = 0$ and 1.880 if $k \ge 2$.</figcaption></figure>

For each $t \in \mathbb{N}$, let $S_t = (S_{t,1}, S_{t,2}, \dots, S_{t,n})$ denote the vector containing the number of jobs at each server right before the arrival of the $(t+1)$th job, and let $A_t \in \{1, 2, \dots, n\}$ denote the server to which this $(t+1)$th job is assigned. (This decision is void if $S_{t,1} + \dots + S_{t,n} = c$ because the job is rejected anyway.) We have $S = \{s \in \mathbb{N}^n : s_1 + s_2 + \dots + s_n \le c\}$ and $A = \{1, 2, \dots, n\}$. The agent obtains a reward of 1 if the job is accepted and 0 otherwise, that is, $R_{t+1} = 1[S_{t,1} + \dots + S_{t,n} \le c - 1]$ for each $t \in \mathbb{N}$.
We consider the following static policy parametrization, with parameter vector $\theta \in \mathbb{R}^n$: irrespective of the system state $s \in S$, an incoming job is assigned to server $i$ with probability

$$ \pi(i|s, \theta) = \pi(i|\theta) = \frac{e^{\theta_i}}{\sum_{j=1}^{n} e^{\theta_j}}, \quad i \in \{1, 2, \dots, n\}. \qquad (53) $$

Assumptions 1 to 3 are satisfied with $n = d$, $\Omega = \mathbb{R}^n$, $\Phi(s) = \prod_{i=1}^{n} (\frac{\lambda}{\mu_i})^{s_i}$ for each $s \in S$, $x_i(s) = s_i$ for $i \in \{1, 2, \dots, n\}$ and $s \in S$, and $\rho_i(\theta) = \pi(i|\theta)$ for $i \in \{1, 2, \dots, n\}$ and $\theta \in \mathbb{R}^n$. Also note that $\nabla \log \rho_i(\theta) = \nabla \log \pi(i|\theta)$ for $i \in \{1, 2, \dots, d\}$ and $\theta \in \mathbb{R}^n$. Except for Assumption 7, the remaining assumptions outlined in Section 5 are also satisfied. We refer to Appendix B.2 for more details. Lastly observe that, in spite of the policy being static and the state space being finite, the function $J$ is still nonconvex for typical system parameters. In fact, our numerical experiments are done in nonconvex scenarios. Furthermore, note that this system can become challenging to optimize if $c$ and $n$ are large.

### 6.2.1 NUMERICAL RESULTS

We study the performance of SAGE and actor-critic under varying numbers of servers and service speed imbalance. Given an integer $n \in \mathbb{N}_{>0}$ multiple of 4 and $\delta > 1$, we consider the following cluster of $n$ servers divided into 4 pools. For each $k \in \{1, 2, 3, 4\}$, pool $k$ consists of the $\frac{n}{4}$ servers indexed from $(k-1)\frac{n}{4} + 1$ to $k\frac{n}{4}$, and each server $i$ in this pool has service rate $\mu_i = \delta^{k-1}$. The total arrival rate is $\lambda = 0.7(\sum_{i=1}^{n} \mu_i)$ and the upper bound on the number of jobs in the system is $c = 10\frac{n}{4}$. Letting $\delta = 1$ gives a system where all servers have the same service speed, while increasing $\delta$ makes the server speeds more and more imbalanced. The initial policy is uniform, i.e., $\pi(i|\Theta_0) = \frac{1}{n}$ for each $i \in \{1, 2, \dots, n\}$, and the initial state is empty, i.e., $S_0 = 0$ with probability 1.

Figure 4 shows performance of SAGE and actor-critic in clusters of $n \in \{4, 20, 100\}$ servers. Solid lines show the evolution of the long-run average reward $J(\Theta_t)$, and dashed lines show the running average reward $\frac{1}{t}\sum_{t'=1}^{t} R_{t'}$. (Recall $J(\Theta_t)$ is the limit of the running average we would see if we ran the system under policy $\pi(\Theta_t)$.) It is defined in (7) and can be computed as shown in Appendix B.2.) As before, transparent areas show the standard deviation around the average. The results under actor-critic are reported only for $n = 4$ servers, as this method already suffers from a combinatorial explosion in the state-action space for $n \in \{20, 100\}$. Indeed, while the memory complexity increases linearly with the number $n$ of servers under SAGE, it increases with the cardinality $\binom{n+c}{c}$ of the state space under actor-critic$^5$, which is already prohibitively large for $n \in \{20, 100\}$.

All four subfigures in Figure 4 show a consistent 2-phase pattern: first the running average reward $\frac{1}{t}\sum_{t'=1}^{t} R_{t'}$ converges to the initial long-run average reward $J(\Theta_0)$, and then the long-run average reward increases to reach the best value, with the running average reward catching up at a slower pace. This suggests that the gradient estimates under both algorithms remain close to zero until the system reaches approximate stationarity. A similar reasoning explains why the algorithms converge at a slower pace when we increase the imbalance factor $\delta$ (as the stationary distribution under the initial uniform policy $\pi(\Theta_0)$ puts mass on states that are further away from the initial empty state) or the number $n$ of servers (as the mixing time increases).

<footnote n="5">5. As shown by applying the stars and bars method in combinatorics.</footnote>
<figure><img src="image_5.png" /><figcaption>Figure 4: Impact of the number of servers and service-rate imbalance on the performance of SAGE and actor-critic in a load-balancing system. Solid lines show the long-run average reward $J(\Theta_t)$, while dashed lines show the running average reward, $\frac{1}{t} \sum_{t'=1}^{t} R_{t'}$. Simulations for $n = 100$ and $\delta = 4$ are omitted because numerical instability of Buzen's algorithm (see Appendix B.2) prevents us from computing $J(\Theta_t)$ in this case.</figcaption></figure>

Focusing on the system with $n = 4$ servers, Figures 4a and 4b show convergence occurs approximately ten times faster under SAGE than under actor-critic. We conjecture that this is again due to the fact that actor-critic relies on estimating the state-value function, so that it needs to estimate $\left(\frac{n+c}{c}\right)$ values. SAGE, on the other hand, exploits the structure of the stationary distribution and only needs to estimate a number of values that grows linearly with $n$ (and is independent of $c$). We also note that actor-critic shows nonmonotonic convergence (i.e., $J(\Theta_t)$ decreases before increasing again between $10^5$ and $10^6$ time steps). We conjecture this is due to a similar phenomenon as described in Section 6.1.

## 6.3 Ising Model and Glauber Dynamics

Consider a system of spin particles spread over a two-dimensional lattice of shape $d_1 \times d_2$, for some $d_1, d_2 \in \{2, 3, 4, \dots\}$. Let $V = \{1, 2, \dots, d_1\} \times \{1, 2, \dots, d_2\}$ denote the set of lattice coordinates. For any two coordinates $v = (v_1, v_2) \in V$ and $w = (w_1, w_2) \in V$, we write $v \sim w$ if and only $v$ and $w$ are neighbors in the lattice, that is, if and only if $|v_1 - w_1| + |v_2 - w_2| = 1$.
A map $\sigma: \mathcal{V} \to \{-1, +1\}$ is called a spin configuration, and the set of all $2^{d_1d_2}$ configurations is denoted by $\Sigma$. Given a configuration $\sigma \in \Sigma$, we refer to $\sigma(v) \in \{-1, +1\}$ as the spin (of the particle located) at $v$. If the system is in some configuration $\sigma \in \Sigma$, we say that the spin at $v \in \mathcal{V}$ is flipped if the system jumps to the configuration $\sigma_{-v} \in \Sigma$ such that $\sigma_{-v}(v) = -\sigma(v)$ and $\sigma_{-v}(w) = \sigma(w)$ for each $w \in \mathcal{V} \setminus \{v\}$. As we will formalize below, the agent's goal is to reach a configuration $\sigma \in \Sigma$ so that the magnetization on the left (resp. right) half of the lattice is close to $\xi_{\text{left}} \in (-1, +1)$ (resp. $\xi_{\text{right}} \in (-1, +1)$), i.e.,

$$ \frac{2}{d_1d_2} M_{\text{left}}(\sigma) \approx \xi_{\text{left}}, $$

$$ \frac{2}{d_1d_2} M_{\text{right}}(\sigma) \approx \xi_{\text{right}}, $$

where

$$ M_{\text{left}}(\sigma) = \sum_{v \in \mathcal{V}: v_2 \le d_2/2} \sigma(v), $$

$$ M_{\text{right}}(\sigma) = \sum_{v \in \mathcal{V}: v_2 > d_2/2} \sigma(v). $$

To each configuration $\sigma \in \Sigma$ is associated an energy $E(\sigma) \triangleq -JI(\sigma) - \mu F(\sigma)$, where $I$ and $F$ are called the interaction and external field terms, respectively, given by

$$ I(\sigma) = \sum_{v,w \in \mathcal{V}: v \sim w} \sigma(v)\sigma(w), $$

$$ F(\sigma) = \sum_{v \in \mathcal{V}} h(v)\sigma(v), $$

where the first sum runs over all pairs of neighboring coordinates (so that each pair appears once). Here, $J \in \mathbb{R}$ is the coupling constant, $\mu \in \mathbb{R}_{\ge 0}$ the magnetic moment, and $h: \mathcal{V} \to \mathbb{R}$ the external magnetic field. Under the dynamics defined below, the probability of a configuration $\sigma \in \Sigma$ will be proportional to $e^{-\beta E(\sigma)}$, where $\beta \in \mathbb{R}_{\ge 0}$ is the inverse temperature. If $J > 0$ (resp. $J < 0$), the interaction term $I$ contributes to increasing the probability of configurations where neighboring spins have the same (resp. opposite) sign. Concurrently, due to the external-field term $F$, the spin at each $v \in \mathcal{V}$ is attracted in the direction pointed by the sign of $h(v)$. The coupling constant $J$ and magnetic moment $\mu$ are fixed and known by the agent (as they depend on the particles), and the agent will fine-tune the inverse temperature $\beta$ and coarse-tune the external magnetic field $h$.

### 6.3.1 GLAUBER DYNAMICS

Given a starting configuration, at every time step, the spin at a coordinate chosen uniformly at random is flipped (or not) with some probability that depends on the current configuration and the parameters set by the agent. This is cast as a Markov decision process as follows. The state and action spaces are given by $S = \Sigma \times \mathcal{V}$ and $A = \{\text{flip, not flip}\}$, respectively. For each $s = (\sigma, v) \in S$ and $a \in A$, the state reached by taking action $a$ in state $s$ is given by $S' = (\sigma', V')$, where $\sigma' = \sigma_{-v}$ if $a = \text{flip}$ and $\sigma' = \sigma$ if $a = \text{not flip}$, and $V'$ is chosen uniformly at random in $\mathcal{V}$, independently of the past states, actions, and rewards. The next reward $r$ is the opposite of the sum of the absolute difference between the next magnetizations and the desired magnetizations, that is,

$$ r = - \left| \xi_{\text{left}} - \frac{2}{d_1d_2} M_{\text{left}}(\sigma') \right| - \left| \xi_{\text{right}} - \frac{2}{d_2d_2} M_{\text{right}}(\sigma') \right|. $$
The agent controls a vector $\theta \in \mathbb{R}^3$ that determines the inverse temperature and the left and right external magnetic fields as follows:

$$ \beta(\theta) = 1 + \tanh(\theta_1), \quad h_{\text{left}}(\theta) = \tanh(\theta_2), \quad h_{\text{right}}(\theta) = \tanh(\theta_3), $$

so that in particular $\beta(\theta) \in (0, 2)$, $h_{\text{left}}(\theta) \in (-1, 1)$, and $h_{\text{right}}(\theta) \in (-1, 1)$. The corresponding external magnetic field and external field term are

$$ h(v|\theta) = h_{\text{left}}(\theta)1[v_2 \le d_2/2] + h_{\text{right}}(\theta)1[v_2 > d_2/2], \quad v \in \mathcal{V}, \\ F(\sigma|\theta) = \sum_{v \in \mathcal{V}} h(v|\theta)\sigma(v) = h_{\text{left}}(\theta)M_{\text{left}}(\sigma) + h_{\text{right}}(\theta)M_{\text{right}}(\sigma), \quad \sigma \in \Sigma. $$

Given $\theta \in \mathbb{R}^3$, for each $s = (\sigma, v) \in \Sigma$, the probability that the spin at the randomly-chosen coordinate $v$ is flipped when the current configuration is $\sigma$ is given by

$$ \pi(\text{flip}|s, \theta) = \frac{1}{1 + e^{\delta(s|\theta)}}, \quad \text{with} \quad \delta(s|\theta) = 2\beta(\theta)\sigma(v) \left( J \sum_{w \in \mathcal{V}: w \sim v} \sigma(w) + \mu h(v|\theta) \right). \quad (54) $$

When $\theta \in \mathbb{R}^3$ is fixed, the dynamics defined by this system are called the Glauber dynamics (Levin and Peres, 2017, Section 3.3). Note that, although we use the word *action* to match the terminology of MDPs, here an *action* should be seen as a random event in the environment, of which only the distribution $\pi$ can be controlled by the agent via the parameter vector $\theta$.

We verify in Appendix B.3 that the stationary distribution of the system state under a particular choice of $\theta \in \mathbb{R}^3$ satisfies

$$ p(s|\theta) \propto e^{\beta(\theta)(JI(\sigma) + \mu F(\sigma|\theta))}, \quad s = (\sigma, v) \in \mathcal{S}, \quad \theta \in \mathbb{R}^3. \quad (55) $$

Assumptions 1 to 3 are satisfied with $n = d = 3$, $\Omega = \mathbb{R}^3$, $\Phi(s) = 1$ for each $s \in \mathcal{S}$, $\log \rho_1(\theta) = \beta(\theta)J$, $\log \rho_2(\theta) = \beta(\theta)\mu h_{\text{left}}(\theta)$, and $\log \rho_3(\theta) = \beta(\theta)\mu h_{\text{right}}(\theta)$ for each $\theta \in \mathbb{R}^3$, and $x_1(s) = I(\sigma)$, $x_2(s) = M_{\text{left}}(\sigma)$, and $x_3(s) = M_{\text{right}}(\sigma)$ for each $s = (\sigma, v) \in \mathcal{S}$. All derivations are given in Appendix B.3.

### 6.3.2 NUMERICAL RESULTS

Figure 5 shows the performance of SAGE in a system with parameters $d_1 = 10$, $d_2 = 20$, $J = \mu = 1$, $\xi_{\text{left}} = -1$, and $\xi_{\text{right}} = 1$. We do not run simulations under actor-critic, as again the state space has size $2^{d_1d_2} = 2^{200}$, which is out of reach for this method. The initial parameter vector is $\Theta_0 = 0$, yielding inverse temperature $\beta(\Theta_0) = 1$ and external fields $h_{\text{left}}(\Theta_0) = h_{\text{right}}(\Theta_0) = 0$. The initial configuration has spins 1 on the left-hand side and -1 on the right-hand side, so that reaching the target configuration requires flipping every spin. In Figure 5a, the reward $R_t$ seems to increase on average monotonically from -4 to 0, which is consistent with the observation that the left (resp. right) magnetization decreases from 1 to -1 (resp. increases from -1 to 1). The increase of the reward is stepwise, with stages where it remains roughly constant for several thousand time steps. Lastly, the standard deviation increases significantly from about $10^4$ to $3 \cdot 10^5$ time steps, and it becomes negligible afterwards.
<figure><img src="image_6.png" /><figcaption>Figure 5: Performance of SAGE in the Ising model.</figcaption></figure>

To help us understand these observations, Figure 5b shows the evolution of the system parameters and of the magnetizations over a particular simulation run. The left magnetization $M_{\text{left}}$ starts decreasing around $10^4$ time steps (bottom plot), approximately when $h_{\text{left}}(\Theta_t)$ and $h_{\text{right}}(\Theta_t)$ become nonzero (top plot), to become roughly -1 around $3 \cdot 10^4$ time steps. At that moment, the system configuration is close to the all-spin-down configuration $\sigma_{-1}$ such that $\sigma_{-1}(v) = -1$ for each $v \in \mathcal{V}$. The right magnetization starts increasing significantly only when the inverse temperature $\beta(\Theta_t)$ has a sudden decrease (top plot). To make sense of this observation, consider $\pi(\text{flip}|s, \theta)$ as given by (54), where $s = (\sigma_{-1}, v)$ for some $v \in \{2, 3, \dots, d_1-1\} \times \{2, 3, \dots, d_2-1\}$. In $\delta(s|\theta)$, the first term $J \sum_{w \in \mathcal{V}: w \sim v} \sigma_{-1}(v)\sigma_{-1}(w)$ is equal to 4, while the absolute value of the second term $\mu\sigma_{-1}(v)h(v|\theta)$ is at most 1; hence, if $\beta(\theta) \simeq 1$ as initially, $\pi(\text{flip}|s, \theta)$ is between $\frac{1}{1+e^{2(4+1)}} \simeq 4.5 \cdot 10^{-5}$ and $\frac{1}{1+e^{2(4-1)}} \simeq 2.5 \cdot 10^{-3}$. The brief decrease of $\beta(\theta)$ is an efficient way of increasing the flipping probability in all states, which allows the system to escape from $\sigma_{-1}$. Other simulation runs are qualitatively similar, but the times at which the qualitative changes occur and the side (left or right) that flips magnetization first vary, which explains the large standard deviation observed earlier.

# 7 Conclusion

In this paper, we incorporated model-specific information about MDPs into the gradient estimator in policy-gradient methods. Specifically, assuming that the stationary distribution is an exponential family, we derived score-aware gradient estimators (SAGEs) that do not require the computation of value functions (Theorem 1). As showcased in Section 6, this
assumption is satisfied by models from stochastic networks, where the stationary distribu-
tion possesses a product-form structure, and by models from statistical mechanics, such as
the Ising model with Glauber dynamics.

The numerical results in Section 6 show that in these systems, policy-gradient algorithms equipped with a SAGE outperform actor-critic. In these examples, the Jacobian of the load function D log ρ(θ) can be computed explicitly in terms of the policy parameter θ. However, SAGE estimators can be harder to compute in more complex cases, for example when D log ρ(θ) depends on some model parameters. Nevertheless, our examples showcase how it is possible to improve the current policy gradient methods by levering information on the MDP, and we expect extensions of SAGEs to cover more challenging cases, for example by combining SAGE with model selection by first estimating the model parameters appearing in D log ρ(θ). We leave such extensions of SAGE for future work.

We have also shown with Theorem 2 that policy gradient with SAGE converges to the
optimal policy under light assumptions, namely, the existence of a local Lyapunov function
close to the optimum, which allows for unstable policies to exist, and a nondegeneracy
property of the Hessian at maxima. The convergence occurs with a probability arbitrarily
close to one provided that the iterates start close enough. Notably, our method proof also
works with other policy-gradients with similar policy-gradient approximation properties.
In Corollary 6, the regret of the algorithm is shown to be $O(T^{2/3+\epsilon} + T\alpha^2/\ell)$, where T is
the number of samples drawn. Unlike most common convergence results, we have gradient
estimation on a countable state space and there is a nonzero probability that such error will
drive the algorithm to an unstable policy. This fact is namely captured by the term $\alpha^2/\ell$.
Remarkably, such instabilities are observed in one of the examples of Section 6. If we had
made stronger assumptions such as the existence of a global Lyapunov function, then such
phenomena would not have been captured by the analysis.
# Appendix A. Policy-Gradient Algorithms

In this section, we give additional details about the two policy-gradient algorithms that are evaluated numerically in Section 6: actor-critic in Appendix A.1, and SAGE-based methods in Appendix A.2.

## A.1 Actor-Critic Algorithm

The actor-critic algorithm is first mentioned in Section 3.4 and compared to our SAGE-based policy-gradient algorithm in Section 6. We focus on the version of actor-critic described in (Sutton and Barto, 2018, Section 13.6) for the average-reward criterion in infinite horizon. The algorithm relies on the following expression for $\nabla J(\theta)$, which is a variant of the policy-gradient theorem (Sutton and Barto, 2018, Chapter 13):

$$ \nabla J(\theta) \propto \mathbb{E}[(R - J(\theta) + v(S') - v(S))\nabla \log \pi(A|S, \theta)], $$

where $(S, A, R, S')$ is a quadruplet of random variables such that $S \sim p(\cdot|\theta)$, $A|S \sim \pi(\cdot|S, \theta)$, and $(R, S')|(S, A) \sim P(\cdot, \cdot|S, A)$ (so that in particular $(S, A, R) \sim \text{STAT}(\theta))$, and $v$ is the state-value function.

**Algorithm 3** Actor-critic algorithm (Sutton and Barto, 2018, Section 13.6) to be called on Line 9 of Algorithm 1, with batch sizes equal to one.

1: **Input:** Positive and differentiable policy parametrization $(s, \theta, a) \mapsto \pi(a|s, \theta)$

2: **Parameters:** Step sizes $\alpha_R > 0$ and $\alpha_v > 0$

3: **Initialization:** $\bullet \bar{R} \leftarrow 0$
$\bullet V[s] \leftarrow 0$ for each $s \in S$

4: **procedure GRADIENT(t)**

5: $\delta \leftarrow R_{t+1} - \bar{R} + V[S_{t+1}] - V[S_t]$

6: Update $\bar{R} \leftarrow \bar{R} + \alpha_R \delta$

7: Update $V[S_t] \leftarrow V[S_t] + \alpha_v \delta$

8: **return** $\delta \nabla \log \pi(A_t|S_t, \Theta_t)$

9: **end procedure**

The pseudocode of the procedure GRADIENT used in the actor-critic algorithm is given in Algorithm 3. This procedure is to be implemented within Algorithm 1 with batch sizes equal to one, meaning that $t_{m+1} = t_m + 1$ for each $m \in \mathbb{N}$. We assume for simplicity that all variables from Algorithm 1 are accessible inside Algorithm 3. The variable $\bar{R}$ updated on Line 6 is a biased estimate of $J(\Theta_m)$, while the table $V$ updated on Line 7 is a biased estimate of the state-value function under policy $\pi(\Theta_m)$. Compared to (Sutton and Barto, 2018, Section 13.6), the value function is encoded by a table $V$ and there are no eligibility traces. If the state space $S$ is infinite, the table $V$ is initialized at zero over a finite subset of $S$ containing the initial state $S_0$ and expanded with zero padding whenever necessary.

## A.2 SAGE-Based Policy-Gradient Method

Algorithm 4 is an extension of Algorithm 2 that allows for batches of size 1. The main advantage of Algorithm 4 over Algorithm 2 is that it estimates $\nabla J(\Theta_m)$ based not only on
**Algorithm 4** SAGE-updated policy-gradient method, to be called on Line 9 of Algorithm 1.

1: **Input:** • Positive and differentiable policy parametrization $(s, \theta, a) \mapsto \pi(a|s, \theta)$

• Jacobian matrix function $\theta \mapsto D \log \rho(\theta)$

• Feature function $s \mapsto x(s)$

2: **Parameters:** Memory factor $\nu \in [0, 1]$

3: **Initialization:** Global variables $N_{-1}, M_{-1}, \bar{X}_{-1}, \bar{R}_{-1}, \bar{C}_{-1}, \bar{E}_{-1} \leftarrow 0$

4: **procedure GRADIENT**($m$)

5: $N_m \leftarrow \nu N_{m-1} + (t_{m+1} - t_m)$

6: $M_m \leftarrow \nu^2 M_{m-1} + (t_{m+1} - t_m)$

7: **return** $D \log \rho(\Theta_m) \text{COVARIANCE}(m) + \text{EXPECTATION}(m)$

8: **end procedure**

9: **procedure COVARIANCE**($m$)

10: Update $\bar{X}_m \leftarrow \nu \bar{X}_{m-1} + \sum_{t=t_m}^{t_{m+1}-1} x(S_t)$

11: Update $\bar{R}_m \leftarrow \nu \bar{R}_{m-1} + \sum_{t=t_m}^{t_{m+1}-1} R_{t+1}$

12: Update $\bar{C}_m \leftarrow \nu \bar{C}_{m-1} + \sum_{t=t_m}^{t_{m+1}-1} (x(S_t) - \frac{1}{N_m} \bar{X}_m)(R_{t+1} - \frac{1}{N_m} \bar{R}_m)$

13: **return** $\frac{N_m}{N_m^2 - M_m} \bar{C}_m$ if $N_m^2 > M_m$ else $\frac{1}{N_m} \bar{C}_m$

14: **end procedure**

15: **procedure EXPECTATION**($m$)

16: Update $\bar{E}_m \leftarrow \nu \bar{E}_{m-1} + \sum_{t=t_m}^{t_{m+1}-1} R_{t+1} \nabla \log \pi(A_t|S_t, \Theta_m)$

17: **return** $\frac{1}{N_m} \bar{E}_m$

18: **end procedure**

batch $\mathcal{D}_m$, but also on previous batches, depending on the memory factor $\nu$ initialized on Line 2. To simplify the signature of procedures in Algorithm 4, we assume variables $N_m$, $M_m$, $\bar{X}_m$, $\bar{R}_m$, $\bar{C}_m$, and $\bar{E}_m$ are global, and that all variables from Algorithm 1 are accessible within Algorithm 4, in particular batch $\mathcal{D}_m$. Line 5 in Algorithm 4 is the counterpart of Line 3 in Algorithm 2. Line 6 in Algorithm 4 is used on Line 13 to compute the counterpart of the quotient $1/(N_m - 1)$ in Line 6 in Algorithm 2. The COVARIANCE($m$) procedure in Algorithm 4 is the counterpart of Lines 4–6 in Algorithm 2. The EXPECTATION($m$) procedure in Algorithm 4 is the counterpart of Line 7 in Algorithm 2. Algorithm 2 can be seen as a special case of Algorithm 4 with memory factor $\nu = 0$. Note that terminology in Algorithm 4 differs slightly compared to Algorithm 2: bar notation refers to cumulative sums instead of averages.

The subroutines COVARIANCE and EXPECTATION compute biased covariance and mean estimates for Cov[$R, x(S)$] and $E[$R $\nabla \log \pi(A|S, \theta)$], where $(S, A, R) \sim \text{STAT}(\Theta_m)$, consistently with Theorem 1. If the memory factor $\nu$ is zero, these procedures return the usual sample mean and covariance estimates taken over the last batch $\mathcal{D}_m$ (as in Algorithm 2), and bias only comes from the fact that the system is not stationary. If $\nu$ is positive, estimates from previous batches are also taken into account, so that the bias is increased in exchange for a (hopefully) lower variance. In this case, the updates on Lines 10–12 and 16 calculate iteratively the weighted sample mean and covariance over the whole history, where observations from epoch $m - \underline{m}$ have weight $\nu^{\underline{m}}$, for each $\underline{m} \in \{0, 1, \dots, m\}$. When $m$ is large, the
mean returned by EXPECTATION is approximately equal to the sample mean over batches $\mathcal{D}_{m-M}$ through $\mathcal{D}_m$, where $M$ is a truncated geometric random variable, independent of all other random variables, such that $\mathbb{P}[M = \underline{m}] \propto \nu^m$ for each $\underline{m} \in \{0, 1, \dots, m\}$; if batches have constant size $c$, then we take into account approximately the last $c(\mathbb{E}[M] + 1) = \frac{c}{1-\nu}$ steps.

# Appendix B. Examples

This appendix provides detailed derivations for the examples of Section 6. We consider the single-server queue with admission control of Section 6.1 in Appendix B.1, the load-balancing example of Section 6.2 in Appendix B.2, and the Ising model of Section 6.3 in Appendix B.3.

## B.1 Single-Server Queue with Admission Control

Consider the example of Section 6.1, where jobs arrive according to a Poisson process with rate $\lambda > 0$, and service times are exponentially distributed with rate $\mu > 0$. Recall the long-run average reward is the difference between an admission reward proportional to the admission probability and a holding cost proportional to the mean queue size. We first verify that Assumptions 1 to 3 are satisfied, then we give a closed-form expression for the objective function, and lastly we discuss the assumptions of Section 5. We consider a random threshold-based policy $\pi_k$ of the form (52) for some $k \in \mathbb{N}$ and some parameter $\theta \in \Omega$, where $\Omega = \{\theta \in \mathbb{R}^{k+1} : \pi_k(\text{admit}|k, \theta) < \frac{\mu}{\lambda}\}$.

### B.1.1 PRODUCT-FORM STATIONARY DISTRIBUTION

The evolution of the number of jobs in the system (either waiting or in service) defines a birth-and-death process with birth rate $\lambda\pi_k(\text{admit}|s, \theta)$ and death rate $\mu\mathbb{1}[s \ge 1]$ in state $s$, for each $s \in \{0, 1, 2, \dots\}$. This process is irreducible because its birth and death rates are positive, and it is positive recurrent because $\lambda\pi_k(\text{admit}|s, \theta) < \mu$ for each $s \in S$ by definition of $\Omega$. This verifies Assumption 1. The stationary distribution is given by

$$ 
\begin{aligned}
p(s|\theta) &= \frac{1}{Z(\theta)} \prod_{q=0}^{s-1} \left(\frac{\lambda}{\mu}\pi(\text{admit}|q, \theta)\right), \\
&= \frac{1}{Z(\theta)} \left(\frac{\lambda}{\mu}\right)^s \left[\prod_{i=0}^{k-1} \pi_k(\text{admit}|i, \theta)^{\mathbb{1}[s \ge i+1]}\right] \pi_k(\text{admit}|k, \theta)^{\max(s-k, 0)}, \quad s \in \mathbb{N}, \quad (56)
\end{aligned}
$$

where the second equality follows by injecting (52), and the value of $Z(\theta)$ follows by normalization. We recognize (10-PF) from Assumption 3, with $n = d = k + 1$, $\Phi(s) = (\frac{\lambda}{\mu})^s$ for each $s \in S$, $x_i(s) = \mathbb{1}[s \ge i + 1]$ for each $i \in \{0, 1, \dots, k-1\}$ and $x_k(s) = \max(s-k, 0)$ for each $s \in S$, and $\rho_i(\theta) = \pi_k(\text{admit}|i, \theta)$ for each $i \in \{0, 1, \dots, k\}$. The function $\rho$ defined in this way is differentiable. Assumption 3 is therefore satisfied, as the distribution of the system seen at arrival times is also (56) according to the PASTA property (Wolff, 1982). For each $s \in \mathbb{N}$ and $a \in \{\text{admit}, \text{reject}\}$, $\nabla \log \pi_k(a|s, \theta)$ is the $(k+1)$-dimensional column vector with value $\mathbb{1}[a = \text{admit}] - \pi_k(\text{admit}|a, \theta)$ in component $i = \min(s, k)$ and zero elsewhere, and $\text{D} \log \rho(\theta)$ is the $(k+1)$-dimensional diagonal matrix with diagonal coefficient
$1 - \pi_k(\text{admit}|i, \theta)$ in position $i$, for each $i \in \{0, 1, \dots, k\}$. This can be used to verify that Assumption 5 is satisfied.

## B.1.2 OBJECTIVE FUNCTION

The objective function is $J(\theta) = \gamma \mathbb{P}[A = \text{admit}] - \eta \lambda \mathbb{E}[S]$, where

$$ \begin{aligned} \mathbb{P}[A = \text{admit}] &= \sum_{i=0}^{k-1} p(i|\theta)\pi_k(\text{admit}|i, \theta) + \left(1 - \sum_{i=0}^{k-1} p(i|\theta)\right)\pi_k(\text{admit}|k, \theta), \\ \mathbb{E}[S] &= \sum_{i=0}^{k-1} ip(i|\theta) + \frac{p(k|\theta)}{1 - \rho_k(\theta)} \left(k + \frac{\frac{\lambda}{\mu}\rho_k(\theta)}{1 - \frac{\lambda}{\mu}\rho_k(\theta)}\right), \\ Z(\theta) &= \sum_{s=0}^{k-1} \left(\prod_{i=0}^{s-1} \frac{\lambda}{\mu}\rho_i(\theta)\right) + \left(\prod_{i=0}^{k-1} \frac{\lambda}{\mu}\rho_i(\theta)\right) \frac{1}{1 - \frac{\lambda}{\mu}\rho_k(\theta)}, \end{aligned} $$

with the convention that empty sums are equal to zero and empty products are equal to one. All calculations remain valid in the limit as $\pi_k(\text{admit}|i, \theta) \to 1$ for some $i \in \{0, 1, \dots, k\}$ (corresponding to $\theta_i \to +\infty$). In the limit as $\pi_k(\text{admit}|i, \theta) \to 0$ for some $i \in \{0, 1, \dots, k\}$, we can study the restriction of the birth-and-death process to the state space $\{0, 1, \dots, c\}$, where $c = \min\{i \in \{0, 1, \dots, k\} : \pi_i(\theta) = 0\}$.

## B.1.3 ASSUMPTIONS OF SECTION 5

For any closed set $U \subset \Omega$, it can be shown that there exists a Lyapunov function $\mathcal{L}$ uniformly over $\theta \in U$ such that $\mathcal{L}(s, a) = \mathcal{L}(s) = \exp/cs$ for some $c > 0$, depending on $U$ and the model parameters. We look at the equivalent geometric ergodicity condition for continuous Markov chains. If $\theta \in U$ we have $\mu - \lambda\pi_k(\theta) > \delta(U) > 0$. Then, for $s > k+1$, the generator of the Markov process $Q_\theta$ satisfies

$$ \begin{aligned} Q_{\theta}\mathcal{L}(s) &= q_{\theta}(s-1|s)\mathcal{L}(s-1) + q_{\theta}(s+1|s)\mathcal{L}(s+1) + q_{\theta}(s|s)\mathcal{L}(s) \\ &= \mu \exp(c(s-1)) + \lambda\pi_k(\theta) \exp(c(s+1)) - (\mu + \lambda\pi_k(\theta)) \exp/cs \\ &= \exp/cs)(\mu \exp(-c) + \lambda\pi_k(\theta)) \exp(c) + (\mu + \lambda\pi_k(\theta)) \\ &= -c(\mu - \lambda\pi_k(\theta) + O(c))\mathcal{L}(s). \end{aligned} \quad (57) $$

For $c$ small enough, from (58) we have that $Q_\theta\mathcal{L}(s) \le -c\delta(U)/2\mathcal{L}(s)$, so that for any $\theta \in U$ the Markov chain corresponding to the policy of $\theta$ is geometrically ergodic. Hence, Assumptions 4, 5 and 6 are satisfied. In general, Assumption 7 does not hold for this example because maxima occur only as $|\theta| \to \infty$. As suggested by Proposition 7, by adding a small regularization term, we can guarantee Assumption 7 while simultaneously ensuring that the maximizer is bounded. In practice, using a regularization term can additionally present some benefits such as avoiding vanishing gradients and saddle points.

## B.1.4 EFFECTIVE STATE SPACE

The effective state space if captured in the term (22). Similarly to the continuous-time Markov chain example from (58), we have that the Lyapunov function is $\mathcal{L}(s, a) = \exp/cs$
for some $c > 0$ small enough. Let $\rho$ be such that $1 > \rho \ge \lambda\pi_k(\theta)/\mu$ for any $\theta \in V$. Let $c$ be small enough such that

$$ \frac{\rho}{1 + \rho} \exp(c) + \frac{1}{1 + \rho} \exp(-c) \le (1 - c \frac{1 - \rho}{2(1 + \rho)}), \quad (59) $$

Then, for any $s \ge k$,

$$ P_{\theta}\mathcal{L}(s) \le \lambda\mathcal{L}(s) + b, \quad (60) $$

where

$$ \lambda = 1 - c \frac{1 - \rho}{2(1 + \rho)}, \quad (61) $$

$$ b = \exp(c_1 k) \quad \text{for some} \quad c_1 > 0. \quad (62) $$

If we let $s_0 \in [k]$, then from the definition of $\mathcal{L}$ in (22) we will have that

$$ \mathcal{L}^* = O(\exp(ck)). \quad (63) $$

We note that in this case $\mathcal{L}^* \sim \text{Volume}(U) \sim \exp(\dim(\theta))$. Hence, it encodes the volume of the optimization space where the algorithm operates. We remark that the Lyapunov function encodes geometric ergodicity and allows to tackle any type of rewards, as long as they satisfy $|r|_{\mathcal{L}} < \infty$. Thus, for a specific reward $r$ better bounds could be attained.

We remark that geometric ergodicity is not equivalent to Foster stability, which is guaranteed by $P_{\theta}\mathcal{L}(s) \le \mathcal{L}(s) - \delta$ for some $\delta > 0$ and all $s \in S$. Foster stability implies stability (i.e., positive recurrence) of the Markov chain, and the existence of $\mathbb{E}[\mathcal{L}(s)]$. In this case, a Lyapunov function of the type $\mathcal{L}(s) \simeq s^2$ would suffice to show positive recurrence.

## B.2 Load-Balancing System

We now consider the load-balancing example of Section 6.2. Recall that jobs arrive according to a Poisson process with rate $\lambda > 0$, there are $n$ servers at which service times are distributed exponentially with rates $\mu_1, \mu_2, \dots, \mu_n$, respectively, and the system can contain at most $c$ jobs, for some $c \in \mathbb{N}_+$. The goal is to choose a static random policy that maximizes the admission probability. We first verify that the system satisfies Assumptions 1 to 3, then we provide an algorithm to evaluate the objective function when the parameters are known; this is used in particular for performance comparison with the optimal policy in the numerical results. Lastly, we discuss the assumptions of Section 5. Throughout this section, we assume that we apply the policy $\pi(\theta)$ defined by (53) for some parameter $\theta \in \mathbb{R}^n$.

### B.2.1 PRODUCT-FORM STATIONARY DISTRIBUTION

That Assumption 1 is satisfied follows from the facts that the rates and probabilities $\lambda, \mu_1, \mu_2, \dots, \mu_n, \pi_1(\theta), \pi_2(\theta), \dots, \pi_n(\theta)$ are positive and that the state space $S$ is finite. Assumption 2 is satisfied because the state space is finite. This system can be modeled either as a loss Jackson network with $n$ queues (one queue for each server in the load-balancing system) or as a closed Jackson network with $n+1$ queues (one queue for each server in the system, plus another queue signaling available positions in the system, with service rate $\lambda$).
Either way, we can verify (for instance by writing the balance equations) that the stationary distribution of the continuous-time Markov chain that describes the evolution of the system state is given by:

$$p(s|\theta) = \frac{1}{Z(\theta)} \prod_{i=1}^{n} \left( \frac{\lambda}{\mu_i} \pi_i(\theta) \right)^{s_i}, \quad s = (s_1, s_2, \dots, s_n) \in \mathcal{S}, \qquad (64)$$

where $Z(\theta)$ follows by normalization. This is exactly (10-PF) from Assumption 3, with $n = d$, $\Omega = \mathbb{R}^n$, $\Phi(s) = \prod_{i=1}^n (\frac{\lambda}{\mu_i})^{s_i}$ for each $s \in \mathcal{S}$, $x_i(s) = s_i$ for each $i \in \{1, 2, \dots, n\}$ and $s \in \mathcal{S}$, and $\rho_i(\theta) = \pi_i(\theta)$ for each $i \in \{1, 2, \dots, n\}$. The function $\rho$ defined in this way is differentiable. Assumption 3 is therefore satisfied, as the distribution of the system seen at arrival times is also (64) according to the PASTA property. Besides the sufficient statistics $x$, the inputs of Algorithm 2 are $\nabla \log \pi(a|s, \theta) = \mathbb{1}_a - \pi(\theta)$, where $\mathbb{1}_a$ is the $n$-dimensional vector with one in component $a$ and zero elsewhere and $\pi(\theta)$ is the policy seen as a (column) vector, and $\mathrm{D} \log \rho(\theta) = \mathrm{Id} - \mathbb{1}\pi(\theta)^{\mathrm{T}}$, where $\mathrm{Id}$ is the $n$-dimensional identity matrix, $\mathbb{1}$ is the $n$-dimensional vector with all-one components, and $\pi(\theta)^{\mathrm{T}}$ is the (row) vector obtained by transposing $\pi(\theta)$. This latter equation can be used to verify Assumption 5.

### B.2.2 OBJECTIVE FUNCTION

When all parameters are known and the number of servers is not too large, the normalizing constant $Z(\theta)$ and admission probability $J(\theta)$ can be computed efficiently using a variant of Buzen's algorithm (Buzen, 1973) for loss networks. Define the array $G = (G_{\underline{c},\underline{n}})_{\underline{c} \in \{0,1,\dots,c\}, \underline{n} \in \{1,2,\dots,n\}}$ by

$$G_{\underline{c},\underline{n}} = \sum_{\substack{s \in \mathbb{N}^{\underline{n}}: \\ |s| \le \underline{c}}} \prod_{i=1}^{\underline{n}} \left( \frac{\lambda}{\mu_i} \rho_i(\theta) \right)^{s_i}, \quad \underline{c} \in \{0,1,\dots,c\}, \quad \underline{n} \in \{1,2,\dots,n\}.$$

The dependency of $G$ on $\theta$ is left implicit to alleviate notation. The normalizing constant and admission probability are given by $Z(\theta) = G_{c,n}$ and $J(\theta) = G_{c-1,n}/G_{c,n}$, respectively. Defining the array $G$ allows us to calculate these metrics more efficiently than by direct calculation, as we have $G_{0,\underline{n}} = 1$ for each $\underline{n} \in \{1,2,\dots,n\}$, and

$$G_{\underline{c},1} = 1 + \frac{\lambda}{\mu_1} \rho_1(\theta) G_{\underline{c}-1,1}, \qquad \underline{c} \in \{1,2,\dots,c\}, \\ G_{\underline{c},\underline{n}} = G_{\underline{c},\underline{n}-1} + \frac{\lambda}{\mu_{\underline{n}}} \rho_{\underline{n}}(\theta) G_{\underline{c}-1,\underline{n}}, \qquad \underline{c} \in \{1,2,\dots,c\}, \quad \underline{n} \in \{2,3,\dots,n\}.$$

### B.2.3 ASSUMPTIONS OF SECTION 5

Assumptions 4, 5, and 6 are automatically satisfied because the state space is finite (with $|\mathcal{S}| = \binom{n+c}{c}$). Verifying Assumption 7 is challenging since it requires computing $\mathrm{Hess}_{\theta^*} J$ at the maximizer $\theta^*$, which depends in an implicit manner on the parameters of the system such as the arrival rate $\lambda$, service rates $\mu_1, \mu_2, \dots, \mu_n$, and policy $\pi(\theta^*)$. However, the nondegeneracy property of the Hessian for smooth functions is a property that is commonly stable in the following sense: if a function satisfies this property, then it will still be satisfied after any small-enough smooth perturbation. In particular, smooth functions with isolated nondegenerate critical points—also known as Morse functions—are dense and form an open
subset in the space of smooth functions (Nicolaescu, 2011, Section 1.2). Thus, unless the example is adversarial or presents symmetries, we can expect Assumption 7 to hold.

## B.3 Ising Model and Glauber Dynamics

Lastly, we focus on the example of Section 6.3. We consider the Markov chain defined by applying the policy (54) parameterized by some vector $\theta \in \mathbb{R}^3$: starting from an arbitrary initial configuration, at every time step, a coordinate is chosen uniformly at random, and the agent flips or not the spin at this coordinate according to the policy.

### B.3.1 PRODUCT-FORM STATIONARY DISTRIBUTION

The Markov chain is irreducible because it has a positive probability of transitioning from any configuration to any other as follows: at every step, choose a coordinate at which the two configurations differ and flip the spin at this coordinate. The Markov chain is positive recurrent because its state space is finite. Hence, Assumptions 1 and 2 are satisfied. We now focus on proving Assumption 3.

Our goal is to verify that the Markov chain that describes the random evolution of the state admits the stationary distribution (55), which we recall here:

$$p(\sigma, v|\theta) = \frac{1}{Z(\theta)} e^{\beta(\theta)JI(\theta) + \beta(\theta)\mu F(\sigma|\theta)}, \quad s = (\sigma, v) \in \mathcal{S}, \quad \theta \in \mathbb{R}^3,$$

Observe that $p(\sigma, v|\theta)$ is independent of $v$, hence we can let $q(\sigma|\theta) \triangleq p(\sigma, v|\theta)$ for each $\sigma \in \Sigma$. The key argument to prove that this is indeed the stationary distribution consists of observing that the policy (54) satisfies, for each $s = (\sigma, v) \in \mathcal{S}$,

$$\pi(\text{flip}|\sigma, v, \theta) = \frac{q(\sigma_{-v}|\theta)}{q(\sigma|\theta) + q(\sigma_{-v}|\theta)}, \quad \pi(\text{not flip}|\sigma, v, \theta) = \frac{q(\sigma|\theta)}{q(\sigma|\theta) + q(\sigma_{-v}|\theta)}, \quad (65)$$

where $\sigma_{-v} \in \Sigma$ is the configuration obtained by flipping the spin at $v$ compared to $\sigma$, that is, $\sigma_{-v}(w) = \sigma(w)$ for each $w \in \mathcal{V} \setminus \{v\}$ and $\sigma_{-v}(v) = -\sigma(v)$.

The balance equation for a particular state $s = (\sigma, v) \in \mathcal{S}$ writes

$$p(\sigma, v|\theta) = \sum_{w \in \mathcal{V}} p(\sigma, w|\theta)\pi(\text{not flip}|\sigma, w, \theta)\frac{1}{d_1d_2} + \sum_{w \in \mathcal{V}} p(\sigma_{-w}, w|\theta)\pi(\text{flip}|\sigma_{-w}, w, \theta)\frac{1}{d_1d_2}.$$

Dropping the dependency on $\theta$ to simplify notation, and injecting (65) into the right-hand side of this balance equation, we obtain successively

$$\sum_{w \in \mathcal{V}} p(\sigma, w) \frac{q(\sigma)}{q(\sigma) + q(\sigma_{-w})} \frac{1}{d_1 d_2} + \sum_{w \in \mathcal{V}} p(\sigma_{-w}, w) \frac{q((\sigma_{-w})_{-w})}{q(\sigma_{-w}) + q((\sigma_{-w})_{-w})} \frac{1}{d_1 d_2} \\ \stackrel{(1)}{=} \sum_{w \in \mathcal{V}} (p(\sigma, w) + p(\sigma_{-w}, w)) \frac{q(\sigma)}{q(\sigma) + q(\sigma_{-w})} \frac{1}{d_1 d_2} \stackrel{(2)}{=} q(\sigma) \sum_{w \in \mathcal{V}} \frac{1}{d_1 d_2} = q(\sigma) \stackrel{(2)}{=} p(\sigma, v),$$

where (1) follows by observing that $(\sigma_{-w})_{-w} = \sigma$ and (2) by recalling that $q(\sigma) = p(\sigma, w)$ for each $(\sigma, v) \in \mathcal{S}$. This proves that the distribution (55) is indeed the stationary distribution of the Markov chain that describes the evolution of the state under the policy (54).
Besides the sufficient statistics $x$, the inputs of Algorithm 2 are given, for each $\theta \in \mathbb{R}^3$ and $s = (\sigma, v) \in S$, by

$$ 
\begin{aligned}
\mathrm{D} \log \rho(\theta) &= \begin{bmatrix} \beta'(\theta)J & 0 & 0 \\ \beta'(\theta)\mu\theta_{\mathrm{left}}(\theta) & \beta(\theta)\mu h_{\mathrm{left}}'(\theta) & 0 \\ \beta'(\theta)\mu\theta_{\mathrm{right}}(\theta) & 0 & \beta(\theta)\mu h_{\mathrm{right}}'(\theta) \end{bmatrix}, \\
\nabla \log \pi(a|\sigma, v, \theta) &= (\mathbb{1}[a = \text{not flip}] - \pi(\text{not flip}|\sigma, v, \theta))\nabla\delta(\sigma, v|\theta), \\
\nabla\delta(\sigma, v|\theta) &= 2 \begin{bmatrix} \beta'(\theta)\sigma(v)(J \sum_{w \in \mathcal{V}: w \sim v} \sigma(w) + \mu h(v|\theta)) \\ \beta(\theta)\sigma(v)\mu h_{\mathrm{left}}'(\theta)\mathbb{1}[v_2 \le d_2/2] \\ \beta(\theta)\sigma(v)\mu h_{\mathrm{right}}'(\theta)\mathbb{1}[v_2 > d_2/2] \end{bmatrix},
\end{aligned}
$$

where $\beta'(\theta)$ (resp. $h'_{\mathrm{left}}(\theta)$, $h'_{\mathrm{right}}(\theta)$) is to be understood as the partial derivative of $\beta$ (resp. $h_{\mathrm{left}}$, $h_{\mathrm{right}}$) with respect to $\theta_1$ (resp. $\theta_2, \theta_3$).

### B.3.2 ASSUMPTIONS OF SECTION 5

In this example, note that since the state space, albeit large, is finite Assumption 4 is satisfied. By inspecting the load function and its derivatives we can also verify that Assumptions 5–6 hold. Finally, Assumption 7 does not hold as the optima in the parameter space is unbounded. Note that as we see in the simulations this assumption is not required in practice for the policy to converge.

## Appendix C. Proof of Theorem 2

In this section, we prove Theorem 2. The outline of this proof is given in Section 5.6.

### C.1 Preliminaries

We are going to use concentration inequalities for Markov chains. Such results are common in the literature (for example, see Karimi et al. 2019), and will be required to get a concentration bound of the plug-in estimators from (16).

Denote by $B_{\epsilon}(\theta)$ the open ball of radius $\epsilon$ centered at $\theta \in \Omega \subseteq \mathbb{R}^n$ and $\mathcal{Y} = S \times A$. Given a function $q: \mathcal{Y} \to \mathbb{R}$ and the Lyapunov function $\mathcal{L}: \mathcal{Y} \to [1, \infty)$ from Assumption 4, define

$$ |q|_{\mathcal{L}} = \sup_{y \in \mathcal{Y}} \frac{|q(y)|}{\mathcal{L}(y)}. \qquad (66) $$

Given a signed measure $\nu$, we also define the seminorm

$$ |\nu|_{\mathcal{L}} = \sup_{|q|_{\mathcal{L}} \le 1} |\nu[q]| = \sup_{|q|_{\mathcal{L}} \le 1} \left| \int q(y)\nu(dy) \right|. \qquad (67) $$

Equations (66) to (67) imply that

$$ |\nu[q]| \le |\nu|_{\mathcal{L}}|q|_{\mathcal{L}}. \qquad (68) $$

Note that we defined $| \cdot |_{\mathcal{L}}$ for a unidimensional function. Given instead $m$ functions $q_i: \mathcal{Y} \to \mathbb{R}$, for the higher-dimensional function $q: \mathcal{Y} \to \mathbb{R}^m$ that satisfies for all $y \in \mathcal{Y}$, $q(y) = (q_1(y), \dots, q_l(y))$, we define $|q|_{\mathcal{L}} = \sqrt{\sum_{i=1}^l |q_i|^2}_{\mathcal{L}}$.

The following lemma yields the concentration inequalities required:
**Lemma 13** Let {$Y_n$}$_{n}$$_{≥}$$_{1}$ be a geometrically ergodic Markov chain with invariant distribution p and transition matrix P(·, ·). Let the Lyapunov function be $\mathcal{L}$ : $\mathcal{Y}$ → $\mathbb{R}$. From geometric ergodicity, there exists $C > 0$ and $\lambda \in (0, 1)$ such that for any $y \in \mathcal{Y}$,

$$|P^m(\cdot|y) - p(\cdot)|_{\mathcal{L}} \le C\lambda^m. \qquad (69)$$

Let $\mathcal{F} = \sigma(Y_1)$ be the $\sigma$-algebra of $Y_1$. Let $q: \mathcal{Y} \to \mathbb{R}^m$ be a measurable function such that $|q|_{\mathcal{L}} < \infty$. For a finite trajectory $Y_1, \dots, Y_M$ of the Markov chain, we define the empirical estimator for $p[q]$ as

$$\hat{p}_M[q] = \frac{1}{M} \sum_{i=1}^{M} q(Y_i). \qquad (70)$$

With these assumptions, there exists $C'$ depending on $C$ and $\lambda$ such that

$$|\mathbb{E}[p[q] - \hat{p}_M[q]|\mathcal{F}]| \le \frac{C'|q|_{\mathcal{L}}}{M} \mathcal{L}(Y_1), \qquad (71)$$

and for $l \in \{1, 2, 4\}$,

$$\mathbb{E}[|p[q] - \hat{p}_M[q]|^l | \mathcal{F}] \le \frac{C' |q|_{\mathcal{L}}^l}{M^{l/2}} \mathcal{L}^l(Y_1). \qquad (72)$$

**Proof** We refer to (Fort and Moulines, 2003, Proposition 12) for a proof of (72). What remains is to prove (71).

Observe that for $y \in \mathcal{Y}$, $P(y) = P(\cdot|y)$ is a distribution over $\mathcal{Y}$. Conditional on $\mathcal{F}$, there exists $C > 0$ such that

$$\begin{aligned} |\mathbb{E}[\frac{1}{M} \sum_{i=1}^{M} q(Y_i) - p[q]|\mathcal{F}]| & \le \frac{1}{M} \sum_{i=1}^{M} |P^i(Y_1)[q] - p[q]| = \frac{1}{M} \sum_{i=1}^{M} |(P^i(Y_1) - p)[q]| \\ & \le \frac{1}{M} \sum_{i=1}^{M} |P^i(\cdot|Y_1) - p(\cdot)|_{\mathcal{L}} |q|_{\mathcal{L}} \mathcal{L}(Y_1) \le \frac{|q|_{\mathcal{L}}}{M} \sum_{i=1}^{M} C \lambda^i \mathcal{L}(Y_1) \\ & \le \frac{C |q|_{\mathcal{L}}}{M(1-\lambda)} \mathcal{L}(Y_1). \end{aligned} \qquad (73)$$

This concludes the proof.

<figure><img src="image_7.png" /></figure>

In epoch $m$, the Markov chain {$S_t$}$_{t}$$_{∈}$$_{[t_m, t_{m+1}]}$ with control parameter $\Theta_m$ has a Lyapunov function $\mathcal{L}_v$. Intuitively, as a consequence of Assumption 4, we can show that the process does not drift to infinity on the event $\mathcal{B}_m$ (despite the changing control parameter $\Theta_m$).

Specifically, for $m > 0$, let {$S_t$}$_{i}$$_{∈}$$_{[t_m, t_{m+1}]}$ be the Markov chain trajectory with transition probabilities $P(\Theta_m)$, where $\Theta_m$ is given by the updates in (3) and (16) and initial state $S_0 \in \mathcal{S}$. Recall that $\mathcal{B}_m$ is defined in (40). We can then prove the following:

**Lemma 14** Suppose Assumption 4 holds. There exists $D < \infty$ such that for $m > 0$, $\mathbb{E}[\mathcal{L}_v(S_{t_{m+1}})1[\mathcal{B}_m]] < D$. In particular, using the notation of Assumption 4 we may choose $D = \mathcal{L}^*$ as defined in (22).
**Proof** We will give an inductive argument. A similar argument can be found in the paper of Atchadé et al. (2017).

First, observe that for $m = 0$, $S_0$ is fixed. Thus, there exists a $D$ such that $\mathcal{L}_v(S_0) \le D$.

Next, assume that $\mathbb{E}[\mathcal{L}_v(S_{t_m})1[\mathcal{B}_{m-1}]] \le D$. On the event $\mathcal{B}_m$, Assumption 4 holds since $\Theta_1, \dots, \Theta_{m-1}, \Theta_m \in V_{\tau,\delta}(\theta^*) \subset U$. Thus, on the event $\mathcal{B}_m$, and when additionally conditioning on $S_{t_{m+1}-1}$ and $\Theta_m$, the following holds true:

$$ \begin{aligned} \mathbb{E}[\mathcal{L}_v(S_{t_{m+1}})1[\mathcal{B}_m]] &\le \mathbb{E}[\mathbb{E}[\mathcal{L}_v(S_{t_{m+1}})1[\mathcal{B}_m]|S_{t_{m+1}-1}]] \\ &= \mathbb{E}[1[\mathcal{B}_m]P_{\Theta_m}\mathcal{L}_v(S_{t_{m+1}-1})] \\ &\le \mathbb{E}[1[\mathcal{B}_m][\lambda\mathcal{L}_v(S_{t_{m+1}-1}) + b]]. \end{aligned} \quad (74) $$

The last step followed from Assumption 4.

Observe finally that the bound in (74) can be iterated by conditioning on $S_{t_{m+1}-2}$; so on and so forth. After $t_{m+1} - t_m$ iterations, one obtains

$$ \mathbb{E}[\mathcal{L}_v(S_{t_{m+1}})1[\mathcal{B}_m]] \le \lambda \mathbb{E}[\mathcal{L}_v(S_{t_m})1[\mathcal{B}_m]] + \frac{b}{1-\lambda}. \quad (75) $$

Noting that $1[\mathcal{B}_m] \le 1[\mathcal{B}_{m-1}]$, the claim follows by induction if we choose $D$ large enough such that $\lambda D + b/(1-\lambda) \le D$, that is $D \ge \mathcal{L}^*$.

## C.2 Proof of Theorem 2

To prove Theorem 2, we more-or-less follow the arguments of (Fehrman et al., 2020, Theorem 25). Modifications are however required because we consider a Markovian setting instead. Specifically, we rely on the bounds in Lemmas 10 to 12 instead of the bounds in (Fehrman et al., 2020, Proposition 20, Proposition 21, Proposition 24), respectively.

Let us begin by bounding

$$ \mathbb{P}[J^* - J(\Theta_m) > \epsilon|\mathcal{B}_0]. \quad (76) $$

Here, $\mathcal{B}_0 = \{\Theta_0 \in V_{\tau,\delta}(\theta^*)\}$—recall (40). Theorem 2 assumes that we initialize in a set $V$ which we will specify later but satisfies $V \subset V_{\tau,\delta}(\theta^*)$. Since we can initialize $\Theta_0$ with positive probability in $V$, we have that $\mathbb{P}[\mathcal{B}_0] \ge \mathbb{P}[\Theta_0 \in V] > 1/c > 0$ for some $c > 0$. Thus, we will focus on finding an upper bound of

$$ \mathbb{P}[J^* - J(\Theta_m) > \epsilon|\mathcal{B}_0] \le c\mathbb{P}[\{J^* - J(\Theta_m) > \epsilon\} \cap \mathcal{B}_0]. \quad (77) $$

Denote the orthogonal projection of $\Theta_m$ onto $M \cap U$ by $\tilde{\Theta}_m = \mathfrak{p}(\Theta_m)$. We can relate the objective gap $J^* - J(\Theta_m)$ to the distance $D_m := \text{dist}(\Theta_m, M \cap U)$ as follows. Since $J$ is twice continuously differentiable with maximum $J^*$ attained at $M \cap U$, the function $J(\theta)$ with $\theta \in V_{\tau,\delta}(\theta^*)$ is locally Lipschitz with constant $\iota_{\tau,\delta}(\theta^*) > 0$. On the event $\mathcal{B}_m$, we have $\Theta_m \in V_{\tau,\delta}(\theta^*)$ and therefore we have the inequality

$$ J^* - J(\Theta_m) = J(\tilde{\Theta}_m) - J(\Theta_m) \le \iota_{\tau,\delta}(\theta^*)|\tilde{\Theta}_m - \Theta_m| = \iota_{\tau,\delta}(\theta^*)D_m. \quad (78) $$

Consequently, we have the bound

$$ \mathbb{P}[\{J^* - J(\Theta_m) > \epsilon\} \cap \mathcal{B}_m] \le \mathbb{P}\left[\left\{D_m \ge \frac{\epsilon}{\iota_{\tau,\delta}(\theta^*)}\right\} \cap \mathcal{B}_m\right]. \quad (79) $$
If we define $\epsilon' = \epsilon/|\mathbf{r}_\delta(\theta^*)|$, the right-hand side of (79) can also be written as

$$ \begin{aligned} \mathbb{P}[\{D_m \ge \epsilon'\} \cap \mathcal{B}_m] &= \mathbb{E}[\mathbf{1}[D_m \ge \epsilon']\mathbf{1}[\mathcal{B}_m]] \\ &= \mathbb{E}[\mathbf{1}[D_m\mathbf{1}[B_m] \ge \epsilon']] = \mathbb{P}[D_m\mathbf{1}[B_m] \ge \epsilon'] \end{aligned} \quad (80) $$

by the positivity of $D_m$.

Next, we use (i) the law of total probability noting that $\mathcal{B}_m \subset \mathcal{B}_0$, (ii) the bound (79) and the inequality $\mathbb{P}[A \cap B] \le \mathbb{P}[A]$ for any two events $A, B$, and finally, (iii) the equality (80). We obtain

$$ \begin{aligned} \mathbb{P}[\{J^* - J(\Theta_m) > \epsilon\} \cap \mathcal{B}_0] &\stackrel{(i)}{\le} \mathbb{P}[\{J^* - J(\Theta_m) > \epsilon\} \cap \mathcal{B}_m] + \mathbb{P}[\{J^* - J(\Theta_m) > \epsilon\} \cap \overline{\mathcal{B}_m}] \\ &\stackrel{(ii)}{\le} \mathbb{P}[\{D_m \ge \epsilon'\} \cap \mathcal{B}_m] + \mathbb{P}[\overline{\mathcal{B}_m}] \\ &\stackrel{(iii)}{\le} \mathbb{P}[D_m\mathbf{1}[\mathcal{B}_m] \ge \epsilon'] + \mathbb{P}[\overline{\mathcal{B}_m}] \\ &\le \mathbb{P}[D_m\mathbf{1}[\mathcal{B}_{m-1}] \ge \epsilon'] + \mathbb{P}[\overline{\mathcal{B}_m}] = \text{Term I} + \text{Term II}. \end{aligned} \quad (81) $$

Term I can be bounded by using Markov's inequality and Lemma 10. This shows that

$$ \text{Term I} \le c\epsilon'^{-2}\mathcal{L}^*m^{-\sigma-\kappa}. \quad (82) $$

Term II can be bounded by Lemma 12. Specifically, one finds that there exists a constant $c > 0$ such that, if $\Theta_0 \in V_{\mathbf{r}/2,\delta}(\Theta^*)$,

$$ \text{Term II} \le 1 - \exp\left(-\frac{c\alpha^2}{\delta^2\ell}\right) + c\delta^{-2}\ell^{-1}m^{1-\sigma-\kappa} + c\alpha \frac{\left(m^{1-3/2\sigma-\kappa/2} + \ell^{-1/2}m^{1-5\sigma/8-\kappa/2}\right)}{(\mathbf{r}/2 - 2\delta)_+} \quad (83) $$

Note next that for any $\alpha \in (0, \alpha_0]$ and $c > 0$ there exists $\delta_0$ such that for any $\delta \in (0, \delta_0]$ there exists $\ell_0$ such that if $\ell \in [\ell_0, \infty)$ there exists a constant $c' > 0$ such that we have the inequality $1 - \exp(-c\alpha^2/\delta^2\ell) \le c'\alpha^2/\delta^2\ell$. We can substitute this bound in (83) to yield

$$ \text{Term II} \le c' \frac{\alpha^2}{\delta^2 \ell} + \text{idem.} \quad (84) $$

Bounding (81) by the sum of (82) and (84), and substituting the bound in (77) reveals that there exists a constant $c'' > 0$ such that if $\Theta_0 \in V_{\mathbf{r}/2,\delta}(\Theta^*)$ then

$$ \begin{aligned} \mathbb{P}[J^* - J(\Theta_m) > \epsilon | \mathcal{B}_0] &\le c''(\epsilon')^{-2}\mathcal{L}^*m^{-\sigma-\kappa} + c''\alpha^2\delta^{-2}\ell^{-1} + c''\delta^{-2}\ell^{-1}m^{1-\sigma-\kappa} \\ &\quad + c''\alpha \frac{\left(m^{1-3/2\sigma-\kappa/2} + \ell^{-1/2}m^{1-5\sigma/8-\kappa/2}\right)}{(\mathbf{r}/2 - 2\delta)_+} \end{aligned} \quad (85) $$

Note that the exponents of $m$ in (85) satisfy that since $\sigma \in (2/3, 1)$, $1 - 3/2\sigma - \kappa/2 \le -\kappa/2$ as well as $1 - 5\sigma/8 - \kappa/2 < 1 - \sigma/2 - \kappa/2$. Finally, let the initialization set be $V = V_{\mathbf{r}/2,\delta}(\Theta^*)$. Note that since $\{\Theta_0 \in V\} \subset \mathcal{B}_0$ there exists a constant $c''' > 0$ such that

$$ \mathbb{P}[J^* - J(\Theta_m) > \epsilon | \Theta_0 \in V] \le c'''\mathbb{P}[J^* - J(\Theta_m) > \epsilon | \mathcal{B}_0]. \quad (86) $$

Substituting the upper bound (85) in (86) concludes the proof. $\square$
### C.3 Proof of Lemma 8

For simplicity, we will denote $t_{m+1} - t_m = T_m$, $X_t = x(S_t)$ throughout this proof. We also temporarily omit the summation indices for the epoch. We note that the policies defined in (6) satisfy that for $(s, a) \in \mathcal{S} \times \mathcal{A}$,

$$ (\nabla \log \pi(a|s, \theta))_{i,a'} = \begin{cases} 1[a = a'] - \pi(a'|s, \theta) & \text{if } i = h(s), \\ 0 & \text{otherwise.} \end{cases} $$

In particular, there exists $c_1 > 0$ such that for any $(s, a) \in \mathcal{S} \times \mathcal{A}$, $|\nabla \log \pi(a|s, \theta)| < c_1$. The proof below, however, can also be extended to other policy classes.

#### C.3.1 PROOF OF (45)

Observe that if the event $B_m$ holds, then the definitions in (16) also imply that

$$ \begin{aligned} \eta_m &= \nabla J(\Theta_m) - H_m = \nabla J(\Theta_m) - (\mathrm{D} \log \rho(\Theta_m)^{\top} \bar{C}_m + \bar{E}_m) \\ &= \nabla J(\Theta_m) - (\mathrm{D} \log \rho(\Theta_m)^{\top} \frac{1}{T_{m+1}} \sum_{t=t_m}^{t_{m+1}-1} (X_t - \bar{X}_m) r(S_t, A_t) \\ &\quad + \frac{1}{T_m} \sum_{t=t_m}^{t_{m+1}-1} r(S_t, A_t) \nabla \log \pi(A_t|S_t, \Theta_m) \\ &= \mathrm{D} \log \rho(\Theta_m)^{\top} (\mathrm{Cov}[R, S] - \frac{1}{T_m} \sum_{t=t_m}^{t_{m+1}-1} (X_t - \bar{X}_m) r(S_t, A_t)) \\ &\quad + (\mathrm{E}[R \nabla \log \pi(A|S, \Theta_m)] - \frac{1}{T_m} \sum_{t=t_m}^{t_{m+1}-1} r(S_t, A_t) \nabla \log \pi(A_t|S_t, \Theta_m)) \\ &= \mathrm{D} \log \rho(\Theta_m)^{\top} \tilde{\eta}_m + \tilde{\zeta}_m. \end{aligned} \quad (87) $$

We will deal with the terms $\tilde{\eta}_m$ in and $\tilde{\zeta}_m$ in (87) one-by-one.

Let us first consider the 1st term, $\tilde{\eta}_m$. Define

$$ \begin{aligned} A &= \mathbb{E}[(X - \mathbb{E}[X])R] - \frac{1}{T_m} \sum_t (X_t - \mathbb{E}[X])r(S_t, A_t), \\ B &= \frac{1}{T_m} \left( \sum_t r(S_t, A_t) \right) (\mathbb{E}[X] - \bar{X}_m), \end{aligned} \quad (88) $$

and observe that

$$ \tilde{\eta}_m = A + B. \quad (89) $$

We look first at A in (88). Recall that $\{Y_t\}_{t>0} = \{(S_t, A_t)\}_{t>0}$ is the chain of state-action pairs (see Section 5.1). Define the function $g: \mathcal{S} \times \mathcal{A} \to \mathbb{R}^n$ as

$$ g(y) = g((s, a)) = (x(s) - \mathbb{E}[x(s)])r(y). \quad (90) $$
Then, we can rewrite

$$A = \mathbb{E}[g(Y)] - \frac{1}{T_m} \sum_t g(Y_t). \qquad (91)$$

We are now almost in position to apply Lemma 13 to A. Observe next that the law of total expectation implies that

$$\mathbb{E}[\eta_m \mathbf{1}[\mathcal{B}_m] | \mathcal{F}_m] = \sum_{a \in \mathcal{A}} \mathbb{E}[\eta_m \mathbf{1}[\mathcal{B}_m] | \mathcal{F}_m, A_{t_m} = a] \pi(a|S_m, \Theta_m), \qquad (92)$$

Without loss of generality, it therefore suffices to consider the case that we have one action $A_{t_m} = a \in \mathcal{A}$. For the first term we have that there exists a constant $c_2 > 0$ such that

$$\begin{aligned} |\mathbb{E}[A \mathbf{1}[\mathcal{B}_m] | \mathcal{F}_m, A_{t_m} = a]| &= |\mathbb{E}[\mathbb{E}[g(Y)] - \frac{1}{T_m} \sum_t g(Y_t) \mathbf{1}[\mathcal{B}_m] | Y_0 = (S_{t_m}, A_{t_m})]| \\ &\stackrel{\text{(Lemma 13)}}{=} \frac{c_2 |g|_{\mathcal{L}}}{T_m} \mathcal{L}((S_{t_m}, a)), \end{aligned} \qquad (93)$$

where we can use that $|g|_{\mathcal{L}} < \infty$ due to Assumption 6.

For the term B in (88). We can add and subtract again the following terms and obtain

$$\begin{aligned} B &= \frac{1}{T_m} \left( \sum_t r(S_t, A_t) \right) \left( \mathbb{E}[X] - \bar{X}_m \right) - \mathbb{E}[R] \left( \mathbb{E}[X] - \bar{X}_m \right) + \mathbb{E}[R] \left( \mathbb{E}[X] - \bar{X}_m \right) \\ &= C + D, \end{aligned} \qquad (94)$$

where

$$\begin{aligned} C &= (\mathbb{E}[X] - \bar{X}_m) \left( \frac{1}{T_m} \sum_t r(S_t, A_t) - \mathbb{E}[R] \right), \\ D &= \mathbb{E}[R] (\mathbb{E}[X] - \bar{X}_m). \end{aligned} \qquad (95)$$

For the term D in (95) we can readily use the concentration of Lemma 13 to obtain

$$\mathbb{E}[\mathbb{E}[R](\mathbb{E}[X] - \bar{X}_m)\mathbf{1}[\mathcal{B}_m] | \mathcal{F}_m, A_{t_m} = a] \le \mathbb{E}[R] \frac{|x(S)|_{\mathcal{L}}}{T_m} \mathcal{L}(S_{t_m}, a), \qquad (96)$$

where we have $|x(S)|_{\mathcal{L}} < \infty$ from Assumption 6 and $\mathbb{E}[R] < J^*$.

For the term C, we use Cauchy-Schwartz together with Lemma 8. In particular, we have

$$\begin{aligned} & \left| \mathbb{E}[(\mathbb{E}[X] - \bar{X}_m) \left( \frac{1}{T_m} \sum_t r(S_t, A_t) - \mathbb{E}[R] \right) \mathbf{1}[\mathcal{B}_m] | \mathcal{F}_m, A_{t_m} = a] \right| \le \\ & \qquad \left| \mathbb{E}[|\mathbb{E}[X] - \bar{X}_m|^2 \mathbf{1}[\mathcal{B}_m] | \mathcal{F}_m, A_{t_m} = a] \right|^{1/2} \times \\ & \qquad \left| \mathbb{E}\left[\left| \frac{1}{T_m} \sum_t r(S_t, A_t) - \mathbb{E}[R] \right|^2 \mathbf{1}[\mathcal{B}_m] | \mathcal{F}_m, A_{t_m} = a \right] \right|^{1/2}. \end{aligned} \qquad (97)$$
For both terms we can repeat the same argument to that in (92) together with Lemma 13 to show that

$$ \begin{aligned} \left|\mathbb{E}\left[\left|\mathbb{E}[X] - \bar{X}_m\right|^2 \mathbf{1}[\mathcal{B}_m]\right| \mathcal{F}_m, A_{t_m} = a\right|\right|^{1/2} &\le c_3 \frac{|X|_{\mathcal{L}}^{1/2}}{T_m^{1/2}} \mathcal{L}(S_{t_m}, a), \\ \left|\mathbb{E}\left[\left|\frac{1}{T_m} \sum_t r(S_t, A_t) - \mathbb{E}[R]\right|^2 \mathbf{1}[\mathcal{B}_m]\right| \mathcal{F}_m, A_{t_m} = a\right|\right|^{1/2} &\le c_4 \frac{|R|_{\mathcal{L}}^{1/2}}{T_m^{1/2}} \mathcal{L}(S_{t_m}, a). \end{aligned} \quad (98) $$

Therefore multiplying both bounds in (98) and using Assumption 6 to bound the $\mathcal{L}$-norms, we obtain that there exists $c_5 > 0$ such that

$$ \left|\mathbb{E}[C|\mathcal{F}_m, A_{t_m} = a]\right| \le \frac{c_5}{T_m} \mathcal{L}(S_{t_m}, a)^2. \quad (99) $$

Adding the bounds (93), (99), and (96) together we have now

$$ \left|\mathbb{E}[\tilde{\eta}_m \mathbf{1}[\mathcal{B}_m] | \mathcal{F}_m, A_{t_m} = a]\right| \le \frac{c_6}{T_m} \mathcal{L}^2(S_{t_m}, a). \quad (100) $$

Finally, averaging this bound over all actions in (92), we obtain

$$ \left|\mathbb{E}[\tilde{\eta}_m \mathbf{1}[\mathcal{B}_m] | \mathcal{F}_m]\right| \le \frac{c_7}{T_m} \left(\sum_a \mathcal{L}(S_{t_m}, a)^2 \pi(a|S_{t_m}, \Theta_m)\right) \le \frac{c_7}{T_m} \mathcal{L}_4(S_{t_m})^{1/2}. \quad (101) $$

Now we use Assumption 5. We can write

$$ \begin{aligned} \left|\mathbb{E}[\nabla \log(\Theta_m) \tilde{\eta}_m \mathbf{1}[\mathcal{B}_m] | \mathcal{F}_m]\right| &= \left|\nabla \log(\Theta_m) \mathbb{E}[\tilde{\eta}_m \mathbf{1}[\mathcal{B}_m] | \mathcal{F}_m]\right| \\ &\le C \left|\mathbb{E}[\tilde{\eta}_m \mathbf{1}[\mathcal{B}_m] | \mathcal{F}_m]\right| \\ &\le \frac{c_8}{T_m} \mathcal{L}(S_{t_m}). \end{aligned} \quad (102) $$

Let us now consider the 2nd term, $\tilde{\zeta}_m$. Define a function of $Y = (S, A)$ as

$$ g(Y) = r(Y) \nabla \log \pi(A|S, \theta), \quad (103) $$

so that

$$ \zeta_m = \mathbb{E}[g(Y)] - \frac{1}{T_m} \sum_t g(Y_t). \quad (104) $$

By combining the argument of (92) with the fact that $|g(Y)|_{\mathcal{L}} < \infty$ by Assumption 6, we find that

$$ \left|\mathbb{E}[\tilde{\zeta}_m \mathbf{1}[\mathcal{B}_m] | \mathcal{F}_m]\right| \le \frac{c_9}{T_m} \mathcal{L}(S_{t_m}). \quad (105) $$

Adding (101) and (105) together with their largest exponents yields

$$ \begin{aligned} \left|\mathbb{E}[\eta_m \mathbf{1}[\mathcal{B}_m] | \mathcal{F}_m]\right| &\le \frac{c_{10}}{T_m} \sum_a \mathcal{L}(S_{t_m}, a)^2 \pi(a|S_{t_m}) \\ &\le \frac{c_{10}}{T_m} \left(\sum_a \mathcal{L}(S_{t_m}, a)^4 \pi(a|S_{t_m})\right)^{1/2} \le \frac{c_{10}}{T_m} \mathcal{L}_4(S_{t_m})^{1/2}. \end{aligned} \quad (106) $$

This concludes the proof of (45).

<figure><img src="image_8.png" /></figure>
### C.3.2 PROOF OF (46)

Note that by using the fact that for a vector-valued random variable $Z$ we have that $\mathbb{E}[|Z|^2] \ge \mathbb{E}[|Z|]^2$, the case for $p=1$ follows from the case $p=2$.

We focus on the case $p=2$. By using the identity $(a+b) \le 2a^2 + b^2$, we estimate

$$ \begin{aligned} & \mathbb{E}[|D \log \rho(\Theta_m)^{\top} \tilde{\eta}_m + \tilde{\zeta}_m|^2 \mathbf{1}[\mathcal{B}_m] | \mathcal{F}_m] \\ & \le 2(\mathbb{E}[|D \log \rho(\Theta_m)^{\top} \tilde{\eta}_m|^2 \mathbf{1}[\mathcal{B}_m] | \mathcal{F}_m] + \mathbb{E}[|\tilde{\zeta}_m|^2 \mathbf{1}[\mathcal{B}_m] | \mathcal{F}_m]) \\ & \stackrel{(5)}{\le} 2c_1^2 \mathbb{E}[|\tilde{\eta}_m|^2 \mathbf{1}[\mathcal{B}_m] | \mathcal{F}_m] + 2\mathbb{E}[|\tilde{\zeta}_m|^2 \mathbf{1}[\mathcal{B}_m] | \mathcal{F}_m]. \end{aligned} \quad (107) $$

We again use the law of total expectation with the action set in (92) and condition on the action $A_m = a$.

For the term involving $\tilde{\zeta}_m$ in (107) we can again use the definition of $g$ in (103). We bound

$$ \begin{aligned} \mathbb{E}[|\tilde{\zeta}_m|^2 \mathbf{1}[\mathcal{B}_m] | \mathcal{F}_m, A_{t_m} = a] &= \mathbb{E}[|\mathbb{E}[g(Y) - \frac{1}{T_m} \sum_t g(Y)|^2 | Y_0 = (S_{t_m}, a)] \\ &\stackrel{\text{(Lemma 13)}}{\le} \frac{c_2}{T_m} \mathcal{L}(S_{t_m}, a)^2. \end{aligned} \quad (108) $$

For the term involving $\tilde{\eta}_m$ in (107), we use the same definition for the terms $A, C$ and $D$ from (93), (99) and (96) as in the proof of (45). We have the bound

$$ \mathbb{E}[|\tilde{\eta}_m|^2 \mathbf{1}[\mathcal{B}_m] | \mathcal{F}_m, A_{t_m} = a] \le 3(\mathbb{E}[|A|^2 \mathbf{1}[\mathcal{B}_m] | \mathcal{F}_m, A_{t_m} = a] + \mathbb{E}[|C|^2 \mathbf{1}[\mathcal{B}_m] | \mathcal{F}_m, A_{t_m} = a] \\ + \mathbb{E}[|D|^2 \mathbf{1}[\mathcal{B}_m] | \mathcal{F}_m, A_{t_m} = a]). \quad (109) $$

For the terms pertaining to $A$ and $D$ in (109) the same argument as those used for $\tilde{\zeta}_m$ in (103) and (108) can be used to show that

$$ \begin{aligned} \mathbb{E}[|A|^2 \mathbf{1}[\mathcal{B}_m] | \mathcal{F}_m, A_{t_m} = a] &\le \frac{c_3}{T_m} \mathcal{L}(S_{t_m}, a)^2 \\ \mathbb{E}[|D|^2 \mathbf{1}[\mathcal{B}_m] | \mathcal{F}_m, A_{t_m} = a] &\le \frac{c_4}{T_m} \mathcal{L}(S_{t_m}, a)^2. \end{aligned} \quad (110) $$

The only remaining term to bound in (109) is $C$. We use again Cauchy-Schwartz's inequality

$$ \begin{aligned} \mathbb{E} \left[ \left| \left( \mathbb{E}[X] - \bar{X}_m \right) \left( \frac{1}{T_m} \sum_t r(S_t, A_t) - \mathbb{E}[R] \right) \mathbf{1}[\mathcal{B}_m] \right|^4 \middle| \mathcal{F}_m, A_{t_m} = a \right] &\le \\ & \mathbb{E} \left[ \left| \left( \mathbb{E}[X] - \bar{X}_m \right)^2 \mathbf{1}[\mathcal{B}_m] \middle| \mathcal{F}_m, A_{t_m} = a \right| \right]^{1/2} \times \\ & \mathbb{E} \left[ \left| \frac{1}{T_m} \sum_t r(S_t, A_t) - \mathbb{E}[R] \right|^4 \mathbf{1}[\mathcal{B}_m] \middle| \mathcal{F}_m, A_{t_m} = a \right]^{1/2} m, \end{aligned} \quad (111) $$

and by Lemma 13 the following hold

$$ \left| \mathbb{E} \left[ \left| \mathbb{E}[X] - \bar{X}_m \right|^4 \mathbf{1}[\mathcal{B}_m] \middle| \mathcal{F}_m, A_{t_m} = a \right] \right|^{1/2} \le c_5 \frac{|X|^2}{T_m} \mathcal{L}(S_{t_m}, a)^2, $$
$$ \left| \mathbb{E} \left[ \left| \frac{1}{T_m} \sum_t r(S_t, A_t) - \mathbb{E}[R] \right|^4 \mathbb{1}[\mathcal{B}_m] \left| \mathcal{F}_m, A_{t_m} = a \right. \right| \right|^{1/2} \le c_6 \frac{|R|_\mathcal{L}^{1/2}}{T_m} \mathcal{L}(S_{t_m}, a)^2. \quad (112) $$

The bound for $C$ thus becomes

$$ \mathbb{E}[|C|^2 | \mathcal{F}_m, A_{t_m} = a] \le \frac{c_7}{T_m^2} \mathcal{L}(S_{t_m}, a)^4. \quad (113) $$

Upper bounding all terms by the largest exponents and adding over the different actions, we finally obtain

$$ \mathbb{E}[|\eta_m|^2 \mathbb{1}[\mathcal{B}_m] | \mathcal{F}_m] \le \frac{c_8}{T_m} \sum_a \mathcal{L}(S_{t_m}, a)^4 \pi(a|S_{t_m}), \quad \Theta_m \le \frac{c_9}{T_m} \mathcal{L}_4(S_{t_m}). \quad (114) $$

That is it.

<figure><img src="image_9.png" /></figure>

### C.4 Proof of Lemma 10

We will again use the notation $t_{m+1} - t_m = T_m$ and without loss of generality we will assume that $T_m = \ell m^{\sigma/2+\kappa}$ instead of $\lfloor \ell m^{\sigma/2+\kappa} \rfloor$. This can be assumed since for $m \ge 1$ there exist constants $c_l, c_u > 0$ such that $c_l \ell m^{\sigma/2+\kappa} \le t_{m+1} - t_m \ge c_u \ell m^{\sigma/2+\kappa}$. The proof of Lemma 10 follows the same steps as in (Fehrman et al., 2020, Proposition 20). However, we have to quickly diverge and adapt the estimates to the case that there the variance of $H_m$ depends on the states of a Markov chain. From the assumptions, it can be shown that there is a unique differentiable orthogonal projection map $\mathfrak{p} : V_{\tau,\delta}(\theta^*) \to \mathcal{M} \cap U$ from $V_{\tau,\delta}(\theta^*) \cap U$ onto $V_{\tau,\delta}(\theta^*) \cap \mathcal{M} \cap U$. The distance of $\Theta_m$ to the set of minima can then be upper bounded by the distance to the projection $\mathfrak{p} : V_{\tau,\delta}(\theta^*) \to \mathcal{M} \cap U$ of $\Theta_{m-1}$ by

$$ \begin{aligned} \mathrm{dist}(\Theta_m, \mathcal{M} \cap U)^2 &\le |\Theta_m - \mathfrak{p}(\Theta_{m-1})|^2 \\ &\le |\Theta_{m-1} - \mathfrak{p}(\Theta_{m-1}) - \alpha_{m-1} \nabla J(\Theta_{m-1}) \\ &\qquad + (\alpha_{m-1} \nabla J(\Theta_{m-1}) - \alpha_{m-1} H_{m-1})|^2. \end{aligned} \quad (115) $$

After expanding (115) and taking expectations, however, the effect of bias already appears, and we must diverge from the analysis from (Fehrman et al., 2020, Equation 44) thereafter. In particular, the effect of the bias of $H_{m-1}$ needs to be handled in the terms

$$ \mathbb{E} \left[ 2 \langle \Theta_{m-1} - \mathfrak{p}(\Theta_{m-1}) - \alpha_{m-1} \nabla J(\Theta_{m-1}), \alpha_{m-1} \nabla J(\Theta_{m-1}) - \alpha_{m-1} H_{m-1} \rangle \mathbb{1}[\mathcal{B}_{m-1}] \right], \quad (116) $$

and

$$ \mathbb{E} \left[ \left| \alpha_{m-1} \nabla J(\Theta_{m-1}) - \alpha_{m-1} H_{m-1} \right|^2 \mathbb{1}[\mathcal{B}_{m-1}] \right] = (\alpha_{m-1})^2 \mathbb{E} \left[ |\eta_{m-1}|^2 \mathbb{1}[\mathcal{B}_{m-1}] \right]. \quad (117) $$

We specifically require bounds of these terms without relying on independence of the iterands.

We focus on (117) first. Recall for $m > 0$, that $\mathcal{F}_m$ is the sigma algebra defined in (41). By using the tower property of the conditional expectation and conditioning on $\mathcal{F}_{m-1}$, from Lemma 8 together with the fact that $T_m < cT_{m-1}$ for some $c > 0$, we obtain directly

$$ (117) = (\alpha_{m-1})^2 \mathbb{E} \left[ \mathbb{E} \left[ |\eta_{m-1}|^2 \mathbb{1}[\mathcal{B}_{m-1}] | \mathcal{F}_{m-1} \right] \right] \stackrel{\text{(Lemma 8)}}{\le} (\alpha_{m-1})^2 \frac{c_1}{T_m} \mathbb{E} \left[ \mathcal{L}_4(S_{t_{m-1}})^2 \mathbb{1}[\mathcal{B}_{m-1}] \right]. \quad (118) $$
Let us next bound (116). Note that this term does not vanish due to dependence of the samples conditional on $\mathcal{F}_{m-1}$. In our case, however, we have a Markov chain trajectory whose kernel will depend on $\Theta_{m-1}$. Let

$$Z_{m-1} = \Theta_{m-1} - \mathbf{p}(\Theta_{m-1}) - \alpha_{m-1} \nabla J(\Theta_{m-1}). \quad (119)$$

We use the law of total expectation again on (116). Note that $Z_{m-1}$ and $\mathcal{B}_{m-1}$ are $\mathcal{F}_{m-1}$-measurable.

$$ (116) \le 2\alpha_{m-1}\mathbb{E}[\langle \mathbf{1}[\mathcal{B}_{m-1}]Z_{m-1}, \mathbb{E}[\eta_{m-1}|\mathcal{F}_{m-1}]\rangle] \\ \le 2\alpha_{m-1}\mathbb{E}[|Z_{m-1}|^2\mathbf{1}[\mathcal{B}_{m-1}]]^{1/2} \mathbb{E}[|\mathbb{E}[\eta_{m-1}\mathbf{1}[\mathcal{B}_{m-1}]|\mathcal{F}_{m-1}]|^2]^{1/2} \\ \le 2\alpha_{m-1}\mathbb{E}[|Z_{m-1}|^2\mathbf{1}[\mathcal{B}_{m-1}]]^{1/2} \mathbb{E}[\mathbf{1}[\mathcal{B}_{m-1}]\mathcal{L}_4(S_{t_{m-1}})^2]^{1/2} \frac{c_2}{T_m}, \quad (120) $$

where (i) have used Cauchy-Schwartz and (ii) Lemma 8 and the fact that for some $c > 0$, $T_m < cT_{m-1}$.

The terms in (118) and (120) containing $\mathcal{L}_4(S_{t_m})$ can be upper bounded as follows. From the definition of (44) and since $v \ge 16$, by a generalized mean inequality and the fact that $\mathcal{L}(s, a) \ge 1$ for any $(s, a) \in \mathcal{S} \times \mathcal{A}$ we have

$$\mathcal{L}_4(s) \le \mathcal{L}_v(s)^{4/v} \le \mathcal{L}_v(s)^{1/4}. \quad (121)$$

Now, by Lemma 14, $\mathcal{L}^*$ is such that for all $m \in \mathbb{N}$

$$\mathbb{E}[\mathbf{1}[\mathcal{B}_{m-1}]\mathcal{L}_4(S_{t_{m-1}})^2] \le \mathbb{E}[\mathbf{1}[\mathcal{B}_{m-2}]\mathcal{L}_4(S_{t_{m-1}})^2] \stackrel{(121)}{\le} \mathbb{E}[\mathbf{1}[\mathcal{B}_{m-2}]\mathcal{L}_v(S_{t_{m-1}})] \le \mathcal{L}^*. \quad (122)$$

For the other term in (120), we can use the same bound used in (Fehrman et al., 2020, Equation 41): There exists constants $y, c > 0$ depending on $J, \theta^*$ and $\mathfrak{r}_0$ such that on the event $\mathcal{B}_{m-1}$ we have

$$|Z_{m-1}|^2 \le (1 - \alpha_{m-1}y)^2 \mathrm{dist}(\Theta_{m-1}, \mathcal{M} \cap U)^2 + c(1 - \alpha_{m-1}y)\alpha_{m-1}\mathrm{dist}(\Theta_{m-1}, \mathcal{M} \cap U)^3 \\ + c(\alpha_{m-1})^2\mathrm{dist}(\Theta_{m-1}, \mathcal{M} \cap U)^4. \quad (123)$$

The bound in (123) characterizes the fact that, close to the manifold of maximizers, the projection is differentiable and can be approximated by an orthogonal expansion of $J$ around the manifold of maximizers. The error terms of this expansion can be bounded depending on the Hessian at $\mathbf{p}(\Theta_{m-1}) \in \mathcal{M} \cap U$, Hess$_{\mathbf{p}(\Theta_{m-1})}$ $J$. We refer to (Fehrman et al., 2020, Proposition 17) for a proof of this fact.

We will now use an induction argument to show the claim of the lemma. Namely, we will assume for the time being that for $m-1$ we have

$$\mathbb{E}[(\mathrm{dist}(\Theta_{m-1}, \mathcal{M} \cap U) \wedge \delta)^2 \mathbf{1}[\mathcal{B}_{m-1}]] \le \delta^2 c(\alpha)(m-1)^{-\sigma-\kappa}, \quad (124)$$

where $c(\alpha) > 0$ is a function of $a$ to be determined. We want to show (124) for $m$. To do so we will use (123) to bound $Z_{m-1}$. Suppose that there exists a sequence $\{b_l\}_{l>0} \subset \mathbb{R}_+$ such that we have

$$\mathbb{E}[|Z_{m-1}|^2 \mathbf{1}[\mathcal{B}_{m-1}]] \le b_{m-1}. \quad (125)$$
Using (125) in (120) yields that for some $c_3 > 0$ we have:

$$ (116) \le 2(b_{m-1})^{1/2} \alpha_{m-1} (\mathcal{L}^*)^{1/2} \frac{c_3}{T_m}. \quad (126) $$

From the expansion of (115) and combining the bounds of (123) and (126) together we obtain

$$ \begin{aligned} \mathbb{E}[\mathrm{dist}(\Theta_m, \mathcal{M} \cap U)^2 \mathbf{1}[\mathcal{B}_{m-1}]] \le & b_{m-1} + 2(b_{m-1})^{1/2} \alpha_{m-1} (\mathcal{L}^*)^{1/2} \frac{c_3}{T_m} \\ & + (\alpha_{m-1})^2 \frac{c_4}{T_m} \mathcal{L}^*. \end{aligned} \quad (127) $$

We show now that from the induction hypothesis, if (124) holds, then we also have the bound

$$ b_{m-1} \le c(\alpha) \delta^2 m^{-\sigma-\kappa} - \delta^2 \frac{\alpha y}{2} c(\alpha) (m-1)^{-\sigma-\kappa} m^{-\sigma}. \quad (128) $$

Indeed, taking expectations in (123) and using the bound (124) yields

$$ \begin{aligned} b_{m-1} \le & (1 - \alpha_{m-1}y)^2 c(\alpha) (m-1)^{-\sigma-\kappa} + c(\alpha)(1 - \alpha_{m-1}y)\alpha_{m-1}\delta c(\alpha)(m-1)^{-\sigma-\kappa} \\ & + c(\alpha)(\alpha_{m-1})^2 \delta^2 c(\alpha)(m-1)^{-\sigma-\kappa}. \end{aligned} \quad (129) $$

Recall that $\alpha_{m-1} = \alpha m^{-\sigma/2-\kappa}$. Adding and subtracting $c(\alpha)m^{-\sigma-\kappa}$ in (129), we obtain that

$$ b_{m-1} \le c(\alpha)m^{-\sigma-\kappa} + c(\alpha)m^{-\sigma}(m-1)^{-\sigma-\kappa} \left( m^{\sigma} - (m-1)^{\sigma+\kappa}m^{-\kappa} - 2\alpha y + \frac{\alpha^2 y}{m^{\sigma}} + (1-\frac{\alpha y}{m^{\sigma}})\alpha\delta + \delta^2 \frac{\alpha^2 y^2}{m^{\sigma}} \right). $$

Note now that there exists $m_0(a) > 0$ such that if $m \ge m_0(a)$, we have

$$ m^{\sigma} - (m-1)^{\sigma+\kappa}m^{-\kappa} - \alpha y + \frac{\alpha^2 y}{m^{\sigma}} < -\frac{\alpha y}{2}. \quad (130) $$

Indeed, note that the latter equation can be satisfied for $m \ge m_0(a)$ since there exists a constant $c > 0$ depending on $\sigma$ and $\kappa$ such that

$$ \begin{aligned} m^{\sigma} - (m-1)^{\sigma+\kappa}m^{-\kappa} \le & m^{-\kappa}(m^{\sigma+\kappa} - (m-1)^{\sigma+\kappa}) \\ \le & m^{-\kappa}(\sigma + \kappa) \max[(m-1)^{\sigma+\kappa-1}, m^{\sigma+\kappa-1}] \\ \le & c_5(\sigma + \kappa)m^{\sigma-1}. \end{aligned} \quad (131) $$

In this case we have that

$$ m_0(\alpha) = \left( \frac{2c_5(\sigma + \kappa)}{y\alpha} \right)^{1-\sigma} > \frac{c'}{\alpha^{1-\sigma}}. \quad (132) $$

Then for $m > m_0(\alpha)$, we will have

$$ b_m \le c(\alpha)m^{-\sigma-\kappa} + c(\alpha)m^{-\sigma}(m-1)^{-\sigma-\kappa} \left( -\frac{\alpha y}{2} + (1-\frac{\alpha y}{m^{\sigma}})\alpha\delta + \delta^2 \frac{\alpha^2 y^2}{m^{\sigma}} \right). $$
Choose $\delta \in (0, \delta_1(\alpha)]$, where $\delta_1(\alpha)$ is a bound that we will choose appropriately, such that for any $m \ge m_0(\alpha)$ we have

$$ (1 - \frac{\alpha y}{m^{\sigma}}) \alpha \delta + \delta^2 \frac{\alpha^2 y^2}{m^{\sigma}} \le \alpha y. \qquad (133) $$

Thus, from (124) we obtain (128). With (128) with an appropriate choice of $c(\alpha)$, we can now show (124) for $m$. We will namely choose $c(\alpha)$ as follows

$$ c(\alpha) = \max\left(\frac{c'}{\alpha^{(1-\sigma)(\sigma+\kappa)}}, \frac{4C^2\mathcal{L}^* + 4yC\mathcal{L}^*\alpha\ell}{\delta^2\ell^2y^2}\right), \qquad (134) $$

where recall that $\delta \in (0, \delta_1(\alpha)]$ and $\delta_1(\alpha)$ were chosen so that (133) holds. Let $L = \ell^{-1}$. Substituting the bound of (128) into (127) and recalling that $T_m = m^{\kappa+\sigma/2}\ell$ yields

$$ \begin{aligned} \mathbb{E}[(\mathrm{dist}(\Theta_m, \mathcal{M} \cap U))^2 \mathbf{1}[\mathcal{B}_{m-1}]] & \le c(\alpha)\delta^2 m^{-\sigma-\kappa} - \frac{\alpha y}{2}c(\alpha)\delta^2(m-1)^{-\sigma-\kappa}m^{-\sigma} \\ & \quad + 2(c(\alpha)\delta^2 m^{-\sigma-\kappa} - \frac{\alpha\lambda}{2}c(\alpha)\delta^2(m-1)^{-\sigma-\kappa}m^{-\sigma})^{1/2}\alpha m^{-\sigma}(\mathcal{L}^*)^{1/2}\frac{c_3}{T_m} + \mathcal{L}^*m^{-2\sigma}\frac{\alpha^2 c_3}{T_m} \\ & \le c(\alpha)\delta^2 m^{-\sigma-\kappa} + m^{-\sigma}(2\sqrt{c(\alpha)}\delta c_3 a(\mathcal{L}^*)^{1/2}Lm^{-\sigma-3\kappa/2} \\ & \quad + c_3\mathcal{L}^*\alpha^2 Lm^{-3\sigma/2-\kappa} - c(\alpha)\delta^2\alpha y(m-1)^{-\sigma-\kappa}) \\ & \le c(\alpha)\delta^2 m^{-\sigma-\kappa} + m^{-\sigma}(m-1)^{-\sigma-\kappa}(2\sqrt{c(\alpha)}\delta c_3 a(\mathcal{L}^*)^{1/2}L + c_3\mathcal{L}^*\alpha^2 L - c(\alpha)\delta^2\alpha y). \end{aligned} \qquad (135) $$

By the choice of $c(\alpha)$ in (134), for any $\kappa \ge 0$ we have the following inequality

$$ 2\sqrt{c(\alpha)}\delta c_5(\mathcal{L}^*)^{1/2}L + c_5\mathcal{L}^*aL - c(\alpha)\delta^2y < 0. \qquad (136) $$

Hence, with this choice of $c(\alpha)$, in (135) the latter term in the right-hand side is negative for any $m \ge 2$ and the induction step follows if $m > m_0(\alpha)$. That is, we have for some $c > 0$ that and when $m > m_0(\alpha)$ that

$$ \mathbb{E}[\mathrm{dist}(\Theta_m, \mathcal{M} \cap U)^2 \mathbf{1}[\mathcal{B}_{m-1}]] \le c \max\left(\frac{\delta^2}{a^{(1-\sigma)(\sigma+\kappa)}}, \frac{\mathcal{L}^*(1+\alpha\ell)}{\ell^2}\right) m^{-\sigma-\kappa}. \qquad (137) $$

We have left to show that the induction hypothesis holds in (124) for some $m$. Recall that $m > m_0(\alpha)$ is the only restriction we needed on the starting point for the induction argument to work—$\delta$ was already chosen depending on $\alpha$ in (133). From the choice

$$ m_0(\alpha) \ge \frac{c'}{\alpha^{1-\sigma}}, \qquad (138) $$

if $m \le m_0(\alpha)$, the following slightly changed version of (124) will hold; namely

$$ \mathbb{E}[(\mathrm{dist}(\Theta_m, \mathcal{M} \cap U)^2 \wedge \delta^2) \mathbf{1}[\mathcal{B}_{m-1}]] \le \delta^2 c(\alpha) m^{-\sigma-\kappa}. \qquad (139) $$

Hence, by same arguments conducted with (139) instead of (124), we have shown by induction that (139) holds for $m > 0$.
For convenience, we will further show that there exists a constant $c_6 > 0$ such that for all $m > 0$ we have

$$ \mathbb{E}[(\mathrm{dist}(\Theta_m, \mathcal{M} \cap U)^2 \wedge \delta^2)1[\mathcal{B}_{m-1}]] \le c_6 \mathcal{L}^* m^{-\sigma-\kappa}. \quad (140) $$

Fix $c_6 > 0$. Choose $\delta_0 \le \delta_1(\alpha)$ depending on $\alpha$ small enough and $\ell_0 > 0$ large enough such that for $\delta \in (0, \delta_0]$ and $\ell \in [\ell_0, \infty)$ we have that

$$ \begin{aligned} \frac{c' \delta^2}{\alpha(1-\sigma)(\sigma+\kappa)} &< c_6 \ge c_6 \mathcal{L}^*, \\ \frac{cD(1+\alpha\ell)}{\ell^2} &< c_6 \mathcal{L}^*, \end{aligned} \quad (141) $$

With the conditions in (141), the proof of the lemma follows noting that $\delta^2 c(\alpha) = \delta^2 c(\alpha, \ell) < c_6 \mathcal{L}^*$.

### C.5 Proof of Lemma 11

We will again use the notation that $t_{m+1} - t_m = T_m$ and without loss of generality assume that $T_m = \ell m^{\sigma/2+\kappa}$ as in Appendix C.4. The proof of Lemma 11 also mainly follows the steps of Fehrman et al. (2020). However, we again need to take care of the terms that the bias and lack of independence generate in the analysis.

The bounding starts noting the inequality

$$ \mathbb{E}\left[\max_{1 \le l \le m} |\Theta_l - \Theta_0| 1[\mathcal{B}_{l-1}]\right] \le \sum_{l=1}^{m} \mathbb{E}[|\Theta_l - \Theta_{l-1}|^2 1[\mathcal{B}_{l-1}]]^{1/2}. \quad (142) $$

We will show that there exists a constant $c > 0$ such that for $l \in [m]$ we have

$$ \mathbb{E}[|\Theta_{l+1} - \Theta_l|^2 1[\mathcal{B}_l]]^{1/2} \le c\alpha \left( l^{-3/2\sigma-\kappa/2} + \sqrt{\frac{1}{\ell}} l^{-5\sigma/8-\kappa/2} \right), \quad (143) $$

where the exponents of $\sigma$ and $\kappa$ already differ from the result of Fehrman et al. (2020), and are required to account for the lack of independence and bias. Following the steps from Fehrman et al. (2020), in the neighborhood $V_{\kappa,\delta}(\theta^*)$, for each $l \le m$ there is a random variable $\epsilon_l : \mathcal{B}_l \to \mathbb{R}^n$ and there exists a constant $c > 0$ such that

$$ |\epsilon_l| < c \mathrm{dist}(\Theta_l, \mathcal{M} \cap U)^2, \quad (144) $$

and such that on the event $\mathcal{B}_l$ we have

$$ \nabla J(\Theta_l) = \mathrm{Hess}_{\mathfrak{p}(\Theta_l)}(\Theta_l - \mathfrak{p}(\Theta_l)) + \epsilon_l. \quad (145) $$

Recalling the definition of $\eta_l$ in (43), we have then the equality

$$ \Theta_{l+1} = \Theta_l - \alpha_l \mathrm{Hess}_{\mathfrak{p}(\Theta_l)}(\Theta_l - \mathfrak{p}(\Theta_l)) - \alpha_l \epsilon_l + \alpha_l \eta_l. \quad (146) $$

Define

$$ \tilde{\Theta}_l = \Theta_l - \alpha_l \mathrm{Hess}_{\mathfrak{p}(\Theta_l)}(\Theta_l - \mathfrak{p}(\Theta_l)). \quad (147) $$
We use the triangle inequality with in (146) separating $\Theta_{l+1} - \Theta_l$ as the summands of $\Theta_{l+1} - \tilde{\Theta}_l$ and $\tilde{\Theta}_l - \Theta_l$.

We estimate first $|\Theta_{l+1} - \tilde{\Theta}_l|^2$. In our case, after expanding $\mathbb{E}[|\Theta_{l+1} - \tilde{\Theta}_l|^2 1[\mathcal{B}_l]|$, we diverge from (Fehrman et al., 2020, Equation 58) and we need to bound

$$ \alpha_l^2 \mathbb{E}[1[\mathcal{B}_l] \langle \epsilon_l, \eta_l \rangle]. \qquad (148) $$

Similar to the proof of Lemma 10, we can condition on $\mathcal{F}_l$ and using that $\epsilon_l$ and $\mathcal{B}_l$ are $\mathcal{F}_l$-measurable together with the Cauchy-Schwartz inequality, we have

$$ \begin{aligned} \alpha_l^2 \mathbb{E}[1[\mathcal{B}_l] \langle \epsilon_l, \eta_l \rangle] &\le \alpha_l^2 \mathbb{E}[\langle 1[\mathcal{B}_l] \epsilon_l, \mathbb{E}[\eta_l 1[\mathcal{B}_l] | \mathcal{F}_l] \rangle] \\ &\le \alpha_l^2 \mathbb{E}[1[\mathcal{B}_l] | \epsilon_l|^2]^{1/2} \mathbb{E}[|\mathbb{E}[\eta_l 1[\mathcal{B}_l] | \mathcal{F}_l] |^2]^{1/2}. \end{aligned} \qquad (149) $$

Since $1[\mathcal{B}_m] \le 1[\mathcal{B}_{m-1}]$, we can bound

$$ \mathbb{E}[|\mathbb{E}[\eta_l 1[\mathcal{B}_l] | \mathcal{F}_l]|^2]^{1/2} \stackrel{\text{Lemma 8}}{\le} \mathbb{E}[1[\mathcal{B}_l] \frac{c_1^2}{T_l^2} \mathcal{L}_v(S_{t_l})]^{1/2} \stackrel{\text{Lemma 14}}{\le} \frac{c_2}{T_l}. \qquad (150) $$

For the remaining term in (149), recall that on the event $\mathcal{B}_l$, since $\Theta_l \in V_{\tau, \delta}(\theta^*)$, we have that $\text{dist}(\Theta_l, \mathcal{M} \cap U) \le \delta$. Hence, we can bound for any $l > 0$ that

$$ \begin{aligned} \mathbb{E}[1[\mathcal{B}_l] | \epsilon_l|^2]^{1/2} &\stackrel{(144)}{\le} (\alpha_l)^2 \mathbb{E}[\text{dist}(\Theta_l, \mathcal{M} \cap U)^4 1[\mathcal{B}_l]]^{1/2} \frac{c_3}{T_{l+1}} \\ &\le (\alpha_l)^2 \delta^2 \mathbb{E}[\text{dist}(\Theta_l, \mathcal{M} \cap U)^2 1[\mathcal{B}_l]]^{1/2} \frac{c_3}{T_{l+1}} \\ &\le (\alpha_l)^2 \delta^2 \mathbb{E}[\text{dist}(\Theta_l, \mathcal{M} \cap U)^2 1[\mathcal{B}_{l-1}]]^{1/2} \frac{c_3}{T_{l+1}} \\ &\stackrel{\text{Lemma 10}}{\le} (\alpha_l)^2 \delta^2 l^{-\sigma/2 - \kappa/2} \frac{c_4}{T_l}. \end{aligned} \qquad (151) $$

The estimation of the remaining terms in the expansion of $\mathbb{E}[|\Theta_l - \tilde{\Theta}_{l-1}|^2 1[\mathcal{B}_{l-1}]]$ can be conducted in the same way as that in Fehrman et al. (2020), to which we refer for the details to the interested reader. Together with the estimate of (151) that accounts for the biases we have that

$$ \begin{aligned} \mathbb{E}[|\Theta_l - \tilde{\Theta}_{l-1}|^2 1[\mathcal{B}_l]] &\le c_5(\alpha_l)^2 \delta^2 \mathbb{E}[\text{dist}(\Theta_l, \mathcal{M} \cap U)^2 1[\mathcal{B}_l]] \\ &\quad + 2\delta \mathbb{E}[\text{dist}(\Theta_l, \mathcal{M} \cap U)^2 1[\mathcal{B}_l]]^{1/2} \frac{c_6}{T_l} + (\alpha_l)^2 \frac{c_7}{T_l} \\ &\le c_8(\alpha_l)^2 \left[ \delta^2 l^{-\sigma-\kappa} + 2\delta l^{-\sigma/2-\kappa/2} \frac{1}{T_l} + \frac{1}{T_l} \right]. \end{aligned} \qquad (152) $$

Substituting $T_l = t_{l+1} - t_l = l^{\kappa+\sigma/2\ell}$ and using $\alpha_l < \alpha_{l-1} = \alpha l^{-\sigma}$ into (152) yields the bound

$$ \mathbb{E}[|\Theta_l - \tilde{\Theta}_{l-1}|^2 1[\mathcal{B}_{l-1}]] \le c_9 \frac{\alpha^2}{l^{2\sigma}} \left( \delta^2 \frac{1}{l^{\sigma+\kappa}} + 2\delta \frac{1}{l^{\sigma+3\kappa/2\ell}} + \frac{1}{l^{\kappa+\sigma/2}} \right) $$
$$ \leq c_{10} \frac{\alpha^2}{l^{5\sigma/4+\kappa\ell}}, \qquad (153) $$

where in the last inequality we have taken the term with the highest order. Using the previous bounds from Lemma 10 we can show that

$$ \mathbb{E}[|\Theta_l - \bar{\Theta}_l|^2 \mathbf{1}[\mathcal{B}_l]] \leq \alpha_l^2 \mathbb{E}[\mathrm{dist}(\Theta_l, \mathcal{M} \cap U) \mathbf{1}[\mathcal{B}_l]] \leq c_{11} \frac{a^2}{l^{3\sigma+\kappa}}, \qquad (154) $$

so that using the triangle inequality and combining the bounds of (153) and (154) we obtain

$$ \mathbb{E}[|\Theta_{l+1} - \Theta_l|^2 \mathbf{1}[\mathcal{B}_l]]^{1/2} \leq c_{12}\alpha(l^{-3/2\sigma-\kappa/2} + \sqrt{\ell}^{-1}l^{-5\sigma/8-\kappa/2}). \qquad (155) $$

Hence, since $\sigma \in (2/3, 1)$ adding the bound (155) in (142) yields

$$ \begin{aligned} \mathbb{E}\left[\max_{1 \le l \le m} |\Theta_l - \Theta_0| \mathbf{1}[\mathcal{B}_{l-1}]\right] &\le \sum_{l=1}^{m} c_{12}\alpha(l^{-3/2\sigma-\kappa/2} + \sqrt{\ell}^{-1}l^{-\sigma-\kappa/2}) \\ &\le c_{13}\alpha(m^{1-3/2\sigma-\kappa/2} + \sqrt{\ell}^{-1}m^{1-5\sigma/8-\kappa/2}). \end{aligned} $$

## C.6 Proof of Lemma 12

The proof mimicks the proof strategy of (Fehrman et al., 2020, Proposition 24), but modifications are required due to our Markovian assumptions and appearances of biases. Specifically, we must carefully consider the adverse effects that these biases could have on the probability that the iterates exit the basin of attraction. Concretely, our effort will go into firstly proving the following sufficiently strong analogue of (Fehrman et al., 2020, Equation 75) that is applicable to our problem:

**Lemma 15** There exist constants $c_1, c_2 > 0$ such that

$$ \mathbb{P}[\mathrm{dist}(\Theta_m, \mathcal{M} \cap U) > \delta, \mathcal{B}_{m-1}] \le \frac{c_1 \alpha^2}{\delta^2 \ell m^{2\sigma}} \mathbb{P}[\mathcal{B}_{m-1}] + \frac{c_2}{\delta^4 \ell m^{\sigma+\kappa}}. \qquad (156) $$

The proof of Lemma 15 can be found in Appendix C.6.2.

Once Lemma 15 has been established, we secondly estimate the combined probability that any of the iterates escape in directions tangential to the manifold. The proof of this fact, which is analogous to (Fehrman et al., 2020, Equation 78-79), can be found in Appendix C.6.3.

**Lemma 16** If $\Theta_0 \in V_{\tau/2, \delta}(\theta^*)$, then

$$ \sum_{l=1}^{m} \mathbb{P}[\mathrm{dist}(\Theta_l, \mathcal{M} \cap U) < \delta, \Theta_l \notin V_{\tau, \delta}(\theta^*), \mathcal{B}_{l-1}] \le \mathbb{P}\left[\max_{1 \le l \le m} |\Theta_l - \Theta_0| \mathbf{1}[\mathcal{B}_{l-1}] > R/2 - 2\delta, \right. \qquad (157) $$
### C.6.1 PROOF THAT LEMMAS 15 AND 16 IMPLY LEMMA 12

First, note that the recursion

$$ \mathbb{P}[\mathcal{B}_m] = \mathbb{P}[\Theta_m \in V_{\mathfrak{r},\delta}(\theta^*), \mathcal{B}_{m-1}] = \mathbb{P}[\mathcal{B}_{m-1}] - \mathbb{P}[\Theta_m \notin V_{\mathfrak{r},\delta}(\theta^*), \mathcal{B}_{m-1}] \quad (158) $$

can be iterated whenever we can control and bound the following probabilities

$$ \mathbb{P}[\Theta_m \notin V_{\mathfrak{r},\delta}(\theta^*), \mathcal{B}_{m-1}] = \mathbb{P}[\mathrm{dist}(\Theta_m, \mathcal{M} \cap U) > \delta, \mathcal{B}_{m-1}] \\ + \mathbb{P}[\mathrm{dist}(\Theta_m, \mathcal{M} \cap U) \le \delta, \Theta_m \notin V_{\mathfrak{r},\delta}(\theta^*), \mathcal{B}_{m-1}]. \quad (159) $$

Using Lemma 15 and induction on (158) and (159), it follows that for some $c > 0$,

$$ \mathbb{P}[\mathcal{B}_m] \ge \prod_{l=1}^{m} \left(1 - \frac{c\alpha^2}{\delta^2 \ell l^2 \sigma}\right)_+ - \sum_{l=1}^{m} \frac{c}{\ell \delta^4 l \sigma + \kappa} - \sum_{l=1}^{m} \mathbb{P}[\mathrm{dist}(\Theta_l, \mathcal{M} \cap U) < \delta, \Theta_l \notin V_{\mathfrak{r},\delta}(\theta^*), \mathcal{B}_{l-1}]. \quad (160) $$

We use Lemma 16 together with Lemma 11 and Markov's inequality to obtain the bound

$$ \sum_{l=1}^{m} \mathbb{P}[\mathrm{dist}(\Theta_l, \mathcal{M} \cap U) < \delta, \Theta_l \notin V_{\mathfrak{r},\delta}(\theta^*), \mathcal{B}_{l-1}] \le c\alpha \frac{(m^{1-3/2\sigma-\kappa/2} + \ell^{-1/2}m^{1-5\sigma/8-\kappa/2})}{(\mathfrak{r}/2 - 2\delta)_+}. \quad (161) $$

Thus, substituting (161) in (160), for some $c > 0$ we have

$$ \mathbb{P}[\mathcal{B}_m] \ge \prod_{l=1}^{m} \left(1 - \frac{c\alpha^2}{\delta^2 \ell l^2 \sigma}\right)_+ - \sum_{l=1}^{m} \frac{c}{\ell \delta^4 l \sigma + \kappa} - c\alpha \frac{(m^{1-3/2\sigma-\kappa/2} + \ell^{-1/2}m^{1-5\sigma/8-\kappa/2})}{(\mathfrak{r}/2 - 2\delta)_+}. \quad (162) $$

Note first that since $\sigma \in (2/3, 1)$ and $\kappa \ge 0$, if $\sigma + \kappa \ne 1$, then there exists a constant $c_1 > 0$ such that

$$ \sum_{l=1}^{m} \frac{c}{\ell \delta^4 l \sigma + \kappa} \le c_1 m^{1-\sigma-\kappa}. \quad (163) $$

Lastly, there also exists a constant $c > 0$, $\alpha_0 > 0$, $\delta_0$ such that if $\alpha \in (0, \alpha_0]$ and $\delta \in (0, \delta_0]$ then there exists $\ell_0 > 0$ such that if $\ell \in [\ell_0, \infty)$ then

$$ \prod_{l=1}^{m} \left(1 - \frac{c\alpha^2}{\delta^2 \ell l^2 \sigma}\right)_+ \ge \exp\left(-\frac{c\alpha^2}{\delta^2 \ell}\right). \quad (164) $$

Lower bounding (162) using (163) and (164) yields Lemma 12. □

### C.6.2 PROOF OF LEMMA 15

We follow first (Fehrman et al., 2020, Equation 69), by fixing $\delta_1$ small enough such that $\delta \in (0, \delta_1]$, on the event $\mathcal{B}_{m-1}$ it is shown by Fehrman et al. (2020) that we have the inequality

$$ \mathrm{dist}(\Theta_m, \mathcal{M} \cap U) \le \left(1 - \frac{\lambda \alpha_{m-1}}{2}\right) \mathrm{dist}(\Theta_{m-1}, \mathcal{M} \cap U) + \alpha_{m-1} |\eta_{m-1}|. \quad (165) $$
We consider now the event {$\text{dist}(\Theta_m, \mathcal{M} \cap U) > \delta$} $\cap \mathcal{B}_{m-1}$. This event occurs when in (165), either $\Theta_{m-1} \in V_{\text{r},\delta/2}(\theta^*)$ and $|\eta_{m-1}| \ge \alpha_{m-1}\delta/2$, or $\Theta_{m-1} \in V_{\text{r},\delta}(\theta^*) \setminus V_{\text{r},\delta/2}(\theta^*)$ and the gradient term can have smaller size. Mathematically, this translates into the inequality

$$ \mathbb{P}[\text{dist}(\Theta_m, \mathcal{M} \cap U) > \delta, \mathcal{B}_{m-1}] \le \mathbb{P}\left[\left|\eta_{m-1}\right| \ge \frac{\delta}{2\alpha_{m-1}}, \Theta_{m-1} \in V_{\text{r},\delta/2}(\theta^*), \mathcal{B}_{m-2}\right] \quad (166) \\ + \mathbb{P}\left[\left|\eta_{m-1}\right| \ge \frac{\delta\lambda}{2}, \Theta_{m-1} \in V_{\text{r},\delta}(\theta^*) \setminus V_{\text{r},\delta/2}(\theta^*), \mathcal{B}_{m-2}\right] =: P_1 + P_2. $$

Contrary to what is done in the proof of (Fehrman et al., 2020, Proposition 24), we cannot use an independence property to estimate the probabilities $P_1$ and $P_2$ in (166). After all, the Markov chain's behavior at epoch $m-1$ depends on $\Theta_{m-1}$.

In order to overcome this issue we will use the characterization of $\eta_{m-1}$ in Lemma 8. Recall Lemma 8, and note that it implies

$$ \begin{aligned} \mathbb{E}[\mathbb{1}[\mathcal{B}_{m-1}]\mathbb{1}[|\eta_{m-1}| \ge \frac{\delta}{2\alpha_{m-1}}]|\mathcal{F}_{m-1}] &= \mathbb{P}[|\eta_{m-1}| \ge \frac{\delta}{2\alpha_{m-1}}, \mathcal{B}_{m-1}|\mathcal{F}_{m-1}] \\ &\le \frac{\mathbb{E}[|\eta_{m-1}|^2\mathbb{1}[\mathcal{B}_{m-1}]|\mathcal{F}_{m-1}]}{\frac{\delta^2}{4(\alpha_{m-1})^2}} \le \frac{4c_2(\alpha_{m-1})^2\mathcal{L}_4(S_{t_{m-1}})}{\delta^2 T_m} \end{aligned} \quad (167) $$

since there exist a constant $c > 0$ such that $T_m < cT_{m-1}$.

Let us first bound $P_1$ in (166). We can write

$$ \begin{aligned} P_1 \stackrel{(i)}{=} & \mathbb{E}[\mathbb{1}[|\eta_{m-1}| \ge \frac{\delta}{2\alpha_{m-1}}]\mathbb{1}[\Theta_{m-1} \in V_{\text{r},\delta/2}(\theta^*)]\mathbb{1}[\mathcal{B}_{m-2}]\mathbb{1}[\mathcal{B}_{m-1}]] \\ &= \mathbb{E}[\mathbb{1}[\Theta_{m-1} \in V_{\text{r},\delta/2}(\theta^*)]\mathbb{1}[\mathcal{B}_{m-2}]\mathbb{E}[\mathbb{1}[\mathcal{B}_{m-1}]\mathbb{1}[|\eta_{m-1}| \ge \frac{\delta}{2\alpha_{m-1}}]|\mathcal{F}_{m-1}]] \\ &\stackrel{(167)}{\le} \frac{4c_2(\alpha_{m-1})^2}{T_m \delta^2} \mathbb{E}[\mathbb{1}[\Theta_{m-1} \in V_{\text{r},\delta/2}(\theta^*)]\mathbb{1}[\mathcal{B}_{m-2}]\mathcal{L}_4(S_{t_{m-1}})], \end{aligned} \quad (168) $$

where for (i) we have used the fact that {$\Theta_{m-1} \in V_{\text{r},\delta/2}(\theta^*)$} $\cap \mathcal{B}_{m-2} \subset \mathcal{B}_{m-1}$.

We deal now with the remaining term in (168). Differently to the independent and unbiased case we need to control the bias and use the tail probability that the Lyapunov function is larger than a certain bound in order to estimate the deviation probability. This step is the crucial different step compared to Fehrman et al. (2020), where we have to explicitly use Assumption 4 and 6. Note that a Cauchy-Schwartz inequality in (168) will not yields an inequality strong enough. See the remark after the proof for further details.

Before bounding the remaining term in (168), we obtain the necessary inequalities. Recall from Lemma 14 that since $\mathbb{E}[\mathcal{L}_4(S_{t_{m-1}})^4\mathbb{1}[\mathcal{B}_{m-2}]] < \mathbb{E}[\mathcal{L}_v(S_{t_{m-1}})\mathbb{1}[\mathcal{B}_{m-2}]] < D < \infty$, then by Markov's inequality we have that there exists $D > 0$ such that for any $m > 0$,

$$ \mathbb{P}[\mathcal{L}(S_{t_{m-1}}) > m^s, \mathcal{B}_{m-2}] \le D^4 m^{-4s}. \quad (169) $$

Note also that under the moment assumptions the following holds

$$ \mathbb{E}[\mathcal{L}(S_{t_{m-1}})\mathbb{1}[\mathcal{B}_{m-2}]\mathbb{1}[\mathcal{L}(S_{t_{m-1}}) > m^s]] = \int_{m^s}^{\infty} \mathbb{P}[\mathcal{L}(S_{t_{m-1}}) > t, \mathcal{B}_{m-2}] dt $$
$$ = \int_{m^s}^{\infty} \frac{D^4}{t^4} dt \leq D^4 m^{-3s+1}. \quad (170) $$

We use the (170) to bound (168) as follows

$$
\begin{align*}
& \mathbb{E}[\mathbf{1}[\Theta_{m-1} \in V_{\mathrm{r},\delta/2}(\theta^*)] \mathcal{L}_4(S_{t_{m-1}}) \mathbf{1}[\mathcal{B}_{m-2}]] \\
& \leq \mathbb{E}[\mathbf{1}[\Theta_{m-1} \in V_{\mathrm{r},\delta/2}(\theta^*)] \mathcal{L}_4(S_{t_{m-1}}) \mathbf{1}[\mathcal{B}_{m-2}](\mathbf{1}[\mathcal{L}_4(S_{t_{m-1}}) > m^s] + \mathbf{1}[\mathcal{L}_4(S_{t_{m-1}}) \leq m^s])] \\
& \leq \mathbb{E}[\mathbf{1}[\Theta_{m-1} \in V_{\mathrm{r},\delta/2}(\theta^*)] m^s \mathbf{1}[\mathcal{B}_{m-2}]] + \mathbb{E}[\mathcal{L}(S_{t_{m-1}}) \mathbf{1}[\mathcal{B}_{m-2}] \mathbf{1}[\mathcal{L}(S_{t_{m-1}}) > m^s]] \tag{170} \\
& \leq m^s \mathbb{P}[\Theta_{m-1} \in V_{\mathrm{r},\delta/2}(\theta^*), \mathcal{B}_{m-2}] + c_3 D m^{-3s+1} \leq m^s \mathbb{P}[\mathcal{B}_{m-1}] + c_3 D m^{-3s+1}. \tag{171}
\end{align*}
$$

Thus, using (171), we can bound $P_1$ in (166). Specifically,

$$ P_1 \leq \frac{4c_4(\alpha_{m-1})^2}{T_m \delta^2} (m^s \mathbb{P}[\mathcal{B}_{m-1}] + m^{-3s+1}). \quad (172) $$

This completes our bound for $P_1$.

We now bound $P_2$ in (166). Repeating the argumentation behind (172), we can show that

$$ P_2 \leq \frac{4c_5}{T_m \lambda^2 \delta^2} \left( m^s \mathbb{P}[\Theta_{m-1} \in V_{\mathrm{r},\delta}(\theta^*) \setminus V_{\mathrm{r},\delta/2}(\theta^*), \mathcal{B}_{m-2}] + m^{-3s+1} \right). \quad (173) $$

Using the facts (i) $\{\Theta_{m-1} \in V_{\mathrm{r},\delta}(\theta^*) \setminus V_{\mathrm{r},\delta/2}(\theta^*)\} \subseteq \{\mathrm{dist}(\Theta_{m-1}, \mathcal{M} \cap U) \geq \delta/2\}$, with (ii) an application of Lemma 10 and Markov's inequality, reveals that

$$ \mathbb{P}[\Theta_{m-1} \in V_{\mathrm{r},\delta}(\theta^*) \setminus V_{\mathrm{r},\delta/2}(\theta^*), \mathcal{B}_{m-2}] \stackrel{(i)}{\leq} \mathbb{P}[\mathrm{dist}(\Theta_{m-1}, \mathcal{M} \cap U) \geq \frac{\delta}{2}, \mathcal{B}_{m-2}] \stackrel{(ii)}{\leq} \frac{4}{\delta^2} c_6 m^{-\sigma-\kappa}. \quad (174) $$

Applying the bound in (173) to (174) yields

$$ P_2 \leq \frac{4c_7}{T_m \lambda^2 \delta^4} \left( m^s m^{-\sigma-\kappa} + m^{-3s+1} \right). \quad (175) $$

This completes the bound for $P_2$ in (166).

Lastly, we return to (166) and select the parameter. Let us now combine (171) and (175) and return to bounding the left-hand side of (166). Specifically, observe that we proved that

$$ \mathbb{P}[\mathrm{dist}(\Theta_m, \mathcal{M} \cap U) > \delta, \mathcal{B}_{m-1}] \leq \frac{4c_8(\alpha_{m-1})^2}{T_m \delta^2} (m^s \mathbb{P}[\mathcal{B}_{m-1}] + m^{-3s+1}) \\ + \frac{4c_9}{T_m \delta^4} (m^{s-\sigma-\kappa} + m^{-3s+1}). \quad (176) $$

We now specify $s = \kappa + \sigma/2$ in (176). Without loss of generality we will again assume that $T_m = \ell m^{\sigma/2+\kappa}$ instead of $\lfloor \ell m^{\sigma/2+\kappa} \rfloor$—there is namely only a constant changed. By choosing the smallest exponents in $m$ in (176) for all $m > 0$ we have

$$ \mathbb{P}[\mathrm{dist}(\Theta_m, \mathcal{M} \cap U) > \delta, \mathcal{B}_{m-1}] \leq c_{10} \frac{a^2}{\delta^2 \ell m^{2\sigma}} \mathbb{P}[\mathcal{B}_{m-1}] + \frac{c_{10}}{\delta^4 \ell} (m^{-3\sigma-4\kappa+1} + m^{-\sigma-\kappa}). \quad (177) $$

Since $\sigma \in (2/3, 1)$, then $-3\sigma - 4\kappa + 1 < -\sigma - \kappa$ for any $\kappa \geq 0$. Upper bounding the leading orders in $m$ completes the proof of Lemma 15. $\square$
**Remark 17** A Cauchy-Schwartz inequality in (168) would only yield a factor $P[\mathcal{B}_{m-1}]^{1/2} > P[\mathcal{B}_{m-1}]$, which would not be sufficient. Similarly, we could have used Lemma 14 directly and obtain a bound on $E[1[\mathcal{B}_{m-2}]\mathcal{L}_4(S_{t_{m-1}})]$. However, this would not give an inequality that can be iterated inductively and is sharp enough. We can directly simplify this term to obtain $P(\mathcal{B}_{m-1})$ in the inequality only when $\mathcal{L}_4(S_{t_{m-1}})$ is bounded.

### C.6.3 PROOF OF LEMMA 16

In the work of Fehrman et al. (2020), it is (Fehrman et al., 2020, Lemma 23) that establishes (Fehrman et al., 2020, Equations 78-79) directly. Since (Fehrman et al., 2020, Lemma 23) is solely a geometric argument, and does not concern the stochastic process, it also applies in our Markovian setting. □

## Appendix D. The Compact Case

In the case that the set of maxima $\mathcal{M}$ is compact, we can improve the convergence rate of Theorem 2. We will namely assume the following

**Assumption 8 (Compactness, Optional)** The open subset $U$ defined in Assumption 7 is such that $\mathcal{M} \cap U$ is compact.

Under this additional assumption we have the following

**Theorem 18 (Compact Case)** Suppose that Assumptions 1 to 8 hold, except that (17) is now relaxed to allow for $\sigma \in (0, 1)$ and $\kappa \in [0, \infty)$. For every maximizer $\theta^* \in \mathcal{M}$, there exist constants $c > 0$ and $\alpha_0 > 0$ such that, for every $\alpha \in (0, \alpha_0]$, there exists a neighborhood $V$ of $\theta^*$ such that there exists $\ell_0 > 0$ such that for any $\ell \in [\ell_0, \infty)$, $m \in \mathbb{N}_+$, and $\epsilon \in (0, 1)$,

$$P[J(\Theta_m) < J^* - \epsilon | \Theta_0 \in V] \le c \left( \epsilon^{-2} m^{-\sigma-\kappa} + \frac{m^{1-\sigma-\kappa}}{\ell} + \frac{\alpha^2}{\ell} \right). \quad (178)$$

The term proportional to $\alpha m^{-\kappa/2} + \alpha m^{1-\sigma/2-\kappa/2} \ell^{-1/2}$ is not in Theorem 18 compared to Theorem 2. This term estimates the probability that the iterates escape $V$ along directions almost parallel to those of $\mathcal{M}$. As it turns out, in the compact case such event cannot occur. The bound in (178) thus holds when the set of maxima is, for example, a singleton $\mathcal{M} \cap U = \{x_0\}$.

### D.1 Proof of Theorem 18

The proof is the same as with Theorem 2, but we can omit the last term in (84) by showing that we can choose $\mathfrak{r}$ arbitrarily large. The argument is as follows. If the manifold $\mathcal{M} \cap U$ is compact, it can be covered by a finite number $k$ of local tubular neighborhoods $V_i = V_{\mathfrak{r}_i, \delta_i}(\theta_i)$ where $\theta_i \in \mathcal{M} \cap U$ and $\mathcal{M} \cap U \subset \cup_{i \in [k]} V_i$. Choose $\delta = \min_{i \in [k]} \delta_i$. Then, any $\theta \in U$ such that $\text{dist}(\theta, \mathcal{M} \cap U) < \delta$ will satisfy that $\mathfrak{p}(\theta) \in \mathcal{M} \cap U$, where $\mathfrak{p}$ is the unique local orthogonal projection on $\mathcal{M} \cap U$ from (39). Now, from compactness, for any $\theta^* \in \mathcal{M} \cap U$ there exists $\tilde{\mathfrak{r}} > 0$ such that $\mathcal{M} \cap U \subset B_{\tilde{\mathfrak{r}}}(\theta^*)$. For any $\mathfrak{r} \ge \tilde{\mathfrak{r}}$ we thus have that $V_{\mathfrak{r}, \delta}(\theta^*) = V_{\tilde{\mathfrak{r}}, \delta}(\theta^*)$ is a tubular neighborhood containing $\mathcal{M} \cap U$. Then, we can choose $\mathfrak{r}$ arbitrarily large and conclude that the last term in the bound for the probability in Theorem 2 vanishes if