We consider the discrete-time infinite-horizon average-reward restless bandit problem.
We propose a novel policy that maintains two dynamic subsets of arms:
one subset of arms has a nearly optimal state distribution and takes actions according to an Optimal Local Control routine; the other subset of arms is driven towards the optimal state distribution and gradually merged into the first subset.
We show that our policy is asymptotically optimal with an $O(\exp(-C N))$ optimality gap for an $N$-armed problem, under the mild assumptions of aperiodic-unichain, non-degeneracy, and local stability.
Our policy is the first to achieve exponential asymptotic optimality under the above set of easy-to-verify assumptions, whereas prior work either requires a strong Global Attractor assumption or only achieves an $O(1/\sqrt{N})$ optimality gap.
We further discuss the fundamental obstacles in significantly weakening our assumptions.
In particular, we prove a lower bound showing that local stability is fundamental for exponential asymptotic optimality.

§ INTRODUCTION

Restless bandits (RBs) [15] is a stochastic sequential decision-making problem consisting of multiple arms coupled by a constraint.
Each arm operates as a Markov decision process (MDP) with two possible actions: activating or keeping passive.
At each time step, the decision maker observes the states of the arms and chooses a fixed fraction of arms to activate, aiming to maximize the expected reward collected from the arms.
The RB problem has a long history and diverse applications. Readers can refer to the recent survey paper [11] for a comprehensive review of the literature.

In this paper, we study the infinite-horizon RBs in discrete time, with the long-run average-reward objective.
Since computing the exact optimal policy for an RB problem is intractable as the number of arms $N$ is large [12], we focus on developing asymptotically optimal policies that can be efficiently computed.
Specifically, the optimality gap of a policy is the difference between the average reward per arm and that of an optimal policy.
A policy is considered to be asymptotically optimal if the optimality gap is $o(1)$ in the asymptotic limit of $N\to\infty$.[In this paper, we use standard Bachmann–Landau notation and focus on the asymptotics with respect to $N$.]

The asymptotic optimality of average-reward RB in the large $N$ regime has been studied for decades, starting with the seminal papers on the renowned Whittle index policy [15, 16].
Since then researchers have been weakening the assumptions for asymptotic optimality [14, 9, 10, 17] or improving the order of optimality gaps [3, 4].
A recent milestone is the exponential asymptotic optimality established in [3, 4], where a class of policies called LP-Priority is proved to achieve an optimality gap  $O(\exp(-CN))$ for some constant $C$, under some assumptions.
Such an exponential gap has significant practical value, since it implies that a policy has nearly-optimal performance even for RBs with small $N$. Moreover, an exponential gap also “beats” the Central Limit Theorem (CLT), an intriguing theoretical property that is not a priori obvious to be possible.
In particular, on a high-level all asymptotic optimality results in the RB literature are based on the concentration of the empirical distribution of the arms' states around a certain optimal distribution. Therefore, an $O(1/\sqrt{N})$ bound owing to CLT was believed to be fundamental.
An important observation made in [3, 4] is that if the system has a certain notion of local linearity,  the optimality gap only depends on the distance between the optimal distribution and the expected empirical state distribution, allowing one to go beyond CLT.

Despite the above intuition,
it remains unclear what is the fundamental mechanism leading to exponential asymptotic optimality. Although a simple assumption called non-degeneracy (or non-singularity) suffices to ensure the local linearity property described above, non-degeneracy alone is not enough to guarantee exponential optimality of the LP-Priority policies considered in [3, 4].
In particular, another crucial assumption is needed, namely the Uniform Global Attractor Property (UGAP). Without UGAP, LP-Priority may even have a constant optimality gap; see [16, 2, 9] for concrete examples.
Unfortunately, UGAP is hard to interpret and verify, as it concerns the global convergence of a certain non-linear difference equation.

In light of the difficulty caused by UGAP, very recent work has studied new policies that are asymptotically optimal without assuming UGAP [9, 10, 17].
Different from LP-Priority and Whittle index policies, these new policies actively control the empirical distribution of the states of the arms, driving the distribution to globally converge towards a certain optimal distribution, without relying on extraneous assumptions.
However, while these results conclude that UGAP is not needed for asymptotic optimality, it remains unclear whether UGAP is fundamental for exponential asymptotic optimality — the best optimality gap proved in this line of work is $O(1/\sqrt{N})$.

Given these developments, it is natural to ask if removing UGAP necessarily results in a sacrifice in the optimality gap.
Is it possible to efficiently find policies that achieve global exponential asymptotic optimality, i.e., an $O(\exp(-CN))$ optimality gap without assuming global attraction?

*Our contributions.
The primary contribution of the paper is to design the first policy, the Two-Set Policy, that achieves exponential asymptotic optimality without UGAP.
The Two-Set Policy maintains two dynamic subsets of arms and applies two different subroutines to the subsets:
one subset of arms has a nearly optimal state distribution and takes actions according to a certain Optimal Local Control; the other subset of arms is driven towards the optimal state distribution and gradually merged into the first subset.
Intuitively, the Optimal Local Control in the first subset induces local linearity, while the second subset achieves global convergence.
While the ideas in each subroutine is not completely new, the Two-Set Policy leverages them in a novel way to achieve the best of both worlds.
The theoretical guarantees of the Two-Set Policy is informally summarized in the theorem below.

[Informal version of Theorem <ref>]
Suppose each arm is unichain, aperiodic, non-degenerate and locally stable. Then we can efficiently compute a $N$-armed policy $\pi$ such that

\begin{align}
\rrel - \rsysn \le O(\exp(-CN)),
\end{align}

where $\rrel$ is an (efficiently computable) upper bound on the optimal long-run average reward, and $\rsysn$ is the long-run average reward under the policy $\pi$ given the vector of initial states $\veS_0$.

We will make precise the definitions of non-degeneracy and local stability. We emphasize that we only need aperiodic unichain under a particular policy rather than all policies.
The assumptions of Theorem <ref> are a subset of, and hence strictly weaker than, the assumptions made in prior work [3, 4] for exponential asymptotic optimality. Moreover, all our assumptions pertain to a linear program defined by problem primitives and are thus much easier to verify than UGAP.

The proof of the above theorem employs a novel multivariate Lyapunov function technique, which generalizes the focus-set approach in [10]. This multivariate Lyapunov function enables us to decouple the complex dynamics under the Two-Set Policy, where the states of the arms and the two dynamic subsets are coupled and change simultaneously.

We next elucidate the roles of our assumptions for exponential optimality.
The importance of unichain, aperiodicity and non-degeneracy is well recognized; we complement the discussion with concrete examples.
We further  establish the fundamental importance of local stability by proving the following.

[Informal version of Theorem <ref>]
For regular unstable RBs, every policy $\pi$ satisfies

\begin{align}
\rrel - \rsysn \ge \Omega(1/\sqrt{N}).
\end{align}

We will formally define the class of regular (locally) unstable RB and provide concrete examples. Note that all known efficient policies rely on the upper bound $\rrel$.
For these policies, <Ref> shows that local stability is crucial for achieving exponential asymptotic optimality for a natural class of RB instances.
The proof of <Ref> uses Lyapunov analysis to reveal an interesting trade-off between two competing needs when local stability fails: to keep the states of the arms near the optimal distribution, and to take the optimal actions when the states are near the optimal distribution.

Other related work. Apart from the average-reward RBs, there are also a rich body of literature on RBs with finite-horizon total-reward [8, 18, 1, 6, 19, 4, 5] or or the infinite-horizon discounted-reward criteria [20, 6].
The approaches in these setting do not directly apply to the average-reward setting, as their computational complexities scale with the (effective) time horizon. However, there are still a large intersection of ideas — in particular, the first $O(1/N)$ optimality gap result was established in the finite-horizon setting by [19].

§ PROBLEM SETUP

Consider the discrete-time average-reward restless bandit problem with $N$ homogeneous arms, referred to as the $N$-armed problem. Index the arms by $[N] \triangleq \{1,2,\dots, N\}.$
To avoid confusion with the Whittle index, We refer to $i$ as  the ID of Arm $i$. Each arm is associated with a Markov decision process (MDP) defined by  $(\sspa, \aspa, P, r),$ which is called the single-armed MDP.
Here $\sspa$ is a finite state space;
$\aspa = \{0, 1\}$ is the action space, where action $1$ is activating/pulling the arm;
$P:\sspa\times\aspa\times\sspa \to [0,1]$ is the transition kernel;
$r: \sspa \times \aspa \to \mathbb{R}$ is the reward function.
Let $\rmax = \max_{s\in\sspa,a\in\aspa} \abs{r(s,a)}$.
The RB problem has a budget constraint, given by a constant $\alpha\in(0,1)$, which requires that $\alpha$ fraction of arms must be pulled at every time step. We assume that $\alpha N$ is an integer.
We focus on the setting where all the model parameters, $\sspa, \aspa, P, r, \alpha$, are known.

A policy $\pi$ for the $N$-armed problem chooses the action for each of the $N$ arms in each time step.
The policy can be randomized and history-dependent.
Under $\pi$, the state vector $\ve{S}_t^\pi \triangleq (S_t^\pi(i))_{i\in[N]} \in \sspa^N$ represents the states of all arms at time $t$, and the action vector  $\ve{A}_t^\pi \triangleq (A_t^\pi(i))_{i\in[N]} \in \aspa^N$ denotes the actions applied to each arm.
For each policy $\pi,$ consider the following two quantities: the limsup average reward
$\rlimsup \triangleq \limsup_{T\to\infty} \frac{1}{T} \sum_{t=0}^{T-1} \frac{1}{N} \sumN \E{r(S_t^\pi(i), A_t^\pi(i))};$ and the liminf average reward $\rliminf \triangleq \liminf_{T\to\infty} \frac{1}{T} \sum_{t=0}^{T-1} \frac{1}{N} \sumN \E{r(S_t^\pi(i), A_t^\pi(i))}$.
When the limsup and liminf coincide, the long-run average reward can be defined as follows:

\[
\rsysn \triangleq \lim_{T\to\infty} \frac{1}{T} \sum_{t=0}^{T-1} \frac{1}{N} \sumN \E{r(S_t^\pi(i), A_t^\pi(i))}.
\]

Note that $\rsysn$ is well-defined when $\pi$ is Markovian w.r.t. the state  $\veS_t$, as $\veS_t$ is a finite-state Markov chain under such $\pi$ <cit.>. Using a similar argument, one can show that $\rsysn$ is well-defined if $\pi$ is Markovian w.r.t. an augmented system state with a finite state space.

Our goal is to solve the following optimization problem.

\begin{align}
\label{eq:N-arm-formulation} \tag{RB}
\underset{\text{policy } \pi}{\text{maximize}} & \quad \rliminf  \qquad
\text{subject to}
\quad  \sumN A_t^\pi(i) = \alpha N,\quad \forall t\ge 0.
\end{align}

Let $\ropt$ be the optimal value of (<ref>), referred to as the optimal reward.
Note that $\ropt = \sup_{\pi'} R^-(\pi', \veS_0) = \sup_{\pi'} R^+(\pi', \veS_0)$ because (<ref>) is an MDP with finite state and action spaces <cit.>.
For any policy $\pi$, we define its optimality gap as $\ropt - \rliminf$. A policy is asymptotically optimal if its optimality gap vanishes as $N\to\infty$, i.e., $\ropt - \rliminf = o(1)$.
This notion of asymptotic optimality is consistent with the literature <cit.>.

We remark that there always exists a stationary Markovian policy that achieves the optimal reward <cit.>.
In later parts of the paper,
although we will not restrict ourselves to stationary Markovian policies,
we focus on policies that induce finite-state Markov chains.
Under each of these policies, $\pi$, the long-run average reward $\rsysn$ is well-defined.
For simplicity, we will refer to $\rsysn$ as the objective function of (<ref>) and write the optimality gap as $\ropt - \rsysn$.

Scaled state-count vector.
We introduce an alternative representation of the $N$-armed system state,
which will be used extensively throughout the paper.
For each subset $D\subseteq[N]$, we define the scaled state-count vector on $D$ as $X_t^\pi(D) = (X_t^\pi(D, s))_{s\in\sspa}$, where
X_t^\pi(D, s) \triangleq \frac{1}{N} \sum_{i\in D} \indibrac{S_t^\pi(i) = s}
$ is the number of arms in $D$ in state $s$ scaled by $1/N.$
When $D = [N]$ is the set of all arms, we simply call $X_t^\pi([N])$ the scaled state-count vector.
Sometimes we view $X_t^\pi(D)$ as a vector-valued function of $D \subseteq [N]$. We refer to this function $X_t$ as the system state at time $t$. Note that the system state $X_t^\pi$ contains the same information as the state vector $\veS_t^\pi$.
We also define the scaled state-action-count vector as $Y_t^\pi = (Y_t^\pi(s, a))_{s\in\sspa, a\in\aspa}$, where
$Y_t^\pi(s,a) \triangleq \frac{1}{N} \sumN \indibrac{S_t^\pi(i)=s, A_t^\pi(i) = a}$ is the fraction of arms with state $s$ and action $a$ at time $t$.

Additional notation.   Let $[0,1]_N = \{0,1/N, 2/N,\dots, 1\}$. For each subset $D\subseteq[N]$, $m(D) = |D| / N$ denotes the fraction of arms contained in $D$, thus $m(D) \in [0,1]_N.$
Let $\Delta(\sspa)$ be probability simplex on  $\sspa$.
We treat each distribution $v\in\Delta(\sspa)$ as a row vector.
Recall that $\pi$ denotes a policy for the $N$-armed problem. When the context is clear, we drop superscript $\pi$ from $\veS_t^\pi$, $\veA_t^\pi$, and $X_t^\pi$.

§ TWO-SET POLICY

In this section, we describe our Two-Set Policy.
We first discuss a linear programming (LP) relaxation of (<ref>) in <Ref> and state assumptions.
The optimal solution of the LP is then used to define two subroutines, Exact Proportional Control in <Ref>, and Optimal Local Control in <Ref>. We then use these two subroutines to construct the Two-Set Policy in <Ref>.
The ideas in the two subroutines are not new; our novelty lies in how we build them into the Two-Set Policy.

§.§ LP relaxation

We consider the linear programming (LP) relaxation of the $N$-armed problem (<ref>), which is crucial for the design and analysis of RB policies.

