Proof. We first consider what happens when Algorithm 2 is used. We proceed by induction on t. This is clearly true for t = 0. Assume that this holds for some t ≥ 0. By Proposition 1, there exists ε_t such that the control y(x) defined in Proposition 1 is optimal for all x ∈ B(x*(t), ε_t). The induction hypothesis and the continuity of y(x) therefore imply that for all ε > 0, there exists C_1, C_2 > 0 such that ||Y^{(N)}(t) - y*(t)||_1 ≤ ε with probability at least 1 - C_1 · e^{-C_2N}, where Y^{(N)}(t) is the control used for X^{(N)}(t). Write

$$ \|\mathbf{X}^{(N)}(t+1) - \mathbf{x}^*(t+1)\|_1 \le \|\mathbf{X}^{(N)}(t+1) - \phi(\mathbf{Y}^{(N)}(t))\|_1 + \|\phi(\mathbf{Y}^{(N)}(t)) - \mathbf{x}^*(t+1)\|_1. $$

Hence, by using Lemma 1 and the union bound, for all ε > 0, there exists C'_1, C'_2 > 0 such that ||X^{(N)}(t+1) - x^*(t+1)||_1 ≤ ε with probability at least 1 - C'_1 · e^{C'_2N}.

To show (28), the only remaining point is to see that if s' is a state such that x_{s'}^*(t+1) = 0, then so is X_{s'}^{(N)}(t+1). By definition of the deterministic evolution x^*(t+1) = φ(y^*(t)) (see e.g. Equation (20)), we have:

$$ x_{s'}^*(t+1) = \sum_{a,s} y_{s,a}^*(t) P_{s,s'}^a. $$

Hence, x_{s'}^*(t+1) equals 0 if for all s, a, either y_{s,a}^*(t) = 0 or P_{s,s'}^a = 0. By construction of Y^{(N)}(t) from y^*(t), y_{s,a}^*(t) = 0 implies that Y_{s,a}^{(N)}(t) = 0, consequently X_{s'}^{(N)}(t+1) = 0.

To study what happens when Algorithm 1 is used, we remark that if the solution to the LP is unique, and X^{(N)}(t) is close enough to x^*(t), then the new LP solution computed by Algorithm 1 is the same y(X^{(N)}(t)) used by Algorithm 2. Hence, the proof for Algorithm 2 also applies in this case. Q.E.D.

### A.4. Proof of the lower bounds

Consider the following two-action restless bandit with two states $S = \{1, 2\}$, with parameters $\mathbf{x}(0) = (0.5, 0.5)$, $T=2$, $A = \{0, 1\}$, $\mathbf{r}^0 = [0, 0]$, $\mathbf{r}^1 = [1, 0]$, and where the transition matrices are

$$ \mathbf{P}^0 = \begin{pmatrix} 0.5 & 0.5 \\ 0.5 & 0.5 \end{pmatrix}, \mathbf{P}^1 = \begin{pmatrix} 0.5 & 0.5 \\ 0.5 & 0.5 \end{pmatrix}. $$

We consider that $D(s, 0) = 0$, $D(s, 1) = 1$ for any state s, and we distinguish the resource constraints $b=0.3$ and $b=0.5$.

For any resource constraint, the solution to the LP is to choose action 1 for as many arms as possible. This gives a reward 2b. If $b=0.3$, the problem is non-degenerate whereas if $b=0.5$ the problem is degenerate.

For the stochastic system with N components, this gives a reward $\frac{1}{N} \lfloor Nb \rfloor$ at time-step 0 and a reward $\min(b, \frac{1}{N} \lfloor NX_1^{(N)}(1) \rfloor)$ at time-step 1. Since $X_1^{(N)}(1)$ follows a binomial distribution of parameter $(N, 0.5)$, the total reward of the LP-update policy is equal to

$$ \frac{1}{N} \lfloor Nb \rfloor + b + \mathbb{E} \left[ \min \left( \frac{1}{N} \lfloor NX_1^{(N)}(1) \rfloor - b, 0 \right) \right]. $$

By the central limit theorem, $\lim_{N \to \infty} \sqrt{N} \mathbb{E} \left[ \min(X_1^{(N)}(1) - 0.5, 0) \right] = \sqrt{2\pi} > 0$. This provides a counter-example for the lower bound of Theorem 1.

If $b=0.3$ and N is not a multiple of 10, the problem does not admit a perfect rounding, and we have $2b - (\frac{1}{N} \lfloor Nb \rfloor + b + \mathbb{E} \left[ \min \left( \frac{1}{N} \lfloor NX_1^{(N)}(1) \rfloor - b, 0 \right) \right]) \ge b - \frac{1}{N} \lfloor Nb \rfloor \ge 0.1/N$. This provides a counter-example for the lower-bound of Theorem 3.

If $b=0.3$ and N is a multiple of 10, then the problem admits a perfect rounding. In this case, classical anti-concentration arguments (Matoušek and Vondrák [23]) show that

$$ \mathbb{P} \left( X_1^{(N)}(1) \le 0.2 \right) \ge \frac{1}{15} \exp(-16N(0.6 - 0.2)^2). $$

This shows that

$$ \mathbb{E} \left[ \min \left( \frac{1}{N} \lfloor NX_1^{(N)}(1) \rfloor - 0.3, 0 \right) \right] \ge -0.1 \mathbb{P} \left( X_1^{(N)}(1) \le 0.2 \right) = -\frac{1}{150} \exp(-2.56N). $$

It provides a counter-example for the lower-bound of Theorem 4. Q.E.D.
**Algorithm 3:** Occupation measure policy for weakly-coupled MDPs.

**Input:** Time horizon $T$ and initial configuration vector $\mathbf{x}(0)$.

Solve the linear program (5) with time horizon $T$ and initial configuration vector $\mathbf{x}(0)$, obtain an optimal solution $\mathbf{y}^*$ and the corresponding $\mathbf{x}^*$;

for $t = 0, 1, 2, \dots, T-1$ do

Compute from the LP solution the occupation measure

$$ \mu_{s,a}^*(t) := \begin{cases} \frac{y_{s,a}^*(t)}{x_s^*(t)}, & \text{if } x_s^*(t) > 0 \\ \mathbf{1}_{\{a=0\}}, & \text{otherwise.} \end{cases} \qquad (29) $$

Observe the current states of the $N$ sub-MDP's $\mathbf{s} = (s_1, s_2, \dots, s_N)$. Initialize $\mathbf{B} := N \cdot \mathbf{b}$, and actions on the $N$ sub-MDP's as $A(\mathbf{s}) = (0, 0, \dots, 0)$;

for $n = 1, 2, \dots, N$ do

Sample an action $a_n$ according to the probability vector $(\mu_{s_n,a}^*(t))_{a \in A(s_n)}$;

if $\mathbf{B} - D(s_n, a_n) \ge 0$ (component-wise) then

$A(s_n) := a_n$;

$\mathbf{B} := \mathbf{B} - D(s_n, a_n)$;

end

end

Apply the actions $A(\mathbf{s})$ to the $N$ sub-MDP's;

end

## Appendix B: THE OCCUPATION MEASURE POLICY

To provide a benchmark for evaluation of the performance of the LP-update policy, we describe in this section a policy that we call the *occupation measure policy*. This policy is a straightforward generalization of the randomized activation control policy in Zayas-Cabán et al. [39] and the occupancy-measured-reward index policy in Xiong et al. [38] to the multi-constraint multi-action case. The occupation measure policy is proven to be asymptotically optimal for the multi-action single-constraint multi-armed bandits in Xiong et al. [38]. Adapting their proof to the multi-constraint case is direct.

The occupation measure policy is a one-pass policy that solves the linear program only once at the very beginning, and constructs the probability vectors $\mu_{s,a}^*(t)$ defined in (29) from the solution. At each decision epoch $t$, the budget left is initialized as $\mathbf{B} := N \cdot \mathbf{b}$, the action on each arm $n$ is initialized as the passive 0. It then samples a new action $a_n$ from the distribution $(\mu_{s_n,a}^*(t))_{a \in A(s_n)}$ on each arm $n$. If choosing action $a_n$ instead of action 0 on arm $n$ does not violate any of the budget constraints, then we apply action $a_n$ on arm $n$, and we decrease the budget left $\mathbf{B}$; otherwise we keep action 0 on arm $n$ and continue to sample an action on the next arm. The detailed implementation is given in Algorithm 3.

## Appendix C: COMPARISON WITH EXISTING NOTIONS OF NON-DEGENERACY

The non-degenerate property as defined in Definition 1 uses a rank condition from the optimal LP solution $\mathbf{y}^*$. In this section, we provide a detailed comparison of our notion to the existing notions presented [41, 14, 6, 40]. To facilitate the comparison, we reformulate our definition.

Let $\mathbf{y}^*$ be an optimal solution to the LP (at time 0) and let us fix a time $t$. For each state $s \in S$, we define by $A_s^+ := \{a \in A | y_{s,a}^*(t) > 0\}$ the number of different actions that are chosen in state $s$ at time $t$. Similar to Brown and Zhang [6], we define:

$$ S^+ := \{s \in S | A_s^+ = 1\}; $$

$$ S^0 := \{s \in S | A_s^+ > 1\}; $$

$$ S^- := \{s \in S | A_s^+ = 0\}. $$

We see that $S^+$, $S^0$, $S^-$ form a partition of $S$. Their interpretations are that (for the LP), a state in $S^+$ has a unique optimal action, a state in $S^0$ has at least two equally optimal actions, and a state in $S^0$ is not reached at time $t$.
**Algorithm 3:** Occupation measure policy for weakly-coupled MDPs.

**Input:** Time horizon $T$ and initial configuration vector $\mathbf{x}(0)$.

Solve the linear program (5) with time horizon $T$ and initial configuration vector $\mathbf{x}(0)$, obtain an optimal solution $\mathbf{y}^*$ and the corresponding $\mathbf{x}^*$;

for $t = 0, 1, 2, \dots, T-1$ do

Compute from the LP solution the occupation measure

$$ \mu_{s,a}^*(t) := \begin{cases} \frac{y_{s,a}^*(t)}{x_s^*(t)}, & \text{if } x_s^*(t) > 0 \\ \mathbf{1}_{\{a=0\}}, & \text{otherwise.} \end{cases} \qquad (29) $$

Observe the current states of the $N$ sub-MDP's $\mathbf{s} = (s_1, s_2, \dots, s_N)$. Initialize $\mathbf{B} := N \cdot \mathbf{b}$, and actions on the $N$ sub-MDP's as $A(\mathbf{s}) = (0, 0, \dots, 0)$;

for $n = 1, 2, \dots, N$ do

Sample an action $a_n$ according to the probability vector $(\mu_{s_n,a}^*(t))_{a \in A(s_n)}$;

if $\mathbf{B} - D(s_n, a_n) \ge \mathbf{0}$ (component-wise) then

$A(s_n) := a_n$;

$\mathbf{B} := \mathbf{B} - D(s_n, a_n)$;

end

end

Apply the actions $A(\mathbf{s})$ to the $N$ sub-MDP's;

end

## Appendix B: THE OCCUPATION MEASURE POLICY

To provide a benchmark for evaluation of the performance of the LP-update policy, we describe in this section a policy that we call the *occupation measure policy*. This policy is a straightforward generalization of the randomized activation control policy in Zayas-Cabán et al. [39] and the occupancy-measured-reward index policy in Xiong et al. [38] to the multi-constraint multi-action case. The occupation measure policy is proven to be asymptotically optimal for the multi-action single-constraint multi-armed bandits in Xiong et al. [38]. Adapting their proof to the multi-constraint case is direct.

The occupation measure policy is a one-pass policy that solves the linear program only once at the very beginning, and constructs the probability vectors $\mu_{s,a}^*(t)$ defined in (29) from the solution. At each decision epoch $t$, the budget left is initialized as $\mathbf{B} := N \cdot \mathbf{b}$, the action on each arm $n$ is initialized as the passive 0. It then samples a new action $a_n$ from the distribution $(\mu_{s_n,a}^*(t))_{a \in A(s_n)}$ on each arm $n$. If choosing action $a_n$ instead of action 0 on arm $n$ does not violate any of the budget constraints, then we apply action $a_n$ on arm $n$, and we decrease the budget left $\mathbf{B}$; otherwise we keep action 0 on arm $n$ and continue to sample an action on the next arm. The detailed implementation is given in Algorithm 3.

## Appendix C: COMPARISON WITH EXISTING NOTIONS OF NON-DEGENERACY

The non-degenerate property as defined in Definition 1 uses a rank condition from the optimal LP solution $\mathbf{y}^*$. In this section, we provide a detailed comparison of our notion to the existing notions presented [41, 14, 6, 40]. To facilitate the comparison, we reformulate our definition.

Let $\mathbf{y}^*$ be an optimal solution to the LP (at time 0) and let us fix a time $t$. For each state $s \in S$, we define by $A_s^+ := \{a \in A | y_{s,a}^*(t) > 0\}$ the number of different actions that are chosen in state $s$ at time $t$. Similar to Brown and Zhang [6], we define:

$$ S^+ := \{s \in S | A_s^+ = 1\}; $$

$$ S^0 := \{s \in S | A_s^+ > 1\}; $$

$$ S^- := \{s \in S | A_s^+ = 0\}. $$

We see that $S^+$, $S^0$, $S^-$ form a partition of $S$. Their interpretations are that (for the LP), a state in $S^+$ has a unique optimal action, a state in $S^0$ has at least two equally optimal actions, and a state in $S^0$ is not reached at time $t$.
As $\mathcal{U}^*(t)$ is the set of indices $(s, a)$ for which $y_{s,a}^*(t) = 0$, we have:

$$ 
\begin{aligned}
|\mathcal{U}^*| &= \sum_{s \in \mathcal{S}} |\mathcal{A}| - A_s^+ \\
&= \sum_{s \in \mathcal{S}^-} A_s + \sum_{s \in \mathcal{S}^+} A_s - |\mathcal{S}^+| + \sum_{s \in \mathcal{S}^0} |\mathcal{A}| - \sum_{s \in \mathcal{S}^0} A_s^+ \\
&= |\mathcal{U}| - |\mathcal{S}^+| - \sum_{s \in \mathcal{S}^0} A_s^+,
\end{aligned}
$$

while $|\mathcal{S}^*| = |\mathcal{S}^0| + |\mathcal{S}^+|$.

As $|\mathcal{J}^*|$ is the number of saturated constraints, we have

$$ |\mathcal{J}^*| + |\mathcal{S}^*| + |\mathcal{U}^*| = |\mathcal{U}| + |\mathcal{J}^*| + |\mathcal{S}^0| - \sum_{s \in \mathcal{S}^0} A_s^+ . $$

The matrix $C^*$ has $|\mathcal{J}^*| + |\mathcal{S}^*| + |\mathcal{U}^*|$ rows and $|\mathcal{U}|$ columns. To simplify the discussion, in this section, we assume that:

$$ \text{All rows of } C^*(t) \text{ are linearly independent.} \qquad (30) $$

So that the rank of $C^*(t)$ boils down to verifying an inequality among dimension. Indeed, under this assumption (30), in order for $C^*$ to have rank $|\mathcal{J}^*| + |\mathcal{S}^*| + |\mathcal{U}^*|$ as in Definition 1, it is necessary and sufficient that

$$ |\mathcal{J}^*| + |\mathcal{S}^0| \leq \sum_{s \in \mathcal{S}^0} A_s^+. \qquad (31) $$

We are now ready to compare the formulation (31) with various existing notions of non-degeneracy in the literature:

1. In Zhang and Frazier [41] and Gast et al. [14] that study two-action bandits with a single equality resource constraint $b = \alpha$. Their notion of non-degeneracy using our notation is to say that $|\mathcal{S}^0| \geq 1$. Since in their model the assumption (30) always holds, and $|\mathcal{J}^*|$ is always 1. In order for (31) to hold, it is necessary and sufficient that $|\mathcal{S}^0| \geq 1$, which shows that the two definitions of non-degeneracy are indeed equivalent.

2. In Brown and Zhang [6], the authors define another notion of non-degeneracy in Definition 5.2 for multi-action single-constraint bandits, which by using the notations in this section can be stated as follows (assuming (30)): the problem is non-degenerate if at least one of the conditions hold: (i) the single resource constraint is not saturated; (ii) $|\mathcal{S}^0| \geq 1$. This is again equivalent to (31). Note that additional care is taken in Definition 5.2 of Brown and Zhang [6], in order to account for the non-satisfaction of (30).

3. In Zhang [40], another definition of non-degeneracy for multi-action multi-constraint bandits is given in Section 3.3. According to this definition, a multi-action multi-constraint bandit is non-degenerate if $|\mathcal{J}^*| \leq 1$ and $|\mathcal{S}^0| \geq 1$. Clearly this implies (31), but under this general model, there are other possible problems that satisfy (31) while not covered by the non-degenerate definition of Zhang [40].

## Appendix D: IMPROVING THE ROUNDING METHOD

When the LP is solved with a given initial condition $\mathbf{X}^{(N)}$, it yields a first-step control $\mathbf{Y} := \mathbf{y}(0) \in \mathcal{Y}(\mathbf{X}^{(N)})$. However, this solution may not be directly transferable to the original problem because the quantities $N \cdot Y_{s,a}$ are not necessarily integers. Throughout our paper, we address this issue by employing a simple truncation method to obtain a vector $\mathbf{Y}^{(N)}$ that is defined as follows:

$$ \mathbf{Y}_{s,a}^{(N)} = \begin{cases} N^{-1} \lfloor N Y_{s,a} \rfloor & \text{if } a \neq 0 \\ X_s^{(N)} - \sum_{a \neq 0} \tilde{Y}_{s,a} & \text{if } a = 0. \end{cases} \qquad (32) $$

This approach results in an admissible solution $\mathbf{Y}^{(N)} \in \mathcal{Y}^{(N)}(\mathbf{X}^{(N)})$, as all elements in $D(s, a)$ and $\mathbf{b}$ are non-negative, and notably, $D(s, 0) = 0$. However, this truncation solution introduces an $\mathcal{O}(1/N)$ sub-optimality gap.

In this section, we delve into potential strategies to enhance performance and reduce this sub-optimality gap. By examining alternative methods to handle the non-integer nature of the solution, we aim to improve the applicability and efficiency of our approach, ensuring closer alignment with the optimal solution of the original problem.
As $\mathcal{U}^*(t)$ is the set of indices $(s, a)$ for which $y_{s,a}^*(t) = 0$, we have:

$$ 
\begin{aligned}
|\mathcal{U}^*| &= \sum_{s \in S} |\mathcal{A}| - A_s^+ \\
&= \sum_{s \in S^-} A_s + \sum_{s \in S^+} A_s - |\mathcal{S}^+| + \sum_{s \in S^0} |\mathcal{A}| - \sum_{s \in S^0} A_s^+ \\
&= |\mathcal{U}| - |\mathcal{S}^+| - \sum_{s \in S^0} A_s^+,
\end{aligned}
$$

while $|\mathcal{S}^*| = |\mathcal{S}^0| + |\mathcal{S}^+|$.

As $|\mathcal{J}^*|$ is the number of saturated constraints, we have

$$ |\mathcal{J}^*| + |\mathcal{S}^*| + |\mathcal{U}^*| = |\mathcal{U}| + |\mathcal{J}^*| + |\mathcal{S}^0| - \sum_{s \in S^0} A_s^+ . $$

The matrix $C^*$ has $|\mathcal{J}^*| + |\mathcal{S}^*| + |\mathcal{U}^*|$ rows and $|\mathcal{U}|$ columns. To simplify the discussion, in this section, we assume that:

All rows of $C^*(t)$ are linearly independent. (30)

So that the rank of $C^*(t)$ boils down to verifying an inequality among dimension. Indeed, under this assumption (30), in order for $C^*$ to have rank $|\mathcal{J}^*| + |\mathcal{S}^*| + |\mathcal{U}^*|$ as in Definition 1, it is necessary and sufficient that

$$ |\mathcal{J}^*| + |\mathcal{S}^0| \leq \sum_{s \in S^0} A_s^+. \qquad (31) $$

We are now ready to compare the formulation (31) with various existing notions of non-degeneracy in the literature:

1. In Zhang and Frazier [41] and Gast et al. [14] that study two-action bandits with a single equality resource constraint $b = \alpha$. Their notion of non-degeneracy using our notation is to say that $|\mathcal{S}^0| \geq 1$. Since in their model the assumption (30) always holds, and $|\mathcal{J}^*|$ is always 1. In order for (31) to hold, it is necessary and sufficient that $|\mathcal{S}^0| \geq 1$, which shows that the two definitions of non-degeneracy are indeed equivalent.

2. In Brown and Zhang [6], the authors define another notion of non-degeneracy in Definition 5.2 for multi-action single-constraint bandits, which by using the notations in this section can be stated as follows (assuming (30)): the problem is non-degenerate if at least one of the conditions hold: (i) the single resource constraint is not saturated; (ii) $|\mathcal{S}^0| \geq 1$. This is again equivalent to (31). Note that additional care is taken in Definition 5.2 of Brown and Zhang [6], in order to account for the non-satisfaction of (30).

3. In Zhang [40], another definition of non-degeneracy for multi-action multi-constraint bandits is given in Section 3.3. According to this definition, a multi-action multi-constraint bandit is non-degenerate if $|\mathcal{J}^*| \leq 1$ and $|\mathcal{S}^0| \geq 1$. Clearly this implies (31), but under this general model, there are other possible problems that satisfy (31) while not covered by the non-degenerate definition of Zhang [40].

## Appendix D: IMPROVING THE ROUNDING METHOD

When the LP is solved with a given initial condition $\mathbf{X}^{(N)}$, it yields a first-step control $\mathbf{Y} := \mathbf{y}(0) \in \mathcal{Y}(\mathbf{X}^{(N)})$. However, this solution may not be directly transferable to the original problem because the quantities $N \cdot Y_{s,a}$ are not necessarily integers. Throughout our paper, we address this issue by employing a simple truncation method to obtain a vector $\mathbf{Y}^{(N)}$ that is defined as follows:

$$ \mathbf{Y}_{s,a}^{(N)} = \begin{cases} N^{-1} \lfloor NY_{s,a} \rfloor & \text{if } a \neq 0 \\ X_s^{(N)} - \sum_{a \neq 0} \tilde{Y}_{s,a} & \text{if } a = 0. \end{cases} \qquad (32) $$

This approach results in an admissible solution $\mathbf{Y}^{(N)} \in \mathcal{Y}^{(N)}(\mathbf{X}^{(N)})$, as all elements in $D(s, a)$ and $\mathbf{b}$ are non-negative, and notably, $D(s, 0) = 0$. However, this truncation solution introduces an $\mathcal{O}(1/N)$ sub-optimality gap.

In this section, we delve into potential strategies to enhance performance and reduce this sub-optimality gap. By examining alternative methods to handle the non-integer nature of the solution, we aim to improve the applicability and efficiency of our approach, ensuring closer alignment with the optimal solution of the original problem.
## D.1. Improvement by using MILP

As previously mentioned, the rounding process we employ in (32) addresses the issue that the first step $y(0)$ from the LP might not yield integer values for $N \cdot y_{s,a}(0)$. To circumvent this, another more refined approach is to explicitly incorporate integer constraints into our LP, leading to the following mixed-integer optimization problem (MILP):

$$ \max_{y \ge 0} \sum_{t=0}^{T-1} \sum_{(s,a) \in \mathcal{U}} r_s^a y_{s,a}(t) \qquad (33a) $$

$$ \text{s.t.} \quad \sum_{a \in \mathcal{A}} y_{s,a}(0) = x_s \quad \forall s, \qquad (33b) $$

$$ \sum_{a \in \mathcal{A}} y_{s,a}(t+1) = (\phi(y(t)))_s \quad \forall s, t, \qquad (33c) $$

$$ Dy(t) \le b \quad \forall t, \qquad (33d) $$

$$ N \cdot y_{s,a}(0) \in \mathbb{N} \quad \forall s, a, \qquad (33e) $$

It is important to note that since the solution primarily utilizes control at the first time-step, we apply the integrity constraint solely at $t = 0$.

We introduce $V_{\text{MILP-update}}^{(N)}(\mathbf{x}, T)$ to denote the expected performance of the policy obtained by implementing Algorithm 1, wherein the solution of the mixed-integer linear program (33) is utilized instead of (5). While this modified problem is more constrained than the original LP, it still represents a relaxation of the initial problem. Therefore, for any admissible policy $\pi \in \text{LP-update}$, MILP-update, the following holds:

$$ V_{\pi}^{(N)}(\mathbf{x}, T) \le V_{\text{rel}}(\mathbf{x}, T) $$

In Figure 5, we present a performance comparison between the LP-update and MILP-update policies. As anticipated, the MILP-update demonstrates superior performance compared to the original LP-update, particularly for smaller values of $N$.

<figure><img src="image_7.png" /><figcaption>Figure 5 Sub-optimality gap of the LP-update policy compared to the MILP-update policy.</figcaption></figure>

From a computational perspective, addressing (33) is more demanding in terms of resources compared to solving the original LP problem. Our experimental findings revealed that solving the MILP required no more than twice the time needed for the original LP for this specific model having 10 states. Notably, the duration for solving the MILP appeared to be unaffected by the size of $N$, suggesting that increases in $N$ do not proportionally complicate the MILP solution process. This consistency in computational time, regardless of the scale of $N$, offers an interesting topic for further investigation, especially to determine if this trend is consistent across various scenarios and applications. It is also important to observe that the zigzag effect in the LP-update policy, resulting from the cutoff rounding, is absent in the MILP-update policy. This suggests that the MILP-update policy might be a more suitable and natural choice for implementation.
## D.1. Improvement by using MILP

As previously mentioned, the rounding process we employ in (32) addresses the issue that the first step $y(0)$ from the LP might not yield integer values for $N \cdot y_{s,a}(0)$. To circumvent this, another more refined approach is to explicitly incorporate integer constraints into our LP, leading to the following mixed-integer optimization problem (MILP):

$$ \max_{y \ge 0} \sum_{t=0}^{T-1} \sum_{(s,a) \in \mathcal{U}} r_s^a y_{s,a}(t) \qquad (33a) $$

$$ \text{s.t.} \quad \sum_{a \in \mathcal{A}} y_{s,a}(0) = x_s \quad \forall s, \qquad (33b) $$

$$ \sum_{a \in \mathcal{A}} y_{s,a}(t+1) = (\phi(y(t)))_s \quad \forall s, t, \qquad (33c) $$

$$ Dy(t) \le b \quad \forall t, \qquad (33d) $$

$$ N \cdot y_{s,a}(0) \in \mathbb{N} \quad \forall s, a, \qquad (33e) $$

It is important to note that since the solution primarily utilizes control at the first time-step, we apply the integrity constraint solely at $t=0$.

We introduce $V_{\text{MILP-update}}^{(N)}(\mathbf{x}, T)$ to denote the expected performance of the policy obtained by implementing Algorithm 1, wherein the solution of the mixed-integer linear program (33) is utilized instead of (5). While this modified problem is more constrained than the original LP, it still represents a relaxation of the initial problem. Therefore, for any admissible policy $\pi \in \text{LP-update}$, MILP-update, the following holds:

$$ V_{\pi}^{(N)}(\mathbf{x}, T) \le V_{\text{rel}}(\mathbf{x}, T) $$

In Figure 5, we present a performance comparison between the LP-update and MILP-update policies. As anticipated, the MILP-update demonstrates superior performance compared to the original LP-update, particularly for smaller values of $N$.

<figure><img src="image_7.png" /><figcaption>Figure 5 Sub-optimality gap of the LP-update policy compared to the MILP-update policy.</figcaption></figure>

From a computational perspective, addressing (33) is more demanding in terms of resources compared to solving the original LP problem. Our experimental findings revealed that solving the MILP required no more than twice the time needed for the original LP for this specific model having 10 states. Notably, the duration for solving the MILP appeared to be unaffected by the size of $N$, suggesting that increases in $N$ do not proportionally complicate the MILP solution process. This consistency in computational time, regardless of the scale of $N$, offers an interesting topic for further investigation, especially to determine if this trend is consistent across various scenarios and applications. It is also important to observe that the zigzag effect in the LP-update policy, resulting from the cutoff rounding, is absent in the MILP-update policy. This suggests that the MILP-update policy might be a more suitable and natural choice for implementation.
## D.2. Perfect rounding

The constructed $\mathbf{Y}^{(N)}$, even using MILP, might still be far from $\mathbf{Y}$ for small values of $N$. To obtain a $\mathbf{Y}^{(N)}$ closer to $\mathbf{Y}$, the approach developed in Gast et al. [14] for two-action bandits is to use the so-called randomized rounding, that consists of using a random vector $\mathbf{Y}^{(N)}$ such that $\mathbb{E}[\mathbf{Y}^{(N)}] = \mathbf{Y}$ and $\mathbf{Y}^{(N)} \in \mathcal{Y}^{(N)}(\mathbf{X}^{(N)})$ almost surely. If such a random vector exists, we say that the problem admits a perfect rounding.

**DEFINITION 2 (PERFECT ROUNDING).** An LP problem admits a perfect rounding for an integer $N$, if for all $\mathbf{x} \in \Delta^{(N)}$ and for all $\mathbf{y} \in \mathcal{Y}(\mathbf{x})$, there exists a random vector $\mathbf{Y}^{(N)}$ such that $\mathbf{Y}^{(N)} \in \mathcal{Y}^{(N)}(\mathbf{X}^{(N)})$ almost surely and $\mathbb{E}[\mathbf{Y}^{(N)}] = \mathbf{y}$.

The next theorem shows that, if there exists a perfect rounding, we obtain the much faster exponential convergence rate by avoiding the rounding errors. The result is summarized in the following convergence theorem:

**THEOREM 4.** Denote by $V_{\text{LP-update}}^{(N)}(\mathbf{x}, T)$ the value of the LP-update policy defined in Algorithm 1 or 2, and by $V_{\text{rel}}(\mathbf{x}, T)$ the value of the linear program (5). Then

* If the weakly-coupled MDP with statistically identical arms is non-degenerate, then there exist constants $C_1, C_2 > 0$ such that for any $N$ for which the problem admits a perfect rounding:

$$ \left| V_{\text{LP-update}}^{(N)}(\mathbf{x}, T) - V_{\text{rel}}(\mathbf{x}, T) \right| \leq C_1 \cdot e^{-C_2 N}. $$

* There exists a non-degenerate weakly-coupled MDP with statistically identical arms that admits a perfect rounding for an infinite number of $N$, and two constants $C'_1, C'_2$ such that for all such $N$:

$$ \left| V_{\text{LP-update}}^{(N)}(\mathbf{x}, T) - V_{\text{rel}}(\mathbf{x}, T) \right| \geq C'_1 \cdot e^{-C'_2 N}. $$

*Proof.* The proof of the first item is identical to the proof of Theorem 3 until the last equation, for which if we have a perfect rounding, then $\mathbb{E}[\mathbf{Y}^{(N)}(t) - \mathbf{Y}(t)] = 0$. The proof of the lower bound in the second item is provided in Section A.4. Q.E.D.

**Comparison with existing literature** The authors of [14] study the single constraint case where $\mathcal{A} = \{0, 1\}$, $D(s, a) = a$ and $b = \alpha$. They show that a perfect rounding exists for all $N$ such that $\alpha N$ is an integer. Another example is when the constraints (12b)-(12c)-(12d) form a totally unimodular matrix (see for instance Hoffman and Kruskal [19]) for all $t$: in that case, the solution to the LP satisfies $\mathbf{y} \in \mathcal{Y}^{(N)}(\mathbf{X}^{(N)})$ automatically and no rounding is needed. There are, however, many situations for which perfect rounding is impossible. In such cases, a method to improve the simple truncation (32) is to find $\mathbf{Y}^{(N)} \in \mathcal{Y}^{(N)}(\mathbf{X}^{(N)})$ that minimizes the distance $\|\mathbf{Y}^{(N)} - \mathbf{Y}\|_1$. This can be computed by solving an integer linear program.

The researchers in [39, 38] implement a randomized approach to derive an admissible $\mathbf{Y}^{(N)}$ from $\mathbf{y}^*$. This method is detailed in our description of their algorithm in Algorithm 3. Initially, their algorithm sets all arms to action 0. Subsequently, it iteratively processes each arm. For an arm in state $s$, the algorithm probabilistically selects a new action $a$ based on the ratio $y_{s,a}^*(t)/x_s^*(t)$. This action $a$ is then assigned to the arm, provided it does not breach any budget constraints. The central limit theorem underpins the rationale for this approach, ensuring that $\mathbf{Y}^{(N)}$ approximates $\mathbf{y}^* + \mathcal{O}(1/\sqrt{N})$. However, it does not extend to achieving a tighter approximation of $\mathbf{y}^* + \mathcal{O}(1/N)$. Consequently, even in the context of non-degenerate problems, the randomized process employed by these algorithms results in a sub-optimality gap of $\mathcal{O}(1/\sqrt{N})$ rather than the refined $\mathcal{O}(1/N)$.
## D.2. Perfect rounding

The constructed $\mathbf{Y}^{(N)}$, even using MILP, might still be far from $\mathbf{Y}$ for small values of $N$. To obtain a $\mathbf{Y}^{(N)}$ closer to $\mathbf{Y}$, the approach developed in Gast et al. [14] for two-action bandits is to use the so-called randomized rounding, that consists of using a random vector $\mathbf{Y}^{(N)}$ such that $\mathbb{E}[\mathbf{Y}^{(N)}] = \mathbf{Y}$ and $\mathbf{Y}^{(N)} \in \mathcal{Y}^{(N)}(\mathbf{X}^{(N)})$ almost surely. If such a random vector exists, we say that the problem admits a perfect rounding.

**DEFINITION 2 (PERFECT ROUNDING).** An LP problem admits a perfect rounding for an integer $N$, if for all $\mathbf{x} \in \Delta^{(N)}$ and for all $\mathbf{y} \in \mathcal{Y}(\mathbf{x})$, there exists a random vector $\mathbf{Y}^{(N)}$ such that $\mathbf{Y}^{(N)} \in \mathcal{Y}^{(N)}(\mathbf{X}^{(N)})$ almost surely and $\mathbb{E}[\mathbf{Y}^{(N)}] = \mathbf{y}$.

The next theorem shows that, if there exists a perfect rounding, we obtain the much faster exponential convergence rate by avoiding the rounding errors. The result is summarized in the following convergence theorem:

**THEOREM 4.** Denote by $V_{\text{LP-update}}^{(N)}(\mathbf{x}, T)$ the value of the LP-update policy defined in Algorithm 1 or 2, and by $V_{\text{rel}}(\mathbf{x}, T)$ the value of the linear program (5). Then

* If the weakly-coupled MDP with statistically identical arms is non-degenerate, then there exist constants $C_1, C_2 > 0$ such that for any $N$ for which the problem admits a perfect rounding:

$$ \left| V_{\text{LP-update}}^{(N)}(\mathbf{x}, T) - V_{\text{rel}}(\mathbf{x}, T) \right| \le C_1 \cdot e^{-C_2 N}. $$

* There exists a non-degenerate weakly-coupled MDP with statistically identical arms that admits a perfect rounding for an infinite number of $N$, and two constants $C'_1, C'_2$ such that for all such $N$:

$$ \left| V_{\text{LP-update}}^{(N)}(\mathbf{x}, T) - V_{\text{rel}}(\mathbf{x}, T) \right| \ge C'_1 \cdot e^{-C'_2 N}. $$

