Model predictive control is almost optimal for restless bandits

## Abstract

Abstract We consider the discrete time infinite horizon average reward restless markovian bandit (RMAB) problem. We propose a model predictive control based non-stationary policy with a rolling computational horizon τ 𝜏 \tau italic_τ . At each time-slot, this policy solves a τ 𝜏 \tau italic_τ horizon linear program whose first control value is kept as a control for the RMAB. Our solution requires minimal assumptions and quantifies the loss in optimality in terms of τ 𝜏 \tau italic_τ and the number of arms, N 𝑁 N italic_N . We show that its sub-optimality gap is O ⁢ ( 1 / N ) 𝑂 1 𝑁 O(1/\sqrt{N}) italic_O ( 1 / square-root start_ARG italic_N end_ARG ) in general, and exp ⁡ ( − Ω ⁢ ( N ) ) Ω 𝑁 \exp(-\Omega(N)) roman_exp ( - roman_Ω ( italic_N ) ) under a local-stability condition. Our proof is based on a framework from dynamic control known as dissipativity . Our solution is easy to implement and performs very well in practice when compared to the state of the art. Further, both our solution and our proof methodology can easily be generalized to more general constrained MDP settings and should thus, be of great interest to the burgeoning RMAB community.

## 1 Introduction

This work investigates the sequential decision making problem of Restless Multi-Armed Bandits (RMAB for short) over an infinite discrete time horizon. In this problem there are $N$ statistically identical arms. At each time step, the decision maker must choose for each arm if they would like to pull the arm or leave it as is. The decision maker has a constraint $\alpha N$ on the maximal number of arms that they may pull at each time instance. Each arm has a known state belonging to a common finite state space and upon choosing an action, produces a known state and action dependent reward. Next, the arms evolve to a new state independently according to a known state-action dependent transition kernel. These arms are only coupled through the budget constraint on the number of arms that maybe pulled at each time instance. The state and reward are both revealed to the decision maker before the next decision needs to be made. The objective of the decision maker is to maximize the long term time average reward.

This problem was first proposed by Whittle (). Over the years RMABs have been used to model a number of practical problems. These applications include web-crawling, queuing, communication systems, scheduling problems and many more, (Veatch & Wein, ), (Dance & Silander, ), (Nino-Mora, ), (O’Meara & Patel, ). The problem of choosing a subset of tasks to perform among a larger collection of tasks under resource constraints shows up time and time again in various resource constrained control problems. For a comprehensive review on RMABs and their applications, the interested reader is directed towards Niño-Mora (). While the existence of an optimal policy for RMABs straightforward, Papadimitriou & Tsitsiklis () showed that the exact solution to this problem is PSPACE-hard. Consequently, most work focused on designing approximate solutions with good performance guarantees.

In the seminal paper, Whittle () suggested that under a condition known as indexability, an index can be associated with each state. This index is now referred to as the Whittle’s index (WI) and it was conjectured that a priority policy based on this index would be an optimal solution for this problem. This setting naturally lends itself to mean field approximations where one may replace the $N$ armed problem with a dynamical system in order to find these approximate solutions. (Weber & Weiss, ) were the first to point out that under a (hard to verify) condition on the dynamical system known as uniform global attractor property (UGAP), the (WI) was asymptotically optimal. Recently, many of the results in the RMAB literature have focused on two major aspects of the problem. Firstly, how quickly do proposed asymptotically optimal policies converge to the optimal solution as a function of the number of arms. Secondly, can the underlying assumptions such as indexability and UGAP be made less restrictive to obtain more general conditions under which optimal solutions can be found. In the former category, under the assumptions of Weber & Weiss (), Gast et al. () showed that the WI policy is exponentially close to the optimal solution. More recently still, several works have been able to show an exponential order of convergence, Gast et al. (); Hong et al. () under less restrictive assumptions. On the other hand, in the latter category, Hong et al. () showed an asymptotic convergence result under a far less restrictive assumption known as the *synchronization assumption*. Similar works following this steering dynamics include Yan () and Hong et al. () to show asymptotically optimal algorithms. Furthermore, Hong et al. () loosened the restrictions further to not only show asymptotic convergence but describe fundamental conditions in order to achieve exponential convergence rates for any algorithm.
A more comprehensive view of recent results can be found in Appendix A with descriptions of both assumptions and the algorithm.

In this work, we return to the dynamical control perspective but do not look to steer our system towards an optimal fixed point. Instead, we focus on the notion of *dissipativity*, to make the case for a well known set of policies known as *model predictive control*. We show that a simple MPC with a finite planning horizon produces an asymptotically optimal solution to the RMAB problem without any need to relax constraints. In doing so, not only do we bring a new perspective to this decades old problem but also describe a simple tractable algorithm that performs exceptionally well in practice with a small planning horizon.

#### Main Contributions

Our main contribution through this paper is to show that a simple computationally tractable algorithm with minimal assumptions is asymptotically optimal in the number of arms. Along the way we make several fresh observations about the nature of RMAB problems:

- 1.
First and foremost, we show that the sub-optimality gap of the natural MPC algorithm (called LP-update in the literature) is $\mathcal{O}(1/\sqrt{N})$, where $N$ is the number of arms. This result uses a minimal set of assumptions.
- 2.
Our proof is of independent interest since it uses a new framework known as *dissipativity*. Dissipativity is a closely studied phenomenon in the model predictive literature and is used to study how a policy drives a system towards optimal fixed points Damm et al. ().
- 3.
Returning to the dynamics around the fixed point, we can tighten the rate of convergence to $e^{-cN}$ under a local stability condition.
- 4.
Perhaps the most helpful portion of our results, for practical purposes; the MPC algorithm works well in practice and is easy to implement. It performs well both in terms of the number of arms $N$ as well as the computational time horizon $T$ beating state of the art algorithms in our bench marks.

#### Roadmap

The rest of the paper is organized as follows. We describe the system model and the corresponding linear relaxation in Section 2. We build the LP-update algorithm in Section 3 and present its performance guarantee in Section 4. We provide the main ingredients of the proof in Section 5 postponing the most technical lemmas to the appendix. We illustrate the performance of the algorithm in Section 6. The appendix contains additional literature review A, details of the algorithm and their extension to multi-constraints MDPs B, additional proofs C, D and details about the numerical experiments F.

## 2 System Model and linear relaxation

### 2.1 System Model

We consider an infinite horizon discrete time restless Markovian bandit problem parameterized by the tuple $\langle\mathcal{S},\mathbf{P}^{0},\mathbf{P}^{1},\mathbf{r}^{0},\mathbf{r}^{1}
;\alpha,N\rangle$. A decision maker is facing $N$ statistically identical arms, and each arm has a state that belongs to the finite state-space $\mathcal{S}$. At each time-instant, the decision maker observes the states of all arms $\textbf{s}=\{s_{1}\dots s_{N}\}$ and chooses a vector of actions $\textbf{a}=\{a_{1}\dots a_{N}\}$, where the action $a_{n}=1$ (respectively $a_{n}=0$) corresponds to pulling the arm $n$ (or to leave it). The decision maker is constrained to pull at most $\alpha N$ arms at each decision epoch.

The matrices $\mathbf{P}^{0}$ and $\mathbf{P}^{1}$ denote the transitions matrices of each arm and the vectors $\mathbf{r}^{0},\mathbf{r}^{1}$ denote the $|\mathcal{S}|$ dimensional vector for the rewards. We assume that all the rewards $\mathbf{r}$ lie between 0 0 and $1$. As the state-space is finite, this assumption can be made without loss of generality by scaling and centering the reward vector.
We assume that the transitions of all arms are Markovian and independent. This means that if the arms are in state $\textbf{s}\in\mathcal{S}^{N}$ and the decision maker takes an action $\textbf{a}\in\{0,1\}^{N}$, then the decision maker earns a reward $\sum_{n}r^{n}_{a_{n}}$ and the next state becomes $\textbf{s}^{\prime}\in\mathcal{S}^{N}$ with probability

$$ $\begin{split}&\mathbb{P}(\textbf{S}(t+1)=\textbf{s}^{\prime}|\textbf{S}(t)= \textbf{s},\textbf{A}(t)=\textbf{a},\dots\textbf{S}(0),\textbf{A}(0))\\ &\qquad=\mathbb{P}(\textbf{S}(t+1)=\textbf{s}^{\prime}|\textbf{S}(t)=\textbf{s },\textbf{A}(t)=\textbf{a})=\Pi_{n=1}^{N}P^{a_{n}}_{s_{n},s_{n}^{\prime}}.\end {split}$ (1) $$

It is important to note that the arms are only coupled through the budget constraint $\alpha$.

Let $\pi$ be a stationary policy mapping each state to a probability of choosing actions, i.e, $\pi:\mathcal{S}^{N}\to\Delta(\{0,1\}^{N})$ subject to the budget constraint, $\alpha$. Let $\Pi^{(N)}$ denote the space of all such policies. Given a policy $\pi$, mapping the joint state of the arms to a realization of the $N$ length joint action vector at time $t$, $\textbf{A}^{N}(t):=\{a_{1}(t),a_{2}(t)\dots a_{N}(t)\}$, and an initial set of arm states $\textbf{S}^{(N)}(0)=\textbf{s}$, we define the average gain of policy $\pi$ as

$$ $V_{\pi}^{(N)}(\textbf{s})=\lim_{T\to\infty}\frac{1}{T}\mathbb{E}_{\pi}\left[ \sum_{t=0}^{T-1}\frac{1}{N}\sum_{n=1}^{N}r^{A_{n}(t)}_{S_{n}(t)}\Bigg{|} \textbf{S}^{(N)}(0)=\textbf{s}\right].$ (2) $$

Here $(S^{N}_{n}(t),A^{N}_{n}(t))$ denotes the state-action pair of the $n$th arm at time $t$ and $r^{A^{N}_{n}(t)}_{S^{N}_{n}(t)}$ denotes the $S^{N}_{n}(t)^{\text{th}}$ entry of the $\mathbf{r}^{A_{n}(t)}$ vector. As the state-space and action space is finite, for any stationary policy $\pi$, the limit is well defined Puterman (). The RMAB problem amounts to computing a policy $\pi$ that maximizes the infinite average reward. We denote the optimal value of the problem as

$$ $\displaystyle V_{\text{opt}}^{N}(\textbf{s})$ $\displaystyle:=\max_{\pi\in\Pi^{(N)}}V_{\pi}^{(N)}(\textbf{s}).$ (3) $$

. The optimal policy exists and the limit is well defined, Puterman (). Under mild conditions (which will be verified in our case), this value does not depend on the initial state in which case we will simply write it as $V_{\text{opt}}^{N}$.

### 2.2 Alternative state representation via empirical distribution

In order to build an approximation of (3), we introduce an alternative representation of the state space, that we will use extensively in the paper. Given any joint state of the arms $\textbf{s}\in\mathcal{S}^{N}$, we denotes the empirical distribution of these arms as $\textbf{x}(\textbf{s})\in\Delta_{S}$, where $\Delta_{\mathcal{S}}$ is the simplex of dimension $|\mathcal{S}|$. $\textbf{x}(\textbf{s})$ is a vector with $|\mathcal{S}|$ dimensions and $x_{i}(\textbf{s})$ is the fraction of arms that are in state $i$. Next, given an action vector a, we denote by $\bm{u}(\textbf{s},\textbf{a})$ the empirical distribution of the state-action pairs ($s,1$). In words, $u_{i}(\textbf{s},\textbf{a})$ is the fraction of arms that are in state $i$ and that are pulled. Since no more than $N\alpha$ arms can be pulled at any time instance and no more than $N\textbf{x}_{s}$ arms can be pulled in state $s$, it follows that when x is fixed, $\bm{u}$ satisfies the following inequalities,

$$ $0\leq\bm{u}\leq\textbf{x}\hskip 21.68121pt\|\bm{u}\|_{1}\leq\alpha\|\textbf{x} \|_{1}=\alpha,$ (4) $$

where $\bm{u}\leq\textbf{x}$ denotes a component-wise inequality and $\|\cdot\|_{1}$ denotes the $l_{1}$ norm. We denote by $\mathcal{U}(\textbf{x})$ the set of feasible actions for a given x, *i.e.*, the set of $\bm{u}$ that satisfy (4).

### 2.3 Linear relaxation

We consider the following linear program:

$$ $\displaystyle V_{\text{opt}}^{\infty}:=\max_{{\bm{x}},\bm{u}\in\Delta_{ \mathcal{S}}}$ $\displaystyle~{}\mathbf{r}^{0}\cdot\textbf{x}+(\mathbf{r}^{1}-\mathbf{r}^{0}) \cdot\bm{u},$ (5a) Subject to: $\displaystyle\bm{u}\in\mathcal{U}(\textbf{x})$ (5b) $\displaystyle\textbf{x}=\textbf{x}\mathbf{P}^{0}+\bm{u}(\mathbf{P}^{1}-\mathbf {P}^{0})$ (5c) $$

This linear program is known to be a relaxation of (3) that is asymptotically tight, that is $V_{\text{opt}}^{N}\leq V_{\text{opt}}^{\infty}$ for all $N$ and $\lim_{N\to\infty}V_{\text{opt}}^{N}=V_{\text{opt}}^{\infty}$, see Gast et al. (); Hong et al. ().

To give some intuition on the relationship between (3) and (5), we remark that if ${\bm{X}(t)}:=\textbf{x}(\textbf{S}^{N}(t))$ is the empirical distribution of states at time $t$ and $\bm{U}(t)=\bm{u}(\textbf{S}(t),\textbf{A}(t))$ is the joint control, then it is shown in Gast et al. () that the Markovian evolution (1) implies

$$ $\mathbb{E}[{\bm{X}(t+1)}\mid{\bm{X}(t)},\bm{U}(t)]={\bm{X}}{(t)}\mathbf{P}^{0} +\bm{U}(t)(\mathbf{P}^{1}-\mathbf{P}^{0}).$ (6) $$