\begin{align}
\label{eq:lp-single} \tag{LP} \underset{\{y(s, a)\}_{s\in\sspa,a\in\aspa}}{\text{maximize}} \mspace{6mu}&\sum_{s\in\sspa,a\in\aspa} r(s, a) y(s, a) \\
\text{subject to}\mspace{6mu}
&\mspace{15mu}\sum_{s\in\sspa} y(s, 1) = \alpha, \label{eq:expect-budget-constraint}\\
& \sum_{s'\in\sspa, a\in\aspa} \!\!\! y(s', a) P(s', a, s) = \sum_{a\in\aspa} y(s,a),\, \forall s;
&\sum_{s\in\sspa, a\in\aspa} \!\!\! y(s,a) = 1;
\;
y(s,a) \geq 0, \; \forall s,a. \label{eq:non-negative-constraint}
\end{align}

To see why (<ref>) is a relaxation of (<ref>), for any stationary Markovian policy $\pi$, consider
y^\pi(s,a) = \lim_{T\to\infty} \frac{1}{T} \sum_{t=0}^{T-1} \Ebig{\frac{1}{N} \sumN \indibrac{S_t^\pi(i) = s, A_t^\pi(i) = a}}, \forall s\in\sspa, a\in\aspa.
It is not hard to see that $\rsysn = \sumsa r(s,a)y^\pi(s,a)$, and $(y^\pi(s,a))_{s\in\sspa, a\in\aspa}$ satisfies the constraints (<ref>)-(<ref>).
Therefore, letting $\rrel$ be the optimal value of (<ref>), it can be shown that $\rrel \geq \ropt$ (See <Ref>).
Note that the LP relaxation can be efficiently solved. We can bound the optimality gap of any policy $\pi$ using the inequality $\ropt - \rliminf \leq \rrel - \rliminf$.
This bound has been widely used in prior work [16, 14, 3, 4, 9].

Optimal single-armed policy. For a fixed optimal solution to (<ref>), $y^* = \{y^*(s, a)\}_{s\in\sspa,a\in\aspa}$, we consider the conditional distribution $\pibs = (\pibs(a|s))_{s\in\sspa, a\in\aspa}$ given by:

\begin{equation}\label{eq:single-arm-opt-def}
\pibs(a | s) =
\begin{cases}
y^*(s, a) / (y^*(s, 0) + y^*(s,1)), & \text{if } y^*(s, 0) + y^*(s,1) > 0,  \\
1/2, &  \text{if } y^*(s, 0) + y^*(s,1) = 0.
\end{cases}
\quad \text{for $s\in\sspa$, $a\in\aspa$.}
\end{equation}

We call $\pibs$ the optimal single-armed policy, due to the fact that $\pibs$ can be viewed as the optimal policy for the single-armed MDP $(\sspa, \aspa, P, r)$ with a relaxed budget constraint <cit.>.

Let $P_{\pibs}$ be the transition matrix induced by the policy $\pibs$ on the single-armed MDP, i.e., $P_{\pibs}(s,s') = \sum\nolimits_{a\in\aspa} \pibs(a|s) P(s,a,s')$.
Since $P_\pibs$ is a stochastic matrix, its largest eigenvalue is equal to $1$. For the other eigenvalues of $P_\pibs$, we consider the following assumption.

[Aperiodic unichain]
The moduli of all eigenvalues of $P_\pibs$ except the largest one are strictly less than $1$;
that is, $P_\pibs$ defines an aperiodic unichain on $\sspa$.

With <Ref>, the transition matrix $P_\pibs$ induces a unique stationary distribution, $\statdist \triangleq (\statdist(s))_{s\in\sspa}$, where $\statdist(s) = y^*(s,1) + y^*(s,0)$, as one can verify using the definition of $P_\pibs$.
We call $\statdist$ the optimal stationary distribution. We define $\sempty \triangleq \{s\in\sspa\colon y^*(s,1) = y^*(s,0) = 0\},$ the set of transient states under policy $\pibs.$

Our next assumption imposes a non-degenerate condition on the optimal solution to (<ref>).

For the fixed optimal solution to (<ref>), $y^*$,
there exists a unique neutral state, i.e., a state $\sneu\in\sspa$ s.t. $y^*(\sneu, 1) > 0$ and $y^*(\sneu, 0) > 0$.

Note that <Ref> stipulates the existence of a neutral state $\sneu$ — one can always find an optimal solution of (<ref>) such that there is at most one neutral state <cit.>.

[Local stability]
Given the non-degenerate condition,
define $\Phi$ as a $|\sspa|$-by-$|\sspa|$ matrix:

\begin{equation}\label{eq:phi-def}
\Phi \triangleq P_\pibs - \vone^\top \statdist - (\costvec-\alpha\vone)^\top  (P_1(\sneu) - P_0(\sneu)),
\end{equation}

where $\costvec \triangleq (\pibs(1|s))_{s\in\sspa}$ and $P_a(\sneu) \triangleq (P(\sneu,a, s))_{s\in\sspa}$ are both row vectors; $\vone$ is the all-one row vector.
We assume that the moduli of all eigenvalues of $\Phi$ are strictly less than $1$.

The implication of <Ref> will become clear after we introduce the subroutine Optimal Local Control in <Ref>, which is constructed from the optimal solution to (<ref>) and has dynamics given by the matrix $\Phi$. We note that [3] also assumes a version of local stability, in addition to UGAP, to achieve exponential asymptotic optimality.
Our local stability assumption is substantially weaker, and it is implied by their version of local stability together with UGAP (see Appendix C.4 of [3]). Importantly, we do not require a global assumption like UGAP.

Assumptions <ref>, <ref> and <ref> are easy to verify as they pertain to either the single-armed  MDP primitives ($P, r, \alpha$) or a linear program define by these problem primitives.
In contrast, verifying UGAP is significantly harder due to the dependency on the non-linear $N$-armed dynamics.

Assumptions <ref>, <ref> and <ref> all depend on an optimal solution  $y^*$ to the (<ref>).
When we state our results with more than one assumption among the three, we consider the same $y^*$.

§.§ Subroutine 1: Exact Proportional Control

We now describe the Exact Proportional Control, as presented in <Ref>. The subroutine takes the optimal single-armed policy $\pibs$ defined in (<ref>) and a subset of arms $D$ from the RB system as the input. The subroutine decides the actions for all arms in $D$ to match with $\pibs.$ Specifically, let the number of arms from $D$ in state $s\in \sspa$ be $z(s).$ The subroutine aims to activate $\pibs(1|s)$ fraction of the $z(s)$ arms. To handle the integer effect, it activates $\floor{\pibs(1|s) z(s)}$ or $\ceil{\pibs(1|s) z(s)}$ arms with properly chosen probability (computed in <Ref>) such that the expected proportion is exactly $\pibs(1|s).$

As shown in <Ref>, if all $N$ arms in the RB system follow Exact Proportional Control, then as $t\to\infty$,
the expected scaled state-count vector $\E{X_t([N])}$ converges to $\statdist$; the expected fraction of arms in state $s$ taking action $a$ converge to $y^*(s,a)$;
the expected budget usage of the arms converges to $\alpha N$. Consequently,
the average expected reward achieves the optimum $\rrel$ of (<ref>).

We remark that <Ref> let the arms follow $\pibs$ in a way different from some prior work [9, 10], where the actions are sampled from $\pibs$ independently across the arms.

Exact Proportional Control

Input: A set of arms with number of arms in each state $(z(s))_{s\in\sspa}$,
optimal single-armed policy $\pibs$

[1]
Compute $p$ such that $p \floor{\pibs(1|s) z(s)} + (1-p) \ceil{\pibs(1|s) z(s)} = \pibs(1|s) z(s)$.
With probability $p$ and $1-p$, activate $\floor{\pibs(1|s) z(s)}$ or $\ceil{\pibs(1|s) z(s)}$ arms in state $s$.

§.§ Subroutine 2: Optimal Local Control

Next, we describe <Ref>, Optimal Local Control, another building block for constructing our policy. The subroutine also takes the optimal single-armed policy $\pibs$ and a subset of arms $D$ as the input.
One key difference from <Ref> is that here <Ref> requires a predetermined budget $B$ as an input, independent of the states of the arms.
Specifically, for the input subset $D$ of $n$ arms, the budget $\bsub $ is set as $\floor{\alpha n}$ or $\ceil{\alpha n}$ with suitable probability such that $\E{\bsub} = \alpha n$, to handle $\alpha n$ being non-integer. When $D = [N]$, we set $\bsub=\alpha N$ as the original budget on the $N$-armed system. For each non-neutral state $s\in \sspa \backslash \{\sneu\},$ the policy also aims to activate $\pibs(1|s)$ fraction of the $z(s)$ arms from $D$ that are in state $s.$ To ensure this is feasible with the budget $B$ constraint,
we consider the following conditions for the input subset $D$:

\begin{align}
\label{eq:lp-priority-condition-1}
\sum\nolimits_{s\neq \sneu} \pibs(1|s) z(s) &\leq  \alpha n - |\sempty|-1, \quad \mbox{and} \quad
\sum\nolimits_{s\neq \sneu} \pibs(0|s) z(s) \leq (1-\alpha)n - |\sempty|-1.
\end{align}

Recall that $\sempty = \{s\in\sspa\colon y^*(s,1) = y^*(s,0) = 0\}$ is the set of transient states under $\pibs.$
By the definition of $\pibs$ in (<ref>), it is not hard to see that under conditions in (<ref>), there hold $\sum_{s\neq \sneu} \ceil{\pibs(1|s) z(s)} \leq \floor{\alpha n}$ and $z(\sneu) + \sum_{s\neq \sneu} \floor{\pibs(1|s) z(s)} \geq  \ceil{\alpha n}$. Therefore, the budget requested by non-neural states $\sspa \backslash \{\sneu\}$ does not exceed $B$ (<Ref>-<ref> of <Ref>), and the neural state ensures exactly $\bsub$ unit of budget is used (<Ref>).

Optimal Local Control

Input: A set of arms with number of arms in each state $(z(s))_{s\in\sspa}$, budget $\bsub$,

optimal single-armed policy $\pibs$, neutral state $\sneu$

[1]
Assert <Ref> and conditions in (<ref>)
$s\in \sspa \backslash \{\sneu\}$
Compute $p$ such that $p \floor{\pibs(1|s) z(s)} + (1-p) \ceil{\pibs(1|s) z(s)} = \pibs(1|s) z(s)$.
With probability $p$ and $1-p$, activate $\floor{\pibs(1|s) z(s)}$ or $\ceil{\pibs(1|s) z(s)}$ arms with state $s$.
Activate the arms with state $\sneu$ using the remaining budget out of $\bsub$.

If all arms in the RB system follow <Ref>, the expected scaled state-count vector $\Eplain{X_t([N])}$ has a linear dynamics, as stated by following the lemma (proof provided in <Ref>).

Under <Ref>, suppose (<ref>) holds for the $N$-armed system at time step $t$. We have

\begin{equation}\label{eq:lp-priority-trans}
\Eplain{X_{t+1}([N]) -\statdist \givenplain X_t, \text{all arms follow \Cref{alg:follow-lp-priority}}} = (X_t([N]) - \statdist) \Phi,
\end{equation}

where $\Phi$ is the $|\sspa|$-by-$|\sspa|$ matrix defined in (<ref>).

Thus under <Ref>,  <Ref> induces locally stable dynamics. Specifically, if the the RB system satisfies (<ref>) for a sufficiently long time and we let all $N$ arms follow <Ref>, then the scaled state-count vector $X_t([N])$ would converge to $\statdist$. Consequently, the expected fraction of arms in state $s$ taking action $a$ would be close to $y^*(s,a)$ and the time-average of expected reward would be close to the upper bound $\rrel$—hence the name “Optimal Local Control" for <Ref>.

§.§ Constructing the Two-Set Policy

We now present our Two-Set Policy for the $N$-armed system, whose pseudo-code is given in <Ref>.
On a high level, the Two-Set Policy works by picking out two disjoint subsets of arms among the $N$ arms in each time step, and letting them follow two subroutines, Optimal Local Control and Exact Proportional Control, respectively.
When the two sets are picked properly, the policy pushes the scaled state-count vector $X_t([N])$ of the $N$-armed system towards $\statdist$. Ultimately, the system will stabilize in the region where all arms can follow <Ref> $\lppriority$.

Specifically, at time step $t$, the Two-Set Policy first picks a subset $\Db_t$, where the arms will follow $\lppriority$.
The policy chooses $\Db_t$ in two steps, specified in <Ref>-<ref> of <Ref>. Roughly speaking, it intends to let $\Db_t$ have a large overlap with $\Db_{t-1}$ as well as a large cardinality, subject the constraint $\slkb(X_t, \Db_t) \geq 0$, where $\slkb(\cdot,\cdot)$ is a function called slack, defined as follows:

\begin{equation}
\label{eq:slkb-def}
\slkb(x, \Db) = \eta m(\Db) -  \norm{x(\Db) - m(\Db) \statdist}_{\umat},
\end{equation}

where $\eta$ is an input parameter with suggested value $|\sspa|^{-1/2}\big(\min\{y^*(\sneu, 0), y^*(\sneu, 1)\} - (|\sempty|+1)/N\big)$, $\umat$ is a matrix defined as $\umat \triangleq \sum_{k=0}^\infty \Phi^k (\Phi^\top)^k$, and  $\norm{v}_\umat = \sqrt{v \umat v^\top}$ for $v\in \R^{|\sspa|}$.
In particular, the criteria that $\slkb(X_{t}, \Db_{t}) \geq 0$ ensures that arms in $\Db_t$ satisfy a “closedness” property: they can continue to follow $\lppriority$ almost infinitely while satisfying condition (<ref>), until some rare event happens due to stochastic transitions. We show this property in <Ref> and <Ref>.

After choosing $\Db_t$, the Two-Set Policy picks a subset $\Da_{t}$ that contains $\beta \triangleq \min(\alpha, 1-\alpha)$ fraction of the remaining arms, rounded to a certain nearby integer.
The set $\Da_{t}$ should be such that $\Da_t \cap \Db_t = \emptyset$ and $\Da_{t} \supseteq \Da_{t-1} \backslash \Db_{t}$.

With $\Db_t$ and $\Da_t$ chosen, the Two-Set Policy let arms in $\Db_t$ follow <Ref> and let the arms in $\Da_t$ follow <Ref> Exact Proportional Control, and chooses the actions of remaining arms to meet the budget constraint.
We will show that this is always feasible in <Ref> in <Ref>.

Two-Set Policy

Input: number of arms $N$, budget $\alpha N$, an optimal solution of LP-relaxation $y^*$, threshold $\eta$,

initial system state $X_0$, initial state vector $\veS_0$, initial focus set $D_{-1}=\emptyset$
[1]
Let $\Dtemp$ be a set with largest cardinality such that $\Dtemp \subseteq \Db_{t-1}$, $\slkb(X_t, \Dtemp) \geq 0$.
Let $\Db_t$ be a set with largest cardinality such that $\Db_t \supseteq \Dtemp$, $\slkb(X_t, \Db_t) \geq 0$.
Let $\Da_{t}$ be s.t. $\Da_{t-1} \backslash \Db_{t} \subseteq \Da_{t} \subseteq [N] \backslash \Db_t$ and $|\Da_{t}| = \big(\big\lfloor \beta (N -|\Db_{t}|)\big\rfloor - 1\big)^+$.
Set $A_t(i)$ for $i\in \Db_{t}$ using <Ref> Optimal Local Control.

Set $A_t(i)$ for $i\in \Da_{t}$ using <Ref> exact propotional policy.

Set $A_t(i)$ for $i\notin \Db_{t} \cup \Da_{t}$ such that $\sumN A_t(i) = \alpha N$.

Apply $A_t(i)$ for $i\in[N]$ and observe the new state vector $\veS_{t+1}$.

§ MAIN RESULTS

In this section, we present our main results.

§.§ Optimality guarantee on the Two-Set Policy

Our first main theorem states an upper bound on the optimality gap of the Two-Set Policy.

Suppose <Ref>, <ref> and <ref> hold. Let $\pi$ be the Two-Set Policy in <Ref> with $\eta =  |\sspa|^{-1/2}\big(\min\{y^*(\sneu, 0), y^*(\sneu, 1)\} - (|\sempty|+1)/N\big)$. Then for some $C > 0$ independent of $N$:

\begin{equation}
\rrel - \rsysn = O\big(\exp(-CN)\big).
\end{equation}

<Ref> establishes the first exponential asymptotic optimality without UGAP.
It implies that even in an RB system with a small number of arms, Two-Set Policy is nearly optimal, under three mild and easy-to-verify assumptions.
In contrast, prior work either requires UGAP to achieve exponential asymptotic optimality or only has $O(1/\sqrt{N})$ optimality gap.
We prove <Ref> in <Ref>.

We discuss the relationship between our Two-Set Policy and the set-expansion policy from prior work [10]. Set-expansion achieves an $O(1/\sqrt{N})$ optimality gap under the aperiodic unichain assumption (<Ref>). It works by letting arms in a specific subset follow $\pibs$ and gradually expanding the subset.
Our Two-Set Policy refines this approach by maintaining two judiciously chosen sets instead of one set.
Furthermore, we argue that when Assumption <ref> or <ref> fails, a variant of the Two-Set Policy can achieve $O(1/\sqrt{N})$ optimality gap, since it becomes almost the same as the set-expansion policy.
To see this, note that when Assumption <ref> (non-degeneracy) fails, it holds that $\eta = -O(1/N)$; when Assumption <ref> (local stability) fails, the matrix $U$ becomes infinity.
In either case, $\Db_t$ will become the maximal set such that $x(\Db_t) = m(\Db_t) \statdist$ if we allow $\slkb(X_t, \Db_t)$ to be slightly negative.
Consequently, $\Db_t$ equivalently follows the Exact Proportional Control and can be combined with $\Da_t$.
Although the Two-Set Policy and the set-expansion policy differ in how the arms follow $\pibs$ and the criteria for choosing the subset, we believe these differences are minor. Using similar arguments as [10], one can prove an $O(1/\sqrt{N})$ bound for the Two-Set Policy.

§.§ Discussion of the assumptions

In this subsection, we discuss the importance of the assumptions in <Ref>. We will show that there are scenarios where, for every policy $\pi$ and certain initial states $\veS_0$, it holds that

\begin{equation}
\label{eq:lower-bound}
\rrel - \rsysn \geq \Omega\Big(\frac{1}{\sqrt{N}}\Big).
\end{equation}

When <Ref> (aperiodic  unichain) or <Ref> (non-degeneracy) is violated, it is not hard to construct examples where (<ref>) holds for every policy $\pi$ and certain initial states $\veS_0$; see <Ref>.
In the non-unichain example, $\ropt$ may depend on the initial states $\veS_0$, in which case the initial-state-independent upper bound $\rrel$ is loose.
In the periodic and degenerate examples, the upper bound $\rrel$ can only be achieved by a policy that satisfies the budget constraint on average over time rather than every time step, again causing a large gap between $\rrel$ and $\ropt$.

We emphsize that although $\rrel$ is just an upper bound of $\ropt$,
all existing bounds on the optimality gap use $\rrel$ or its variants. Hence, the above examples reveal the limitations of existing technology.
Specifically, while the unichain assumption may be circumvented by a tighter upper bound than $\rrel$, the aperiodic and non-degenerate assumptions appear more fundamental — they are closely tied to the per-step constraint of RB problems.
Therefore, relaxing the aperiodic and non-degeneracy assumptions requires fundamentally new approaches.

§.§.§ A lower bound when <Ref> is violated

When local stability (<Ref>) fails, instead of finding one counterexample, we take a step further to establish the lower bound (<ref>) for a large class of instances.
To establish this lower bound, we need a fundamental lemma that relates the gap $\rrel-\ropt$ to our <Ref>.

Suppose that the LP relaxation (<ref>) has a unique optimal solution $y^*$, and $y^*(s,1) + y^*(s,0) > 0,$ $\forall s\in\sspa$. Then

\begin{equation}
\rrel - \rsysn \geq \rewardgap \lim_{T\to\infty} \frac{1}{T} \sum\nolimits_{t=0}^{T-1} \Ebig{\disy(Y_t^{\pi})},
\end{equation}

where $\rewardgap$ is a positive constant; $\disy$ is a function of state-action distribution given by

\[
\disy(y) = \sum\nolimits_{s\in S^+} y(s,0)+ \sum\nolimits_{s\in S^-}y(s,1),
\]

where $S^+ = \{s\in\sspa\colon y^*(s,1) > 0, y(s,0) = 0\}$, $S^- = \{s\in\sspa\colon y^*(s,1) = 0, y^*(s,0) > 0\}$.

Note that $\disy(Y_t)$ is zero if and only if $Y_t$ is generated by <Ref> (Optimal Local Control). Therefore, <Ref> shows that the gap $\rrel - \rsysn$ is determined by how closely a policy $\pi$  follows <Ref>.
<Ref> is proved  in <Ref> using a primal-dual argument.

Because $\disy(Y_t) \geq \absplain{\sum_{s\in S^+} X_t([N],s) - \alpha}$, <Ref> implies that for an RB satisfying conditions of the lemma, to achieve exponential asymptotic optimality, one must keep the long-run average of $\absplain{\sum_{s\in S^+} X_t([N],s) - \alpha}$ exponentially small. Doing so could be hard when the RB is degenerate, as there is no neutral state to absorb the randomness of $X_t([N])$.

<Ref> can also be used to show a lower bound on $\rrel - \rsysn$ for a class of RB instances with an unstable $\Phi$ matrix, which we call regular unstable RB, as defined below.

For an RB instance satisfying <Ref> and <Ref>, we call it regular unstable if it has the following properties:

(<ref>) has a unique optimal solution $y^*$, and $y^*(s,1) + y^*(s,0) > 0$ for all $s\in\sspa$.

The matrix $\Phi$ defined in (<ref>) has at least one eigenvalue whose modulus is strictly larger than $1$, and the modulus of any eigenvalue of $\Phi$ does not equal $1$.

For each $\xi \in \R^{|\sspa|}$ with $\normplain{\xi}_2 = 1$ and $\statdist \xi^\top = 0$, for each large-enough $N$  and each time step $t$,

\begin{equation*}
\Ebig{\absbig{\big(X_{t+1}([N]) - \statdist - (X_t([N]) - \statdist)\Phi \big)  \xi^\top } \givenbig X_t, \veA_t} \geq \frac{\noiseortho}{\sqrt{N}} \quad \text{ if } \norm{X_t([N]) - \statdist}_1 \leq \radiusus,
\end{equation*}

where $\radiusus, \noiseortho > 0$ are some constants; $\veA_t$ follows the Optimal Local Control (<Ref>), which is always feasible if $\norm{X_t([N]) - \statdist}_1 \leq \min\{y^*(\sneu, 0), y^*(\sneu, 1)\}-1/N$.

Our next main theorem is a lower bound for regular unstable RBs; it is proved in <Ref>.

For every regular unstable RB, every policy $\pi$ and any initial state vector $\veS_0$, we have

\[
\rrel - \rsysn \geq \Omega\Big(\frac{1}{\sqrt{N}}\Big).
\]

<Ref> implies that an exponential optimality gap cannot be achieved for a significant subset of  RB instances violating local stability (<Ref>).
Therefore, in addition to the widely accepted non-degenerate condition, local stability is also fundamental for exponential  optimality.

We briefly comment on Property <ref> of regular unstable RB: this property requires the state-count vector to have sufficient randomness in most directions under the Optimal Local Control when the system is near the optimal state distribution $\statdist$.
In <Ref>, we provide an intuitive sufficient condition for <ref>, and a weaker version of the condition under which <Ref> remains valid.

The proof of <Ref> is based on the following idea: By <Ref>, arms need to closely follow the Optimal Local Control to maintain a small optimality gap; however, if Optimal Local Control is locally unstable, the system will quickly drift away from the region of the state space where Optimal Local Control is applicable.
In the proof, we use Lyapunov functions to quantify how much the system drifts away from $\statdist$ in each step, based on how closely it follows the Optimal Local Control.

We prove the following in <Ref> by giving a concrete example of regular unstable RB.

There exists an RB satisfying Assumptions <ref> and <ref> and violating <Ref> such that
\rrel - \rsysn \geq \Omega\Big(\frac{1}{\sqrt{N}}\Big).

§ PROOF SKETCH OF THEOREM <REF>

The Two-Set Policy maintains two dynamics subsets, each running a different subroutine. The main challenge of the proof lies in analyzing the complex dynamics of the two subsets.
To tackle this challenge, we use a novel multivariate Lyapunov function that takes not only the system state $X_t$ but also the subsets $\Db_t$, $\Da_t$ as variables.

Step 1. We first show that the optimality gap can be bounded as follows:
\rrel - \rsysn \leq O(1)\cdot \Probbig{\Db_\infty \neq [N]},
where the subscript $\infty$ indicates the steady-state distribution. The proof of this bound is similar to those of Theorems 3 and 8 of [3] and is provided in <Ref>.

Step 2. Next, by the definition of $\Db_t$, it is not hard to see that if $\norm{X_t([N]) - \statdist}_\umat \leq \eta$, then $\Db_t = [N]$. Therefore, it suffices to show that $\Prob{\norm{X_\infty([N]) - \statdist}_\umat > \eta}$ is exponentially small.
To bound the distance $\normplain{X_\infty([N]) - \statdist}_\umat$, an important step is to construct proper bivariate functions that capture the dynamics of the arms in the subsets $\Db_t$ and $\Da_t$. We define

\begin{align*}
\hw(x, D) &= \norm{x(D) - m(D)\statdist}_\wmat, \quad
\hu(x, D) = \norm{x(D) - m(D)\statdist}_\umat,
\end{align*}

where $\wmat$ is a $|\sspa|$-by-$|\sspa|$ weight matrix and ${v}_\wmat = \sqrt{v \wmat v^\top}$ for $v\in \R^{|\sspa|}$.
The matrices $\wmat$ and $\umat$ are designed such that $P_\pibs$ and $\Phi$ are pseudo-contractions on $\simplex(\sspa)$ with fixed point $\statdist$ under the $\umat$-weighted norm and the $\wmat$-weighted norm, respectively (<Ref>).
The bivariate functions $\hw$ and $\hu$ are referred to as subset Lyapunov functions, because they witness the convergence of state distributions on a fixed subset of arms when that subset continues to follow either the Exact Proportional Control or the Optimal Local Control (<Ref>).

Step 3. We use subset Lyapunov functions $\hw$ and $\hu$ as building block for constructing our global Lyapunov function $V(x, \Db, \Da)$, defined as

\begin{equation*}
V(x, \Db, \Da) = \hu(x, \Db) + \hw(x, \Da) + \Ltemp (1 - m(\Db)),
\end{equation*}

where $L$ is a constant ensuring that $V(x, \Db, \Da)$ bounds $\hu(x, [N]) = \normplain{X_\infty([N]) - \statdist}_\umat$.
We show that $V(X_t, \Db_t, \Da_t)$ decreases by a large amount every time step with high probability,  unless $V(X_t, \Db_t, \Da_t) = O(1/\sqrt{N})$, which implies that the Lyapunov function is small in steady-state.

A key idea in analyzing the changes of $V(X_t, \Db_t, \Da_t)$ is to analyze those changes caused by state transitions and those caused by set-updates separately, leading to the key inequality below:

\begin{align}
\nonumber
& V(X_{t+1}, \Db_{t+1}, \Da_{t+1}) - V(X_t, \Db_t, \Da_t) \\
\nonumber
&\leq \big(\hu(X_{t+1}, \Db_t) - \hu(X_{t}, \Db_t)\big)+ \big(\hw(X_{t+1}, \Da_t) - \hw(X_{t}, \Da_t)\big)\\
\nonumber
&\mspace{20mu} + 4\big(\lamu^{1/2} + \lamw^{1/2}\big) m(\Db_t \backslash \Db_{t+1}) + O\Big(\frac{1}{N}\Big).
\end{align}

To upper bound the RHS of this inequality, we prove two lemmas to show that $m(\Db_t \backslash \Db_{t+1})$ is small (<Ref>) and that $\big(\hu(X_{t+1}, \Db_t) - \hu(X_{t}, \Db_t)\big)+ \big(\hw(X_{t+1}, \Da_t) - \hw(X_{t}, \Da_t)\big)$ is sufficiently negative (<Ref>). See <Ref> for the full details of the proof.

[1]
David B. Brown and James E. Smith.
Index policies and performance bounds for dynamic selection problems.
Management Science, 66(7):3029–3050, 2020.

[2]
Nicolas Gast, Bruno Gaujal, and Chen Yan.
Exponential convergence rate for the asymptotic optimality of whittle
index policy.
arXiv:2012.09064 [cs.PF], 2020.

[3]
Nicolas Gast, Bruno Gaujal, and Chen Yan.
Exponential asymptotic optimality of whittle index policy.
Queueing Systems, 104(1):107–150, 2023.

[4]
Nicolas Gast, Bruno Gaujal, and Chen Yan.
Linear program-based policies for restless bandits: Necessary and
sufficient conditions for (exponentially fast) asymptotic optimality.
Math. Oper. Res., 2023.

[5]
Nicolas Gast, Bruno Gaujal, and Chen Yan.
Reoptimization nearly solves weakly coupled markov decision
arXiv preprint arXiv:2211.01961, 2024.

[6]
Abheek Ghosh, Dheeraj Nagaraj, Manish Jain, and Milind Tambe.
Indexability is not enough for whittle: Improved, near-optimal
algorithms for restless bandits.
In Proc. Int. Conf. Autonomous Agents and Multiagent Syst.,
page 1294–1302, Richland, SC, 2023.

[7]
Alan J Goldman and Albert W Tucker.
Theory of linear programming.
Linear inequalities and related systems, 38:53–97, 1956.

[8]
Weici Hu and Peter Frazier.
An asymptotically optimal index policy for finite-horizon restless
arXiv:1707.00205 [math.OC], 2017.

[9]
Yige Hong, Qiaomin Xie, Yudong Chen, and Weina Wang.
Restless bandits with average reward: Breaking the uniform global
attractor assumption.
In A. Oh, T. Neumann, A. Globerson, K. Saenko, M. Hardt, and
S. Levine, editors, Advances in Neural Information Processing Systems,
volume 36, pages 12810–12844. Curran Associates, Inc., 2023.

[10]
Yige Hong, Qiaomin Xie, Yudong Chen, and Weina Wang.
Unichain and aperiodicity are sufficient for asymptotic optimality of
average-reward restless bandits.
arXiv preprint arXiv:2402.05689, 2024.

[11]
José Niño-Mora.
Markovian restless bandits and index policies: A review.
Mathematics, 11(7), 2023.

[12]
Christos H. Papadimitriou and John N. Tsitsiklis.
The complexity of optimal queuing network control.
Math. Oper. Res., 24(2):293–305, 1999.

[13]
Martin L Puterman.
Markov decision processes: Discrete stochastic dynamic
John Wiley & Sons, 2005.

[14]
I. M. Verloop.
Asymptotically optimal priority policies for indexable and
nonindexable restless bandits.
Ann. Appl. Probab., 26(4):1947–1995, 2016.

[15]
Peter Whittle.
Restless bandits: activity allocation in a changing world.
J. Appl. Probab., 25:287 – 298, 1988.

[16]
Richard R. Weber and Gideon Weiss.
On an index policy for restless bandits.
J. Appl. Probab., 27(3):637–648, 1990.

[17]
Chen Yan.
An optimal-control approach to infinite-horizon restless bandits:
Achieving asymptotic optimality with minimal assumptions.
arXiv preprint arXiv:2403.11913, 2024.

[18]
Gabriel Zayas-Cabán, Stefanus Jasin, and Guihua Wang.
An asymptotically optimal heuristic for general nonstationary
finite-horizon restless multi-armed, multi-action bandits.
Advances in Applied Probability, 51:745–772, 2019.

[19]
Xiangyu Zhang and Peter I. Frazier.
Restless bandits with many arms: Beating the central limit theorem.
arXiv:2107.11911 [math.OC], 2021.

[20]
Xiangyu Zhang and Peter I. Frazier.
Near-optimality for infinite-horizon restless bandits with many arms.
arXiv:2203.15853 [cs.LG], 2022.

§ PRELIMINARY LEMMAS AND PROOFS

In this section, we state and prove several preliminary lemmas.
In <Ref>, we show that the linear program (<ref>) is a relaxation of the $N$-armed problem.
In <Ref>, we prove two lemmas about the transition dynamics under the two subroutines, Exact Proportional Control (<Ref>) and Optimal Local Control (<Ref>).
In <Ref>, we define two weight matrices, $\wmat$ and $\umat$, and show the pseudo-contraction properties of $P_\pibs$ and $\Phi$ under the $L_2$ norms weighted by $\wmat$ and $\umat$, respectively.

§.§ Justifying the LP relaxtion

In this section, we show that linear program (<ref>) is indeed a relaxation of the $N$-armed problem, i.e., $\rrel \geq \ropt$.
While this fact is used throughout the literature, our paper assumes weaker assumptions than most prior work, so for the sake of rigor, we prove it in our setting for completeness.

Let $\rrel$ be the optimal value of the linear program (<ref>), and let $\ropt$ be the optimal reward of the $N$-armed restless bandit problem (<ref>). Then we have

\begin{equation}
\rrel \geq \ropt.
\end{equation}

For any stationary Markovian policy $\pi$, define

\[
y^{\pi}(s,a) \triangleq \lim_{T\to\infty} \frac{1}{T} \sum_{t=0}^{T-1} \EBig{\frac{1}{N} \sumN \indibrac{S_t^\pi(i) = s, A_t^\pi(i) = a}} \quad \forall s\in\sspa, a\in\aspa.
\]

We first show that $\rsysn = \sumsa r(s,a)y^\pi(s,a)$.

\begin{align}
\sumsa r(s,a)y^\pi(s,a)
\nonumber
&= \sumsa r(s,a)\lim_{T\to\infty} \frac{1}{T} \sum_{t=0}^{T-1} \EBig{\frac{1}{N} \sumN \indibrac{S_t^\pi(i) = s, A_t^\pi(i) = a}} \\
\nonumber %\label{eq:upper-bound:interchange-sum-and-limit}
&= \lim_{T\to\infty} \frac{1}{NT}\sum_{t=0}^{T-1} \sumN \EBig{\sumsa r(s,a) \indibrac{S_t^\pi(i) = s, A_t^\pi(i) = a}} \\
\nonumber
&= \lim_{T\to\infty} \frac{1}{NT}\sum_{t=0}^{T-1} \sumN \E{r(S_t^\pi(i), A_t^\pi(i)} \\
\nonumber
&= \rsysn.
\end{align}

Next, we show that $(y^\pi(s,a))_{s\in\sspa, a\in\aspa}$ satisfies the constraints of (<ref>).
We first consider (<ref>):

\begin{align*}
\sums y^\pi(s,1)
&= \sums \lim_{T\to\infty} \frac{1}{T} \sum_{t=0}^{T-1} \EBig{\frac{1}{N} \sumN \indibrac{S_t^\pi(i) = s, A_t^\pi(i) = 1}}  \\
&= \lim_{T\to\infty} \frac{1}{NT}\sum_{t=0}^{T-1} \sumN \EBig{\sums \indibrac{S_t^\pi(i) = s, A_t^\pi(i) = 1}} \\
&= \lim_{T\to\infty} \frac{1}{NT}\sum_{t=0}^{T-1} \alpha N \\
&= \alpha.
\end{align*}

Next, we look at the constraint (<ref>):

\begin{align*}
\sum_{s'\in\sspa, a\in\aspa} y^\pi(s',a) P(s',a,s)
&= \lim_{T\to\infty} \frac{1}{NT}\sum_{t=0}^{T-1} \sumN \sumsa P(s',a,s) \Prob{S_t^\pi(i) = s', A_t^\pi(i) = a} \\
&=  \lim_{T\to\infty} \frac{1}{NT}\sum_{t=0}^{T-1} \sumN \Prob{S_{t+1}^\pi(i) = s} \\
&= \lim_{T\to\infty} \frac{1}{NT}\sum_{t=1}^{T} \sumN \Prob{S_t^\pi(i) = s} \\
&= \sum_{a\in\aspa} y^\pi(s,a) \\
\sumsa y^\pi(s,a)
&= \frac{1}{NT}\sum_{t=0}^{T-1} \sumN \EBig{\sumsa \indibrac{S_t^\pi(i) = s, A_t^\pi(i) = a}}
= 1,
\end{align*}

and it is obvious that $\sumsa y^\pi(s,a) \geq 0$.

Combining the above argument, $(y^\pi(s,a))_{s\in\sspa, a\in\aspa}$ is a feasibile solution to (<ref>), so
$\rsysn = \sumsa r(s,a)y^\pi(s,a) \leq \rrel$.

By standard results for MDP with finite state and action spaces, there always exists a stationary Markovian policy whose long-run average reward achieves the optimal reward <cit.>. Letting $\pi$ be this optimal stationary Markovian policy, then $\ropt = \rsysn \leq \rrel$.

§.§ Transition dynamics under subroutines

The lemma below characterizes the dynamics of a $N$-armed RB system if all arms follow Exact Proportional Control (<Ref>).

For any time step $t$,

\begin{equation}
\E{X_{t+1}([N]) \givenplain X_t, \textup{all arms follow \Cref{alg:follow-pibs}}} = X_t([N]) P_{\pibs} \quad a.s.,
\end{equation}

where expectation is taken entry-wise for the random vector $X_{t+1}([N])$; $P_{\pibs}$ is the transition matrix induced by $\pibs$ on the single-armed MDP, i.e., $P_{\pibs}(s,s') = \sum_{a\in\aspa} \pibs(a|s) P(s,a,s')$ for $s,s'\in\sspa$.

By definition, if all $N$ arms follow Exact Proportional Control at time $t$, the expected number of arms in state $s$ taking action $a$ is equal to $X_t([N], s)\pibs(a|s)$. Among these arms, the fraction of arms transitioning to state $s'$ in the next time step equals $P(s,a,s')$. Since the state transition of each arm is Markovian given the state and action, the expected number of arms going from state $s$ to state $s'$ and taking action $a$ equals $X_t([N], s)\pibs(a|s)P(s,a,s')$. Summing over $s$ and $a$, we get

\begin{align*}
&\mspace{23mu} \E{X_{t+1}([N], s') \givenplain X_t, \textup{all arms follow \Cref{alg:follow-pibs}}} \\
&= \sumsa X_t([N], s)\pibs(a|s)P(s,a,s') \quad \forall s'\in\sspa,
\end{align*}

which completes the proof of the lemma.

If all arms in the RB system follow Optimal Local Control (<Ref>), the expectation of the scaled state-count vector $X_t([N])$ has a linear dynamics, as stated in the lemma below. The lemma is basically a restatement of Lemma 2 in [3], although in [3], the LP-Priority policy was considered, and it was assumed that $y^*(s,1) + y^*(s,0) > 0$ for all $s\in\sspa$.

Throughout this proof, we use the shorthand $x$ for $x([N])$, and omit “$\veA_t\sim \lppriority$” in the conditionings.

We fixed an arbitrary $x$ and let $X_t = x$.
For any $s\neq \sneu$ and $a\in\aspa$, the expected number of arms taking action $a$ equals $x(s) \pibs(a|s)$; for $s=\sneu$, the expected number of arms taking action $1$ equals $\alpha -\sum_{s\neq\sneu }x(s)\pibs(1|s)$. and the expected number of arms taking action $0$ equals $1-\alpha- \sum_{s\neq\sneu }x(s)\pibs(0|s)$.
for any $s'\in\sspa$ such that $s'\neq \sneu$,

\begin{align}
&\E{X_{t+1}(s') \givenplain X_t=x} \nonumber \\
\nonumber
&= \sum_{s\neq \sneu, a\in\aspa} x(s) \pibs(a|s) P(s,a,s') \\
\nonumber
&\mspace{20mu} + \Big(\alpha - \sum_{s\neq \sneu} x(s)\pibs(1|s)\Big) P(\sneu, 1, s')
\nonumber
+ \Big(1 - \alpha - \sum_{s\neq \sneu} x(s) \pibs(0|s) \Big) P(\sneu, 0, s')  \\
\nonumber
&= \sum_{s\neq \sneu} x(s)P_\pibs(s, s')
+ \Big(\alpha - \sum_{s\neq \sneu} x(s) \pibs(1|s)\Big)P_1(\sneu, s')
+ \Big(1 - \alpha - \sum_{s\neq \sneu} x(s)\pibs(0|s) \Big)P_0(\sneu, s') \\
\label{eq:derive-phi-1}
&= \sums x(s)P_\pibs(s,s')  - x(\sneu)P_\pibs(\sneu, s')  \\
\nonumber
&\mspace{20mu} + \Big(\alpha - \sum_{s\neq \sneu} x(s) \pibs(1|s)\Big)P_1(\sneu, s')
+ \Big(1 - \alpha - \sum_{s\neq \sneu} x(s)\pibs(0|s) \Big)P_0(\sneu, s').
\end{align}

Because $P_\pibs(\sneu,s') =  \pibs(1|\sneu) P_1(\sneu,s')+ \pibs(0|\sneu)P_0(\sneu, s')$, we expand the second term in (<ref>) and rearrange the terms to get

\begin{align}
&\mspace{20mu} \E{X_{t+1}(s') \givenplain X_t=x} \nonumber\\
&=  \sums x(s)P_\pibs(s,s')
+  \Big(\alpha - \sums x(s) \pibs(1|s) \Big) P_1(\sneu, s')
+ \Big(1 - \alpha - \sums x(s) \pibs(0|s) \Big) P_0(\sneu, s')\nonumber \\
&= \sums x(s) P_\pibs(s,s')
+ \Big(\alpha - \sums x(s)\pibs(1|s) \Big) P_1(\sneu, s')
-\Big(  \alpha - \sums x(s)\pibs(1|s) \Big) P_0(\sneu, s') \label{eq:derive-phi-2}\\
&= \sums x(s)\Big(P_\pibs(s,s') -  \pibs(1|s)\big(P_1(\sneu, s') - P_0(\sneu, s')\big) \Big) + \alpha (P_1(\sneu, s') - P_0(\sneu, s')),  \label{eq:derive-phi-3}
\end{align}

where (<ref>) uses the facts that $\pibs(0|s) = 1-\pibs(1|s)$ and $\sums x(s) = 1$.

Next we calculate $\E{X_{t+1}(\sneu) \givenplain X_t=x}$.
Observe that $\E{X_t(\sneu)} = 1 - \sum_{s'\neq \sneu} \E{X_t(s')}$, so

\begin{align}
\nonumber
&\mspace{20mu} \E{X_{t+1}(\sneu) \givenplain X_t=x}\\
\nonumber
&= 1 - \sum_{s'\neq \sneu} \E{X_t(s')} \nonumber \\
\nonumber
&= \sums  x(s)\Big(1 - \sum_{s'\neq \sneu} P_\pibs(s,s') + \sum_{s'\neq\sneu}  \pibs(1|s) \big(P_1(\sneu, s') - P_0(\sneu, s')\big) \Big) - \alpha \sum_{s'\neq\sneu}(P_1(\sneu, s') - P_0(\sneu, s')) \\
\label{eq:derive-phi-3-sneu}
&= \sums x(s) \Big(P_\pibs(s,\sneu) -  \pibs(1|s) \big(P_1(\sneu, \sneu) - P_0(\sneu, \sneu)\big)\Big)  + \alpha (P_1(\sneu, \sneu) - P_0(\sneu, \sneu)).
\end{align}

Note that (<ref>) has the same form as (<ref>), with $s'$ replaced by $\sneu$.

We rewrite (<ref>) and (<ref>) in matrix form to get

\begin{align}\label{eq:derive-phi-4}
\E{X_{t+1} \givenplain X_t=x}
&= x\Big(P_\pibs - \costvec^\top \big(P_1(\sneu) - P_0(\sneu)\big) \Big) + \alpha (P_1(\sneu) - P_0(\sneu)).
\end{align}

Recall that we have let $\Phi = P_\pibs - \vone^\top \statdist - (\costvec - \alpha \vone)^\top \big(P_1(\sneu) - P_0(\sneu)\big)$, we can verify that the RHS of (<ref>) is the same as $(x-\statdist)\Phi + \statdist$. Specifically,

\begin{align}
(x-\statdist)\Phi + \statdist
\nonumber
&= (x-\statdist)\Big(P_\pibs - \vone^\top \statdist - (\costvec - \alpha \vone)^\top \big(P_1(\sneu) - P_0(\sneu)\big)\Big) + \statdist \\
\nonumber
&= x \Big(P_\pibs - \costvec^\top \big(P_1(\sneu) - P_0(\sneu)\big) \Big) + x\Big(- \vone^\top \statdist -  \alpha \vone^\top \big(P_1(\sneu) - P_0(\sneu)\big)\Big) \\
\nonumber
&\mspace{20mu} - \statdist \Big(P_\pibs - \costvec^\top \big(P_1(\sneu) - P_0(\sneu)\big) \Big) - \statdist \Big(- \vone^\top \statdist -  \alpha \vone^\top \big(P_1(\sneu) - P_0(\sneu)\big)\Big) + \statdist\\
\nonumber
&= x \Big(P_\pibs - \costvec^\top \big(P_1(\sneu) - P_0(\sneu)\big) \Big) - \statdist - \alpha \big(P_1(\sneu) - P_0(\sneu)\big) \\
\label{eq:derive-phi-5}
&\mspace{20mu} - \statdist  + \alpha \big(P_1(\sneu) - P_0(\sneu)\big) +  \statdist +  \alpha \big(P_1(\sneu) - P_0(\sneu)\big) + \statdist \\
\nonumber
&= x\Big(P_\pibs - \costvec^\top \big(P_1(\sneu) - P_0(\sneu)\big) \Big) + \alpha (P_1(\sneu) - P_0(\sneu)).
\end{align}

where to get (<ref>) we use the facts that $x\vone^\top = \statdist \vone^\top = 1$ and $\statdist \costvec^\top = \alpha$.
Substituting the RHS of (<ref>)  with $(x-\statdist)\Phi + \statdist$, we have proved that under Optimal Local Control, $\E{X_{t+1} \givenplain X_t=x} = (x-\statdist)\Phi + \statdist$.
This completes the proof of <Ref>.

§.§ Weighted $L_2$ norms for quantifying the convergence of distributions

As discussed in <Ref> and <Ref>, if a set of arms follow Exact Proportional Control or Optimal Local Control for a long time, the distribution of their states converge to the optimal stationary distribution $\statdist$.
In our analysis, it is convenient to define a proper metric to witness this convergence in each time step and quantify the convergence rate.
For this purpose, we introduce two matrices, $\wmat$ and $\umat$, and study the $\wmat$-weighted ($\umat$-weighted) $L_2$ norm on $\sspa$.
Note that the same matrix $\wmat$ and $\wmat$-weighted norm have also been used in [10].

The definition of $\wmat$ and $\umat$ are given as follows.

Assume Assumption <ref>, <ref> and <ref>.
Let $\wmat$ and $\umat$ be $|\sspa|$-by-$|\sspa|$ matrices given by

\begin{align}
\label{eq:w-def}
\wmat &= \sum_{k=0}^\infty (P_\pibs - \vone^\top \statdist)^k ((P_\pibs - \vone^\top \statdist)^\top)^k \\
\umat &= \sum_{k=0}^\infty \Phi^k (\Phi^\top)^k,
\end{align}

Next, we show that the matrices $\wmat$ and $\umat$ are well-defined and positive definite, and lower bound their eigenvalues by $1$. Then if we let $\lamw$ and $\lamu$ denote maximal eigenvalues of $\wmat$ and $\umat$, the eigenvalues of $\wmat$ and  $\umat$ are in the range $[1, \lamw]$ and $[1, \lamu]$, respectively.

Assume <Ref>, <ref> and <ref>, the matrices $\wmat$ and $\umat$ given in <Ref> are well-defined. Moreover, their eigenvalues are lower bounded by $1$.

The well-definedness and eigenvalue lower bound for $\wmat$ has been proved in [10]. The main fact used in the proof is that all eigenvalues of $P_\pibs - \vone^\top \statdist$ have moduli strictly less than $1$ under <Ref>.
Since we have assumed the same things for $\Phi$ in <Ref>, the proof for $\umat$ is analogous with $P_\pibs - \vone^\top \statdist$ replaced by  $\Phi$.

With matrices $\wmat$ and $\umat$, we can define some weighted $L_2$ norms.
Let $\norm{\cdot}_\wmat$ denote the $\wmat$-weighted $L_2$ norm, given by $\norm{u}_\wmat = \sqrt{u \wmat u^\top}$ for any row vector $u\in \R^{|\sspa|}$; let $\norm{\cdot}_\umat$ denote the $\umat$-weighted $L_2$ norm, defined similarly.

<Ref> below shows that $P_\pibs$ and $\Phi$ are pseudo-contractions on the probability simplex $\simplex(\sspa)$, under the $\wmat$-weighted $L_2$ norm and the $\umat$-weighted $L_2$ norm, respectively.

[Pseudo-contraction under the weighted $L_2$ norms]lemmawucontraction
Assume Assumption <ref>, <ref> and <ref>.
For any distribution $v\in \Delta(\sspa)$, we have

\begin{align}
\label{eq:pibar-contraction}
\norm{(v  - \statdist)P_\pibs}_\wmat &\leq \rhow \norm{v - \statdist}_\wmat \\
\label{eq:lp-contraction}
\norm{(v  - \statdist)\Phi}_\umat &\leq \rhou \norm{v - \statdist}_\umat.
\end{align}

where $\norm{\cdot}_\wmat$ is the $\wmat$-weighted $L_2$ norm, i.e., $\norm{u}_\wmat = \sqrt{u \wmat u^\top}$ for any row vector $u\in \R^{|\sspa|}$; $\norm{\cdot}_\umat$ is the $\umat$-weighted $L_2$ norm defined similarly; $\rhow = 1 - 1/(2\lamw)$; $\rhou = 1 - 1/(2\lamu)$.

The inequality (<ref>) has been proved in [10]. The proof of (<ref>) is analogous.

§ PROVING THE EXPONENTIAL ASYMPTOTIC OPTIMALITY OF TWO-SET POLICY (THEOREM <REF>)

In this section, we prove <Ref>, which states that Two-Set Policy achieve exponential asymptotic optimality under Assumptions <ref>, <ref> and <ref>.
In <Ref>, we state all the lemmas needed for proving <Ref>.
Then in <Ref>, we prove <Ref> using these lemmas.
Finally, we prove the lemmas. In <Ref>, we first prove the core lemmas that are closely related to the properties of Two-Set Policy; then in <Ref>, we prove the rest of the lemmas, which are technical and relatively routine.

[node distance=2cm and 3cm]

(lemma7) [lemma] at (-4, 0) <Ref>;
(lemma8) [lemma] at (0, 0) <Ref>;
(lemma9) [lemma] at (4, 0) <Ref>;
(lemma10) [lemma] at (2, -3) <Ref>;
(lemma11) [lemma] at (6, -3) <Ref>;
(theorem1) [theorem] at (0, -6) <Ref>;

[arrow] (lemma9) – (lemma10);
[arrow] (lemma9) – (lemma11);
[arrow] (lemma7) – (theorem1);
[arrow] (lemma8) – (theorem1);
[arrow] (lemma9.south) to[out=-90, in=0] (theorem1.east);
[arrow] (lemma10) – (theorem1);
[arrow] (lemma11) – (theorem1);

The flowchart for the proof of <Ref>. We omit <Ref> from the flowchart for readability since it is used in the proofs of all other lemmas and the theorem.

§.§ Lemma statements

In this section, we state a sequence of lemmas to prepare for proving <Ref>. The logical relationship between the lemmas and the theorem is illustrated in the flowchart <Ref>.

We first introduce a few notational conventions.
Note that the Two-Set Policy in <Ref> induces a finite-state Markov chain whose state at time $t$ is $(X_t, \Db_t, \Da_t)$. We use $\Fullstate_t$ as a shorthand for $(X_t, \Db_t, \Da_t)$.
We use $X_\infty, \Db_\infty, \Da_\infty, \Fullstate_\infty$ to represent variables following long-run average distribution: For example, for any function $f$ of $\Fullstate_t$,

\[
\E{f(\Fullstate_\infty)} \triangleq \lim_{t\to\infty} \frac{1}{T}\sum_{t=0}^{T-1} \E{f(\Fullstate_t)},
\]

where the limit always exists because $\Fullstate_t$ is a finite-state Markov chain, though it could depend on initial state $(X_0, \Db_0, \Da_0)$.

Our first lemma, <Ref>, shows that the arms in $\Db$ and $\Da$ can indeed follow the intended subroutine when $\eta$ is small enough. This lemma is proved in <Ref>.

[Subroutine conformity]lemmasubroutineconformity
Assume Assumptions <ref>, <ref> and <ref> hold.
Consider the Two-Set Policy described in <Ref> with $\eta =  |\sspa|^{-1/2}\big(\min\{y^*(\sneu, 0), y^*(\sneu, 1)\} - (|\sempty|+1)/N\big)$:
The steps in <Ref>-<ref> are well-defined. In other words, we can always let the arms in $\Db_t$ follow Optimal Local Control, let the arms in $\Da_t$ follow Exact Proportional Control, and choose actions for the arms not in $\Db_t\cup\Da_t$ to satisfy the budget constraint $\sumN A_t(i) = \btotal$.

Our next lemma, <Ref>, relates the optimality gap with the probability of $\Db_\infty \neq [N]$.
We prove it in <Ref> with a similar argument as the proofs of Theorem 3 and Theorem 8 in [3]; however, we circumvent the use of UGAP in our proof.

Assume Assumptions <ref>, <ref> and <ref> hold.
The optimality gap under the Two-Set Policy (<Ref>) is bounded as

\begin{equation}
\rrel - \rsysn \leq O(1) \cdot \Prob{\Db_\infty \neq [N]}.
\end{equation}

<Ref> implies that it suffices to show that $\Prob{\Db_\infty \neq [N]}$ is exponentially small, which we will prove via a Lyapunov argument, with the help of <Ref> below. <Ref> is proved in <Ref>.

Suppose there exists a function $V(x, \Db, \Da)$ and constants $V_{\max}$ and $\threshbar > 0$ such that for any $t\geq 0$,

\begin{equation}
\label{eq:two-set:drift-lem:bounded}
0 \leq V(\Fullstate_t) \leq V_{\max} \quad a.s.;
\end{equation}

if $V(\Fullstate_{t}) > \threshbar$,

\begin{align}
\label{eq:two-set:drift-lem:drift-1}
\Ebig{V(\Fullstate_{t+1}) - V(\Fullstate_t) \givenbig \Fullstate_t} &\leq -\gamma \\
\label{eq:two-set:drift-lem:jump-1}
\Ebig{(V(\Fullstate_{t+1}) - V(\Fullstate_t))^+  \givenbig \Fullstate_t} &= O(1/\sqrt{N});
\end{align}

if $\threshbar/2 < V(\Fullstate_t) \leq \threshbar$,

\begin{equation}
\label{eq:two-set:drift-lem:jump-2}
\Ebig{(V(\Fullstate_{t+1}) - V(\Fullstate_t))^+  \givenbig \Fullstate_t} \leq O(\exp(-C N));
\end{equation}

if $V(\Fullstate_{t}) \leq \threshbar/2$,

\begin{equation}
\label{eq:two-set:drift-lem:jump-3}
\Probbig{V(\Fullstate_{t+1}) > \threshbar/2  \givenbig \Fullstate_t} \leq O(\exp(-C' N)),
\end{equation}

where $\gamma, C, C'$ are positive constants.
Then we have

\begin{align}
\label{eq:two-set:drift-lem:conclusion}
\Prob{V(\Fullstate_{\infty}) > \threshbar} = O(\exp(-C'' N)),
\end{align}

where $C''$ is another positive constant.

The remaining lemmas prepare for constructing the Lyapunov function $V(x, \Db, \Da)$ satisfying the conditions in <Ref>.
We first define two classes of functions, $\hw=\{\hw(x, D)\}_{D\subseteq [N]}$ and $\hu=\{\hu(x, D)\}_{D\subseteq [N]}$:
For any system state $x$ and subset $D\subseteq[N]$, define

\begin{align}
\label{eq:hw-def}
\hw(x, D) &= \norm{x(D) - m(D)\statdist}_\wmat \\
\label{eq:hu-def}
\hu(x, D) &= \norm{x(D) - m(D)\statdist}_\umat.
\end{align}

These two classes of functions quantify the convergence of state distributions for subsets of arms that follow Exact Proportional Control or Optimal Local Control. Such functions are referred to as subset Lyapunov functions and are used as building blocks in the Lyapunov analysis of the prior work [10].

Next, we show several properties of $\hw$ and $\hu$. They are proved in <Ref> using properties of the two weighted $L_2$ norms.

[Properties of $\hw$ and $\hu$]lemmahwhuproperties
Assume Assumptions <ref>, <ref> and <ref> hold.
For any time step $t$ and $D\subseteq [N]$, let $X_{t+1}'$ be the system state at time $t+1$ if the arms in $D$ follow Exact Proportional Control at time $t$, then we have

\begin{align}
\label{eq:hw:drift}
\Ebig{(\hw(X_{t+1}', D) - \rhow \hw(X_t, D))^+ \givenbig X_t} &\leq O(1/\sqrt{N}) \quad a.s. \\
\label{eq:hw:high-prob}
\Probbig{\hw(X_{t+1}', D) > \rhow \hw(X_t, D) + r   \givenbig X_t} &\leq O(\exp(-C N r^2)) \quad a.s. \quad \forall r \geq 0,
\end{align}

where $\rhow = 1 - 1/(2\lamw)$, and $C$ is a positive constant.
In addition, for any $D, D'\subseteq [N]$,

\begin{align}
\label{eq:hw:strength}
\hw(x, D) &\geq \frac{1}{|\sspa|^{1/2}} \norm{x(D) - m(D)\statdist}_1 \\
\label{eq:hw:lipschitz}
\abs{\hw(x, D) - \hw(x, D')} &\leq 2\lamw^{1/2} (m(D'\backslash D) + m(D\backslash D')),
\end{align}

The same inequalities hold for $\hu$ with Exact Proportional Control and $W$ replaced by the Optimal Local Control and $U$, which we omit here.

Finally, we state two important lemmas characterizing the dynamics of the set $\Db_t$. In <Ref>, we show that the arms in $\Db_t$ are most likely to remain in $\Db_{t+1}$, so that $\Db_t$ is almost non-shrinking in the set-inclusive sense.
In <Ref>, we utilize the maximality of $\Db_t$ to show that the fraction of arms not in $\Db_t$ is small when $X_t(\Da_t)$ is close to $m(\Da_t)\statdist$ or when $X_t([N])$ is close to $\statdist$.
These two lemmas imply that $\Db_t$ tends to gradually expand and cover all arms as the two subroutines drive the states of the arms in $\Db_t$ and $\Da_t$ to the optimal distribution.
Both lemmas are proved in <Ref>.

[Almost non-shrinking]lemmaalmostnonshrinking
Assume Assumptions <ref>, <ref> and <ref> hold.
Under the Two-Set Policy (<Ref>), for any $t\geq0$,

\begin{align}
\label{eq:two-set:non-shrink-drift}
\Ebig{m(\Db_t\backslash \Db_{t+1}) \givenbig \Fullstate_t} &\leq O(m(\Db_t)\exp(-C N m(\Db_t)^2)) \quad a.s. \\
\label{eq:two-set:non-shrink-high-prob}
\Probbig{m(\Db_{t}\backslash \Db_{t+1}) > 0 \givenbig \Fullstate_t} &\leq O(\exp(-C N m(\Db_t)^2)) \quad a.s.,
\end{align}

where $C$ is a positive constant.

[Sufficient coverage]lemmasufficientcoverage
Assume Assumptions <ref>, <ref> and <ref> hold.
Consider the Two-Set Policy (<Ref>) with $\eta =  |\sspa|^{-1/2}\big(\min\{y^*(\sneu, 0), y^*(\sneu, 1)\} - (|\sempty|+1)/N\big)$. For any $t\geq0$,

\begin{align}
\label{eq:two-set:sufficient-coverage:goal-1}
1-m(\Db_t) &\leq  \frac{\lamu^{1/2}}{\eta\beta}  \hw(X_t, \Da_t) + \frac{2}{\beta N} \quad a.s.\\
\label{eq:two-set:sufficient-coverage:goal-2}
m(\Db_t) &= 1 \quad \text{ if } \hu(X_t, [N]) \leq \eta.
\end{align}

§.§ Proof of Theorem <ref>

Define the Lyapunov function $V(x, \Db, \Da)$ as

\begin{equation}\label{eq:two-set:full-lyapunov}
V(x, \Db, \Da) = \hu(x, \Db) + \hw(x, \Da) + \Ltemp (1 - m(\Db)),
\end{equation}

where $L = 2\lamu^{1/2} + 4\lamw^{1/2} - 2\lamw^{1/2}\beta$.
We claim that

\begin{equation}
\label{eq:two-set:thm:interm-goal-1}
\Probbig{V(X_\infty, \Db_\infty, \Da_\infty) > \eta} = O(\exp(-CN)),
\end{equation}

for some positive constant $C$.
To see why (<ref>) helps, note that by <Ref>,

\[
V(x, \Db, \Da) \geq  \hu(x, \Db) + 2\lamu^{1/2} (1 - m(\Db)) \geq \hu(x, [N]),
\]

so $V(X_\infty, \Db_\infty, \Da_\infty) \leq \eta$ implies $\hu(X_\infty, [N]) \leq \eta$, which further implies $\Db_\infty =[N]$ by <Ref>. Therefore, suppose we can show (<ref>), with <Ref> and the above argument, we have

\begin{align*}
\rrel - \rsysn
&\leq O(1)\cdot \Probbig{\Db_\infty \neq [N]} \\
&\leq O(1)\cdot \Probbig{\hu(X_\infty, [N]) > \eta} \\
&\leq O(1)\cdot \Probbig{V(X_\infty, \Db_\infty, \Da_\infty) > \eta} \\
&= O(\exp(-CN)),
\end{align*}

which will finish the proof.

In the rest of the proof, we will show (<ref>) by invoking <Ref> with $\threshbar = \eta$.
First, $V(x, \Db, \Da)$ is non-negative and upper bounded by $\lamu^{1/2}+\lamw^{1/2} + \Ltemp$, so the condition (<ref>) holds. To show the other conditions, we first examine the difference $V(\Fullstate_{t+1}) - V(\Fullstate_t)$.
Consider the decomposition

\begin{align}
V(\Fullstate_{t+1}) - V(\Fullstate_t)
\label{eq:two-set:thm:set-change-diff-term}
&= \big(V(X_{t+1}, \Db_{t+1}, \Da_{t+1}) - V(X_{t+1}, \Db_{t}, \Da_{t}))  \\
\label{eq:two-set:thm:state-change-diff-term}
&\mspace{20mu} + \big(V(X_{t+1}, \Db_{t}, \Da_{t}) -  V(X_{t}, \Db_{t}, \Da_{t})\big).
\end{align}

We call the difference term in RHS of (<ref>) the set-update term, and the difference term in (<ref>) the state-transition term.
We bound them separately.

For the set-update term in (<ref>), we can derive a bound fully in terms of the sizes of sets, using the Lipschitz continuity of $\hw(x, D)$ and $\hu(x, D)$ with respect to $D$ given in (<ref>) of <Ref>:

\begin{align}
\nonumber
&\mspace{20mu} V(X_{t+1}, \Db_{t+1}, \Da_{t+1}) - V(X_{t+1}, \Db_{t}, \Da_{t})  \\
\nonumber
&= \hu(X_{t+1}, \Db_{t+1}) - \hu(X_{t+1}, \Db_{t}) + \hw(X_{t+1}, \Da_{t+1}) - \hw(X_{t+1}, \Da_{t}) \\
\nonumber
&\mspace{20mu} - \Ltemp (m(\Db_{t+1}) - m(\Db_{t})) \\
\label{eq:two-set:thm:set-update-term-1}
&\leq 2\lamu^{1/2} \big(m(\Db_{t+1}\backslash \Db_t) + m(\Db_t \backslash \Db_{t+1})\big) + 2\lamw^{1/2} \big(m(\Da_{t+1}\backslash \Da_t) + m(\Da_t \backslash \Da_{t+1})\big)\\
&\mspace{20mu} - \Ltemp (m(\Db_{t+1}\backslash \Db_t) - m(\Db_{t}\backslash \Db_{t+1})).
\end{align}

We can bound $m(\Da_{t+1}\backslash \Da_t)$ and $m(\Da_t \backslash \Da_{t+1})$ in the last expression in terms of $\Db_{t+1}$ and $\Db_t$.
By definition, $\Da_{t+1} \supseteq \Da_{t} \backslash \Db_{t+1}$, and $\Da_t$ is disjoint from $\Db_t$, so it is not hard to see that

\begin{equation}
\label{eq:two-set:thm:set-update-term-2}
\Da_{t}\backslash \Da_{t+1} \subseteq \Db_{t+1}\backslash \Db_t,  \quad m(\Da_{t}\backslash \Da_{t+1}) \leq m(\Db_{t+1}\backslash \Db_t).
\end{equation}

Moreover, by definition $m(\Da_t) \geq \beta (1-m(\Db_{t})) - \frac{2}{N}$ and  $m(\Da_{t+1}) \leq \beta (1-m(\Db_{t+1}))$, so

\begin{align}
\nonumber
m(\Da_{t+1}\backslash \Da_{t}) &= m(\Da_{t}\backslash \Da_{t+1}) + m(\Da_{t+1}) - m(\Da_{t}) \\
\label{eq:two-set:thm:set-update-term-3}
&\leq m(\Db_{t+1}\backslash \Db_t) - \beta(m(\Db_{t+1}) - m(\Db_t)) + \frac{2}{N}.
\end{align}

Plugging (<ref>) and (<ref>) into (<ref>)
and rearranging the terms, we get

\begin{align}
V(X_{t+1}, \Db_{t+1}, \Da_{t+1}) - V(X_{t+1}, \Db_{t}, \Da_{t})
\label{eq:two-set:thm:set-update-term-4}
\leq 4\big(\lamu^{1/2} + \lamw^{1/2}\big) m(\Db_t \backslash \Db_{t+1}) + \frac{4\lamw^{1/2}}{N},
\end{align}

where the terms related to $m(\Db_{t+1}\backslash \Db_t)$ have been canceled out.

For the state-transition term in (<ref>), it follows directly from the definition of $V$ that

\begin{align}
\nonumber
&\mspace{20mu} V(X_{t+1}, \Db_{t}, \Da_{t}) -  V(X_{t}, \Db_{t}, \Da_{t})\\
\label{eq:two-set:thm:state-transition-term-1}
&\leq \big(\hu(X_{t+1}, \Db_t) - \hu(X_{t}, \Db_t)\big)+ \big(\hw(X_{t+1}, \Da_t) - \hw(X_{t}, \Da_t)\big)
\end{align}

Combine the above calculations, we get the following key inequality:

\begin{align}
\nonumber
& V(\Fullstate_{t+1}) - V(\Fullstate_t) \\
\label{eq:two-set:thm:v-diff-term-1}
&\leq \big(\hu(X_{t+1}, \Db_t) - \hu(X_{t}, \Db_t)\big)+ \big(\hw(X_{t+1}, \Da_t) - \hw(X_{t}, \Da_t)\big)\\
\label{eq:two-set:thm:v-diff-term-2}
&\mspace{20mu} + 4\big(\lamu^{1/2} + \lamw^{1/2}\big) m(\Db_t \backslash \Db_{t+1}) + \frac{4\lamw^{1/2}}{N}.
\end{align}

Next, we verify the drift condition (<ref>) of <Ref>, restated below: when $V(\Fullstate_t)\geq \eta$,

\begin{equation}
\label{eq:two-set:thm:interm-goal-2}
\Ebig{V(\Fullstate_{t+1}) - V(\Fullstate_t)\givenbig \Fullstate_t} \leq -\gamma + O(1/\sqrt{N}) \quad a.s.,
\end{equation}

for a constant $\gamma > 0$.
We take expectation in the bound (<ref>)-(<ref>). Because
$\Eplain{\hu(X_{t+1}, \Db_t)\givenplain \Fullstate_t} \leq \rhou \hu(X_{t}, \Db_t) + O(1/\sqrt{N})$,
$\Eplain{\hw(X_{t+1}, \Da_t)\givenplain \Fullstate_t} \leq \rhow \hw(X_{t}, \Da_t) + O(1/\sqrt{N})$ (<Ref>),
$\Eplain{m(\Db_t \backslash \Db_{t+1}) \givenplain \Fullstate_t} = O(1/\sqrt{N})$ (<Ref>), we get

\begin{equation}
\label{eq:two-set:thm:full-drift-term-1}
\Ebig{V(\Fullstate_{t+1}) -  V(\Fullstate_t) \givenbig \Fullstate_t} \leq - \big((1-\rhou) \hu(X_{t}, \Db_t) + (1-\rhow) \hw(X_{t}, \Da_t)\big) + O(1/\sqrt{N}).
\end{equation}

To show (<ref>), it remains to lower bound the term $(1-\rhou) \hu(X_{t}, \Db_t) + (1-\rhow) \hw(X_{t}, \Da_t)$ when $V(\Fullstate_t) > \eta$.
We invoke <Ref> to get

\begin{align}
\nonumber
&\leq  \hu(X_{t}, \Db_t) +  \hw(X_{t}, \Da_t) + L\Big( \frac{\lamu^{1/2}}{\eta\beta}  \hw(X_t, \Da_t) + \frac{2}{\beta N}\Big) \\
\label{eq:v-h-bound}
&= K_{Vh} \cdot \big((1-\rhou) \hu(X_{t}, \Db_t) + (1-\rhow) \hw(X_{t}, \Da_t)\big) + O(1/N),
\end{align}

for some positive constant $K_{Vh}$. So $ \Ebig{V(\Fullstate_{t+1}) -  V(\Fullstate_t) \givenbig \Fullstate_t} \leq - \eta / K_{Vh} + O(1/\sqrt{N})$ when $V(\Fullstate_t) > \eta$, which implies (<ref>).

Next, we verify the condition (<ref>) of <Ref>, restated below: when $V(\Fullstate_t) > \eta$,

\begin{equation}
\Ebig{(V(\Fullstate_{t+1}) - V(\Fullstate_t))^+  \givenbig \Fullstate_t} \leq O(1/\sqrt{N}) \quad a.s.
\end{equation}

Again, using the bound in (<ref>)-(<ref>), we have

\begin{align}
\nonumber
& \Ebig{(V(\Fullstate_{t+1}) - V(\Fullstate_t))^+ \givenbig \Fullstate_t}\\
\nonumber
&\leq \Ebig{\big(\hu(X_{t+1}, \Db_t) - \hu(X_{t}, \Db_t)\big)^+ \givenbig \Fullstate_t} + \Ebig{\big(\hw(X_{t+1}, \Da_t) - \hw(X_{t}, \Da_t)\big)^+ \givenbig \Fullstate_t} \\
\nonumber
&\mspace{20mu} + 4\big(\lamu^{1/2} + \lamw^{1/2}\big) \Ebig{m(\Db_t \backslash \Db_{t+1})\givenbig \Fullstate_t} + \frac{4\lamw^{1/2}}{N}.
\end{align}

Each term on the RHS above is $O(1/\sqrt{N})$: In particular, by <Ref>. $\E{m(\Db_t \backslash \Db_{t+1}) \givenbig \Fullstate_t} = O(m(\Db_{t})O(\exp(-CN m(\Db_{t})))) = O(1/\sqrt{N})$.

Finally, we show that when $V(\Fullstate_t) \leq \eta$,

\begin{equation}
\label{eq:two-set:thm:interm-goal-4}
\Probbig{V(\Fullstate_{t+1}) > \rhoFinal V(\Fullstate_t) + r \givenbig \Fullstate_t} = O(\exp(-CN)) \quad a.s.,
\end{equation}

where $r = (1-\rhoFinal)\eta / 2$, $\rhoFinal$ is a constant in $(0,1)$ and $C$ is a positive constant. It is not hard to see that (<ref>) implies the conditions (<ref>) and (<ref>) of <Ref>,
We rearrange the bound in (<ref>)-(<ref>) and apply (<ref>) to get

\begin{align*}
\nonumber
& V(\Fullstate_{t+1}) - \big(1-\frac{1}{K_{Vh}}\big) V(\Fullstate_t) \\
&\leq \big(\hu(X_{t+1}, \Db_t) - \rhou \hu(X_{t}, \Db_t)\big)+ \big(\hw(X_{t+1}, \Da_t) - \rhow\hw(X_{t}, \Da_t)\big)\\
&\mspace{20mu} + 4\big(\lamu^{1/2} + \lamw^{1/2}\big) m(\Db_t \backslash \Db_{t+1}) + O(1/N).
\end{align*}

For $N$ larger than a certain constant threshold, the $O(1/N)$ term in the last expression is less than $r/4$.
Using a union bound and invoking <Ref>, we get

\begin{align*}
&\mspace{20mu} \ProbBig{V(\Fullstate_{t+1}) > \big(1-\frac{1}{K_{Vh}}\big) V(\Fullstate_t) + r\givenBig \Fullstate_t} \\
&\leq \Probbig{\hu(X_{t+1}, \Db_t) > \rhou \hu(X_{t}, \Db_t) + r/4 \givenbig \Fullstate_t} \\
&\mspace{20mu} + \Probbig{\hw(X_{t+1}, \Da_t) > \rhow \hw(X_{t}, \Da_t) + r/4 \givenbig \Fullstate_t} \\
&\mspace{20mu} + \Probbig{4\big(\lamu^{1/2} + \lamw^{1/2}\big) m(\Db_t \backslash \Db_{t+1})  > r/4 \givenbig \Fullstate_t} + \indibrac{O(1/N) > r/4} \\
&\leq O(\exp(-CN)) + O(\exp(-CN m(\Db_t))).
\end{align*}

By <Ref>, $V(\Fullstate_t) < \eta$ implies $m(\Db_t) =1$, so $O(\exp(-CN m(\Db_t))) = O(\exp(-CN))$.
Therefore,  (<ref>) holds with $\rhoFinal = 1 - 1/(K_{Vh})$.

§.§ Proofs of the core lemmas

In this proof, we need to verify the condition of following the $\lppriority$, (<ref>), for the arms in $\Db_t$.
Then we need to show that $A_t(i)$ for $i\notin \Db_t \cup \Da_t$ can be chosen such that $\sumN A_t(i) = \btotal$, if the arms in $\Db_t$ and $\Da_t$ follow Optimal Local Control and Exact Proportional Control, respectively.

The condition (<ref>) specialized to $\Db_t$ are stated as follows.

\begin{align}
\label{eq:two-set:db-condition-1}
\sum_{s\neq \sneu} \pibs(1|s) X_t(\Db_t, s) &\leq  \alpha m(\Db_t) - (|\sempty|+1)/ N, \\
\label{eq:two-set:db-condition-2}
\sum_{s\neq \sneu} \pibs(0|s) X_t(\Db_t, s) &\leq (1-\alpha)m(\Db_t) - (|\sempty|+1)/N.
\end{align}

We first show (<ref>). Observe that $\sums \pibs(1|s) \statdist(s) = \sumsa y^*(s,a) = \alpha$, so

\begin{align}
\nonumber
\alpha m(\Db_t) - \sum_{s\neq \sneu} \pibs(1|s) X_t(\Db_t, s)
&= y^*(\sneu, 1)m(\Db_t)  - \sum_{s\neq\sneu} \pibs(1|s) \big(m(\Db_t)\statdist(s) - X_t(\Db_t, s)\big) \\
\nonumber
&\geq y^*(\sneu, 1)m(\Db_t) - \sum_{s\neq\sneu} \abs{m(\Db_t)\statdist(s) - X_t(\Db_t, s)} \\
\nonumber
&\geq y^*(\sneu, 1)m(\Db_t) - \norm{m(\Db_t)\statdist - X_t(\Db_t, s)}_1 \\
\nonumber
&\geq y^*(\sneu, 1)m(\Db_t) - |\sspa|^{1/2}\norm{m(\Db_t)\statdist - X_t(\Db_t, s)}_\umat,
\end{align}

where to get the last expression, we use the argument that $U$'s eigenvalues are at least $1$, so $\norm{v}_1\leq |\sspa|^{1/2} \norm{v}_2 \leq  |\sspa|^{1/2} \norm{v}_U$ for any $|\sspa|$-dim vector $v$.
To bound $\norm{m(\Db_t)\statdist - X_t(\Db_t, s)}_\umat$, recall that $\Db_t$ is chosen such that $\slkb(X_t, \Db_t) \geq 0$, so by definition, we have

\[
\norm{X_t(\Db_t) - m(\Db_t) \statdist}_U \leq \eta m(\Db_t) \leq |\sspa|^{-1/2} \big(y^*(\sneu, 1) - (|\sempty|+1)/ N\big) m(\Db_t).
\]

Therefore, combining the above calculations,

\[
\alpha m(\Db_t) - \sum_{s\neq \sneu} \pibs(1|s) X_t(\Db_t, s) \geq (|\sempty|+1)/ N,
\]

which implies (<ref>). (<ref>) can be similarly proved.

It remains to show that $\sumN A_t(i) = \btotal$ can be satisfied.
If $\Db = [N]$, by definition we have $\sum_{i\in \Db_t} A_t(i) = \btotal$. Thus we only need to consider the case when $\Db \subsetneq [N]$.
By definition, $\sum_{i\in \Db_t} A_t(i) \in \{\lfloor\alpha |\Db_t|\rfloor, \lceil\alpha |\Db_t|\rceil\}$.
Because $A_t(i)\in \{0,1\}$ for $i\in \Da_t$, we have

\[
\floor{\alpha |\Db_t|} \leq \sum_{i\in \Db_t\cup\Da_t} A_t(i)\leq \ceil{\alpha |\Db_t|} + |\Da_t|.
\]

Since $\sum_{i\notin \Db_t\cup\Da_t} A_t(i)$ can be anything integer between $0$ and $N - |\Db_t| - |\Da_t|$, it suffices if

\begin{align*}
\floor{\alpha |\Db_t|} + (N - |\Db_t| - |\Da_t|) &\geq \btotal  \quad a.s.  \\
\ceil{\alpha |\Db_t|} + |\Da_t| &\leq \btotal \quad a.s.
\end{align*}

These two inequalities can be verified using the facts that $|\Da_t| = \big(\big\lfloor \beta (N -|\Db_{t}|)\big\rfloor - 1\big)^+$, $\btotal \in \{\floor{\alpha N}, \ceil{\alpha N}\}$ and $\Db_t\subsetneq [N]$.

By the definition of $\Db_{t+1}$, it suffices to construct a set $\overline{D}\subseteq \Db_t$ based on the realization of $X_{t+1}$ such that

\begin{align}
\label{eq:two-set:non-shrinking:interm-req-1}
\slkb(X_{t+1}, \overline{D}) &\geq 0 \quad a.s.\\
\label{eq:two-set:non-shrinking:interm-req-2}
\Ebig{(m(\Db_t \backslash \overline{D}) \givenbig \Fullstate_t} &\leq O(\exp(-C N m(\Db_t)^2)) \quad a.s. \\
\label{eq:two-set:non-shrinking:interm-req-3}
\Probbig{m(\Db_{t}\backslash\overline{D}) > 0 \givenbig \Fullstate_t} &\leq O(\exp(-C N m(\Db_t)^2)) \quad a.s.
\end{align}

Consider taking $\overline{D}$ as:

\begin{equation}\label{eq:two-set:non-shrinking:dbar-def}
\overline{D} =
\begin{cases}
\Db_t \quad &\text{ if } \slkb(X_{t+1}, \Db_t) \geq 0\\
\emptyset \quad &\text{ if } \slkb(X_{t+1}, \Db_t) < 0.
\end{cases}
\end{equation}

With $\overline{D}$ defined in (<ref>), (<ref>) is obviously true. To show (<ref>) and (<ref>), observe that

\begin{align}
\nonumber
\Ebig{(m(\Db_t\backslash\overline{D}))^+ \givenbig X_t, \Db_t, \Da_t}
&= m(\Db_t) \Prob{\slkb(X_{t+1}, \Db_t) < 0} \\
\nonumber
\Probbig{m(\Db_{t}\backslash\overline{D}) > 0 \givenbig \Fullstate_t}
&= \Prob{\slkb(X_{t+1}, \Db_t) < 0}.
\end{align}

Therefore, it remains to bound $\Prob{\slkb(X_{t+1}, \Db_t) < 0}$. By definition, $\slkb(x, D) = \eta m(D) -   \normplain{x(D)- m(D) \statdist}_\umat = \eta m(D) - \hu(x, D)$, so

\begin{align}
\nonumber
\Prob{\slkb(X_{t+1}, \Db_t) < 0}
&= \Prob{\hu(X_{t+1},\Db_t) > \eta m(\Db_t)},
\end{align}

Because $\slkb(X_t, \Db_t) \geq 0$, we have $\hu(X_t, \Db_t) \leq \eta m(\Db_t)$.

\begin{align}
\nonumber
\hu(X_{t+1},\Db_t) &\leq  \hu(X_{t+1},\Db_t) - \rhou \hu(X_t, \Db_t) + \rhou \eta m(\Db_t)  \\
\nonumber
\Prob{\hu(X_{t+1}, \Db_t) \geq \eta m(\Db_t)}
&\leq \Prob{\hu(X_{t+1}, \Db_t) - \hu(X_t, \Db_t) \geq (1-\rhou) \eta m(\Db_t)} \\
\nonumber
&= O(\exp(-C' N (1-\rhou)^2 \eta^2 m(\Db_t)^2)),
\end{align}

where the last inequality comes from <Ref> and $C'$ is a positive constant. Because $(1-\rhou)^2 \eta^2 > 0$ and is a constant,  $\Prob{\slkb(X_{t+1}, \Db_t) < 0} = O(\exp(-C N m(\Db_t)^2))$ with $C = C'(1-\rhou)^2 \eta^2$.

First, we show that

\begin{equation}
\label{eq:two-set:sufficient-coverage:interm-goal-1}
\hw(X_t, \Da_t) \geq  \lamu^{-1/2} \eta m(\Da_t).
\end{equation}

If we can show (<ref>), combined with the fact that $m(\Da_{t})= \big(\big\lfloor \beta (N -|\Db_{t}|)\big\rfloor - 1\big)^+ \big/ N$, we get

\[
\hw(X_t, \Da_t) \geq  \lamu^{-1/2} \eta \big(\beta(N -|\Db_{t}|)-2\big)/N =  \lamu^{-1/2} \eta \beta(1-m(\Db_t)) -   2\lamu^{-1/2} \eta/N,
\]

which implies (<ref>).

To show (<ref>), we only need to consider the case where $m(\Da) \neq 0$. Suppose, to get a contradiction, that $ \hw(X_t, \Da_t) <  \lamu^{-1/2} \eta m(\Da_t)$. Because $\norm{v}_\umat \leq \lamu^{1/2}\norm{v}_2 \leq \lamu^{1/2}\norm{v}_\wmat$ for any $v\in \R^{|\sspa|}$,

\begin{align}
\nonumber
\normbig{X_t(\Da_t) - m(\Da_t) \statdist}_\umat
&\leq \lamu^{1/2}\normbig{X_t(\Da_t) - m(\Da_t) \statdist}_\wmat
<  \eta m(\Da_t).
\end{align}

Combined with the triangular inequality and the definition of $\Db_t$, we have

\begin{align*}
\normbig{X_t(\Db_t\cup\Da_t) - m(\Db_t\cup\Da_t) \statdist}_\umat &\leq \normbig{X_t(\Db_t) - m(\Db_t) \statdist}_\umat  + \normbig{X_t(\Da_t) - m(\Da_t) \statdist}_\umat \\
&< \eta m(\Db_t) + \eta m(\Da_t).
\end{align*}

Consequently, $\Db_t\cup\Da_t$ is a strict superset of $\Db_t$ (thus also a superset of $\Dtemp$) such that $\slkb(X_t, \Db_t\cup\Da_t) \geq 0$. So the arms in $\Da_t$ should have instead been taken into $\Db_t$ when the algorithm chooses $\Db_t$ to be the maximal set such that $\Db_t\supseteq \Dtemp$ and $\slkb(X_t, \Db_t) \geq 0$, which causes a contradiction. Therefore, (<ref>) holds.

Finally, we show (<ref>).
When $\hu(X_t, [N]) \leq \eta$, we have $\slkb(X_t, [N])\geq 0$. Again, because $\Db_t$ is chosen as the largest-cardinality set containing $\Dtemp$, we must have $\Db_t = [N]$.

§.§ Proofs of the technical lemmas

We will first bound the norm of the long-run expectation of the $X_t([N]) - \statdist$ using $\Probbig{\Db_\infty \neq [N]}$, and then derive the optimality gap bound.

Let $Q = (I - \Phi)^{-1}$, which is well-defined because all eigenvalues of $\Phi$ have moduli strictly less than $1$ (<Ref>).
By direct calculation, we have that for any $v\in\Delta(\sspa)$,

\begin{equation}\label{eq:phi-stein-equation}
(v - \statdist) Q - (v - \statdist)\Phi Q = v - \statdist.
\end{equation}

We define the positive constant $\lamq$ as

\[
\lamq = \max\Big\{\max_{v\in\Delta(\sspa)} \norm{(v-\statdist) Q}_1,  \max_{v\in\Delta(\sspa)} \norm{(v-\statdist) \Phi Q}_1\Big\},
\]

which must be finite because $\Delta(\sspa)$ is a compact set.

Letting $v = X_t([N])$ and taking the expectation (for random vectors) in (<ref>), we get

\begin{align}
\Ebig{X_t([N]) - \statdist}
\nonumber
&= \Ebig{(X_t([N])-\statdist)Q - (X_t([N])-\statdist)\Phi Q} \\
\nonumber
&= \Ebig{(X_{t+1}([N])-\statdist)Q - (X_t([N])-\statdist)\Phi Q} \\
\nonumber
&\mspace{20mu} - \Ebig{(X_{t+1}([N])-\statdist)Q - (X_t([N])-\statdist)Q}.
\end{align}

Using the fact that

\[
\lim_{T\to\infty} \frac{1}{T} \sum_{t=0}^{T-1} \Ebig{(X_{t+1}([N])-\statdist)Q - (X_t([N])-\statdist)Q} = 0,
\]

we get

\begin{equation}\label{eq:two-set:state-diff-steins-diff}
\lim_{T\to\infty} \frac{1}{T} \sum_{t=0}^{T-1} \Ebig{X_t([N]) - \statdist} = \lim_{T\to\infty} \frac{1}{T} \sum_{t=0}^{T-1} \Ebig{(X_{t+1}([N])-\statdist)Q - (X_t([N])-\statdist)\Phi Q}.
\end{equation}

Theerefore, it remains to bound $\Ebig{(X_{t+1}([N])-\statdist)Q - (X_t([N])-\statdist)\Phi Q}$ for every time step $t$.

Let $\Fullstate_t = (X_t, \Db_t, \Da_t)$.
We will bound $\Ebig{(X_{t+1}([N])-\statdist)Q - (X_t([N]-\statdist)\Phi Q \givenbig \Fullstate_t}$ by decomposing it into two terms:

\begin{align}
\nonumber
&\mspace{20mu} \Ebig{(X_{t+1}([N])-\statdist)Q - (X_t([N])-\statdist)\Phi Q \givenbig \Fullstate_t} \\
\label{eq:steins-difference-lp-is-N}
&= \Ebig{(X_{t+1}([N])-\statdist)Q - (X_t([N])-\statdist)\Phi Q  \givenbig \Fullstate_t} \indibrac{\Db_t = [N]}\\
\label{eq:steins-difference-lp-not-N}
&\mspace{20mu} + \Ebig{(X_{t+1}([N])-\statdist)Q - (X_t([N])-\statdist)\Phi Q  \givenbig \Fullstate_t} \indibrac{\Db_t \neq [N]}.
\end{align}

For the term in (<ref>), note that by <Ref>, $\Ebig{X_{t+1}([N]) -\statdist \givenbig \Fullstate_t} = (X_t([N])-\statdist)\Phi$ when $\Db_t = [N]$, so

\begin{equation*}
\Ebig{(X_{t+1}([N])-\statdist)Q \givenbig \Fullstate_t}  \indibrac{\Db_t = [N]}
= (X_t([N])-\statdist)\Phi Q \indibrac{\Db_t = [N]},
\end{equation*}

which implies that the term in (<ref>) is zero. For the term in (<ref>), its $L_1$ norm can be bounded as

\begin{align}
\nonumber
&\mspace{20mu} \norm{\Ebig{(X_{t+1}([N])-\statdist)Q - (X_t([N])-\statdist)\Phi Q  \givenbig \Fullstate_t}}_1 \indibrac{\Db_t \neq [N]} \\
\nonumber
&\leq  \Ebig{\norm{(X_{t+1}([N])-\statdist)Q - (X_t([N])-\statdist)\Phi Q}_1  \givenbig \Fullstate_t} \indibrac{\Db_t \neq [N]} \\
\nonumber
&\leq \Ebig{\norm{(X_{t+1}([N])-\statdist)Q}_1 + \norm{(X_t([N])-\statdist)\Phi Q}_1  \givenbig \Fullstate_t} \indibrac{\Db_t \neq [N]}  \\
\label{eq:steins-difference-bound-by-lamq}
&\leq 2\lamq \indibrac{\Db_t \neq [N]},
\end{align}

where (<ref>) follows from the definition of $\lamq$ and the fact that $X_{t+1}([N]), X_t([N]) \in \Delta(\sspa)$.
Combining the calculations of the two terms in (<ref>) and (<ref>), we have

\begin{align}
\nonumber
&\mspace{20mu} \norm{\Ebig{(X_{t+1}([N])-\statdist)Q - (X_t([N])-\statdist)\Phi Q}}_1 \\
\nonumber
&=   \EBig{ \norm{\Ebig{(X_{t+1}([N])-\statdist)Q - (X_t([N])-\statdist)\Phi Q  \givenbig \Fullstate_t}}_1 \indibrac{\Db_t \neq [N]} }  \\
\label{eq:two-set:steins-diff-expect-ell1-bound}
&\leq 2\lamq \Prob{\Db_t \neq [N]}.
\end{align}

Plugging (<ref>) into (<ref>), we get

\begin{align}
\nonumber
\normbig{\Ebig{X_\infty([N]) - \statdist} }_1
&\leq \lim_{T\to\infty} \frac{1}{T} \sum_{t=0}^{T-1} \norm{\Ebig{(X_{t+1}([N])-\statdist)Q - (X_t([N])-\statdist)\Phi Q}} \\
\label{eq:two-set:state-diff-expect-prob-bound}
&\leq 2\lamq \Prob{\Db_\infty \neq [N]}.
\end{align}

where the expectations or probabilities of variables or events with subscript $\infty$ are shorthands for the long-run time averages.

Now we start to bound the optimality gap. We first rewrite the long-run average reward $\rsysn$ in terms of an expectation with respect to the states:

\begin{align}
\nonumber
\rsysn
&= \Ebig{r(\veS_\infty, \veA_\infty)}
=  \E{r^\pi(X_\infty, \Db_\infty, \Da_\infty)},
\end{align}

where $r^\pi(x, \Db, \Da)$ is the expected reward at a given state under the Two-Set Policy:

\[
r^\pi(x, \Db, \Da) \triangleq \Ebig{r(\veS_t, \veA_t) \givenbig X_t = x, \Db_t = \Db, \Da_t = \Da}.
\]

Note that the RHS of $r^\pi(x, \Db, \Da)$'s definition does not depend on $t$ because the Two-Set Policy is stationary and Markovian policy with state $(X_t, \Db_t, \Da_t)$.

Next, we examine the value of $r^\pi(x, \Da, \Db)$ when $\Db = [N]$ and $\Da = \emptyset$. Let $Y_t(s,a)$ be the random variable representing the fraction of arms in state $s$ taking action $a$ at time $t$. Then

\begin{align*}
r^\pi(X_t, [N], \emptyset)
&= \sumsa r(s,a) \Ebig{Y_t(s,a) \givenbig X_t, \Db_t = [N], \Da_t = \emptyset}.
\end{align*}

By the definition of Optimal Local Control,

\begin{align}
\Ebig{Y_t(s,a) \givenbig X_t, \Db_t = [N], \Da_t = \emptyset} &= \pibs(a|s) X_t([N], s) \quad s\neq \sneu \\
\Ebig{Y_t(\sneu,1) \givenbig X_t, \Db_t = [N], \Da_t = \emptyset, } &= \alpha - \sum_{s\neq \sneu} \pibs(1|s) X_t([N], s)  \\
\Ebig{Y_t(\sneu,0) \givenbig X_t, \Db_t = [N], \Da_t = \emptyset, } &= 1 - \alpha - \sum_{s\neq \sneu} \pibs(0|s) X_t([N], s).
\end{align}

\begin{equation}\label{eq:lp-priority-inst-reward}
\begin{aligned}
r^\pi(X_t, [N], \emptyset)
&= \sum_{s\neq \sneu, a\in\aspa} \big(r(s,a) - r(\sneu,a)\big) \pibs(a|s) X_t([N],s) \\
&\mspace{20mu} + \alpha r(\sneu,1) + (1-\alpha) r(\sneu,0),
\end{aligned}
\end{equation}

Let $\hat{r}(v)$ be the affine function on $\simplex(\sspa)$ such that $\hat{r}(X_t([N])) = r^\pi(X_t, \emptyset, [N])$:

\[
\hat{r}(v) \triangleq \sum_{s\neq \sneu, a\in\aspa} \big(r(s,a) - r(\sneu,a)\big) \pibs(a|s) v(s) + \alpha r(\sneu,1) + (1-\alpha) r(\sneu,0).
\]

It is easy to see that $\hat{r}(v)$ is $4\rmax$-Lipschitz continuous with respect to $L_1$-norm, and $\hat{r}(\statdist) = \rrel$.

With the above definitions, we can lower bound $\rsysn$ as

\begin{align*}
\rsysn
&= \Ebig{r^\pi(X_\infty, \Db_\infty, \Da_\infty)}   \\
&\geq \Ebig{ \hat{r}(X_\infty([N])) \indibrac{\Db_\infty = [N]}} - \rmax \Probbig{\Db_\infty \neq [N]} \\
&= \Ebig{ \hat{r}(X_\infty([N]))}  - 2\rmax \Probbig{\Db_\infty \neq [N]} \\
&= \hat{r}\big(\Ebig{X_\infty([N])}\big) - 2\rmax \Probbig{\Db_\infty \neq [N]} \\
&= \hat{r}(\statdist) + \hat{r}\big(\Ebig{X_\infty([N])}\big) - \hat{r}\big(\statdist\big)   - 2\rmax \Probbig{\Db_\infty \neq [N]} \\
&\geq  \rrel - 4\rmax \norm{\Ebig{X_\infty([N]) - \statdist}}_1 - 2\rmax \Probbig{\Db_\infty \neq [N]} \\
&\geq \rrel - 2\rmax(4 \lamq + 1) \Probbig{\Db_t \neq [N]}.
\end{align*}

Because $2\rmax(4 \lamq + 1) = O(1)$, we have $\rrel - \rsysn \leq O(1) \cdot  \Probbig{\Db_t \neq [N]}$.

Let $f(a) = (a-\threshbar/2)^+$.
To obtain an bound on $\Prob{V(\Fullstate_\infty) > \threshbar}$, we will calculate $\E{f(V(\Fullstate_{t+1}))} - \E{f(V(\Fullstate_{t}))}$ for each time step $t$ and invoking the fact that

\begin{equation}\label{eq:bar}
\lim_{T\to\infty} \frac{1}{T} \sum_{t=0}^{T-1} \big(\E{f(V(\Fullstate_{t+1}))} - \E{f(V(\Fullstate_{t}))}\big) = 0.
\end{equation}

We will bound $\E{f(V(\Fullstate_{t+1}))} - \E{f(V(\Fullstate_{t}))}$ in three cases,  $V(\Fullstate_{t}) > \threshbar$, $\threshbar/2 < V(\Fullstate_{t}) \leq \threshbar$, and  $V(\Fullstate_{t}) \leq \threshbar/2$.

First, if $V(\Fullstate_{t}) > \threshbar$,

\begin{align}
\nonumber
&\mspace{20mu} \Ebig{f(V(\Fullstate_{t+1})) \givenbig \Fullstate_t} - f(V(\Fullstate_t)) \\
\nonumber
&= \Ebig{\big(V(\Fullstate_{t+1}) - V(\Fullstate_t)\big)\indibrac{V(\Fullstate_{t+1}) > V(\Fullstate_t)} \givenbig \Fullstate_t}  \\
\nonumber
&\mspace{20mu} + \Ebig{\big(V(\Fullstate_{t+1}) - V(\Fullstate_t)\big) \indibrac{\threshbar/2 < V(\Fullstate_{t+1}) \leq V(\Fullstate_t)} \givenbig \Fullstate_t}  \\
\nonumber
&\mspace{20mu} + \Ebig{\big(\threshbar/2 - V(\Fullstate_t)\big) \indibrac{V(\Fullstate_{t+1}) \leq \threshbar/2} \givenbig \Fullstate_t}  \\
\nonumber
&\leq \Ebig{\big(V(\Fullstate_{t+1}) - V(\Fullstate_t)\big)\indibrac{V(\Fullstate_{t+1}) > V(\Fullstate_t)} \givenbig \Fullstate_t}  \\
\label{eq:two-set:drift-lem:interm-1}
&\mspace{20mu} + \frac{1}{2} \Ebig{\big(V(\Fullstate_{t+1}) - V(\Fullstate_t)\big) \indibrac{ V(\Fullstate_{t+1}) \leq V(\Fullstate_t)} \givenbig \Fullstate_t}  \\
&= \frac{1}{2}\Ebig{\big(V(\Fullstate_{t+1}) - V(\Fullstate_t)\big)^+ \givenbig \Fullstate_t}  + \frac{1}{2} \Ebig{V(\Fullstate_{t+1}) - V(\Fullstate_t) \givenbig \Fullstate_t} \\
\label{eq:two-set:drift-lem:interm-2}
&\leq -\frac{1}{2}\gamma  + O(1/\sqrt{N}).
\end{align}

where the inequality in (<ref>) follows uses the fact that when $V(\Fullstate_{t}) \geq \threshbar$, $V(\Fullstate_{t+1}) \geq 0$,

\[
\frac{1}{2}\big(V(\Fullstate_{t+1}) - V(\Fullstate_t)\big) - (\threshbar/2 - V(\Fullstate_t)) = \frac{1}{2}\big(V(\Fullstate_{t+1}) + V(\Fullstate_t)\big) - \threshbar/2 \geq 0;
\]

the inequality in (<ref>) applies the properties of $V$ assumed in (<ref>) and (<ref>).

Next, consider the case of $\threshbar/2 < V(\Fullstate_{t}) \leq \threshbar$.

\begin{align}
\nonumber
&\mspace{20mu} \Ebig{f(V(\Fullstate_{t+1})) \givenbig \Fullstate_t} - f(V(\Fullstate_t)) \\
\nonumber
&\leq  \Ebig{\big(f(V(\Fullstate_{t+1})) - f(V(\Fullstate_t))\big) \indibrac{V(\Fullstate_{t+1}) > V(\Fullstate_t)} \givenbig \Fullstate_t}  \\
\nonumber
&=  \Ebig{\big(V(\Fullstate_{t+1}) - V(\Fullstate_t)\big)^+ \givenbig \Fullstate_t} \\
\label{eq:two-set:drift-lem:interm-3}
&\leq O(\exp(-C N)),
\end{align}

where (<ref>) follows from (<ref>).

Finally, consider the case of $V(\Fullstate_{t}) \leq \threshbar/2$.

\begin{align}
\nonumber
&\mspace{20mu} \Ebig{f(V(\Fullstate_{t+1})) \givenbig \Fullstate_t} - f(V(\Fullstate_t)) \\
\nonumber
&= \Ebig{(V(\Fullstate_{t+1})-\threshbar/2)^+ \givenbig \Fullstate_t} \\
\nonumber
&\leq V_{\max} \Probbig{V(\Fullstate_{t+1})-\threshbar/2 \givenbig \Fullstate_t} \\
\label{eq:two-set:drift-lem:interm-4}
&= O(\exp(-C' N)),
\end{align}

where (<ref>) follows from (<ref>).

Combining the above calculations,

\begin{equation*}
\Ebig{f(V(\Fullstate_{t+1}))} - \E{f(V(\Fullstate_t))}
\leq \Big(-\frac{1}{2}\gamma + O(1/\sqrt{N}) \Big) \Probbig{V(\Fullstate_{t}) > \threshbar} + O(\exp(-\min(C, C') N)).
\end{equation*}

Rearranging the terms, we get

\begin{equation}
\label{eq:two-set:drift-lem:interm-5}
\Probbig{V(\Fullstate_{t}) > \threshbar} \leq \frac{1}{\gamma / 2 - O(1/\sqrt{N})}\Big( \Ebig{f(V(\Fullstate_{t+1}))} - \E{f(V(\Fullstate_t))} + O(\exp(-\min(C, C') N))\Big),
\end{equation}

for any time step $t$.
Plugging the bound (<ref>) into (<ref>), we get

\[
\Probbig{V(\Fullstate_{\infty}) > \threshbar} = O(\exp(-\min(C, C') N)).
\qedhere
\]

Proving (<ref>) and (<ref>)
By the pseudo-contraction properties of $\wmat$-weighted norms given in <Ref>,

\begin{align*}
\rhow \hw(X_t, D)
&=  \rhow \norm{X_{t}(D) - m(D)\statdist}_\wmat \\
&\geq \norm{(X_{t}(D) - m(D)\statdist)P_\pibs }_\wmat \\
& = \norm{X_{t}(D)P_\pibs - m(D)\statdist}_\wmat.
\end{align*}

\begin{align*}
\hw(X_{t+1}', D) - \rhow \hw(X_t, D)
&\leq  \norm{X_{t+1}'(D)- m(D)\statdist}_\wmat -  \norm{X_{t}(D)P_\pibs - m(D)\statdist}_\wmat \\
&\leq  \norm{X_{t+1}'(D)- X_t(D) P_{\pibs}}_\wmat \\
&\leq \lamw^{1/2} \norm{X_{t+1}'(D)- X_t(D) P_{\pibs}}_2 \\
&\leq \lamw^{1/2} \norm{X_{t+1}'(D)- X_t(D) P_{\pibs}}_1.
\end{align*}

Therefore, it suffices to show that

\begin{align}
\label{eq:hw:interm-goal-1}
\Ebig{ \norm{X_{t+1}'(D)- X_t(D) P_{\pibs}}_1 \givenbig X_t}
&= O(1/\sqrt{N}) \quad a.s.  \\
\label{eq:hw:interm-goal-2}
\Prob{\norm{X_{t+1}'(D)- X_t(D) P_{\pibs}}_1 > r \givenbig X_t}
&= O(\exp(-C\lamw N r^2))   \quad a.s. \quad \forall r\geq 0.
\end{align}

To show (<ref>) and (<ref>), we examine the distribution of $\normbig{X_{t+1}'(D)- X_t(D) P_{\pibs}}_1$ given $X_t$.
Note that by <Ref>, $\Eplain{X_{t+1}'(D) \givenplain X_t} = X_t(D) P_{\pibs}$.

\begin{align}
\nonumber
\normbig{X_{t+1}'(D)- X_t(D) P_{\pibs}}_1
&= \normbig{X_{t+1}'(D) - \E{X_{t+1}'(D) \givenplain X_t} }_1   \\
\label{eq:hw:interm-3}
&\leq \normbig{X_{t+1}'(D) - \E{X_{t+1}'(D) \givenplain X_t, (A_t(i))_{i\in D}} }_1 \\
\label{eq:hw:interm-4}
&\mspace{20mu}  +  \normbig{ \E{X_{t+1}'(D) \givenplain X_t, (A_t(i))_{i\in D}} - \E{X_{t+1}'(D) \givenplain X_t} }_1.
\end{align}

We first derive an expectation bound and a high-probability bound for (<ref>).
Note that for each $s\in\sspa$,

\begin{align*}
X_{t+1}'(D, s) = \frac{1}{N} \sum_{i\in D} \indibrac{S_{t+1}'(i) = s},
\end{align*}

where $S_{t+1}'(i)$ denote the state of arm $i$ at time $t+1$ if the arms in $D$ follow Exact Proportional Control at time $t$.
Given $X_t$ and $(A_t(i))_{i\in D}$, $\indibrac{S_{t+1}'(i) = s}$ for $i\in D$ are independent Bernoulli random variables, each with mean $P(S_t(i), A_t(i), s)$.
By Cauchy-Schwartz,

\begin{align}
\nonumber
&\mspace{20mu} \EBig{\absBig{ X_{t+1}'(D, s) - \Eplain{X_{t+1}'(D, s) \givenplain X_t, (A_t(i))_{i\in D}}} \givenBig X_t, (A_t(i))_{i\in D}}  \\
\nonumber
&\leq  \Var{X_{t+1}'(D, s) \givenBig X_t, (A_t(i))_{i\in D}}^{1/2} \\
\nonumber
&= \frac{1}{N} \left(\sum_{i\in D} \Var{\indibrac{S_{t+1}'(i) = s}, \givenBig X_t, (A_t(i))_{i\in D}} \right)^2 \\
\nonumber
&\leq \frac{1}{\sqrt{N}}.
\end{align}

By Hoeffding inequality, for all $r\geq 0$,

\begin{align}
\nonumber
\ProbBig{\absBig{ X_{t+1}'(D, s) - \Eplain{X_{t+1}'(D, s) \givenplain X_t, (A_t(i))_{i\in D}}} > r \givenBig X_t, (A_t(i))_{i\in D}} \leq 2\exp\big(-2N r^2\big).
\end{align}

Therefore, summing over all $s\in\sspa$ and taking expectations over $(A_t(i))_{i\in D}$, we get

\begin{align}
\label{eq:hw:interm-5}
\EBig{\normBig{ X_{t+1}'(D, s) - \Eplain{X_{t+1}'(D, s) \givenplain X_t, (A_t(i))_{i\in D}}}_1 \givenBig X_t} &\leq \frac{|\sspa|}{\sqrt{N}} \\
\label{eq:hw:interm-6}
\ProbBig{\normBig{ X_{t+1}'(D, s) - \Eplain{X_{t+1}'(D, s) \givenplain X_t, (A_t(i))_{i\in D}}}_1 > r \givenBig X_t} &\leq 2|\sspa|\exp\Big(-\frac{2N r^2}{|\sspa|^2}\Big),
\end{align}

where (<ref>) comes from applying a union bound.

Next, we bound the term in (<ref>).
For any $s\in\sspa$, we have

\begin{align*}
\E{X_{t+1}'(D, s) \givenplain X_t, (A_t(i))_{i\in D}} &= \sum_{s'\in\sspa, a\in\aspa} \pibs(a|s) X_t(D, s') P(s', a, s) \\
\E{X_{t+1}'(D, s) \givenplain X_t}
&= \sum_{s'\in\sspa, a\in\aspa} \Gamma_t(s',a) P(s', a, s),
\end{align*}

where $N\cdot\Gamma_t(s,a)$ represents the number of arms in $D$ with state $s$ taking action $a$ at time $t$. By the definition of Exact Proportional Control (<Ref>),

\[
\Gamma_t(s',1) \in \{\ceil{N\pibs(1|s') X_t(D, s')}/N, \floor{N\pibs(1|s') X_t(D, s')}/N \},
\]

and $\Gamma_t(s,0) = X_t(D, s) - \Gamma_t(s,1)$.

\begin{align}
\nonumber
&\mspace{20mu} \absbig{\E{X_{t+1}'(D, s) \givenplain X_t, (A_t(i))_{i\in D}} - \E{X_{t+1}'(D, s) \givenplain X_t} }  \\
\nonumber
&\leq \sum_{s'\in\sspa, a\in\aspa} \abs{\pibs(a|s) X_t(D, s')  -  \Gamma_t(s',a)} P(s', a, s) \\
\label{eq:hw:interm-7}
&\leq \frac{1}{N} \sum_{s'\in\sspa, a\in\aspa} P(s', a, s) \quad a.s.
\end{align}

Summing over all $s\in\sspa$ in (<ref>), we get

\begin{align}
\label{eq:hw:interm-8}
\normbig{\E{X_{t+1}'(D, s) \givenplain X_t, (A_t(i))_{i\in D}} - \E{X_{t+1}'(D, s) \givenplain X_t}}_1
&\leq \frac{2|\sspa|}{N} \quad a.s.
\end{align}

Combining the above calculations, we get

\begin{align}
\nonumber
&\mspace{20mu} \Ebig{ \norm{X_{t+1}'(D)- X_t(D) P_{\pibs}}_1 \givenbig X_t} \\
\nonumber
&\leq   \Ebig{ \norm{X_{t+1}'(D) - \E{X_{t+1}'(D, s) \givenplain X_t, (A_t(i))_{i\in D}}}_1 \givenbig X_t}  \\
\nonumber
&\mspace{20mu} +  \Ebig{\normbig{\E{X_{t+1}'(D, s) \givenplain X_t, (A_t(i))_{i\in D}} - \E{X_{t+1}'(D, s) \givenplain X_t}}_1 \givenbig X_t}  \\
&\leq \frac{1}{\sqrt{N}} + \frac{2|\sspa|}{N},
\end{align}

which implies (<ref>).
We also have that for any $r\geq 2|\sspa| / N$,

\begin{align}
\nonumber
&\mspace{20mu} \Prob{\norm{X_{t+1}'(D)- X_t(D) P_{\pibs}}_1 > r \givenbig X_t} \\
\nonumber
&\leq  \Probbig{ \norm{X_{t+1}'(D) - \E{X_{t+1}'(D, s) \givenplain X_t, (A_t(i))_{i\in D}}}_1 > r - \frac{2|\sspa|}{N} \givenbig X_t} \\
\nonumber
&\leq 2|\sspa| \exp\Big(-\frac{2N (r-2|\sspa|/N)^2}{|\sspa|^2}\Big) \\
\label{eq:hw:interm-9}
&\leq 2|\sspa| \exp\Big(-\frac{N r^2}{|\sspa|^2} + \frac{8}{N}\Big)
\end{align}

where we use the fact that $(r - 2|\sspa|/N)^2 \geq r^2 / 2 - (2|\sspa|/N)^2$ in the (<ref>).
Note that the above upper bound also holds for $r< 2|\sspa| / N$ because the term in (<ref>) is greater than $1$ for $r < 2|\sspa| / N$. Because $\exp(8/N) = O(1)$, we have proved (<ref>):

\[
\Prob{\norm{X_{t+1}'(D)- X_t(D) P_{\pibs}}_1 > r \givenbig X_t} = O(\exp(-C\lamw Nr^2)),
\]

where $C = 1/ (\lamw |\sspa|^2)$. This finishes the proof that

\begin{align}
\tag{\ref{eq:hw:drift}}
\Ebig{(\hw(X_{t+1}', D) - \rhow \hw(X_t, D))^+ \givenbig X_t} &\leq O(1/\sqrt{N}) \quad a.s. \\
\tag{\ref{eq:hw:high-prob}}
\Probbig{\hw(X_{t+1}', D) > \rhow \hw(X_t, D) + r   \givenbig X_t} &\leq O(\exp(-C N r^2)) \quad a.s. \quad \forall r \geq 0.
\end{align}

Proving the analogs of (<ref>) and (<ref>) for $\hu$.
Now let $X_{t+1}'$ denote the system at time $t+1$ given that the arms in $D$ follow the $\lppriority$ at time $t$.
The proof is almost verbatim: We first show

\begin{align}
\label{eq:hu:interm-1}
\EBig{\normBig{ X_{t+1}'(D, s) - \Eplain{X_{t+1}'(D, s) \givenplain X_t, (A_t(i))_{i\in D}}}_1 \givenBig X_t} &\leq \frac{|\sspa|}{\sqrt{N}} \\
\label{eq:hu:interm-2}
\ProbBig{\normBig{ X_{t+1}'(D, s) - \Eplain{X_{t+1}'(D, s) \givenplain X_t, (A_t(i))_{i\in D}}}_1 > r \givenBig X_t} &\leq 2|\sspa|\exp\Big(-\frac{2N r^2}{|\sspa|^2}\Big),
\end{align}

using identical arguments. Then we can show that

\begin{align}
\label{eq:hu:interm-3}
\normbig{\E{X_{t+1}'(D, s) \givenplain X_t, (A_t(i))_{i\in D}} - \E{X_{t+1}'(D, s) \givenplain X_t}}_1
&\leq \frac{4|\sspa|}{N} \quad a.s.,
\end{align}

using the fact that under Optimal Local Control (<Ref>), for all $s'\in\sspa$ and $a\in\aspa$,

\begin{align*}
\abs{\pibs(a|s) X_t(D, s')  -  \Gamma_t(s',a)} &\leq \frac{1}{N} \quad a.s. \quad \forall s'\neq \sneu \\
\abs{\pibs(a|\sneu) X_t(D, s')  -  \Gamma_t(\sneu,a)} &\leq \frac{|\sspa|}{N} \quad a.s.
\end{align*}

Finally, we combine (<ref>),  (<ref>) and (<ref>) to get that under $\hu$,

\begin{align}
\label{eq:hu:drift}
\Ebig{(\hw(X_{t+1}', D) - \rhow \hw(X_t, D))^+ \givenbig X_t} &\leq O(1/\sqrt{N}) \quad a.s. \\
\label{eq:hu:high-prob}
\Probbig{\hw(X_{t+1}', D) > \rhow \hw(X_t, D) + r   \givenbig X_t} &\leq O(\exp(-C_U N r^2)) \quad a.s. \quad \forall r \geq 0,
\end{align}

for some constant $C_U > 0$.

The rest of the inequalities in this lemma.
The other two inequalities, (<ref>) and (<ref>), and their analogs for $\hu$ are relatively straightforward from the definitions of the weighted $L_2$ norms, so we skip them here. Similar calculations have also been done in [10] for $\hw$, so we refer the readers to [10] for the detailed calculations.

§ EXAMPLES AND PROOFS OF VIOLATING ASSUMPTIONS

§.§ Examples of violating Assumptions <ref>-<ref>

An aperiodic but non-unichain example.
We first give an example of RBs where $P_\pibs$ is aperiodic but not unichain.
Consider the RB problem whose single-armed MDP has two states, $A$ and $B$, with no transitions between these two states.
The reward function is given by $r(A, 1) = r(B, 0) = 1$ and $r(A, 0) = r(B, 1) = 0$.
Let $\alpha = 1/2$ in the budget constraint.
Since there is no transition between these two states, the Markov chain induced by the optimal single-armed policy, $\pibs$, $P_\pibs$, is not unichain. However, each of the recurrent classes $\{A\}$ and $\{B\}$ is aperiodic in this Markov chain.

It is easy to verify that the optimal solution to the LP relaxation (<ref>) is  $y^*(A,1) = y^*(B, 0) = 1/2$ and $y^*(s,a) = 0$ for the rest of $s\in\sspa, a\in\aspa$; the optimal value of  (<ref>) is $\rrel = 1$.
Now consider the $N$-armed problem where all arms are initialized in state $A$. In this case, the long-run average reward of any policy is $1/2$. Therefore, $\rrel - \rsysn \geq 1/2 = \Omega(1/\sqrt{N})$ for any $\pi$.

A unichain but periodic example.
Next, we give an example of RBs where  $P_\pibs$ is unichain but periodic, and $\rrel - \rsysn = \Theta(1)$. This example was included in [10], and we paraphrase it here for completeness.

Consider the RB problem whose single-armed MDP has two states, $A$ and $B$.
The state of the single-armed MDP deterministically transitions to the other state in each time step, regardless of the action applied.
Let the reward function be $r(A,0)=r(B,1)=1$ and $r(A,1)=r(B,0)=0$; let $\alpha = 1/2$ in the budget constraint.
Obviously, the Markov chain induced by any single-armed policy is unichain and periodic with a period $2$.

One can easily see the optimal solution to the LP relaxation (<ref>) is by $y^*(A, 0) = y^*(B, 1) = 1/2$ and $y^*(A, 1) = y^*(B, 0) = 0$; the optimal value of (<ref>) is $\rrel = 1$.
Now consider the $N$-armed problem where all arms are initialized in state $A$. Then at any time $t$, either all arms are in state $A$ or all arms are in state $B$. In this case, under any policy $\pi$, the long-run average $\rsysn$ is $1/2$, so $\rrel - \rsysn= 1/2 = \Omega(1/\sqrt{N})$.

An example satisfying <Ref> but not <Ref>.

Next, we give an example of RB satisfying <Ref> but violating <Ref>.
Suppose each arm in the $N$-armed problem has state $S_t(i) = \text{Ber}(1/2)$ i.i.d. across $i$ and $t$. Let $\alpha = 1/2$ in the budget constraint. The reward function $r(s,1) = s$ and $r(s,0) = 0$ for $s = 0, 1$.
Intuitively, the decision maker observes the state of each arm, and pulls $N/2$ of them. Each arm pulled by the decision maker generates an amount of reward equal to its state.
The single-armed MDP is clearly aperiodic unichain under any policy because its state is i.i.d.

The optimal solution of the LP relaxation (<ref>) is $y^*(1,1) = y^*(0,0) = 1/2$, and $y^*(s,a) = 0$ for other $s\in\sspa, a\in\aspa$; the optimal value $\rrel = 1/2$. Intuitively, the optimal single-armed policy pulls the arm if and only if it is in state $1$.

In the $N$-armed problem, the optimal policy is obviously pulling $N/2$ arms with the largest states. Then the optimal reward $\ropt$ is given by

\[
\ropt = \frac{1}{N}\EBig{\min\Big(\sumN S_i(t), \frac{N}{2}\Big)},
\]

for an arbitrary $t$.
Since $S_t(i) = \text{Ber}(1/2)$ and are i.i.d. across $i$, it can be shown that $\ropt = 1/2 - \Theta(1/\sqrt{N})$. Therefore, $\rrel - \ropt = \Theta(1/\sqrt{N})$.

§.§ Proofs of Lemma <ref> and Theorem <ref>

In this section, we prove our second main result, <Ref>, which is a lower bound on $\rrel - \rsysn$ for regular unstable RBs (<Ref>); meanwhile we will also prove <Ref> stated in <Ref> as one of the preliminary steps for proving <Ref>.
This section is organized linearly, with a sequence of lemma statements and proofs; <Ref> will be proved at the end of this section.

We first restate and prove <Ref>.

We first investigate the dual of (<ref>), whose properties help us prove the lemma:

\begin{align}
\tag{Dual-LP}\label{eq:lp-dual} \underset{\lambda, \mu, \{f(s)\}_{s\in\sspa}}{\text{minimize}} &  \mu \\
\text{subject to}\mspace{12mu}
\label{eq:dual-lp-constraint}
&  r(s,a) +  (\alpha - \indibrac{a=1}) \lambda - \mu  + \sum_{s'} P(s,a,s')f(s') \leq f(s)  \quad \forall s\in\sspa, a\in\aspa \\
\nonumber
& \lambda, \mu, f(s) \in \R \quad \forall s\in\sspa
\end{align}

Note that the primal problem (<ref>) is always feasible and bounded, so (<ref>) is also feasible and bounded.
By strong duality of linear programs, the optimal value of the dual problem is $\mu^* = \rrel$.
Moreover, there always exists an optimal solution of the dual problem $(\lambda^*, \vlam, \mu^*)$ that satisfies strict complementarity with the unique primal optimal solution $y^*$ [7]:

\begin{equation}
\label{eq:c-s}
r(s,a) - \lambda^* \indibrac{a=1} + \lambda^* \alpha - \mu  + \sum_{s'} P(s,a,s')\vlam(s') - \vlam(s) = 0 \quad \text{ iff } \quad y^*(s,a) = 0.
\end{equation}

Let $\rlam(s,a) = r(s,a) -  \lambda^* \indibrac{a=1} + \lambda^* \alpha$.
By (<ref>) and the definition of $S^+$ and $S^-$,
there exists a constant $\rewardgap > 0$ such that

\begin{align}
\label{eq:lb:scs-1}
\rlam(s, 0) - \rrel + \sum_{s'} P(s,0,s')\vlam(s') &\leq \vlam(s) - \rewardgap   \quad \forall s\in S^+ \\
\label{eq:lb:scs-2}
\rlam(s, 1) - \rrel + \sum_{s'} P(s,1,s')\vlam(s') &\leq \vlam(s) - \rewardgap \quad \forall s\in S^-.
\end{align}

Define the function $\qlam\colon \sspa\times \aspa \to \R$ as

\[
\qlam(s,a) = \rlam(s,a) - \rrel + \sum_{s'} P(s,a,s')\vlam(s').
\]

Then (<ref>) implies that $\vlam(s) \geq \qlam(s,a)$ for any $s\in\sspa,a\in\aspa$;
(<ref>) and (<ref>) imply that $\qlam(s,1) - \qlam(s,0)$ is at least $\rewardgap$ for $s\in S^+$ and at most $-\rewardgap$ for $s\in S^-$.

Now, we are ready to bound $\rrel - \rsysn$: for any policy $\pi$ and initial state vector $\veS_0$,

\begin{align}
\nonumber
&\rrel - \rsysn\\
\nonumber
&=  \EBig{\rrel - \frac{1}{N} \sumN \rlam(S_\infty(i), A_\infty(i))} \\
\nonumber
&= \frac{1}{N} \sumN \EBig{\rrel - \rlam(S_\infty(i), A_\infty(i))} \\
\label{eq:lb:apply-rate-conservation}
&=  \frac{1}{N}\sumN \EBig{\rrel -  \rlam(S_\infty(i), A_\infty(i)) - \sum_{s'\in\sspa}P(S_\infty(i),A_\infty(i),s')\vlam(s') + \vlam(S_\infty(i))} \\
\nonumber
&= \frac{1}{N}\sumN \EBig{-\qlam(S_\infty(i), A_\infty(i)) + \vlam(S_\infty(i))}\\
\label{eq:lb:apply-action-gap}
&\geq  \frac{1}{N}\sumN \rewardgap \EBig{\indibrac{S_\infty(i)\in S^+, A_\infty(i) = 0} + \indibrac{S_\infty(i)\in S^-, A_\infty(i) = 1}} \\
&= \rewardgap \EBig{\sum_{s\in S^+} Y_\infty(s,0) + \sum_{s\in S^-} Y_\infty(s,1)}.
\end{align}

where to get (<ref>), we use the fact that  $\Eplain{\sum_{s'\in\sspa}P(S_\infty(i),A_\infty(i),s')f(s')}=\Eplain{f(S_\infty(i))}$ for any function $f\colon\sspa\to\R$.

Next, we show <Ref>, which implies that under any actions for the $N$-arms, one can find a subset of arms (can be empty) such that these almost follow the Optimal Local Control; the size of the subset is determined by $\disy(Y_t)$.
Essentially, this lemma shows that under any policy, if $\disy(Y_t)$ is small, the transition dynamics will be close to the transition dynamics under Optimal Local Control.

Consider an instance of regular unstable RBs (<Ref>). Under any policy $\pi$, at any time $t$, there exists a set $D_t\subseteq[N]$ depending on $X_t$ and $\veA_t$ with size $m(D_t) \geq 1 - 2\disy(Y_t)/\beta - 1/(\beta N)$, where $\beta = \min(\alpha,1-\alpha)$,
such that if we resample actions for the arms in $D_t$ using Optimal Local Control, we get almost the same actions as $A_t$ up to at most $1$ arm.
For such $D_t$, we have

* The expectation of $X_{t+1}(D_t)$ satisfies

\begin{equation}
\label{eq:subset-lp-priority}
\normbig{\Ebig{X_{t+1}(D_t) - m(D_t) \statdist \givenbig X_t, \veA_t} - (X_t(D_t)- m(D_t)\statdist)\Phi}_1 \leq \frac{2}{N}.
\end{equation}

* Let $X_{t+1}'(D_t)$ be the state distribution at time $t+1$ when the actions for the arms in $D_t$ are resampled using Optimal Local Control, then there exists a coupling between $X_{t+1}'$ and $X_{t+1}$ such that

\begin{equation}
\label{eq:subset-lp-priority-as}
\norm{X_{t+1}'(D_t) - X_{t+1}(D_t)}_1 \leq \frac{2}{N} \quad a.s.
\end{equation}

We construct the subset $D_t$ by removing a suitable subset of arms from $[N]$.
We first remove each arm $i$ with state $S_t(i)\in S^+$ and $A_t(i) = 0$ or $S_t(i) \in S^-$ and $A_t(i) = 1$. After this step, we get a subset $D'$ with $N(1-\disy(Y_t))$ arms. The budget usage of $D'$ equals $\alpha N - N \sum_{s\in S^-} Y_t(s,1)$.
Then $D_t$ is obtained by removing arms that take action $1$ or action $0$ from $D'$ such that the budget usage of $D_t$ equals $\floor{\alpha |D_t|}$.
If $\alpha N - N \sum_{s\in S^-} Y_t(s,1)  > \alpha |D'|$, it is not hard to show that the number of arms removed from $D'$ is bounded by

\[
|D'| - |D_t|  \leq \frac{1}{1-\alpha} \Big(\absbig{\alpha N - N \sum_{s\in S^-} Y_t(s,1) - \alpha |D'| }  + 1\Big) \leq  \frac{1}{1-\alpha} \big(N\disy(Y_t) + 1\big);
\]

If $\alpha N - N \sum_{s\in S^-} Y_t(s,1)  > \alpha |D'|$, it is not hard to show that the number of arms removed from $D'$ is bounded by

\[
|D'| - |D_t| \leq \frac{1}{\alpha} \Big(\absbig{\alpha N - N \sum_{s\in S^-} Y_t(s,1) - \alpha |D'|}  + 1\Big) \leq  \frac{1}{\alpha} \big(N\disy(Y_t) + 1\big).
\]

In both cases, we have

\begin{equation}
|D'| - |D_t|  \leq \frac{1}{\beta}\big(N\disy(Y_t) + 1\big),
\end{equation}

where recall that $\beta = \min(\alpha, 1-\alpha)$.

Now we get a subset $D_t$ such that $m(D_t) \geq 1 - \disy(Y_t)- (\disy(Y_t)+ 1/N) / \beta \geq 1-2\disy(Y_t)/\beta - 1/(\beta N)$. Moreover, for each $i\in D_t$, $A_t(i) = 1$ if $S_t(i) \in S^+$; $A_t(i) = 0$ if $S_t(i) \in S^-$; the total budget usage of $D_t$ equals $\floor{\alpha |D_t|}$.

Following similar calculations as the proof of <Ref>, we can show that

\[
\E{X_{t+1}(D_t) \givenplain X_t, \veA_t}
= X_t(D_t) \Big(P_\pibs - \costvec^\top \big(P_1(\sneu) - P_0(\sneu)\big) \Big) + m(D_t) \frac{\floor{\alpha N}}{N} (P_1(\sneu) - P_0(\sneu)).
\]

whereas by the definition of $\Phi$,

\[
(X_t(D_t) - m(D_t)\statdist) \Phi + m(D_t) \statdist = X_t(D_t) \Big(P_\pibs - \costvec^\top \big(P_1(\sneu) - P_0(\sneu)\big) \Big) + m(D_t) \alpha (P_1(\sneu) - P_0(\sneu)).
\]

Comparing the above two equations, we get (<ref>).

To show (<ref>), note that if we resample the actions using the Optimal Local Control for the arms in $D_t$,
the actions of the arms in $S^+$ and $S^-$ are unchanged; at most one arm with state $\sneu$ gets a different action after resampling.
Therefore, we can couple $X_{t+1}'$ and $X_{t+1}$ such that for each arm that gets the same actions after resampling, the next states of the arm are the same. In this way, we get $\normplain{X_{t+1}'(D_t) - X_{t+1}(D_t)}_1 \leq 2/N$.

Next, we conduct Lyapunov analysis with the help of <Ref>.  To construct Lyapunov functions, we first define two seminorms.
Consider the Jordan decomposition of $\Phi$, $\Phi = V J V^{-1}$. We split $J$ into $J = \Jus + \Jst$, where $\Jus$ contains the Jordan blocks of $J$ with eigenvalues' moduli greater than $1$, and $\Jst$ contains the Jordan blocks of $J$ with eigenvalues' moduli smaller than $1$.
Let $\Ius$ be the diagonal matrix with $1$ on the entries where $\Jus$ is non-zero; let $\Ist$ be the identity matrix subtracting $\Ius$.
We define the matrices $\Uus$ and $\Ust$ via the infinite series:

\begin{align}
\label{eq:Uus-def}
\Uus &= \sum_{k=1}^\infty \big(V \Jus^{\dagger} V^{-1}\big)^k \big(V \Jus^{\dagger} V^{-1}\big)^{\top k}  \\
\label{eq:Ust-def}
\Ust &= V \Ist V^{-1}  \big(V \Ist V^{-1}\big)^{\top} + \sum_{k=1}^\infty \big(V \Jst V^{-1} \big)^k \big(V \Jst V^{-1}\big)^{\top k}
\end{align}

where the superscript $\dagger$ denotes pseudo-inverse.

In the next lemma, we show that $\Uus$ and $\Ust$ are well-defined, and prove several useful facts about the semi-norms weighted by $\Uus$ and $\Ust$.

For any regular unstable RBs (<Ref>),
the matrices $\Uus$ and $\Ust$ given in (<ref>) and (<ref>) are well-defined and positive-semi-definite.
Further, if we let $\norm{v}_{\Uus} = \sqrt{v \Uus v^\top}$ and $\norm{v}_{\Ust} = \sqrt{v \Ust v^\top}$ for any $v\in \R^{|\sspa|}$, then

\begin{align}
\label{eq:Uus-expansion}
\norm{v \Phi}_{\Uus} &\geq (1 + \gapus) \norm{v}_{\Uus} \\
\label{eq:Uus-expansion-max}
\norm{v \Phi}_{\Uus} &\leq (1 + \gapusmax) \norm{v }_{\Uus} \\
\label{eq:Ust-contraction}
\norm{v \Phi}_{\Ust} &\leq (1-\gapst) \norm{v }_{\Ust},
\end{align}

where $\gapus, \gapusmax, \gapst$ are positive constants.
Moreover, there exists a positive constant $\Knorm$ such that for any $v\in \R^{|\sspa|}$,

\begin{equation}
\label{eq:Uus-plus-Ust-strength}
\norm{v}_{\Uus} +  \norm{v}_{\Uus}  \geq \Knorm \norm{v}_1.
\end{equation}

To show the well-definedness $\Uus$ and $\Ust$, note that the moduli of all eigenvalues of $V \Jus^{\dagger} V^{-1}$ and $V \Jst V^{-1}$ are strictly less than $1$, so their $k$-th powers must converge exponentially as $k\to\infty$, implying the convergence of the series in (<ref>) and (<ref>). Moreover, since each term in the two series is positive semi-definite, both series converge to positive semi-definite matrices.

Next, we show (<ref>). Observe that

\begin{align*}
\Phi \Uus \Phi^\top
&= \sum_{k=1}^\infty \big(V J \Jus^{\dagger k} V^{-1}\big)\big(V J \Jus^{\dagger k} V^{-1}\big)^\top \\
&= \big(V J \Jus^{\dagger} V^{-1}\big)\big(V J \Jus^{\dagger} V^{-1}\big)^\top +
\sum_{k=2}^\infty \big(V J \Jus^{\dagger k} V^{-1}\big)\big(V J \Jus^{\dagger k} V^{-1}\big)^\top \\
&= \big(V \Ius V^{-1}\big)\big(V \Ius V^{-1}\big)^\top +
\sum_{k=2}^\infty \big(V \Jus^{\dagger (k-1)} V^{-1}\big)\big(V \Jus^{\dagger (k-1)} V^{-1}\big)^\top  \\
&= \big(V \Ius V^{-1}\big)\big(V \Ius V^{-1}\big)^\top + \Uus.
\end{align*}

We claim that there exists a positive constant $\gapus'$ such that for any $v \in \R^{|\sspa|}$,

\begin{equation}
\label{eq:ius-bound-uus}
v\big(V \Ius V^{-1}\big)\big(V \Ius V^{-1}\big)^\top v^\top \geq \gapus' v \Uus v^\top,
\end{equation}

Without loss of generality, we only need to show (<ref>) for $v$ not being the zero vector. By normalization, it suffices to show that the function $v\mapsto v \Uus v^\top$ is bounded on the set  $\{v\colon \norm{v\big(V \Ius V^{-1}\big)}_2 = 1\}$.
Let $\spus$ be the linear subspace spanned by rows of $V$ corresponding to the non-zero entries on $\Ius$'s diagonal, and let $\spst$ be the linear subspace spanned by rows of $V$ corresponding to the zero entries on $\Ius$'s diagonal.
Observe that $\R^{|\sspa|}$ can be represented as the direct sum $\spus \oplus \spst$, which implies that any $v\in\R^{|\sspa|}$ can be decomposed as $v = v^\parallel + v^\perp$, where $v^\parallel \in \spus$ and $v^\perp \in \spst$.
Moreover, $\spst$ is the kernel of $V \Ius V^{-1}$, and $V \Ius V^{-1}$ is the identity mapping when restricted to $\spus$. Therefore, if $\normplain{v\big(V \Ius V^{-1}\big)}_2 = 1$, then $\normplain{v^\parallel}_2 = 1$.
Observe that $\spst$ is also the kernel of $\Ust$, so

\[
v \Uus v^\top =  (v^\parallel + v^\perp) \Uus (v^\parallel+v^\perp)^\top = v^\parallel \Uus (v^\parallel)^\top.
\]

Because $\{v\in\spus \colon \norm{v}_2 = 1\}$ is a compact set, $v \Uus v^\top$ is bounded when $\normplain{v\big(V \Ius V^{-1}\big)}_2 = 1$.
After getting (<ref>), we have that

\begin{align*}
\norm{v \Phi}_{\Uus}
&= \sqrt{v \Phi \Uus \Phi^\top v^\top} \\
&\geq \sqrt{(1 + \gapus') v \Uus v^\top} \\
&= \sqrt{1+\gapus'} \norm{v}_{\Uus}.
\end{align*}

Letting $\gapus = \sqrt{1+\gapus'} - 1$, we obtain (<ref>).

Next, we prove (<ref>).
For any $v\in \R^{|\sspa|}$, we have $v \Uus v^\top \geq \normbig{v V \Jus^\dagger V^{-1}}_2^2$.
We claim that there exists a positive constant $\gapusmax'$ such that

\begin{equation}
\label{eq:ius-bound-uus-reverse}
\norm{v V \Ius V^{-1}}_2 \leq \gapusmax' \normbig{v V \Jus^\dagger V^{-1}}_2.
\end{equation}

To show (<ref>),
without loss of generality, it suffices to show that
it suffices to show that the function $v\mapsto \normplain{v V \Ius V^{-1}}_2$ is bounded on the set $\{v\colon \normplain{v V \Jus^\dagger V^{-1}}_2 = 1\}$.
Consider the decomposition $\R^{|\sspa|} = \spus \oplus \spst$.
Because the kernel of $V \Jus^\dagger V^{-1}$ is $\spst$, and $V \Jus^\dagger V^{-1}$ is invertible when restricted to $\spus$, the set $\{v\colon \normplain{v V \Jus^\dagger V^{-1}}_2 = 1\}$ is the sum of $\spst$ and a compact subset of $\spus$. Because $\spst$ is also the kernel of $v V \Ius V^{-1}$, $\norm{v V \Ius V^{-1}}_2$ is bounded. This proves (<ref>). Because $\Phi \Uus \Phi^\top = \big(V \Ius V^{-1}\big)\big(V \Ius V^{-1}\big)^\top + \Uus$, (<ref>) implies that

\begin{align*}
\norm{v \Phi}_{\Uus}
&\leq \sqrt{v\Phi \Uus \Phi v^\top} \\
&= \sqrt{v\big(V \Ius V^{-1}\big)\big(V \Ius V^{-1}\big)^\top v^\top + v\Uus v^\top} \\
&\leq \sqrt{(\gapusmax')^2 +1 } \norm{v}_{\Uus}.
\end{align*}

Taking $\gapusmax = \sqrt{(\gapusmax')^2 +1 } - 1$ finishes the proof.

Next, we show (<ref>). Observe that

\begin{align*}
\Phi \Ust \Phi^\top
&= V J \Ist V^{-1}  \big(V J \Ist V^{-1}\big)^{\top} + \sum_{k=1}^\infty \big(V J \Jst^k V^{-1} \big) \big(V J \Jst^k V^{-1}\big)^{\top} \\
&= \sum_{k=1}^\infty \big(V \Jst^k V^{-1} \big) \big(V \Jst^k V^{-1}\big)^{\top} \\
&= \Ust -  V \Ist V^{-1}  \big(V \Ist V^{-1}\big)^{\top}.
\end{align*}

We claim that there exists a positive constant $\gapst'$ such that for any $v\in \R^{|\sspa|}$,

\begin{equation}
\label{eq:ist-bound-ust}
v\big(V \Ist V^{-1}\big)\big(V \Ist V^{-1}\big)^\top v^\top \geq \gapst' v \Ust v^\top.
\end{equation}

The argument for (<ref>) is similar to the argument for (<ref>), so we omit it here.

\begin{align*}
\norm{v\Phi}_{\Ust} &= \sqrt{v\Phi \Ust \Phi^\top v^\top} \\
&\leq \sqrt{(1-\gapst')v\Ust v^\top} \\
&= \sqrt{1-\gapst'} \norm{v}_{\Ust}.
\end{align*}

Taking $\gapst = 1 - \sqrt{1-\gapst'}$, we obtain (<ref>).

To show (<ref>), note that

\begin{align*}
\norm{v}_{\Uus} +  \norm{v}_{\Uus}
&= \norm{v \Uus v^\top} + \norm{v \Ust v^\top}  \\
&\geq \norm{v (\Uus+\Ust) v^\top} \\
&= \norm{v}_{\Uus+\Ust}.
\end{align*}

Since all norms in a finite-dimensional Euclidean space are equivalent, it suffices to show $\norm{\cdot}_{\Uus+\Ust}$ is a norm, or equivalently, $\Uus+\Ust$ is positive definite. For any $v\in \R^{|\sspa|}$, we have

\begin{align*}
v (\Uus+\Ust) v^\top
&\geq v \Big(\big(V \Jus^{\dagger} V^{-1}\big) \big(V \Jus^{\dagger} V^{-1}\big)^{\top} + V \Ist V^{-1}  \big(V \Ist V^{-1}\big)^{\top} \Big) v^\top \\
&= vV \big(\Jus^\dagger V^{-1} V^{-\top} \Jus^\dagger + \Ist V^{-1} V^{-\top} \Ist  \big)   V^\top v^\top.
\end{align*}

By the definitions of $\Jus$ and $\Ist$, it is not hard to see that we always have $v (\Uus+\Ust) v^\top \geq 0$, and $v$ is the zero vector whenever $v (\Uus+\Ust) v^\top = 0$. Therefore, $\Uus+\Ust$ is positive definite, which implies (<ref>).

Now we are ready to conduct the Lyapunov analysis.
We define the Lyapunov functions $\hst$ and $\hus$: for any system state $x$ and $D\subseteq[N]$, let

\begin{align}
\hst(x, D) &=  \norm{x(D) - m(D)\statdist}_{\Ust}  \\
\hus(x, D) &=  \norm{x(D) - m(D)\statdist}_{\Uus}.
\end{align}

Our next lemma gives a lower bound for the drift of $\hst(x, D)$ under any policy.

Consider an instance of regular unstable RBs (<Ref>). Let $K_1$ be any positive constant. Under any policy $\pi$, for any time step $t$, we have

\begin{align}
\nonumber
&\mspace{23mu} \Ebig{\hus(X_{t+1}, [N]) \givenbig \givenbig X_t, \veA_t} - \hus(X_t, [N])\\
\frac{\gapus K_1- (2+\gapus + \gapusmax) K_1 \indibrac{\goodevent_t^c}}{\sqrt{N}} - \frac{6\lamus^{1/2}}{N} -  8 \lamus^{1/2} \disy(Y_t) \quad a.s.,
\end{align}

where the good event $\goodevent_t$ is given by

\[
\goodevent_t = \Big\{\hus(X_t, [N]) > \frac{K_1}{\sqrt{N}} + \frac{2\lamus^{1/2} \disy(Y_t)}{\beta} + \frac{\lamus^{1/2}}{\beta N}, \, \norm{X_t([N]) - \statdist}_1 \leq \radiusus -  \frac{6\disy(Y_t)}{\beta} - \frac{3}{\beta N}\Big\},
\]

and $K_1 = \lamus^{1/2} \noiseortho / (2+\gapus+\gapusmax)$.

Let $D_t$ be the set given in <Ref>. By the definition of $\hus$ and the fact that $m(D_t) \geq 1-2\disy(Y_t) - 1/N$, we have

\begin{align}
\nonumber
&\mspace{23mu} \Ebig{\hus(X_{t+1}, [N]) \givenbig X_t, \veA_t} - \hus(X_t, [N])  \\
\label{eq:pf-hus-drift:interm-0}
&\geq  \Ebig{\hus(X_{t+1}, D_t)  \givenbig X_t, \veA_t} - \hus(X_t, D_t) - 4\lamus^{1/2} (1-m(D_t))  \\
\nonumber
&\geq  \Ebig{\hus(X_{t+1}, D_t)  \givenbig X_t, \veA_t} - \hus(X_t, D_t) - 8 \lamus^{1/2} \disy(Y_t) / \beta - 4\lamus^{1/2} /(\beta N).
\end{align}

Note that (<ref>) uses the fact that $\norm{v}_{\Uus} \leq \lamus^{1/2} \norm{v}_{1}$.

We first consider the case with $\hus(X_t, D_t) > K_1 / \sqrt{N}$. We have

\begin{align}
\nonumber
\Ebig{\hus(X_{t+1}, D_t) \givenbig \givenbig X_t, \veA_t}
&= \Ebig{\norm{X_{t+1}(D_t) - m(D_t)\statdist}_{\Uus} \givenbig X_t, \veA_t} \\
\label{eq:pf-hus-drift:interm-1}
&\geq \normbig{\Ebig{X_{t+1}(D_t) - m(D_t)\statdist \givenbig X_t, \veA_t}}_{\Uus} \\
\label{eq:pf-hus-drift:interm-2}
&\geq  \normbig{\big(X_{t}(D_t) - m(D_t)\statdist\big)\Phi}_{\Uus}  - 2\lamus^{1/2}/(\beta N),
\end{align}

where (<ref>) uses Jensen's inequality;
(<ref>) invokes <Ref> and the fact that $\norm{v}_{\Uus} \leq \lamus^{1/2} \norm{v}_{1}$ for any $v\in\R^{|\sspa|}$.
By <Ref>, we have

\begin{equation*}
\normplain{v\Phi}_{\Uus} \geq (1 + \gapus) \normplain{v}_{\Uus}.
\end{equation*}

Letting $v = X_{t}(D_t) - m(D_t)\statdist$, combined with (<ref>), we get

\begin{align*}
\Ebig{\hus(X_{t+1}, D_t) \givenbig \givenbig X_t, \veA_t}
(1+\gapus) \hus(X_t, D_t) - \frac{2\lamus^{1/2}}{\beta N} \\
&\geq \hus(X_t, D_t) + \frac{\gapus K_1}{\sqrt{N}} - \frac{2\lamus^{1/2}}{\beta N},
\end{align*}

which implies that

\begin{equation}
\label{eq:pf-hus-drift:goal-1}
\Ebig{\hus(X_{t+1}, [N]) \givenbig \givenbig X_t, \veA_t} - \hus(X_t, [N]) \geq \frac{\gapus K_1}{\sqrt{N}} - \frac{6\lamus^{1/2}}{\beta N} -  \frac{8 \lamus^{1/2} \disy(Y_t)}{\beta}.
\end{equation}

Next, we consider the case of $\hus(X_t, D_t) \leq K_1 / \sqrt{N}$:

\begin{align}
\nonumber
&\mspace{23mu} \Ebig{\hus(X_{t+1}, D_t) \givenbig \givenbig X_t, \veA_t}  \\
\nonumber
&=  \Ebig{\norm{X_{t+1}(D_t) - m(D_t) \statdist}_{\Uus}  \givenbig X_t, \veA_t} \\
\nonumber
&\geq \Ebig{\norm{X_{t+1}(D_t) - m(D_t) \statdist - (X_t(D_t)-m(D_t) \statdist)\Phi}_{\Uus}  \givenbig X_t, \veA_t} - \normbig{\big(X_t(D_t) - m(D_t)\statdist\big)\Phi}_{\Uus} \\
\label{eq:pf-hus-drift:interm-3}
&\geq \Ebig{\norm{X_{t+1}(D_t) - m(D_t) \statdist - (X_t(D_t) - m(D_t)\statdist)\Phi}_{\Uus}  \givenbig X_t, \veA_t} - (1+\gapusmax) \hus(X_t, D_t),
\end{align}

where (<ref>) uses the fact that $\normplain{v}_{\Phi\Uus \Phi^\top} \leq (1+\gapusmax) \normplain{v}_{\Uus}$ for all $v\in \R^{|\sspa|}$, proved in <Ref>.
We plan to invoke Property <ref> of regular unstable RBs to bound the first term in (<ref>).
Let $X_{t+1}'(D_t)$ be the state distribution at time $t+1$ when the arms in $D_t$ follow $\lppriority$. By <Ref>, there exists a coupling such that $\norm{X_{t+1}'(D_t) - X_{t+1}(D_t)}_1 \leq 2/N$ almost surely, so

\begin{align}
\nonumber
&\mspace{23mu} \Ebig{\norm{X_{t+1}(D_t) - m(D_t) \statdist - (X_t(D_t) - m(D_t)\statdist)\Phi}_{\Uus}  \givenbig X_t, \veA_t} \\
\label{eq:pf-hus-drift:interm-4}
&\geq \Ebig{\norm{X_{t+1}'(D_t) - m(D_t) \statdist - (X_t(D_t) - m(D_t)\statdist)\Phi}_{\Uus}  \givenbig X_t, \veA_t} - \frac{2\lamus^{1/2}}{N}.
\end{align}

Let $\xi^\top$ be the right eigenvector of $\Uus$ corresponding to the eigenvalue $\lamus$, with norm $\norm{\xi}_2 = 1$.
Because $\Uus$ is positive semidefinite, $\Uus \geq \lamus\xi^\top \xi$, so

\[
\norm{v}_{\Uus} \geq \sqrt{\lamus v\xi^\top \xi v^\top} = \lamus^{1/2} \absplain{v\xi^\top},
\]

for any $v\in \R^{|\sspa|}$.
Moreover, we argue that $\xi$ satisfies the conditions in Property <ref> of regular unstable RBs: the column space of $\Uus$ is a subspace of the column space of $\Phi$. Because $\statdist \Phi = 0$, we have $\statdist \Uus = 0$. Therefore, $\statdist \xi^\top = \lamus^{-1} \statdist \Uus \xi^\top = 0$.
Therefore, by Property <ref> of regular unstable RBs,

\begin{align}
\nonumber
&\mspace{23mu} \Ebig{\norm{X_{t+1}'(D_t) - m(D_t) \statdist - (X_t(D_t) - m(D_t)\statdist)\Phi}_{\Uus}  \givenbig X_t, \veA_t} \\
\nonumber
&\geq \lamus^{1/2} \Ebig{\absbig{\big(X_{t+1}'(D_t) - m(D_t) \statdist - (X_t(D_t) - m(D_t)\statdist)\Phi\big)\xi^\top}  \givenbig X_t, \veA_t} \\
\label{eq:pf-hus-drift:interm-5}
&\geq \frac{\lamus^{1/2}\noiseortho}{\sqrt{N}} \indibracbig{\norm{X_t(D_t) - m(D_t)\statdist}_1 \leq m(D_t) \radiusus}.
\end{align}

Combining  (<ref>),  (<ref>) and (<ref>), we get

\begin{align}
\nonumber
&\mspace{23mu} \Ebig{\hus(X_{t+1}, D_t) \givenbig \givenbig X_t, \veA_t} - \hus(X_t, D_t)  \\
\nonumber
&\geq \frac{\lamus^{1/2}\noiseortho}{\sqrt{N}} \indibracbig{\norm{X_t(D_t) - m(D_t)\statdist}_1 \leq m(D_t) \radiusus} -  (2+\gapusmax)\hus(X_t, D_t) - \frac{2\lamus^{1/2}}{N} \\
\nonumber
&\geq \frac{\lamus^{1/2}\noiseortho}{\sqrt{N}} \indibracbig{\norm{X_t(D_t) - m(D_t)\statdist}_1 \leq m(D_t) \radiusus} -  \frac{(2+\gapusmax)K_1}{\sqrt{N}} - \frac{2\lamus^{1/2}}{N}.
\end{align}

Therefore, when $\hus(X_t, D_t) \leq K_1 / \sqrt{N}$,

\begin{align}
\nonumber
&\mspace{23mu} \Ebig{\hus(X_{t+1}, [N]) \givenbig \givenbig X_t, \veA_t} - \hus(X_t, [N]) \\
\label{eq:pf-hus-drift:goal-2}
&\geq \frac{\lamus^{1/2}\noiseortho}{\sqrt{N}} \indibracbig{\norm{X_t(D_t) - m(D_t)\statdist}_1 \leq m(D_t) \radiusus} -  \frac{(2+\gapusmax)K_1}{\sqrt{N}} - \frac{6\lamus^{1/2}}{\beta N} - \frac{8 \lamus^{1/2} \disy(Y_t)}{\beta}.
\end{align}

Summarizing the bounds (<ref>) and (<ref>) and substituting $K_1 = \lamus^{1/2} \noiseortho / (2+\gapus+\gapusmax)$ in the bounds, we get

* When $\hus(X_t, D_t) > K_1 / \sqrt{N}$ or $\norm{X_t(D_t) - m(D_t)\statdist}_1 \leq m(D_t) \radiusus$, we have

\begin{equation}
\label{eq:pf-hus-drift:case-1-2}
\Ebig{\hus(X_{t+1}, [N]) \givenbig \givenbig X_t, \veA_t} - \hus(X_t, [N]) \geq \frac{\gapus K_1}{\sqrt{N}} - \frac{6\lamus^{1/2}}{\beta N} -  \frac{8 \lamus^{1/2} \disy(Y_t)}{\beta}.
\end{equation}

* When $\hus(X_t, D_t) \leq K_1 / \sqrt{N}$ and $\norm{X_t(D_t) - m(D_t)\statdist}_1 > m(D_t) \radiusus$, we have

\begin{align}
\Ebig{\hus(X_{t+1}, [N]) \givenbig \givenbig X_t, \veA_t} - \hus(X_t, [N])
\label{eq:pf-hus-drift:case-3}
&\geq -  \frac{(2+\gapusmax)K_1}{\sqrt{N}} - \frac{6\lamus^{1/2}}{\beta N} -  \frac{8 \lamus^{1/2} \disy(Y_t)}{\beta}.
\end{align}

Note that the conditions separating the two cases above are in terms of $X_t(D_t)$. To complete the proof, we need to state the conditions in terms of $X_t([N])$ instead, observe that

\begin{align*}
\hus(X_t, D_t)
&\geq \hus(X_t, [N]) - \lamus^{1/2} (1-m(D_t)) \\
&\geq  \hus(X_t, [N]) - 2\lamus^{1/2} \disy(Y_t) / \beta - \lamus^{1/2}/(\beta N) \\
\norm{X_t(D_t) - m(D_t)\statdist}_1
&\leq \norm{X_t([N]) - \statdist}_1 + 2(1-m(D_t)) \\
&\leq \norm{X_t([N]) - \statdist}_1 + 4\disy(Y_t)/\beta + 2/(\beta N).
\end{align*}

Therefore, $\hus(X_t, [N]) > K_1 / \sqrt{N} + 2\lamus^{1/2} \disy(Y_t) / \beta + \lamus^{1/2}/(\beta N)$ implies $\hus(X_t, D_t) > K_1 / \sqrt{N}$; $\norm{X_t([N]) - \statdist}_1 \leq \radiusus -  6\disy(Y_t)/\beta - 3/(\beta N)$ implies that

\[
\norm{X_t(D_t) - m(D_t)\statdist}_1 \leq \radiusus - \frac{2\disy(Y_t)}{\beta} - \frac{1}{\beta N}  \leq \radiusus - (1-m(D_t)) \leq m(D_t) \radiusus,
\]

which implies the bounds in the lemma statement.

In the next lemma, we upper bound the drift of $\hst(x, D)$ under any policy.

Consider an instance of regular unstable RBs (<Ref>). Under any policy $\pi$, for any time step $t$, we have

\begin{equation}
\Ebig{\hst(X_{t+1}, [N]) \givenbig X_t, \veA_t} \leq (1-\gapst) \hst(X_t, [N]) + \frac{\lamst^{1/2}|\sspa|}{\sqrt{N}} + \frac{10\lamst^{1/2}|\sspa|}{\beta N} +  \frac{8\lamst^{1/2} \disy(Y_t)}{\beta} \quad a.s.,
\end{equation}

where $\lamst$ is the largest eigenvalue of $\Ust$.

Let $D_t$ be the set given in <Ref>.
Let $X_{t+1}'(D_t)$ be the state distribution at time $t+1$ when the actions for the arms in $D_t$ are resampled using $\lppriority$.