*Proof.* The proof of the first item is identical to the proof of Theorem 3 until the last equation, for which if we have a perfect rounding, then $\mathbb{E}[\mathbf{Y}^{(N)}(t) - \mathbf{Y}(t)] = 0$. The proof of the lower bound in the second item is provided in Section A.4. Q.E.D.

**Comparison with existing literature** The authors of [14] study the single constraint case where $\mathcal{A} = \{0, 1\}$, $D(s, a) = a$ and $b = \alpha$. They show that a perfect rounding exists for all $N$ such that $\alpha N$ is an integer. Another example is when the constraints (12b)-(12c)-(12d) form a totally unimodular matrix (see for instance Hoffman and Kruskal [19]) for all $t$: in that case, the solution to the LP satisfies $\mathbf{y} \in \mathcal{Y}^{(N)}(\mathbf{X}^{(N)})$ automatically and no rounding is needed. There are, however, many situations for which perfect rounding is impossible. In such cases, a method to improve the simple truncation (32) is to find $\mathbf{Y}^{(N)} \in \mathcal{Y}^{(N)}(\mathbf{X}^{(N)})$ that minimizes the distance $\|\mathbf{Y}^{(N)} - \mathbf{Y}\|_1$. This can be computed by solving an integer linear program.

The researchers in [39, 38] implement a randomized approach to derive an admissible $\mathbf{Y}^{(N)}$ from $\mathbf{y}^*$. This method is detailed in our description of their algorithm in Algorithm 3. Initially, their algorithm sets all arms to action 0. Subsequently, it iteratively processes each arm. For an arm in state $s$, the algorithm probabilistically selects a new action $a$ based on the ratio $y_{s,a}^*(t)/x_s^*(t)$. This action $a$ is then assigned to the arm, provided it does not breach any budget constraints. The central limit theorem underpins the rationale for this approach, ensuring that $\mathbf{Y}^{(N)}$ approximates $\mathbf{y}^* + \mathcal{O}(1/\sqrt{N})$. However, it does not extend to achieving a tighter approximation of $\mathbf{y}^* + \mathcal{O}(1/N)$. Consequently, even in the context of non-degenerate problems, the randomized process employed by these algorithms results in a sub-optimality gap of $\mathcal{O}(1/\sqrt{N})$ rather than the refined $\mathcal{O}(1/N)$.