In (5), the variable $x_{i}$ corresponds to the time-averaged fraction of arms in state $i$ and similarly the variable $u_{i}$ corresponds to the time-averaged fraction of arms in state $i$ that are pulled. The constraint (5b) imposes that *on average*, no more than $\alpha N$ arms are pulled. This is in contrast with the condition imposed for problem (3) that enforces this condition at each time step.

## 3 Construction of the LP-update policy

### 3.1 The Finite-Horizon Mean Field Control Problem

To build the LP-update policy, we consider a *controlled dynamical system*, also called the *mean field model*, that is a finite-time equivalent of (5). For a given initial condition ${\bm{x}}(0)$ and a time-horizon $\tau$, the states and actions of this dynamical system are constrained by the evolution equations

$$ $\displaystyle\bm{u}(t)$ $\displaystyle\in\mathcal{U}({\bm{x}}(t))$ (7a) $\displaystyle{\bm{x}}({t+1})$ $\displaystyle={\bm{x}}(t)\mathbf{P}^{0}+\bm{u}(t)(\mathbf{P}^{1}-\mathbf{P}^{0 }),$ (7b) $$

$\forall t\in\{0\dots\tau-1\}$. In the above equation, (7b) should be compared with (6) and indicates that ${\bm{x}}(t)$ and $\bm{u}(t)$ correspond to the quantities $\mathbb{E}[{\bm{X}(t)}]$ and $\mathbb{E}[\bm{U}(t)]$ of the original stochastic system. As the constraint (7a) must be ensured by ${\bm{x}}(t)$ and $\bm{u}(t)$, this constraint (4) must be satisfied for the expectations: $\mathbb{E}[\bm{U}(t)]\in\mathcal{U}(\mathbb{E}[{\bm{X}(t)}])$.

The reward collected at time $t$ for this dynamical system is $\mathbf{r}^{0}\cdot{\bm{x}}(t)+(\mathbf{r}^{1}-\mathbf{r}^{0})\cdot\bm{u}(t)$. Let $\lambda$ be the dual multiplier of the constraint (5c) of an optimal solution of (5). We define a *deterministic* finite-horizon optimal control problem as:

$$ $\displaystyle W_{\tau}{(\textbf{x}(0))}$ $\displaystyle=\max_{{\bm{x}},\bm{u}}\sum_{t=0}^{\tau-1}\left(\mathbf{r}^{0} \cdot\textbf{x}(t)+(\mathbf{r}^{1}-\mathbf{r}^{0})\cdot\bm{u}(t)\right)+ \lambda\cdot{\bm{x}}(\tau),$ (8a) Subject to: ${\bm{x}}$ and $\bm{u}$ satisfy ( 7 ) for all $t\in\{0,\tau-1\}$ , (8b) $$

Before moving forward, the above equation deserves some remarks. First, for any finite $\tau$, the objective and the constraints (7) of the optimization problem (8) are linear in the variables $({\bm{x}}(t),\bm{u}(t))$. This means that this optimization problem is computationally easy to solve. In what follows, we denote by $\mu_{\tau}(\textbf{x})$ the value of $\bm{u}(0)$ of an optimal solution to (8).

Second, the definition of (8) imposes that the constraint $\|\bm{u}\|_{1}\leq\alpha\|\textbf{x}\|_{1}=\alpha$ holds for each time $t$. This is in contrast to the way this constraint is typically treated in RMAB problems, in which case (4) is replaced with the time-averaged constraint $\frac{1}{T}\sum_{t=0}^{T-1}\|\bm{u}(t)\|_{1}\leq\alpha$. The latter relaxation was introduced in Whittle () and is often referred to as Whittle’s relaxation (Avrachenkov & Borkar, ; Avrachenkov et al., ). This is the constraint that we use to write (5). Gast et al. () showed that for any finite $T$, the finite-horizon equivalent of (3) converges to (8) as $N$ goes to infinity. The purpose of this paper is to show that the solution of the finite $T-$horizon LP (8) provides an almost-optimal solution to the original $N-$arm average reward problem (3).

Last, as we will discuss later, taking $\lambda$ as the dual multiplier of the constraint (5c) helps to make a connection between the finite and the infinite-horizon problems (8) and (5). Our proofs will hold with minor modification by replacing $\lambda$ by 0 0 and in practise we do not use this multiplier for our algorithm.

### 3.2 The Model Predictive Control Algorithm

Figure: Algorithm 1 Evaluation of the LP-Update policy

The pseudo-code of the LP-update policy is presented in Algorithm 1. The LP-update policy takes as an input a time-horizon $\tau$. At each time-slot, the policy solves the finite horizon linear program (8) to obtain $\mu_{\tau}(\textbf{x})$ that is the value of $\bm{u}(0)$ of an optimal solution to (8). Note that such a policy may not immediately translate to an applicable policy as we do not require that $N\mu_{\tau}(\textbf{x})$ to be integers. We therefore use *randomized rounding* to obtain a feasible policy for our $N$ armed problem, $\textbf{A}^{N}(t)$. Applying these actions to each arm gives an instantaneous reward and a next state. This form of control has been referred to as *rolling horizon* Puterman () but more commonly referred to as *model predictive control* Damm et al. (). Our algorithm maybe summarized more succinctly as:

$$ $\textbf{S}^{(N)}(t)\xrightarrow[]{\text{Solve LP \eqref{EQ::VOPTDYN}}}\mu_{ \tau}(\textbf{x}(\textbf{S}^{(N)}(t)))\xrightarrow[\text{rounding}]{\text{ Randomized}}\textbf{A}^{N}(t)\xrightarrow[\text{new state}]{\text{Observe}} \textbf{S}^{(N)}(t+1).$ (9) $$

The randomized rounding procedure that we use is similar to the one described in Gast et al. (). We discuss it in Appendix B.1.

## 4 Main theoretical results

The main result of our paper is to show that a finite-horizon model predictive control algorithm, that we call the LP-update policy, is asymptotically optimal for the infinite horizon bandit problem. Note that this LP-update policy is introduced in Gast et al. (; ); Ghosh et al. () for finite-horizon restless bandit.

### 4.1 First result: LP-update is asymptotic optimal

Our result will show that the LP-update policy is asymptotically optimal as the number of arms $N$ goes to infinity, under an easily verifiable mixing assumption on the transition matrices. To express this condition, for a fixed integer and a sequence of action $\textbf{a}=(a_{1}\dots a_{k})\in\{0,1\}^{k}$, we denote by $P^{\textbf{a}}_{i,j}$ the $(i,j)$th entry of the matrix $\prod_{t=1}^{k}\mathbf{P}^{a_{k}}$. We then denote by $\rho_{k}$ the following quantity(^1^11The definition (10) is related to the notion of ergodic coefficient defined in (Puterman, ) that used a related constant to prove convergence of value iteration algorithms in span norm for the unconstrained discounted Markov Decision Process.):

$$ $\displaystyle\rho_{k}\triangleq\min_{s,s^{\prime}\in\mathcal{S},\textbf{a}\in \{0,1\}^{k}}\sum_{s^{*}\in\mathcal{S}}\min\{P^{a_{1},a_{2},\dots a_{k}}_{s,s^{ *}},P^{0,0,\dots 0}_{s^{\prime},s^{*}}\}$ (10) $$

In the above equation, the minimum is taken over all possible initial state $s,s^{\prime}$ and all possible sequence of action. The quantity $\rho_{k}$ can be viewed as the probability (under the best coupling) that two arms starting in states $s$ and $s^{\prime}$ reach the same state after $k$ iterations, if the sequence $a_{1}\dots a_{k}$ is used for the first arm while the second arm only uses the action 0 0. The assumption that we use for our result is that $\rho_{k}>0$ for some integer $k$.

###### Assumption 1 .

There exists a finite $k$ such that $\rho_{k}>0$.

While the assumption may look abstract, note that when the $\mathbf{P}^{0}$ matrix is ergodic, it ensures that assumption 1 holds. Indeed in this case, there exists a $k>0$ such that $P^{0\dots 0}_{ij}>0$ for all $i,j$ which would imply that $\rho_{k}>0$. Related assumptions and their relationship to ergodicity can be found in Hernandez-Lerma & Lasserre (). Assumption 1 is similar to the unichain and aperiodic condition imposed in Hong et al. (). Note that as this quantity involves the best coupling and not a specific coupling, it is more general than the synchronization assumption used in Hong et al. ().
We are now ready to state our first theorem, where we provide a performance bound of the average reward of the LP-update policy, that we denote by $V_{\text{LP},\tau}^{(N)}(\infty)$.

###### Theorem 4.1 .

Assume 1. There exist constants $C_{\lambda},C_{\Phi}>0$ such that for any $\epsilon>0$, there exists $\tau(\epsilon)$ such that, Algorithm 1 has the following guarantee of performance:

$$ $V_{\text{LP},\tau(\epsilon)}^{(N)}(\infty)\geq V_{\text{opt}}^{N}(\infty)- \epsilon-\left(\frac{kC_{\Phi}}{\rho}+C_{\lambda}C_{\Phi}+1\right)\frac{( \alpha N-\lfloor{\alpha N}\rfloor)}{N}-\left(\frac{k}{\rho_{k}}+C_{\lambda} \right)\left(\sqrt{\frac{|\mathcal{S}|}{N}}\right)$ (11) $$

This result shows that the LP-update policy becomes optimal as $\tau$ and $N$ go to infinity. The sub-optimality gap of LP-update decomposes in three terms. The first term corresponds to an upper-bound on the sub-optimality of using the finite-horizon $\tau$ when solving the LP-problem (8). Note, our proof shows that one can take $\tau(\epsilon)=\mathcal{O}(\frac{1}{\epsilon})$. Moreover, in the numerical section, we will show that choosing a small value like $\tau=10$ is sufficient for most problems. The second term corresponds to a rounding error for the budget $\alpha$: for all $\alpha$ such that $N\alpha$ is an integer, this term equals 0 0. The dominating error term is the last term, $O(1/\sqrt{N})$. It corresponds to the difference between the $N$-arm problem and the LP-problem with $\tau+\infty$.

###### Assumption 1 .

###### Theorem 4.1 .

### 4.2 Second result: Exponentially small gap under a stability condition

Theorem 4.1 shows that the sub-optimality gap of the LP-update policy is of order $O(1/\sqrt{N})$ under quite general conditions. While one could not hope for a better convergence rate in general, there are known cases for which one can construct policies that become optimal exponentially fast when $N$ goes to infinity. This is the case for Whittle index under the conditions of indexability, uniform global attractor property (UGAP), non-degeneracy and global exponential stability Gast et al. (). More details can be found in Appendix A. In this section, we show that LP-update also becomes optimal exponentially fast under essentially same conditions as the ones presented in Hong et al. (). The first condition that we impose is that the solution of the above LP-problem is non-degenerate (as defined in Gast et al. (); Hong et al. ()).

###### Assumption 2 (Non-degenerate) .

We assume that the solution (${\bm{x}}^{*},\bm{u}^{*})$ to the linear program (5) is unique and satisfies that $x^{*}_{i}>0$ for all $i\in\mathcal{S}$ and that there exists a (unique) state $i^{*}\in\mathcal{S}$ such that $0<u^{*}_{i^{*}}<x^{*}_{i^{*}}$.

The second condition concerns the local stability of a map around the fixed point.

###### Assumption 3 .

Assume 2 and let $P^{*}$ be the $|\mathcal{S}|\times|\mathcal{S}|$ matrix such that

$$ $\displaystyle P^{*}_{ij}=\left\{\begin{array}[]{ll}P^{0}_{ij}&\text{ if $i$ is such that $u^{*}_{i}<x^{*}_{i}$.}\\ P^{1}_{ij}-P^{1}_{i^{*}j}+P^{0}_{i^{*}j}&\text{ if $i$ is such that $u^{*}_{i} =x^{*}_{i}$.}\end{array}\right.$ $$

We assume that the matrix $P^{*}$ is stable, *i.e.*, that all of the norm of all the eigenvalues of $P^{*}$ are strictly smaller than $1$, except for one eigenvalue that equals $1$.

The two above conditions are equivalent to the assumption of non-degeneracy and local stability defined in Hong et al. ().

The last condition that we impose is a technical assumption that simplify the proofs.

###### Assumption 4 (Unicity) .

We assume that for all ${\bm{x}}\in\Delta_{\mathcal{S}}$, the LP program (8) has a unique solution.

It guarantees that the LP-update policy is uniquely defined. Note that the assumptions of unicity of the fixed point are often made implicitly in papers when authors talk about “the” optimal solution instead of “an” optimal solution.

###### Theorem 4.2 .

Assume 1, 2, 3, and 4 then there exist constants $C^{\prime},C^{\prime\prime}>0$ (independent on $N$ and $\tau$) such that for all $\epsilon>0$, with $\tau(\epsilon)$ (set according to Theorem  4.1) and $N$ such that $\alpha N$ is an integer, Algorithm 1 has the following guarantee of performance

$$ $V_{\text{LP},\tau(\epsilon)}^{(N)}(\infty)\geq V_{\text{opt}}^{N}(\infty)- \epsilon-C^{\prime}e^{-C^{\prime\prime}N}.$ (12) $$

The first term of the bound above is identical to the one of Theorem 4.1. What is more important is that the last term decays exponentially with $N$.

###### Assumption 2 (Non-degenerate) .

###### Assumption 3 .

###### Assumption 4 (Unicity) .

###### Theorem 4.2 .

## 5 Proofs: Main Ideas

In this section, we provide the major ingredients of the proofs of the two main theorems. We provide more details for the proof of Theorem 4.1 because this is the more original of the two. The proofs of all lemmas and some details of computation are deferred to Appendix  C.

### 5.1 Sketch for Theorem 4.1 :

Three major components are required in-order to complete the proof.

#### Part 1, Properties of the dynamical control problem ( 8 )

For ${\bm{x}},\bm{u}$, we denote by $\Phi({\bm{x}},\bm{u}):={\bm{x}}\mathbf{P}^{0}+\bm{u}(\mathbf{P}^{1}-\mathbf{P}
^{0})$ the deterministic transition kernel, and we recall that the instantaneous reward is $R({\bm{x}},\bm{u}):=\mathbf{r}^{0}\cdot{\bm{x}}+(\mathbf{r}^{1}-\mathbf{r}^{0})$. In Lemma C.2, we establish several properties that relate the finite-horizon problem (8) and the finite-horizon problem (5) that hold under Assumption 1. First we show that the average gain of the finite-time horizon problem (8) converges to the average gain of the infinite-horizon problem, that is $\lim_{\tau\to\infty}W_{\tau}({\bm{x}})/\tau=g^{*}$ for all ${\bm{x}}$. Second, we also show that the *bias function* $h^{\star}(\cdot):\Delta_{\mathcal{S}}\to\mathbb{R}$ given by :

$$ $h^{\star}({\bm{x}}):=\lim_{\tau\to\infty}W_{\tau}{({\bm{x}})}-\tau g^{\star}$ (13) $$

is well defined and Lipschitz-continuous with constant $k/\rho_{k}$ and that the convergence in (13) is uniform in ${\bm{x}}$. Moreover, the gain and the bias function satisfy

$$ $h^{\star}({\bm{x}})+g^{\star}=\max_{\bm{u}\in\mathcal{U}({\bm{x}})}R({\bm{x}}, \bm{u})+h^{\star}(\Phi(\textbf{x},\bm{u})).$ (14) $$

While both these definitions are well known in the average reward *unichain MDP* setting *without constraints*, Lemma C.2 establishes these definitions for the *constrained, deterministic problem*.

#### Part 2, Dissipativity and rotated cost

Let $({\bm{x}}^{*},\bm{u}^{*})$ be the(^2^22For clarity we will present the proof as if this point is unique although the proof can be rendered without this requirement.) optimal solution of the infinite-horizon problem (5), and let $l(\textbf{x},\bm{u}):=g^{\star}-\mathbf{r}^{0}\cdot\textbf{x}-(\mathbf{r}^{1}-
\mathbf{r}^{0})\cdot\bm{u}$. Following Damm et al. (), an optimal control problem with stage cost $l({\bm{x}},\bm{u})$ and dynamic ${\bm{x}}(t+1):=\Phi({\bm{x}},\bm{u})$ is called *dissipative* if there exists a *storage function* $\lambda:\Delta_{\mathcal{S}}\to\mathbb{R}$ that satisfy the following equation:

$$ $\tilde{l}({\bm{x}},\bm{u}):=l(\textbf{x},\bm{u})+\lambda(\textbf{x})-\lambda( \Phi(\textbf{x},\bm{u}))\geq l({\bm{x}}^{*},\bm{u}^{*})=\tilde{l}({\bm{x}}^{*} ,\bm{u}^{*})=0.$ $$

The cost, $\tilde{l}(\textbf{x},\bm{u})$ is called the *rotated cost function*.

In Lemma C.3, we show that our problem is dissipative by setting the storage function $\lambda(x):=\lambda\cdot x$, where $\lambda$ is the optimal dual multiplier of the constraint (8a). It is important to note, the rotated cost so defined is always non-negative.

#### Part 3, MPC is optimal for the deterministic control problem

By using our definition of rotated cost, we define the following minimization problem

$$ $\displaystyle L_{\tau}(\textbf{x})$ $\displaystyle:=\min\sum_{t=0}^{\tau-1}\tilde{l}(\textbf{x}(t),\bm{u}(t)).$ Subject to: ${\bm{x}}(t)$ and $\bm{u}(t)$ satisfy ( 7 ) for all $t\in\{0,T-1\}$ . (15) $$

By dissipativity, $L_{\tau}({\bm{x}})$ is monotone increasing. Moreover, we have that $L_{T}(\textbf{x})\geq 0=L_{T}({\bm{x}}^{*})$. Hence, the problem is *operated optimally at the fixed point*, $(x^{*},u^{*})$. The optimal operation at a fixed point is a key observation, made by several works, (Goldsztajn & Avrachenkov, ; Gast et al., ; Yan, ; Hong et al., ), we recover this result as a natural consequence of dissipativity.

By Lemma C.4, $L_{\tau}(\textbf{x})=\tau g^{\star}-W_{\tau}{({\bm{x}})}+\lambda\cdot{\bm{x}}$ because(^3^33In fact, the two control problems are equivalent: $\bm{u}^{*}$ is optimal for (8) if and only if it is for (15).) of a telescopic sum of the terms $\lambda\cdot{\bm{x}}(t)$. Combining this with (13) implies that $\lim_{\tau\to\infty}L_{\tau}({\bm{x}})=h^{\star}({\bm{x}})+\lambda\cdot{\bm{x}}$. As $L_{\tau}({\bm{x}})$ is monotone, it follows that, for any $\epsilon>0$, there exists a $\tau(\epsilon)$ such that for $L_{\tau(\epsilon)}({\bm{x}})-L_{\tau(\epsilon)-1}({\bm{x}})<\epsilon$. Note that the fact that $L_{\tau}({\bm{x}})$ is monotone is crucial to obtain this property, leading us to use dissipativity. Putting the three components together, we can now prove the final steps.

###### Proof of Theorem 4.1 .

We begin by considering the difference between the optimal value and the model predictive value function, we let ${\bm{X}(t)}$ denote the empirical distribution of the arms at time $t$ and $\bm{U}(t)$ be the corresponding empirical distribution of the actions. We drop the superscript $N$ for convenience in the proof below, but it should be noted that $\bm{U}(t)$ is always obtained from a randomized rounding procedure and hence, dependent on $N$.

$$ $\displaystyle V_{\text{opt}}^{N}(\textbf{x})-V_{\text{LP},T}^{(N)}(\textbf{x})$ $\displaystyle\leq\lim_{T\to\infty}\frac{1}{\tau}\sum_{t=0}^{T-1}\mathbb{E} \left[g^{\star}-R({\bm{X}}(t),{\bm{U}}(t))\right]$ $\displaystyle=\lim_{T\to\infty}\frac{1}{T}\sum_{t=0}^{T-1}\mathbb{E}\left[g^{ \star}-\mathbf{r}^{0}\cdot{\bm{X}}(t)-(\mathbf{r}^{1}-\mathbf{r}^{0})\cdot\bm{ u}(t)+(\mathbf{r}^{1}-\mathbf{r}^{0})\cdot\left(\bm{u}(t)-\bm{U}(t)\right)\right]$ $\displaystyle\leq\lim_{T\to\infty}\frac{1}{T}\sum_{t=0}^{T-1}\mathbb{E}\left[l ({\bm{X}}(t),\bm{u}(t))\right]+\frac{\alpha N-\lfloor{\alpha N}\rfloor}{N}$ (16) $$

The first inequality follows from the well known result $V_{\text{opt}}^{N}(\textbf{x})\leq g^{\star}$, (Yan, ; Gast et al., ; Hong et al., ). The last inequality follows from randomized rounding and Lemma B.1. Let $(A):=\lim_{T\to\infty}\frac{1}{\tau}\sum_{t=0}^{T-1}\mathbb{E}\left[l({\bm{X}}
(t),\bm{u}(t))\right]$ denotes the first term of (16). Adding and subtracting the storage cost we have:

$$ $\displaystyle(A)=\lim_{T\to\infty}\frac{1}{T}\sum_{t=0}^{T-1}\mathbb{E}\left[ \tilde{l}({\bm{X}}(t),\bm{u}(t))-\lambda\cdot{\bm{X}}(t)+\lambda\cdot\Phi({\bm {X}}(t),\bm{u}(t))\right].$ $$

Now note, by the dynamic programming principle $\tilde{l}({\bm{x}},\bm{u})=L_{\tau}({\bm{x}})-L_{\tau-1}(\Phi({\bm{x}},\bm{u}))$. Further, for any state ${\bm{x}}$ and its corresponding control $\bm{u}$ consider:

$$ $L_{\tau}({\bm{x}})-L_{\tau-1}(\Phi({\bm{x}},\bm{u}))=L_{\tau}({\bm{x}})-L_{ \tau}(\Phi({\bm{x}},\bm{u}))+L_{\tau}(\Phi({\bm{x}},\bm{u}))-L_{\tau-1}(\Phi({ \bm{x}},\bm{u})).$ $$

By choosing $\tau=\tau(\epsilon)$ so that, $L_{\tau}({\bm{x}})-L_{\tau-1}({\bm{x}})<\epsilon$, we can ensure,

$$ $\displaystyle L_{\tau}({\bm{x}})-L_{\tau-1}(\Phi({\bm{x}},\bm{u}))\leq L_{\tau }({\bm{x}})-L_{\tau}(\Phi({\bm{x}},\bm{u}))+\epsilon$ (17) $$

Plugging these inequalities together, introducing a telescopic sum, and manipulating the order of variables slightly (shown in Appendix C.4), we have the following relation:

$$ $\displaystyle(A)\leq\epsilon+\lim_{T\to\infty}$ $\displaystyle\frac{1}{T}\sum_{t=0}^{T-1}\mathbb{E}\left[L_{\tau}({\bm{X}}(t+1) )-L_{\tau}(\Phi({\bm{X}}(t),\bm{u}(t)))-\lambda\cdot{\bm{X}}(t+1)+\lambda\cdot \Phi({\bm{X}}(t),\bm{u}(t))\right].$ (18) $$

$L_{\tau}(\textbf{x})$ is Lipschitz in x for any $\tau$. Since, $\mathbb{E}[{\bm{X}}(t+1)|{\bm{X}}(t),\bm{U}(t)]=\Phi({\bm{X}}(t),\bm{U}(t))$, Lemma 1 of Gast et al. () implies:

$$ $\|\mathbb{E}[{\bm{X}}(t+1)|{\bm{X}}(t),\bm{U}(t)]-\Phi({\bm{X}}(t),\bm{u}(t)) \|_{1}\leq\frac{\sqrt{|\mathcal{S}|}}{\sqrt{N}}+\frac{C_{\Phi}(\alpha N- \lfloor{\alpha N}\rfloor)}{N},$ $$

where $C_{\Phi}$ is the Lipschitz constant of the map $\Phi$. The Lipschitz constant of $L_{\tau}(\cdot)$ in the $l1$ norm is $\frac{k}{\rho_{k}}$ and $C_{\lambda}$ is equal to $\|\lambda\|_{\infty}$. Passing to the limit $\tau\to\infty$ we have:

$$ $V_{\text{opt}}^{N}(\textbf{x})-V_{\text{LP},\tau}^{(N)}(\textbf{x})\leq \epsilon+\left(\frac{k}{\rho}+C_{\lambda}\right)\left(\frac{\sqrt{|\mathcal{S} |}}{\sqrt{N}}+\frac{C_{\Phi}(\alpha N-\lfloor{\alpha N}\rfloor)}{N}\right)+ \frac{(\alpha N-\lfloor{\alpha N}\rfloor)}{N}$ (19) $$

∎

###### Proof of Theorem 4.1 .

### 5.2 Theorem 4.2

The proof Theorem 4.2 is more classical and follows the same line as the proof of the exponential asymptotic optimality of Gast et al. (); Hong et al. (). The first ingredient is Lemma D.1 that shows that, by non-degeneracy, there exists a neighborhood of ${\bm{x}}^{*}$ such that $\mu_{\tau}({\bm{x}}^{*})$ is locally linear around ${\bm{x}}^{*}$. The second ingredient is Lemma D.2 that shows that ${\bm{X}(t)}$ is concentrated around ${\bm{x}}^{*}$ as $t\to\infty$, *i.e.*, for all $\epsilon>0$, there exists $C>0$ such that $\lim_{t\to\infty}\Pr[\|{\bm{X}(t)}-{\bm{x}}^{*}\|\geq\varepsilon]\leq e^{-CN}$. Combining the two points imply the result. For more details, see Appendix D.

## 6 Numerical illustrations

In this section, we illustrate numerically the performance of the LP-update policy. We show that it performs very well in practice and outperform classical heuristics in most cases. We choose to compare against two heuristics of the literature: the LP-index policy of Gast et al. () and the FTVA of Hong et al. (). We choose these two heuristics because they are natural and simple to implement and do not rely on hard-to-tune hyperparameters (further discussion for our choice can be found in Appendix  F). All parameters of the examples are provided in Appendix F. All codes are provided in supplementary material. In the future, a Github link will be provided.

#### Comparison on representative examples

For our first comparison, we consider three representative examples: a randomly generated example (in dimension 8), the main example used in Hong et al. () and described in their Appendix G.2 (abbreviated Example Hong et al. () in the following), and Example 2 of Figure 7.4 of Chen () (abbreviated Example Chen ()). We report the result in Figure 1. We observe that in all three cases, both the LP-update and the FTVA policies are asymptotically optimal, but the LP-update policy outperforms FTVA substantially. The situation for LP-priority is quite different: as shown in Gast et al. (), in dimension 8, the LP-priority policy is asymptotically optimal for most of the randomly generated examples. This is the case for example (a) of Figure 1 for which LP-update and LP-priority give an essentially equivalent performance and are essentially almost optimal.
Note that the authors of Chen () provide the numerical value of the optimal policy $N$ for Example Chen (), which can be done because $|\mathcal{S}|=3$. We observe that the LP-update is very close to this value.

Figure: Figure 1: Performance as a function of $N$
Refer to caption: x1.png

#### System dynamics:

To explore the difference of performance between the LP-update policy and FTVA, we study the dynamics as a function of time of the different policies for the first two examples of Figure 1. We plot in Figure 2 the evolution of the system as a function of time (for a single trajectory). We observe that if the distance between $X_{t}$ and $x^{*}$ are similar for LP-update and FTVA, the rotated cost is much smaller for LP-update, which explains its good performance.

Figure: Figure 2: Distance $\|X_{t}-x^{*}\|$ and rotated cost as a function of time for $N=100$.
Refer to caption: x4.png

To explore the system dynamics in more details, we focus on Example Chen () and study the behavior of the LP-update, FTVA and LP-priority policy. In Figures 3(a,b,c), we present a trajectory of each of the three policies: each orange point corresponds to a value of $X_{t}\in\Delta_{\mathcal{S}}$ (as this example is in dimension $|\mathcal{S}|=3$, the simplex $\Delta_{\mathcal{S}}$ can be represented as a triangle). This example has an unstable fixed point (Assumption 3 is not satisfied). Both LP-update and FTVA concentrate their behavior around this fixed point but the LP-priority policy exhibits two modes. When concentrating on the rotated cost (Figure 3(d)), we observe that it is much smaller for LP-update when compared to FTVA. This explains why LP-update performs better as shown in Figure 1(c).

Figure: Figure 3: Example 2 from Gast et al. 22. Simulation for $N=100$.
Refer to caption: x8.png

#### Influence of parameters

In Figure 4, we study the influence of different parameters on the performance of LP-update. In each case, we generated $20$ examples and plot the average “normalized” performance among these $20$ examples. By “normalized”, we mean that for each example, we divide the performance of the policy by the value of the LP problem (5). The first plot 4(a) shows that the influence of the parameter $\tau$ is marginal (the curves are not distinguishable for $\tau\in\{3,5,10\}$). In the second plot 4(b), we observe that the sub-optimality gap of LP-update does not seem to depend much on the state space size whereas the one of FTVA does: FTVA performs worse when $|\mathcal{S}|$ is large probably because there are fewer synchronized arms. In the last plot 4(c), we study the performance as a function of the budget $\alpha$.

Figure: Figure 4: Comparison of the gains as a function of some parameters.
Refer to caption: x12.png

## 7 Conclusion

In this paper, we study the problem of constructing an efficient policy for restless multi-armed bandit for the average reward case. We show that under quite general condition, a simple model-predictive control algorithm provides a solution that achieves both the best known upper bound up to know ($O(1/\sqrt{N})$ or $\exp(-\Omega{N})$ for stable non-degenerate problems), but also works very efficiently in practice. Our paper provides the first analysis of this policy for the average reward criterion.

Our paper uses a novel framework based on dissipativity that helps up to make a crisp connection between the finite- and the infinite-horizon problems and we are the first to us it in this context. We believe that this framework is what makes our analysis simple and easily generalizable to other settings. As an example, we discuss in Appendix B.2 potential generalization to the multi-action multi-constraint bandits, that is an interesting direction for future work.

## Appendix A Additional Literature Review

The study of bandit problems dates back at least to the 1930’s in the form of sequential design of experiments, (Thompson, ), (Wald, ). However, the more modern perspective of Markovian bandits in terms of MDPs is largely credited to the monograph by Bellman et al. (). This monograph brought to light the combinatorial difficulty of this problem even over a finite horizon. The *rested* variant of the multi-armed bandit problem was famously resolved through Gittins () which proposed a remarkably simple structure to the solution by allotting indices to the states of the bandit arms.

In his seminal work Whittle () generalized the problem to the *restless multi-armed bandit problem* in its modern form and conjectured that under a condition known as indexability, an index policy is asymptotically optimal (in the number of arms) for the RMAB problem. Weber & Weiss () showed a counter example to this conjecture, further going on to show in the same work that the conjecture did hold under a *global attractor* and ergodicity condition of a single-arm. This result was generalized to the multi-action setting by Hodge & Glazebrook (). On the other hand, Verloop (); Gast et al. () generalized the set of index policies by constructing a set of priority policies known as LP priority index. More surprisingly, Gast et al. (; ) showed these priority policies are in fact exponentially close to being optimal, under three conditions: uniform global attractor property (UGAP), *non-degeneracy* of the fixed point and global exponential stability Gast et al. (). The restriction of indexability is removed in Gast et al. () and the condition of UGAP was first removed Hong et al. ().

Gast et al. (; ) observed that under the conditions of non-degeneracy, there exists a local neighborhood around the fixed point where the priority policy set is affine. This is the key observation which leads to the exponentially small error bounds under simple local stability conditions. This observation around the fixed point drove several works (Hong et al., ; ; ; Goldsztajn & Avrachenkov, ; Yan, ) to shift their perspective towards policies that look to steer the dynamical system towards the optimal fixed point. Hong et al. () created a virtual system that was driven to the fixed point. To transform this into an asymptotically optimal policy, they required the real system to synchronize with the virtual system, hence, they developed the *synchronization assumption*. Unlike Hong et al. () which passively allowed arms to align with each other eventually, Hong et al. () used a *two set policy* to actively align non aligned arms to the fixed point. In contrast to this form of control, Yan () designed the *align and steer* policy to steer the mean field dynamical control system towards the fixed point. Yan ()’s policy is asymptotically optimal under the assumptions of *controllability* of the dynamical system to more directly move the state of the system towards the fixed point over time. In parallel to this work, Goldsztajn & Avrachenkov () used a similar idea of mean field dynamical control for weakly coupled systems under a *relaxed but generalized constraint set* if the single arm process is unichain and aperiodic.

In contrast to these approaches our algorithm is not explicitly designed at the outset to steer the corresponding dynamical system towards a fixed point. Rather, we show that the model predictive control algorithm is designed to solve an equivalent *rotated cost minimization* problem. In doing so, it produces policies that are close to optimal to the infinite horizon average reward problem. We thus produce a policy that is easily generalizable to the multi-action, general constraint setting *without relaxing the constraints* using ideas from Whittle’s relaxation. Further, when the fixed point is unique, minimizing the rotated cost problem will coincide with steering the dynamical control system towards the fixed point, allowing us to recover exponentially close asymptotic optimality bounds.

There are also a lot of recent papers on RMABs for the finite-horzon reward criteria
Hu & Frazier (); Zayas-Caban et al. (); Brown & Smith (); Ghosh et al. (); Zhang & Frazier (); Gast et al. (), or the infinite-horizon discounted reward Zhang & Frazier (); Ghosh et al. (). Our results are close in spirit to those since our policy shows how to use a finite-horizon policy to construct a policy that is asymptotically optimal for the infinite-horizon rewards.

## Appendix B Algorithm details

### B.1 Rounding procedure

The rounding procedure is composed of two steps:

- 1.
First, given $\bm{u}=\mu_{\tau}({\bm{x}}(\textbf{S}(t)))$, we need to construct a vector $\bm{U}^{N}$ that is *as close as possible* to $\bm{u}$ while satisfying the constraints that $\sum_{i}u_{i}<=\alpha$ and $N\bm{U}^{N}$ is an integer.
- 2.
Second, we use $\bm{U}^{N}$ to construct a feasible sequence of actions $\textbf{A}(t)\in\{0,1\}^{N}$ that we can apply to $\textbf{S}(t)$.

The second step is quite easy: as $Nx_{i}$ and $NU^{N}_{i}$ are integers for all $i$, we can construct a feasible sequence of action by applying the action $1$ to $NU^{N}_{i}$ arms that in state $i$ and action 0 0 to the $N(x_{i}-U^{N}_{i})$ others. The first step is more complicated because we want to construct a $\bm{U}^{i}$ that is *as close as possible* as $u_{i}$. Before giving our procedure, let us illustrate it through two examples of $(N{\bm{x}},N\bm{u})$ for $\alpha=0.5$:

| $N{\bm{x}}$ | $N\bm{u}$ | Difficulty | Possible $\bm{U}^{N}$ by Algorithm 2 |
| --- | --- | --- | --- |
| $(10,10,10,9)$ | $(10,9.5,0,0)$ | $\|u\|_{1}\geq\alpha$ | $(10,9,0,0)$ |
|  |  | We need to randomize | $(10,6,0,0)$ with proba $0.7$ |
| $(10,10,10,9)$ | $(10,5.7,0.2,0)$ | to be as close as | $(10,5,1,0)$ with proba $0.2$ |
|  |  | possible to $\bm{u}$ | $(10,5,0,0)$ with proba $0.1$ |
| $(10,10,10,9)$ | $(10,4.9,4.6,0)$ | $\|u\|_{1}\geq\alpha$ and we | $(10,5,4,0)$ with proba $0.4$ |
|  |  | need to randomize | $(10,4,5,0)$ with proba $0.6$ |

This leads us to construct Algorithm 2. This algorithm first construct a $\bm{v}$ that is as close as possible to $\bm{u}$ while satisfying $\|\bm{v}\|\leq\alpha$. There might be multiple choices for this first step and the algorithm can choose any. For instance, applied to the example such that ${\bm{x}}=(10,10,10,9)$ and $\bm{u}=(10,4.9,4.6,0)$, this algorithm would first construct a vector $\bm{v}$ that can be any convex combination of $(10,4.4,4.6,0,0)$ and $(10,4.9,4.1,0,0)$. For instance, if the algorithm chose the vector $\bm{v}=(10,4.4,4.6,0)$, then the algorithm would produce $(10,5,4,0)$ with probability $0.4$ and $(10,4,5,0)$ with probability $0.6$. Once this is done, the algorithm outputs a $\bm{U}^{N}$ that satisfies the constraints and such that $\mathbb{E}[\bm{U}^{N}]=\bm{v}$. In some cases, there might be multiple distributions that work. The algorithm may output any. An efficient procedure to implemented this last step is provided in Section 5.2.3 of Ioannidis & Yeh ().

Figure: Algorithm 2 Randomized rounding

###### Lemma B.1 .

Algorithm 2 outputs a random vector $\bm{U}^{N}$ such that $\bm{U}^{N}\in\mathcal{U}({\bm{x}})$ and such that

$$ $\displaystyle\|\mathbb{E}[\bm{U}^{N}]-\bm{u}\|_{1}\leq\frac{\lfloor\alpha N \rfloor-\alpha N}{N}.$ $$

###### Lemma B.1 .

### B.2 Generalization to multi-constraints MDPs

In this section, we give some clues on how and why our algorithm and results can be generalized to multi-action multi-constraints MDPs. Specifically, we make the following modifications to our model:

- •
Instead of restricting an action to be in $\{0,1\}$, we can consider any finite action set $\mathcal{A}=\mathcal{A}=\{0\dots A-1\}$.
- •
Instead of having a single constraint $\sum_{n}A_{n}(t)\leq\alpha N$, we consider that there are $K$ types of constraints and that when taking action $a$ in state $i$, this action consumes $D^{a}_{i,k}\geq 0$ of resource $k$. The action 0 0 is special and does not consume resources. We impose $K$ resource constraints:
$\displaystyle\sum_{n=1}^{N}D^{A_{n}(t)}_{S_{n}(t),k}\leq C_{k}$
$\displaystyle\forall k\in{1\dots K}.$

As these new constraints are linear, one can define LP relaxation equivalence of (5) and (8), which leads to an LP-update adapted to this new setting. The only major modification of the algorithm concerns the randomized rounding and is where one would need to have an action 0 0 in order to guarantee the feasibility of the solution.

- •
The generalization of Assumption 1 to the multi-action case is straightforward. Hence, one can prove an equivalent of Theorem 4.1 for multi-action multi-constrained MDPs with a rate of convergence of $O(1/\sqrt{N})$.
- •
To obtain a generalization of Theorem 4.2, there are two main difficulties
1.
The first is to redefine the notion of a non-degenerate fixed-point in order to replace Assumption 2 by one adapted to the multi-action multi-constraint case. To that end, we believe that the most appropriate notion is the one presented in Gast et al. (). It provides the notion of non-degeneracy that uses a linear map.
2.
The second is to provide an equivalent of the stability Assumption 3. To do so, we can again use the notion of non-degeneracy defined in Gast et al. () that defines a linear map. The stability of this linear map should suffice to prove the theorem.

## Appendix C Proof of Theorem 4.1

### C.1 Part 1, Properties of the dynamical control problem 8 :

We begin by formally defining the space of stationary policies $\bar{\Pi}$ as the set of all policies that map ${\bm{x}}\in\Delta_{\mathcal{S}}$ to an action $\bm{u}\in\mathcal{U}({\bm{x}})$. Note, any policy in $\Pi^{(N)}$ must satisfy this condition. Now, for an initial distribution ${\bm{x}}\in\Delta_{\mathcal{S}}$ define the discounted infinite horizon reward problem,

$$ $\displaystyle V^{\infty}({\bm{x}}):$ $\displaystyle\max_{\pi\in\bar{\Pi}}\sum_{t=0}^{\infty}\beta^{t}R({\bm{x}}(t), \bm{u}(t))$ (20) Subject to the dynamics ( 7 ) (21) $$

where $x(t)$ denotes the trajectory induced by the policy $\pi$. Note, this limit is always well defined for $\beta<1$ and the optimal policy exists (Puterman, ). The main idea of this section will be an exercise in taking appropriate limit sequences of the discount factor $\beta$ to define the gain and bias for the constrained average reward problem.

Before we proceed we introduce a little notation for convenience in writing down our proof. For any pair $({\bm{x}},\bm{u})$ define an equivalent $\{y_{s,a}\}_{s,a}$ with:

$$ $\displaystyle y_{s,1}$ $\displaystyle:=\bm{u}_{s}$ $\displaystyle y_{s,0}$ $\displaystyle:={\bm{x}}_{s}-\bm{u}_{s}$ $$

Here $y_{s,a}(t)$ represents the fraction of arms in state $s$, taking action $a$ at time $t$. It is also convenient to introduce a concatenated reward vector $\mathbf{r}:=[\mathbf{r}^{0},\mathbf{r}^{1}]$ and a concatenated transition kernel $\mathbf{P}:=[\mathbf{P}^{0},\mathbf{P}^{1}]$. If $\textbf{y}:=\{y_{s,a}\}$ represents the $\mathcal{S}\times 2$ vector, it is not hard to check that $R({\bm{x}},\bm{u}):=\mathbf{r}^{0}\cdot{\bm{x}}+(\mathbf{r}^{1}-\mathbf{r}^{0}
)\cdot\bm{u}=\mathbf{r}\cdot\textbf{y}$ and $\Phi({\bm{x}},\bm{u}):=\mathbf{P}^{0}\cdot{\bm{x}}+(\mathbf{P}^{1}-\mathbf{P}^
{0})\cdot\bm{u}=\mathbf{P}\cdot\textbf{y}$. We can rewrite the discounted problem as follows:

$$ $\displaystyle V^{\infty}({\bm{x}}):$ $\displaystyle\max_{\pi\in\bar{\Pi}}\sum_{t=0}^{\infty}\beta^{t}\mathbf{r}\cdot \textbf{y}(t)$ (22) $\displaystyle y_{s,0}(t)+y_{s,1}(t)=$ $\displaystyle{\bm{x}}_{s}(t)$ (23) $\displaystyle\sum_{a}y_{s,a}(t+1)=$ $\displaystyle\sum_{s^{\prime},a^{\prime}}P^{a^{\prime}}_{s^{\prime},s}y_{s^{ \prime},a^{\prime}}(t)\hskip 14.45377pt\forall\hskip 3.61371ptt$ (24) s.t. $\displaystyle\sum_{s}y_{s,1}(t)\leq$ $\displaystyle\alpha\hskip 14.45377pt\forall\hskip 3.61371ptt$ (25) $$

We now state the following lemma on the continuity of the value function in the state:

###### Lemma C.1 .

Under assumption 1, for any $\beta\in(0,1]$ we have,

$$ $V^{\infty}({\bm{x}})-V^{\infty}({\bm{x}}^{\prime})\leq\frac{k}{\rho_{k}}\|{\bm {x}}-{\bm{x}}^{\prime}\|_{1}$ (26) $$

Note, this continuity result is independent of the discount factor $\beta$. We postpone the proof to Appendix E. This will be critical to proving the main result of the subsection below.

###### Lemma C.2 .

Consider the infinite horizon average reward problem under the synchronization assumption. There exists a constant gain and a bias function from $\Delta_{\mathcal{S}}\to\mathbb{R}$ denoted by $g^{\star}$ and $h^{\star}(\cdot)$ respectively defined by:

$$ $g^{\star}:=\lim_{i\to\infty}(1-\beta_{i})V^{\infty}_{\beta_{i}}({\bm{x}})$ $$

for an appropriate sequence $\beta_{i}\to 1$ as $i\to\infty$. For the same sequence of $\beta_{i}$ we can define the Lipschitz continuous bias function,

$$ $h^{\star}(\textbf{x}):=\lim_{i\to\infty}\sum_{t=0}^{\infty}\beta_{i}^{t}R({\bm {x}}(t),\bm{u}(t))-(1-\beta_{i})^{-1}g^{\star}$ $$

Furthermore, they satisfy the following fixed point equations,

$$ $g^{\star}+h^{\star}(\textbf{x})=\max_{\bm{u}\in\mathcal{U}({\bm{x}})}R({\bm{x} }(t),\bm{u}(t))+h^{\star}(\Phi({\bm{x}},\bm{u}))$ (27) $$

###### Proof.

Now note, $\textbf{y}\cdot\textbf{1}=1$, where 1 is the all $1^{\prime}$s vector. We will overload the notation slightly by letting $\textbf{y}(t)$ be the optimal trajectory taken by the optimal policy $\pi^{*}_{\beta}$ for 20. Hence, for any constant $g$ we have:

$$ $V^{\infty}({\bm{x}})=(1-\beta)^{-1}g+\sum_{t=0}^{\infty}\beta^{t}[\mathbf{r}-g \textbf{1}]\cdot\textbf{y}(t)$ (28) $$

One can check that $V^{\infty}({\bm{x}})$ satisfies the following *Bellman equation*,

$$ $V^{\infty}({\bm{x}})=\mathbf{r}\cdot\textbf{y}(t)+\beta V^{\infty}(\mathbf{P} \cdot\textbf{y})$ (29) $$

Now note, given $\textbf{y}(t)$ at time $t$, ${\bm{x}}(t+1)$ is given by a linear (hence, continuous) map ${\bm{x}}(t+1):=\mathbf{P}\cdot\textbf{y}(t)$. Further note, given $\pi^{*}_{\beta}$, $\textbf{y}(t)$ is upper hemicontinuous with respect to ${\bm{x}}(t)$. Let ${\bm{H}}_{\beta}$ denote the map induced by $\pi^{*}_{\beta}$ from ${\bm{x}}$ to y, then ${\bm{H}}_{\beta}({\bm{x}}(t))=\textbf{y}(t)$. Combining the two maps, we have $\mathbf{P}{\bm{H}}_{\beta}:\Delta_{\mathcal{S}}\to\Delta_{\mathcal{S}}$. Since, $\Delta_{\mathcal{S}}$ is closed and bounded, by Kakutani’s fixed point theorem there exists atleast one fixed point ${\bm{x}}_{\beta}:=\mathbf{P}{\bm{H}}_{\beta}\cdot{\bm{x}}_{\beta}$ if $\textbf{y}_{\beta}$ is the corresponding value for ${\bm{x}}_{\beta}$ we have ${\bm{x}}_{\beta}:=\mathbf{P}\cdot\textbf{y}_{\beta}$. Combining these observations allows us to write,

$$ $V^{\infty}({\bm{x}}_{\beta})=\mathbf{r}\cdot\textbf{y}_{\beta}+\beta V^{\infty }({\bm{x}}_{\beta})$ $$

We will set $g_{\beta}:=\mathbf{r}\cdot\textbf{y}_{\beta}$ and designate this value as the *gain* of our problem. It then follows that $V^{\infty}({\bm{x}}_{\beta})=(1-\beta)^{-1}g_{\beta}$. Now, we can define a *bias function* $h_{\beta}(\cdot)$ as the remaining terms of the equation (28), concretely,

$$ $h_{\beta}({\bm{x}})=\sum_{t=0}^{\infty}\beta^{t}[\mathbf{r}-g_{\beta}\textbf{1 }]\cdot\textbf{y}(t)$ (30) $$

and this definition of $h_{\beta}(\cdot)$ is well defined since the sum is bounded. It now follows that,

$$ $V^{\infty}({\bm{x}})=(1-\beta)^{-1}g_{\beta}+h_{\beta}({\bm{x}})$ (31) $$

More importantly, by plugging these definitions into (29) and making the dependence of y on the maximal policy, we obtain the following recursion equation:

$$ $g_{\beta}+h_{\beta}({\bm{x}})=\max_{\textbf{y}:\{\textbf{y}_{1}\in\mathcal{U}( {\bm{x}})\}}g_{\beta}+\beta h_{\beta}(\mathbf{P}\cdot\textbf{y})$ (32) $$

Further, due to the Lipschitz continuity, Lemma C.1 of $V^{\infty}({\bm{x}})$ with respect to ${\bm{x}}$ we have, for any discount factor $\beta_{i}<1$,

$$ $\displaystyle(1-\beta_{i})\|V^{\infty}_{\beta_{i}}({\bm{x}})-V^{\infty}_{\beta _{i}}({\bm{x}}^{\prime})\|_{1}\leq\frac{k(1-\beta_{i})}{\rho_{k}}\|{\bm{x}}-{ \bm{x}}^{\prime}\|$ $$

. Hence, by choosing a sequence $\beta_{i}\uparrow 1$ as $i\to\infty$ we obtain:

$$ $\lim_{i\to\infty}(1-\beta_{i})\|V^{\infty}_{\beta_{i}}({\bm{x}})-V^{\infty}_{ \beta_{i}}({\bm{x}}^{\prime})\|_{1}\to 0$ $$

. From Puterman () we know there exists a sequence $i\to\infty$ such that $(1-\beta_{i})V^{\infty}_{\beta_{i}}({\bm{x}})\to\lim_{T\uparrow\infty}\frac{W_
{T}{({\bm{x}})}}{T}$. In particular, this implies that $(1-\beta_{i})V^{\infty}_{\beta_{i}}({\bm{x}})$ converges to a constant value which we shall denote by $g^{\star}$. This gives us the first result.
Critically, by noting that $g_{\beta_{i}}:=(1-\beta_{i})V^{\infty}_{\beta_{i}}({\bm{x}}_{\beta_{i}})$ we see that $g^{\star}$ must be the limit point of $R({\bm{x}}_{\beta_{i}},\bm{u}_{\beta_{i}})=:\mathbf{r}\cdot\textbf{y}_{\beta_{
i}}=:g_{\beta_{i}}$, hence, it is the solution to (5).

Next, once again leveraging Lemma C.1 we have, for any $\beta<1$,

$$ $\|V^{\infty}({\bm{x}})-(1-\beta)^{-1}g_{\beta}\|_{1}:=\|V^{\infty}({\bm{x}})-V ^{\infty}({\bm{x}}_{\beta})\|_{1}\leq\frac{k}{\rho_{k}}\|{\bm{x}}-{\bm{x}}_{ \beta}\|_{1}$ $$

Plugging these results into (31) we see that,

$$ $\|h_{\beta}({\bm{x}})\|_{1}\leq\frac{k}{\rho_{k}}\|{\bm{x}}-{\bm{x}}_{\beta}\| _{1}$ (33) $$

is bounded by a constant for all $\beta<1$. By the same token $h_{\beta}$ is Lipschitz with constant $\frac{k}{\rho_{k}}$. Choosing the same sequence $\beta_{i}$ in (30) and passing to the limit we obtain :

$$ $h^{\star}({\bm{x}}):=\lim_{i\to\infty}h_{\beta_{i}}({\bm{x}})=\lim_{i\to\infty }\sum_{t=0}^{\infty}\beta_{i}^{t}[\mathbf{r}-g_{\beta_{i}}\textbf{1}]\cdot \textbf{y}(t)=\lim_{i\to\infty}\sum_{t=0}^{\infty}\beta_{i}^{t}\mathbf{r}\cdot \textbf{y}(t)-(1-\beta_{i})^{-1}g^{\star}$ $$

Thus, we have defined both the gain and bias $g^{\star}$ and $h^{\star}(\cdot)$ for our *deterministic average reward problem*. Recall, $R({\bm{x}},\bm{u}):=\mathbf{r}^{0}\cdot{\bm{x}}+(\mathbf{r}^{1}-\mathbf{r}^{0}
)\cdot\bm{u}=\mathbf{r}\cdot\textbf{y}$ giving us the second result:

$$ $h^{\star}({\bm{x}})=\lim_{i\to\infty}\sum_{t=0}^{\infty}\beta_{i}^{t}R({\bm{x} }(t),\bm{u}(t))-(1-\beta_{i})^{-1}g^{\star}$ (34) $$

Further, this sequence $\beta_{i}\uparrow 1$ in (32) yields the following fixed point equation:

$$ $g^{\star}+h^{\star}(\textbf{x})=\max_{\bm{u}\in\mathcal{U}({\bm{x}})}R({\bm{x} }(t),\bm{u}(t))+h^{\star}(\Phi({\bm{x}}(t),\bm{u}(t)))$ $$

which completes the proof.
∎

###### Lemma C.1 .

###### Lemma C.2 .

###### Proof.

### C.2 Part 2: Dissipativity

###### Lemma C.3 .

Let $\tilde{l}({\bm{x}},\bm{u})=g^{*}-R({\bm{x}},\bm{u})+\lambda\cdot{\bm{x}}-
\lambda\cdot\Phi({\bm{x}},\bm{u})$. Then:

- •
$\tilde{l}({\bm{x}},\bm{u})\geq 0$ for all ${\bm{x}}\in\Delta_{\mathcal{S}}$ and $\bm{u}\in\mathcal{U}({\bm{x}})$.
- •
If $({\bm{x}}^{*},\bm{u}^{*})$ is a solution to (5), then $\tilde{l}({\bm{x}}^{*},\bm{u}^{*})=0$.

This shows that our problem is dissipative.

###### Proof.

Recall that $\lambda$ is the optimal dual variable of the constraint $\textbf{x}=\mathbf{P}^{0}\cdot\textbf{x}+(\mathbf{P}^{1}-\mathbf{P}^{0})\cdot
\bm{u}$ of problem (5). By strong duality and Lemma  C.2, this implies that

$$ $\displaystyle g^{*}$ $\displaystyle=\max_{{\bm{x}},\bm{u}}~{}R({\bm{x}},\bm{u})+\lambda\cdot(\mathbf {P}^{0}\cdot{\bm{x}}+(\mathbf{P}^{1}-\mathbf{P}^{0})\cdot\bm{u})-\lambda\cdot{ \bm{x}}$ $\displaystyle\text{Subject to: }\quad\bm{u}\in\mathcal{U}(\textbf{x}).$ $$

Recall that $\Phi({\bm{x}},\bm{u})=\mathbf{P}^{0}\cdot{\bm{x}}+(\mathbf{P}^{1}-\mathbf{P}^{
0})\cdot\bm{u}$. This implies that

$$ $\displaystyle 0$ $\displaystyle=\min_{{\bm{x}},\bm{u}}~{}\tilde{l}({\bm{x}},\bm{u})$ $\displaystyle\text{Subject to: }\bm{u}\in\mathcal{U}({\bm{x}}),$ $$

which implies the results of the lemma.
∎

###### Lemma C.3 .

###### Proof.

### C.3 Part 3: MPC is optimal for the deterministic problem

Recall that $L_{\tau}(\textbf{x}):=\min\sum_{t=0}^{\tau-1}\tilde{l}({\bm{x}}(t),\bm{u}(t))$ with ${\bm{x}}(0)=\textbf{x}$.

###### Lemma C.4 .

For all ${\bm{x}}\in\Delta_{\mathcal{S}}$, for any $\epsilon>0$, there exists $t$ such that:

$$ $\displaystyle L_{t}({\bm{x}})-L_{t-1}(\Phi({\bm{x}},\bm{u}))<\epsilon.$ $$

where $\bm{u}\in\mathcal{U}({\bm{x}})$ is $\mu_{t}({\bm{x}})$.

###### Proof.

Recall that the objective function of the optimization problems $W_{\tau}$ and $L_{\tau}$ are:

$$ $\displaystyle\sum_{t=0}^{T-1}\tilde{l}({\bm{x}}(t),\bm{u}(t))$ For $L_{\tau}$ $\displaystyle\sum_{t=0}^{T-1}R({\bm{x}}(t),\bm{u}(t))+\lambda({\bm{x}}(\tau)).$ For $W_{\tau}$ $$

By Lemma  C.3, the rotated cost is $\tilde{l}(\textbf{x},\bm{u})=g^{*}-R({\bm{x}},\bm{u})+\lambda({\bm{x}})-
\lambda(\Phi({\bm{x}},\bm{u}))$. As any valid control of $L_{\tau}$ and $W_{\tau}$ satisfy ${\bm{x}}(t+1)=\Phi({\bm{x}}(t),\bm{u}(t))$, the objective of $L_{\tau}$ can be rewritten as:

$$ $\displaystyle\sum_{t=0}^{\tau-1}\tilde{l}({\bm{x}}(t),\bm{u}(t))$ $\displaystyle=\sum_{t=0}^{\tau-1}g^{*}-R({\bm{x}}(t),\bm{u}(t))+\lambda({\bm{x }}(t))-\lambda(\Phi({\bm{x}}(t),\bm{u}(t)))$ $\displaystyle=\lambda({\bm{x}}(0))+\tau g^{*}-\sum_{t=0}^{T-1}R({\bm{x}}(t), \bm{u}(t))-\lambda({\bm{x}}(\tau)),$ $$

As the last two terms correspond to the objective function for $W_{\tau}$, this shows that $L_{\tau}({\bm{x}})=\tau g^{\star}+\lambda\cdot{\bm{x}}-W_{\tau}{({\bm{x}})}$. Clearly, the two objectives are equivalent and $\mu_{\tau}({\bm{x}})=\bm{u}(0)$. Now note, by Lemma  C.2 we have,

$$ $\lim_{\tau\to\infty}L_{\tau}({\bm{x}})=\lambda\cdot{\bm{x}}+h^{\star}({\bm{x}} )<\infty$ $$

is bounded. Due to dissipativity, $L_{\tau}(\textbf{x})$ is monotone increasing in $\tau$. It follows that

$$ $L_{\tau}({\bm{x}})-L_{\tau-1}(\Phi({\bm{x}},\bm{u}))<\epsilon.$ $$

This completes the proof.
∎

###### Lemma C.4 .

###### Proof.

### C.4 Proof of Theorem 4.1 : Computation details

Here we detail the analysis of $(A)$ and how we go from (16) to (18) in the proof of Theorem 4.1. This term is equal to

$$ $\displaystyle\frac{1}{T}\sum_{t=0}^{T-1}$ $\displaystyle\mathbb{E}\left[l({\bm{X}(t)},\bm{u}(t))\right]=\frac{1}{T}\sum_{ t=0}^{T-1}\mathbb{E}\left[\tilde{l}({\bm{X}}(t),\bm{u}(t))-\lambda\cdot{\bm{X} }(t)+\lambda\cdot\Phi({\bm{X}}(t),\bm{u}(t))\right]$ $\displaystyle=\frac{1}{T}\sum_{t=0}^{T-1}\mathbb{E}\left[L_{\tau}({\bm{X}(t)}) -L_{\tau-1}(\Phi({\bm{X}(t)},\bm{u}(t))-\lambda\cdot{\bm{X}}(t)+\lambda\cdot \Phi({\bm{X}}(t),\bm{u}(t))\right]$ $\displaystyle\leq\frac{1}{T}\sum_{t=0}^{T-1}\mathbb{E}\left[L_{\tau}({\bm{X}(t )})-L_{\tau}(\Phi({\bm{X}(t)},\bm{u}(t)))+\epsilon-\lambda\cdot{\bm{X}}(t)+ \lambda\cdot\Phi({\bm{X}}(t),\bm{u}(t))\right],$ (35) $$

where we use the definition of the rotated cost for the first equality, the dynamic principle for the second line and the identity (17) for the last time.

$$ $\displaystyle\frac{1}{T}\sum_{t=0}^{T-1}[L_{\tau}({\bm{X}(t)})-\lambda\cdot{ \bm{X}(t)}]$ $\displaystyle=\frac{1}{T}\sum_{t=0}^{T-1}[L_{\tau}({\bm{X}(t+1)})-\lambda\cdot {\bm{X}(t+1)}]$ $\displaystyle\quad+\frac{1}{T}\left(L_{\tau}({\bm{X}(0)})-L_{\tau}({\bm{X}(T)} )+\lambda\cdot{\bm{X}(0)}-\lambda\cdot{\bm{X}(T)}\right)$ $$

When $T$ goes to infinity, the second line of the above result goes to 0 0. This shows that when we take the limit as $T$ goes to infinity, (35) is equal to

$$ $\displaystyle\epsilon+\lim_{T\to\infty}\frac{1}{T}\sum_{t=0}^{T-1}\left[L_{ \tau}({\bm{X}}(t+1))-L_{\tau}(\Phi({\bm{X}}(t),\bm{u}(t)))-\lambda\cdot{\bm{X} }(t+1)+\lambda\cdot\Phi({\bm{X}}(t),\bm{u}(t))\right],$ $$

which shows (18)

## Appendix D Proof of Theorem 4.2

###### Lemma D.1 .

Assume 2 and 3, then there exists a neighborhood $\mathcal{N}$ of ${\bm{x}}^{*}$ and a matrix $A$ such that $\mu_{\tau}({\bm{x}})=\mu_{\tau}({\bm{x}}^{*})+A({\bm{x}}-{\bm{x}}^{*})$ for all ${\bm{x}}\in\mathcal{N}$.

###### Proof.

Let us denote by $S^{+}:=\{i:u^{*}_{i}=x^{*}_{i}\}$ the set of states for which all the action $1$ is taken for all arms in those state and by $S-+:=\{i:u^{*}_{i}=x^{*}_{i}\}$ the set of states for which all the action 0 0 is taken for all arms in those state.
Recall that $i^{*}$ is the unique state such that $0<u^{*}_{i}<x^{*}_{i}$ (which exists and is unique by Assumption 2. We define the function $f:\Delta_{\mathcal{S}}\to\mathcal{S}$ by:

$$ $\displaystyle f_{i}(x)=\left\{\begin{array}[]{ll}x_{i}&\text{ For all $i\in \mathcal{S}^{+}$}\\ \alpha-\sum_{i\in\mathcal{S}^{+}}x_{i}&\text{ For $i=i^{*}$}\\ 0&\text{ For all $i\in\mathcal{S}^{-}$.}\end{array}\right.$ $$

We claim that there exists a neighborhood $\mathcal{N}$ such that

- 1.
For all ${\bm{x}}\in\mathcal{N}$, we have $f({\bm{x}})\in\mathcal{U}({\bm{x}})$, which means that $f({\bm{x}})$ is a feasible control for ${\bm{x}}$.
*Proof*. We remark that by construction, one has $\sum_{i}f_{i}({\bm{x}})=\alpha$. Hence, $f({\bm{x}})$ is a valid control if and only if $0\leq\bm{u}\leq{\bm{x}}$. This is clearly true for all $i\neq i^{*}$. For $i^{*}$, it is true if $0\leq\alpha-\sum_{i\mathcal{S}^{+}}x_{i}\leq x^{*}_{i}$ which holds in a neighborhood of ${\bm{x}}$ because of non-degeneracy (Assumption 2) that implies that $0<u^{*}_{i}=\alpha-\sum_{i\mathcal{S}^{+}}x^{*}_{i}<x^{*}_{i^{*}}$.
- 2.
There exists a neighborhood $\mathcal{N^{\prime}}$ of ${\bm{x}}^{*}$ such that if ${\bm{x}}(0)\in\mathcal{N^{\prime}}$, and we construct the sequence ${\bm{x}}(t)$ by setting ${\bm{x}}(t+1)=\Phi({\bm{x}}(t),f({\bm{x}}(t)))$, then ${\bm{x}}(t)\in\mathcal{N}$ for all $t$.
*Proof*. We remark that ${\bm{x}}(t+1)={\bm{x}}^{*}+({\bm{x}}(t)-{\bm{x}}^{*})P^{*}$, where $P^{*}$ is the matrix defined in Assumption 3. Indeed, one has:
$\displaystyle\Phi_{j}({\bm{x}},f({\bm{x}}))$
$\displaystyle=({\bm{x}}P^{0})_{j}+(f({\bm{x}})(P^{1}-P^{0}))_{j}$
$\displaystyle=\sum_{i\in\mathcal{S}^{+}}x_{i}P^{1}_{ij}+\sum_{i\in\mathcal{S}^
{-}}x_{i}P^{0}_{ij}+x_{i^{*}}P^{0}_{i^{*},j}+(\alpha-\sum_{i\in\mathcal{S}^{+}
}x_{i})(P^{1}_{i^{*},j}-P^{0}_{i^{*},j})$
$\displaystyle=\sum_{i\in\mathcal{S}^{-}\cup\{i^{*}\}}x_{i}P^{0}_{ij}+\sum_{i
\in\mathcal{S}^{+}}x_{i}(P^{1}_{i,j}-P^{1}_{i^{*},j}+P^{0}_{i^{*},j})+\alpha(P
^{1}_{i^{*},j}-P^{0}_{i^{*},j})$
$\displaystyle=({\bm{x}}P^{*})_{j}+\alpha(P^{1}_{i^{*},j}-P^{0}_{i^{*},j}),$
where $P^{*}$ is the matrix defined in Assumption 3.
Note that by construction, one has $\Phi({\bm{x}}^{*},f({\bm{x}}^{*}))={\bm{x}}^{*}$. This implies that $\Phi_{j}({\bm{x}},f({\bm{x}}))-\Phi_{j}({\bm{x}}^{*},f({\bm{x}}^{*}))=({\bm{x}
}-{\bm{x}}^{*})P^{*}$. The result follows by Assumption 3 that imposes that the matrix $P^{*}$ is stable.
- 3.
Let $Z({\bm{x}})$ be the reward function collected by the control $f({\bm{x}})$. $Z({\bm{x}})$ is a linear function of ${\bm{x}}$ and $Z({\bm{x}}^{*})=\tau g^{*}+\lambda\cdot x^{*}=W_{\tau}{({\bm{x}}^{*})}$. As a result, $Z({\bm{x}})=W_{\tau}{({\bm{x}})}$ for all ${\bm{x}}\in\mathcal{N^{\prime}}$.
*Proof*. The linearity of $Z$ is a direct consequence of the linearity of $f$ and $\Phi$. The fact that $Z({\bm{x}}^{*})=\tau g^{*}+\lambda\cdot x^{*}$ is because $f({\bm{x}}^{*})=\bm{u}^{*}$. It is also the optimal value of $W_{\tau}({\bm{x}}^{*})$ by dissipativity. Last, as the function $W_{t}()$ is the solution of an LP, it is concave. Since $Z$ is linear and $0<x_{i}<1$ for all state $i$, it follows that the two functions must coincide.

∎

###### Lemma D.2 .

Assume 1, and 4. Assume that $\alpha N$ is an integer. For any neighborhood $\mathcal{N}$ of ${\bm{x}}^{*}$ There exists a $C>0$ such that if ${\bm{X}(t)}$ is a trajectory of the optimal control problem, then:

$$ $\displaystyle\lim_{t\to\infty}\Pr[{\bm{X}(t)}\in\mathcal{N}]\leq e^{-CN}.$ $$

###### Proof.

We use two main steps:

- 1.
For any initial condition ${\bm{x}}(0)$, let ${\bm{x}}_{\tau}(t+1)=\Phi({\bm{x}}_{\tau}(t),\mu_{\tau}({\bm{x}}_{\tau}(t)))$, then there exists $\tau$ and $T$ such that for any initial condition ${\bm{x}}(0)$ and any $t\geq T$, one has ${\bm{x}}_{\tau}(t)\in\mathcal{N}$.
*Proof*. The result follows by dissipativity (which holds because of Assumption 1). Indeed, as ${\bm{x}}_{\tau}(\cdot)$ is an optmal trajectory, it must hold that $\lim_{\tau\to\infty,T\to\infty}\tilde{l}({\bm{x}}_{\tau}(t),\mu_{\tau}({\bm{x}
}_{\tau}(t)))=0$.
- 2.
The map ${\bm{x}}\mapsto\mu_{\tau}({\bm{x}})$ is Lipschitz-continuous in ${\bm{x}}$.
*Proof*. By assumption 4, the control $\mu_{\tau}({\bm{x}})$ is unique. The result then follows because $\mu_{\tau}({\bm{x}})$ is the solution of a linear program parametrized by ${\bm{x}}$.

The lemma can then be proven by adapting the proof of Theorem 3 of Gast et al. () and in particular their Equation (17). This is proven by using Lemma 1 Gast et al. () that implies that $\Pr[\|{\bm{X}(t+1)}-\Phi({\bm{X}(t)},\bm{U}(t))\|\geq\epsilon]\leq e^{-C^{
\prime}N}$.
∎

We are now ready to prove Theorem 4.2. We can follow the proof of Theorem 4.1 up to
(18) that shows that $V_{\text{opt}}^{N}(\textbf{x})-V_{\text{LP},\tau}^{(N)}(\textbf{x})$ is bounded(^4^44In this expression, the term $(\alpha N-\lfloor{\alpha N}\rfloor)/N$ is equal to 0 0 here because we assumed that $\alpha N$ is an integer.) by

$$ $\displaystyle\epsilon+\lim_{\tau\to\infty}\frac{1}{\tau}\sum_{t=0}^{\tau-1} \mathbb{E}\left[L_{\tau}({\bm{X}}(t+1))-L_{\tau}(\Phi({\bm{X}}(t),\bm{u}(t)))- \lambda{\bm{X}}(t+1)+\lambda\Phi({\bm{X}}(t),\bm{u}(t))\right]$ $$

Let $g({\bm{x}})=L_{\tau}({\bm{x}})-\lambda{\bm{x}}$. By Lemma D.1, this function is linear on $\mathcal{N}$. Let us denote by $E(t)$ the event

$$ $\displaystyle E(t):=\{{\bm{X}(t+1)}\in\mathcal{N}\land\Phi({\bm{X}}(t),\bm{u}( t))\in\mathcal{N}\}.$ $$

Hence, this shows that:

$$ $\displaystyle\mathbb{E}$ $\displaystyle\left[L_{\tau}({\bm{X}(t+1)})-L_{\tau}(\Phi({\bm{X}(t)},\bm{u}(t) ))-\lambda{\bm{X}(t+1)}+\lambda\Phi({\bm{X}(t)},\bm{u}(t))\right]$ $\displaystyle=\mathbb{E}\left[g({\bm{X}(t+1)})-g(\Phi({\bm{X}(t)},\bm{u}(t)))\right]$ $\displaystyle=\mathbb{E}\left[(g({\bm{X}(t+1)})-g(\Phi({\bm{X}(t)},\bm{u}(t))) )\mathbf{1}_{E(t)}\right]+\mathbb{E}\left[(g({\bm{X}(t+1)})-g(\Phi({\bm{X}(t)} ,\bm{u}(t))))\mathbf{1}_{\not E(t)}\right]$ $$

By linearity of the $g$ when $E(t)$ is true and by the fact that $\mathbb{E}[\bm{U}(t)]=\bm{u}(t)$ (Lemma B.1), the first term is equal to 0 0. Moreover, as $E(t)$ is true with probability at least $1-2e^{-CN}$, the second term is bounded by $C^{\prime}e^{-CN}$ when $t$ is large.

This concludes the proof of Theorem 4.2.

###### Lemma D.1 .

###### Proof.

###### Lemma D.2 .

###### Proof.

## Appendix E Proof of Lemma C.1

### Proof Outline

This section is dedicated to proving the Lipschitz property for the value function $V^{\infty}({\bm{x}})$ in ${\bm{x}}$ under the ergodicity Assumption 1. For simplicity of exposition we will set $k=1$, although it is not hard to extrapolate the proof for $k>1$ and is left as an exercise for the reader.

The key (rather counter-intuitive) idea behind this proof is to rewrite the problem in terms of an $M$ component vector $\textbf{s}\in\mathcal{S}^{(M)}$. As $M$ tends to infinity we will use the dense nature of rational numbers in the real line and continuity arguments to argue that the proof holds for any ${\bm{x}}\in\Delta_{\mathcal{S}}$.

Hence, to start with, we will assume that ${\bm{x}}\in\Delta_{\mathcal{S}}^{(M)}$, the set of all points on the simplex that can be represented by an $M$ component vector $\mathcal{S}^{(M)}$. In order to complete this proof we will need an intermediate result that verifies the Lipschitz property for all ${\bm{x}},{\bm{x}}^{\prime}\in\Delta_{\mathcal{S}}^{M}$ such that $\|{\bm{x}}-{\bm{x}}^{\prime}\|\leq\frac{2}{M}$ i.e, they differ on at most two components.

### Proof of Lemma

Following the proof outline above let ${\bm{x}},{\bm{x}}^{\prime}\in\Delta_{\mathcal{S}}^{(M)}$. Hence, there exist unique (up to permutation) $M$ component vectors $\textbf{s}({\bm{x}}):=\{s_{0},s_{1}\dots s_{M-1}\}$ and $\textbf{s}^{\prime}({\bm{x}}):=\{s^{\prime}_{0},s^{\prime}_{1}\dots s^{\prime}
_{M-1}\}$.

###### Lemma E.1 .

Under assumption 1, for any $\beta\in(0,1]$, let ${\bm{x}},{\bm{x}}^{\prime}\in\Delta_{\mathcal{S}}^{(M)}$ with $\|{\bm{x}}-{\bm{x}}^{\prime}\|_{1}\leq\frac{2}{M}$, then,

$$ $V^{\infty}({\bm{x}})-V^{\infty}({\bm{x}}^{\prime})\leq\frac{1}{\rho}\|{\bm{x}} -{\bm{x}}^{\prime}\|_{1}$ (36) $$

###### Proof.

Let $\textbf{s}({\bm{x}}):=\{s_{0},s_{1}\dots s_{M-1}\}$, WLOG we can assign $\textbf{s}^{\prime}({\bm{x}}^{\prime}):=\{s_{0}^{\prime},s_{1},s_{2}\dots s_{M
-1}\}$. Note, in this case if $s_{0}$ is the $1^{\text{st}}$ component ($s_{0}=1$) then one can look at the fraction of arms in state $1$, $\textbf{x}(\textbf{s})_{1}:=x_{1}$. Now if $s_{0}^{\prime}=i\neq 1$, we have $x^{\prime}_{1}=x_{1}-\frac{1}{M}$ and $x^{\prime}_{i}=x_{i}+\frac{1}{M}$, allowing us to conclude $\|{\bm{x}}-{\bm{x}}^{\prime}\|_{1}=\frac{2}{M}$. Note, for any ${\bm{x}}$, let a be an action vector i.e. $\textbf{a}=\{a_{0},a_{1}\dots a_{M-1}\}$. Since, $V^{\infty}({\bm{x}})$ is always well defined for any $\beta<1$, one can write the $Q$ function for state s and action a as follows:

$$ $Q(\textbf{s},\textbf{a}):=\frac{1}{M}\sum_{n=0}^{M-1}R^{a_{n}}_{s_{n}}+\beta \sum_{\textbf{s}^{\prime\prime}\in\Delta_{\mathcal{S}}^{(M)}}V^{\infty}( \textbf{x}(\textbf{s}^{\prime\prime}))\Pi_{n=0}^{M-1}P^{a_{n}}_{s_{n},s_{n}^{ \prime\prime}}$ $$

Now, there exists $\textbf{a}^{*}:=\{a_{0},a_{1}\dots a_{M-1}\}$ such that $Q(\textbf{s},\textbf{a}^{*})=V^{\infty}(\textbf{x}(\textbf{s}))$. Suppose we pick $\textbf{a}^{\prime}:=\{0,a_{1},a_{2}\dots a_{M-1}\}$. Note, $\textbf{a}^{\prime}$ always satisfies the constraint on the action space (if needed by pulling one less arm). Hence,

Note, $0\leq R^{a_{0}}_{s_{0}}\leq 1$, so the first term is bounded above by $\frac{1}{M}$. We will focus on the second term, to this end let $\rho(i)=\min\{P^{a_{0}}_{s_{0},i},P^{0}_{s_{0}^{\prime},i}\}$.

Note, by assumption (1), $\rho:=\sum_{i\in\mathcal{S}}\rho(i)>0$. Further, let

$$ $\sigma:=\max_{{\bm{x}},{\bm{x}}^{\prime}\in\Delta_{\mathcal{S}}^{(M)},\|{\bm{x }}-{\bm{x}}^{\prime}\|_{1}\leq 2/M}|V^{\infty}({\bm{x}})-V^{\infty}({\bm{x}}^{ \prime})|$ $$

Let us denote by $\tilde{{\bm{x}}}:=\textbf{x}(\{i,s_{1}^{\prime\prime}\dots s_{M-1}^{\prime
\prime}\})$ and $\tilde{{\bm{x}}}^{\prime}:=\textbf{x}(\{i^{\prime},s_{1}^{\prime\prime}\dots s
_{M-1}^{\prime\prime}\})$ we have: $\|\tilde{{\bm{x}}}-\tilde{{\bm{x}}}^{\prime}\|_{1}\leq 2/M$. This leads us to the following result, for any ${\bm{x}},{\bm{x}}^{\prime}$ such that $\|{\bm{x}}-{\bm{x}}^{\prime}\|_{1}\leq 2/M$ we have:

$$ $\displaystyle|V^{\infty}({\bm{x}})-V^{\infty}({\bm{x}}^{\prime})|\leq\frac{1}{ M}+(1-\rho)$ $$

In particular, choose ${\bm{x}},{\bm{x}}^{\prime}\in\Delta_{\mathcal{S}}^{(M)}$ that maximize the value of $V^{\infty}({\bm{x}})-V^{\infty}({\bm{x}}^{\prime})$, we then have,

$$ $\sigma\leq\frac{1}{M}+(1-\rho)\sigma$ $$

Finally, putting the results together we have for any ${\bm{x}},{\bm{x}}^{\prime}\in\Delta_{\mathcal{S}}^{(M)}$ with $\|{\bm{x}}-{\bm{x}}^{\prime}\|_{1}\leq 2/M$:

$$ $|V^{\infty}({\bm{x}})-V^{\infty}({\bm{x}}^{\prime})|\leq\sigma\leq\frac{1}{2 \rho}\|{\bm{x}}-{\bm{x}}^{\prime}\|_{1}$ $$

∎

We are now ready for the final steps of the proof.

###### Proof of Lemma C.1 .

The result then follows by noting that for any ${\bm{x}},{\bm{x}}^{\prime}\in\Delta_{\mathcal{S}}^{(M)}$ there exists a shortest path from ${\bm{x}}$ to ${\bm{x}}^{\prime}$, say $\{{\bm{x}}={\bm{x}}_{0},{\bm{x}}_{1},{\bm{x}}_{2}\dots,{\bm{x}}_{P}={\bm{x}}^{
\prime}\}$ such that between any two sequential components exactly one component is changed at a time $\|{\bm{x}}_{i}-{\bm{x}}_{i+1}\|=\frac{2}{M}$. Clearly, no more than $2M$ changes can occur along this path, it follows that $P<2M$. Summing along this path allows us to see that for all ${\bm{x}},{\bm{x}}^{\prime}\in\Delta_{\mathcal{S}}^{(M)}$ we have,

$$ $|V^{\infty}({\bm{x}})-V^{\infty}({\bm{x}}^{\prime})|\leq\frac{1}{\rho}\|{\bm{x }}-{\bm{x}}^{\prime}\|_{1}$ $$

We can now complete the proof by noting that $\cup_{M>0}\Delta_{\mathcal{S}}^{(M)}$ is dense in $\Delta_{\mathcal{S}}$. Hence, the result must hold for any ${\bm{x}},{\bm{x}}^{\prime}\in\Delta_{\mathcal{S}}$. This completes the proof.
∎

###### Lemma E.1 .

###### Proof.

###### Proof of Lemma C.1 .

## Appendix F Details about the experiments

### F.1 Generalities

All our simulations are implemented in Python, by using numpy for the random generators and array manipulation and pulp to solve the linear programs. All simulations are run on a personal laptop (macbook pro from 2018). To ensure reproducibility, we will make the code and the Python notebook publicly available (this not done now for double-blind reasons).

#### Value of τ 𝜏 \tau italic_τ

– Except specified otherwise, we use the value $\tau=10$ in all of our example except for the example Chen () where we use $\tau=50$.

#### Number of simulations and confidence intervals

– To obtain a estimate of the steady-state average performance, we simulate the system up to time $T=1000$ and estimate the average performance by computing the average over values from $t=200$ to $T=1000$. All reported confidence intervals correspond to 95% interval on the mean computed by using $\hat{\mu}\pm 2\frac{\hat{\sigma}}{\sqrt{K-1}}$, where $\hat{\mu}$ and $\hat{\sigma}$ are the empirical mean and standard deviation and $k$ is the number of samples (the number of independent trajectories and/or number of randomly generated examples).

#### A note on choice of algorithms for comparison

As specified in the main body of the paper, we choose to restrict our comarison to two algorithms: LP-update and FTVA. We note here that Hong et al. () is theoretically an exponentially optimal solution to the RMAB problem and its absence in our comparison might seem conspicuous to the reader. Hong et al. () is a two set policy which requires the *identification of a largest set* that can be aligned with the fixed point and actively steering the remaining arms towards this set. Theoretically, we know that for a suitable definition such a set always exists for any finite joint state, however, from a practical perspective the problem of choosing a “largest” set that aligns with the fixed point seems unclear to us. Unfortunately, the authors do not provide an implementation for their algorithm and hence, we refrain from comparing our work in order to avoid misleading our readers about the efficacy of their algorithm in comparison to our own as the performance of such an algorithm might heavily depend on the choice of hyperparameter used to implement the alignment procedure. On the contrary, the LP-update does not have any hyperparameters (except for $\tau$ that we set to $\tau=10$ in all our simulations), and so have FTVA and LP-update.

### F.2 Example Hong et al. ( 2023 )

This example has $8$ dimensions and its parameters are:

$$ $\displaystyle P^{0}=\left(\begin{array}[]{cccccccc}1&&&&&&&\\ 1&&&&&&&\\ &0.48&0.52&&&&&\\ &&0.47&0.53&&&&\\ &&&&0.9&0.1&&\\ &&&&&0.9&0.1&\\ &&&&&&0.9&0.1\\ 0.1&&&&&&&0.9\\ \end{array}\right)$ $$

$$ $\displaystyle P^{1}=\left(\begin{array}[]{cccccccc}0.9&0.1&&&&&&\\ &0.9&0.1&&&&&\\ &&0.9&0.1&&&&\\ &&&0.9&0.1&&&\\ &&&0.46&0.54&&&\\ &&&&0.45&0.55&&\\ &&&&&0.44&0.56&\\ &&&&&&0.43&0.57\\ \end{array}\right)$ $$

$R^{1}_{i}=0$ for all state and $R^{0}=(0,0,0,0,0,0,0,0.1)$.

For $\alpha=0.5$ (which is the parameter used in Hong et al. ()), the value of the relaxed LP is $0.0125$ and the LP index are $[0.025,0.025,0.025,0.025,0,-0.113,-0.110,-0.108]$. Note that these numbers differ from the “Lagrangian optimal indices” given in Appendix G.2 of Hong et al. (). This is acknowledge in Hong et al. () when the authors write “In [their] setting, because the optimal solution $y$ remains optimal even without the budget constraint, we can simply remove the budget constraint to get the Lagrangian relaxation […] A nuance is that the optimal Lagrange multiplier for the budget constraint is not unique in this setting, so there can be different Lagrange relaxations”. As a result, the priorities given by their indices and the true LP-index are different. In our simulation, we use the “true” LP-index and not their values. We also tested their values and obtained results that are qualitatively equivalent for the LP-priority (*i.e.*, for their order or ours, the LP-priority essentially gives no reward).

### F.3 Example Chen ( 2022 )

This example has $|\mathcal{S}|=3$ dimensions and its parameters are:

$$ $\displaystyle P^{0}=\left(\begin{array}[]{ccc}0.022&0.102&0.875\\ 0.034&0.172&0.794\\ 0.523&0.455&0.022\\ \end{array}\right)\qquad P^{1}=\left(\begin{array}[]{ccc}0.149&0.304&0.547\\ 0.568&0.411&0.020\\ 0.253&0.273&0.474\\ \end{array}\right)$ $$

$R^{1}=(0.374,0.117,0.079)$ and $R^{0}_{i}=0$ for all $i\in\mathcal{S}$.

For $\alpha=0.4$ (which was the parameter used in Chen ()), the value of the relaxed LP is $0.1238$ and its LP-index are $(0.199,-0.000,-0.133)$.

### F.4 Random example

To generate random example, we use functions from the library numpy of Python:

- •
To generate the transition matrices, We generate an array of size $S\times 2\times S$ by using the function np.random.exponential(size=(S, 2, S)) from the library numpy of python and we then normalize each line so that the sum to $1$.
- •
We generate reward by using np.random.exponential(size=(S, 2))

The example used in Figure 1(a) corresponds to setting the seed of the random generator to $3$ by using np.random.seed(3) before calling the random functions. Its parameters (rounded to $3$ digits) are:

$$ $\displaystyle P^{0}=\left(\begin{array}[]{cccccccc}0.101&0.155&0.043&0.090&0.2 81&0.285&0.017&0.029\\ 0.006&0.207&0.076&0.136&0.085&0.299&0.147&0.043\\ 0.317&0.254&0.065&0.013&0.144&0.111&0.061&0.035\\ 0.098&0.183&0.069&0.068&0.218&0.028&0.200&0.136\\ 0.053&0.080&0.009&0.038&0.483&0.036&0.159&0.143\\ 0.018&0.105&0.027&0.397&0.150&0.102&0.161&0.040\\ 0.110&0.050&0.088&0.024&0.023&0.142&0.169&0.393\\ 0.055&0.043&0.017&0.494&0.227&0.034&0.119&0.011\\ \end{array}\right)$ $$

$R^{0}=(0.073,0.087,0.778,0.186,1.178,0.417,1.996,1.351)$

$$ $\displaystyle P^{1}=\left(\begin{array}[]{cccccccc}0.011&0.124&0.006&0.131&0.2 24&0.070&0.241&0.191\\ 0.071&0.138&0.033&0.023&0.045&0.250&0.339&0.101\\ 0.093&0.113&0.056&0.061&0.109&0.351&0.157&0.059\\ 0.158&0.176&0.151&0.150&0.060&0.142&0.053&0.109\\ 0.370&0.185&0.261&0.020&0.022&0.064&0.047&0.030\\ 0.199&0.139&0.099&0.050&0.141&0.104&0.082&0.187\\ 0.214&0.088&0.011&0.075&0.295&0.174&0.075&0.068\\ 0.028&0.157&0.126&0.078&0.039&0.127&0.376&0.069\\ \end{array}\right)$ $$

$R^{1}=(0.059,3.212,1.817,0.302,2.259,0.067,0.344,0.172)$.

For $\alpha=0.5$, the value of the relaxed problem is $1.3885$ and the LP-index are $0.377,3.273,0.846,-0.116,0.802,,-1.230,-0.562$.

For the other randomly generated examples, we compute the average performance over 50 examples by varying the seed between 0 0 and $49$ for most of the figures except for Figures 4(b) where we only use $20$ examples (and vary the seed between 0 0 and $19$) to improve computation time.

## Peer Reviews

confidence: 3
contribution: 3
initial_rating: 6
presentation: 3
questions: 1.	In Section 6, the algorithm LP-priority is not formally introduced. It is confusing to distinguish between the LP-update and LP-priority algorithms due to the lack of a clear definition.
2.	What are the main technical contributions of the theoretical analysis compared to existing works? I suggest highlighting the novelty and primary contributions of the theoretical analysis more clearly.
3.	In Line 88, the paper states, "It performs well both in terms of the number of arms N as well as the computational time horizon T, beating state-of-the-art algorithms in our benchmarks." However, in the numerical experiments section, the authors did not compare the computational efficiency of the proposed MPC approach with existing algorithms. Moreover, it would be beneficial to provide a more rigorous discussion on why the MPC approach reduces the computational burden.
soundness: 4
strengths: 1.	The paper introduces a novel application of MPC to the RMAB problem which achieve suboptimal gap $O(1/\sqrt{N})$ with a minimal set of assumptions and exponential convergence rate under local stability conditions.
2.	The proposed MPC approach reduces the computational burden associated with solving RMAB problems and perform well in numerical experiments.
3.	This paper presents an interesting framework based on dissipativity, and provides theoretical analysis in this paper.
summary: The paper addresses the discrete time infinite horizon average reward Restless Markovian Bandit (RMAB) problem with a Model Predictive Control (MPC) approach. The proposed MPC algorithm achieve the suboptimality gap $O(1/\sqrt{N})$ with a minimal set of assumptions, and can achieve exponential convergence under local stability conditions. Moreover, the MPC algorithm works well in practice with SOTA computational complexity.
weaknesses: 1.	In Section 6, the algorithm LP-priority is not formally introduced. It is confusing to distinguish between the LP-update and LP-priority algorithms due to the lack of a clear definition.
2.	What are the main technical contributions of the theoretical analysis compared to existing works? I suggest highlighting the novelty and primary contributions of the theoretical analysis more clearly.
3.	In Line 88, the paper states, "It performs well both in terms of the number of arms N as well as the computational time horizon T, beating state-of-the-art algorithms in our benchmarks." However, in the numerical experiments section, the authors did not compare the computational efficiency of the proposed MPC approach with existing algorithms. Moreover, it would be beneficial to provide a more rigorous discussion on why the MPC approach reduces the computational burden.

---

confidence: 4
contribution: 2
initial_rating: 3
presentation: 1
questions: - On one hand, the LP based relaxation has been extensively used in the RMAB literature, such as Verloop 2016, Zhang and Frazier 2022. On the other hand, the randomized rounding procedure is almost the same as that in Gast et al. 2023a and a finite-horizon MPC algorithm (LP-update policy) was proposed in Gast et al. 2023a,b. It is more like a straightforward extension. From the algorithmic perspective, can you more explicitly state what you see as the key novel aspects of the MPC based algorithm compared to prior work, particularly Gast et al. 2023a,b.? 
- The first result in Section 4.1. This is not surprising and it is commonly known in the RMAB literature that LP-based method for RMAB is provably asymptotically optimal. Indeed, the LP based method has been leveraged to design index policy for RMAB problem without the hard-to-verify indexability condition as required by the Whittle index policy, and such a LP-based method to design index policy is provably asymptotically optimal as shown in the literature, e.g., Verloop 2016 is one of the first works in this domain.  Can you discuss how your result in Section 4.1 advances the state-of-the-art beyond what was already known from works like Verloop 2016? Are there any aspects of your analysis or bounds that are novel or improved? 
- The second result in Section 4.2. Likewise, the proof is directly from Gast et al. 2023a and Hong et al. 2024a. For both results, can your clarify exactly how the use of dissipativity differs from or improves upon previous approaches? can you discuss any limitations of previous methods that your approach overcomes? 
- In practice, how to determine how many arms to pull at each time, given that $\alpha N$ may not be an integer? If it is not an integer (may consider as an average constraint, there is no need to design an index policy). 
-Equation (3) itself is not a RMAB problem. It should be properly defined with the budget constraint to be satisfied at teach time. It may be better to rigorously define (3)  
- Many typos in the paper, just to name a few here: line 120, "the budget constraint, $\alpha$", line 146, $\boldsymbol{x}_s$ is not defined. line 144, since $u(s,a)$ is denoted as the empirical distribution of the state-action pairs $(s,1)$, why not just express it as $u(s,1)$ since the action is fixed.
- In Section 3.1 after introducing the optimization problem in (8), the authors claimed that this problem is computationally easy to solve. On one hand, indeed this is a linear problem and it is "relatively" easy to solve when the state space and the number of arms is small, given many LP solvers. On the other hand, this claim is not precise, since solving a LP with large-scale parameters/spaces can still be very computationally expensive and take a lot of time in practice. This may be one of the limitations of LP based method for RMAB compared to the Whittle index policy although LP based methods do not require the indexability condition. 
- Can you elaborate why the definition of (8) imposes the constraint to be satisfied for each time $t$ as claimed in lines 198-199?
soundness: 2
strengths: - RMAB has been extensively studied in recent years, from both the offline and online settings. This paper investigates the fundamental property of RMAB (e.g., the asymptotic optimality, optimality gap) in the offline setting. It is a challenging and interesting problem. 
- A model predictive control (MPC) based non-stationary policy has been developed and theoretically analyzed, i.e., its optimality gap is characterized in terms of the number of arms $N$. 
- Some experimental results were presented to validate the performance of this MPC based policy and its comparison with baselines.
summary: This paper studied the discrete time infinite horizon average reward restless markovian bandit (RMAB) problem and focused on the asymptotic optimality in this setting. In particular, a model predictive control (MPC) based non-stationary policy with a rolling computational horizon $\tau$ is proposed and its sub-optimality gap is presented. The performance of this policy is also evaluated via simulations.
weaknesses: - The asymptotic optimality performance of RMAB has raised many attentions in recent years. Despite that this paper proposed a policy by leveraging the MPC idea, it is hard to identify the technical novelty from the perspectives of algorithm design and technical proofs and results (See questions below). This paper heavily relies on the previous works such as Gast et al. 2023 a,b. 
- The simulations are rather weak in both the settings and the baseline methods considered. 
- This paper in general is poorly written with many typos, broken sentences and abused notations.

---

confidence: 4
contribution: 2
initial_rating: 5
presentation: 3
questions: Please refer to **Weaknesses**.
soundness: 3
strengths: The paper presents a novel use of dissipativity for analyzing RMABs, which offers fresh insights into this field. The suboptimality bounds are rigorously derived, showcasing that the MPC-based policy approaches optimality as the number of arms, 𝑁, increases.
summary: This paper addresses the infinite-horizon average reward RMAB problem. It proposed a model predictive control (MPC) policy using a rolling computational horizon of length $\tau$, which achieves a suboptimality gap with order $\mathcal{O}(1/\sqrt{N})$ with $N$ being the number of arms.
weaknesses: Though this paper conveys a solid algorithm design and theoretical proof, we must admit that the current submission has significant weaknesses.  

1. My main concern is that the order of $\mathcal{O}(1/\sqrt{N})$ gap has been a well-known result for a long time, which various types of algorithms can achieve. In particular, the LP-based algorithms in many related works cited by this submission without an additional MPC layer can also achieve this optimality gap.  Though the reviewer admits that the proposed MPC layer can bring some advantages, the key reason that this reviewer does not champion this paper is that there is a new work (https://arxiv.org/pdf/2410.15003) recently has pushed the gap to the order of $\mathcal{O}(1/N)$ by using a diffusion approximation technique. I understand that the authors submitted their work earlier than this recent work and may not be aware of this new work, the main concern always holds as the proposed MPC-based algorithm does not improve the well-known optimality gap.

2.  Resolving the LP at each time step is not really a novel idea, which can be found in multiple works in related work cited by this submission. Though controlling $\tau$ is new, resolving LP at each time step causes a significant amount of computational complexity.  Aligned with my first concern, why do we need such a complex algorithm that does not even improve the theoretical guarantee, which is claimed as a main contribution in this work?

3. In many real-world applications that can be modeled as an RMAB framework, the underlying MDP for each arm may not be known in advance. My question is how can we leverage the proposed MPC-based algorithm for such a setting? 

4. This paper considers a homogeneous arm setting. The reviewer quite doubts the scalability of the proposed algorithm when arms are heterogeneous and the number of arms is very large.

---

confidence: 3
contribution: 3
initial_rating: 8
presentation: 3
questions: - Are there any comparison on the running time of different algorithms, and your algorithm with different $\tau$? 

- Are there any experiment on real-data?

- In which part of the proof shows that $\tau(\epsilon) = O({1\over \epsilon})$?


=========After Rebuttal=========

Thanks for your reply, I do not have other questions.
soundness: 3
strengths: - The analysis shows a good suboptimality gap.

- The algorithm is quite efficient to implement.

- The experiment results are also convincing.
summary: This paper proposes a new way to solve restless bandit: putting both the state information and the action information into a continuous vector space, and do linear programming to maximize the reward in the next $\tau$ rounds. When we look forward $\tau$ steps with a large enough $\tau$, this algorithm can find near optimal solutions. The authors show that normally the suboptimality gap is $O(\sqrt{1/N})$, but with certain conditions, it can be reduced to be $O(\exp(-cN))$, where $N$ is the number of arms and $c$ is some constant. As for the gap corresponding to $\tau$, they also show it is about $O({1\over \tau})$. In experiments, it is shown that choosing $\tau$ as a small value (e.g., 10) leads to very good performance, indicating the efficiency of model predictive control.
weaknesses: - The author assumes that all the arms are statistically identical. Is this a common assumption in restless bandits? I think there are many cases that the arms are not identical. It seems that your algorithm cannot be adapted to this setting easily (e.g., if all the arms have distinct transitions and rewards)?

Meta-review: This is a borderline paper. Overall, the reviewers are critical of various aspects of the paper. Most notably about the somewhat incremental novelty and computational scalability issues. I therefore believe that a major revision of the paper is necessary. The reviews contain a variety of suggests of how to improve and revise the paper.

Decision: Reject