2.1. Restless bandit process

A discrete-time restless bandit process (RB) is a controlled Markov process $(\mathcal{X}, \{0, 1\}, \{P(a)\}_{a \in \{0, 1\}}, c, x_0)$ where $\mathcal{X}$ denotes the state space, which is a finite or countable set; $\{0, 1\}$ denotes the action space, where the action 0 is called the passive action and the action 1 is the active action; P(a), $a \in \{0, 1\}$ , denotes the transition matrix when the action a is chosen; $c\,:\, \mathcal{X} \times \{0, 1\} \to \mathbb{R}$ denotes the cost function; and $x_0$ denotes the initial state. We use $X_t$ and $A_t$ to denote the action of the process at time t. The process evolves in a controlled Markov manner; i.e., for any realization $x_{0\,:\,t+1}$ of $X_{0\,:\,t+1}$ and $a_{0\,:\,t+1}$ of $A_{0\,:\,t+1}$ , we have $\mathbb{P}\big( {X}_{t+1} = {x}_{t+1} | {X}_{0\,:\,t} = {x}_{0\,:\,t}, {A}_{0\,:\,t} = {a}_{0\,:\,t}\big)=\mathbb{P} \big( X_{t+1} = x_{t+1} | X_{t} = x_{t}, A_{t} = a_t\big)$ , which we denote by $P_{x_t x_{t+1}}(a_t)$ .

2.2. Restless multi-armed bandit problem

A restless multi-armed bandit is a collection of n independent RBs $\big(\mathcal{X}^i, \{0, 1\}, \{P^i(a)\}_{a \in \{0, 1\}}, c^i, x^i_0\big)$ , $i \in {\mathcal N} \,:\!=\, \{1, \ldots, n\}$ . A decision-maker observes the state of all RBs, may choose to activate only $m < n$ of them, and incurs a cost equal to the sum of the costs incurred by the individual RBs.

Let $\boldsymbol{\mathcal{X}} \,:\!=\, \prod_{i \in \mathcal N} \mathcal{X}^i$ and $\boldsymbol{\mathcal A}(m) \,:\!=\, \bigl\{ \boldsymbol{a} = \big(a^1, \ldots, a^n\big) \in {\mathcal A}^{n} \,:\, \sum_{i \in {\mathcal N}} a^i = m\bigr\}$ denote the joint state space and the feasible action space, respectively. Let ${\boldsymbol{X}}_t \,:\!=\, \big(X^1_t, \dots X^n_t\big)$ and ${\boldsymbol{A}}_t = \big(A^1_t, \dots, A^n_t\big)$ denote the joint state and actions at time t. As the RBs evolve independently, for any realization ${\boldsymbol{x}}_{0\,:\,t}$ of ${\boldsymbol{X}}_{0\,:\,t}$ and $\boldsymbol{a}_{0\,:\,t}$ of ${\boldsymbol{A}}_{0\,:\,t}$ , we have

\begin{align*}\mathbb{P} \big( {\boldsymbol{X}}_{t+1} = {\boldsymbol{x}}_{t+1} | {\boldsymbol{X}}_{0\,:\,t} = {\boldsymbol{x}}_{0\,:\,t}, {\boldsymbol{A}}_{0\,:\,t} = \boldsymbol{a}_{0\,:\,t} \big) = \prod_{i = 1}^{n} \mathbb{P} \Big( X^{i}_{t+1} = x^{i}_{t+1} | X^{i}_{t} = x^i_{t}, A^{i}_{t} = a^i_t \Big).\end{align*}

When the system is in the state ${\boldsymbol{x}}_t = \big(x^1_t, \ldots, x^n_t\big)$ and the decision-maker chooses action $\boldsymbol{a}_t = \big(a^1_t, \ldots, a^n_t\big)$ , the system incurs a cost $\bar{c}({\boldsymbol{x}}_t, \boldsymbol{a}_t) \,:\!=\, \sum_{i \in {\mathcal N}} c^i\big(x^i_t, a^i_t\big)$ . The decision-maker chooses his actions using a time-homogeneous Markov policy ${\boldsymbol{g}}\,:\, \boldsymbol{\mathcal X} \to \boldsymbol{\mathcal A}(m)$ , i.e., chooses ${\boldsymbol{A}}_t = {\boldsymbol{g}}({\boldsymbol{X}}_t)$ . The performance of any Markov policy ${\boldsymbol{g}}$ is given by

\begin{equation*}J^{({\boldsymbol{g}})}({\boldsymbol{x}}_0) \,:\!=\, (1-\beta) \mathbb{E}\Biggl[ \sum_{t = 0}^{\infty}\beta^t \bar{c}({\boldsymbol{X}}_t, {\boldsymbol{g}}({\boldsymbol{X}}_t)) \bigg| {\boldsymbol{X}}_0 = {\boldsymbol{x}}_0\Biggr],\end{equation*}

where $\beta \in (0, 1)$ is the discount factor and ${\boldsymbol{x}}_0$ is the initial state of the system.

We are interested in the following optimization problem.

Problem 1 is a multi-stage stochastic control problem and one can obtain an optimal solution using dynamic programming. However, the dynamic programming solution is intractable for large n since the cardinality of the state space is $\prod_{i \in {\mathcal N}}|\mathcal{X}^i|$ , which grows exponentially with n. In the next section, we describe a heuristic known as the Whittle index to efficiently obtain a suboptimal solution of the problem.

2.3. Indexability and the Whittle index

Consider an RB $(\mathcal{X}, \{0, 1\}, \{P(a)\}_{a \in \{0, 1\}}, c, x_0)$ . For any $\lambda \in \mathbb{R}$ , we consider a Markov decision process $\{ \mathcal{X}, \{0, 1\}, \{P(a)\}_{a \in \{0, 1\}}, c_\lambda, x_0 \}$ , where

(1) \begin{equation} c_\lambda(x, a) \,:\!=\, c(x, a) + \lambda a, \quad \forall x \in \mathcal{X}, \forall a \in \{0, 1\}.\end{equation}

The parameter $\lambda$ may be viewed as a penalty for taking the active action. The performance of any time-homogeneous policy $g\,:\, \mathcal{X} \to \{0, 1\}$ is

(2) \begin{equation}J^{(g)}_{\lambda}(x_0) \,:\!=\, (1-\beta) \mathbb{E}\Biggl[ \sum_{t = 0}^{\infty} \beta^t c_\lambda(X_t, g(X_t)) \bigg| X_0 = x_0 \Biggr].\end{equation}

Consider the following optimization problem.

Problem 2 is also a Markov decision process and one can obtain an optimal solution using dynamic programming. Let $V_\lambda\,:\, \mathcal{X} \to \mathbb{R}$ be the unique fixed point of the following:

(3) \begin{equation}V_\lambda(x) = \min \bigl\{ H_\lambda(x, 0) , H_\lambda(x, 1) \bigr\},\quad\forall x \in \mathcal{X},\end{equation}

where

(4) \begin{equation} H_\lambda(x, a) = (1-\beta) c_\lambda(x, a) + \beta \sum_{y \in \mathcal{X}} P_{xy}(a) V_\lambda(y), \quad a \in \{0, 1\}.\end{equation}

Let $g_\lambda(x)$ denote the minimizer of the right-hand side of $(3)$ , where we set $g_\lambda(x) = 1$ if $H_\lambda(x, 0) = H_\lambda(x, 1)$ . Then, from Markov decision theory [Reference Puterman28], we know that the time-homogeneous policy $g_\lambda$ is optimal for Problem 2.

Define the passive set ${\Pi}_\lambda$ to be the set of states where the passive action is optimal, i.e.,

(5) \begin{equation}{\Pi}_\lambda \,:\!=\, \left\{ x \in \mathcal{X}\,:\, g_\lambda(x) = 0 \right\}\!.\end{equation}

Alternatively, the Whittle index w(x) is a value of the penalty $\lambda$ for which the optimal policy is indifferent between taking the active and the passive action when the RB is in state x.

2.4. Whittle index heuristic

A restless multi-armed bandit problem is said to be indexable if all RBs are indexable. For indexable problems, the Whittle index heuristic is as follows: Compute the Whittle indices of all arms offline. Then, at each time, obtain the Whittle indices of the current state of all arms and play arms with the m largest Whittle indices.

As mentioned earlier, the Whittle index policy is a popular approach for restless bandits because (i) its complexity is linear in the number of alternatives and (ii) its performance is often close to optimal in practice [Reference Ansell, Glazebrook, Niño-Mora and O’Keeffe5, Reference Glazebrook and Mitchell16, Reference Glazebrook, Ruiz-Hernandez and Kirkbride19]. However, there are only a few general conditions to check indexability for general models.

2.5. Alternative characterizations of the passive set

We now present alternative characterizations of the passive set, which are important for the sufficient conditions for indexability that we provide later.

Let $\Sigma$ denote the family of all stopping times with respect to the natural filtration of $\{X_t\}_{t \ge 0}$ . For any state $x \in \mathcal{X}$ , penalty $\lambda \in \mathbb{R}$ , and stopping time $\tau \in \Sigma$ , define

(6) \begin{align} M(x, \tau) &\,:\!=\, \mathbb{E}\bigl[ \beta^{\tau} | X_0 = x, \{ A_t = 0 \}_{t = 0}^{\tau-1} \bigr], \notag \\ L(x, \tau) &\,:\!=\, \mathbb{E}\Bigg[ \sum_{t = 0}^{\tau-1} \beta^{t} c(X_t, 0) + \beta^{\tau} c(X_{\tau}, 1) \Bigm| X_0 = x, \{ A_t = 0 \}_{t = 0}^{\tau-1} \Bigg], \notag \\ W_\lambda(x) &\,:\!=\, (1-\beta)\lambda + \beta \sum_{y \in \mathcal{X}} P_{xy}(1) V_\lambda(x).\end{align}

Let $h_{\tau, \lambda}$ denote the (history-dependent) policy that takes the passive action up to time $\tau-1$ , takes the active action at time $\tau$ , and then follows the optimal policy $g_\lambda$ (for Problem 2). We now present different characterizations of the passive set.

See Appendix A for the proof.

3.1. Preliminary results

Consider an RB $(\mathcal{X}, \{0, 1\}, \{P(a)\}_{a \in \{0, 1\}}, c, x_0)$ . For any Markov policy $g \,:\, \mathcal{X} \to \{0,1\}$ and $\lambda \in \mathbb{R}$ , we can write

(7) \begin{equation} J^{(g)}_\lambda(x) = D^{(g)}(x) + \lambda N^{(g)}(x),\end{equation}

where

\begin{align*} D^{(g)}(x) &\,:\!=\, (1-\beta) \mathbb{E} \Bigg[ \sum_{t = 0}^{\infty} \beta^t c(X_t, g(X_t)) \bigg| X_0 = x \Bigg] \end{align*}

and

\begin{align*} N^{(g)}(x) &\,:\!=\, (1-\beta) \mathbb{E} \Bigg[ \sum_{t = 0}^{\infty} \beta^t g(X_t) \bigg| X_0 = x \Bigg] \end{align*}

are the expected discounted total cost and the expected number of activations under the policy g starting at the initial state x. $D^{(g)}({\cdot})$ and $N^{(g)}({\cdot})$ can be computed using policy evaluation formulas. In particular, define $P^{(g)} \,:\, \mathcal{X} \times\mathcal{X} \to \mathbb{R}$ and $c^{(g)} \,:\, \mathcal{X} \to \mathbb{R}$ as follows: $P^{(g)}_{xy} = P_{xy}(g(x)) \text{ and } c^{(g)}_\lambda(x) =c_\lambda(x, g(x)) = c^{(g)}(x, g(x)) + \lambda g(x)$ for any $x \in\mathcal{X}$ . We also view g as an element in $\{0,1\}^{|\mathcal{X}|}$ . Then, using the policy evaluation formula for infinite-horizon Markov decision processes [Reference Puterman28], we obtain

(8) \begin{equation} D^{(g)}(x) = (1-\beta)\bigl[(I - \beta P^{(g)})^{-1} c^{(g)}\bigr](x) \text{ and } N^{(g)}(x) =(1-\beta) \bigl[(I - \beta P^{(g)})^{-1} g\bigr](x).\end{equation}

We now provide a geometric interpretation of the value function $V_\lambda(x)$ as a function of $\lambda$ . For any $g \in {\mathcal G}$ , the plot of $J^{(g)}_\lambda(x) = D^{(g)}(x) + \lambda N^{(g)}(x)$ as a function of $\lambda$ is a straight line with y-intercept $D^{(g)}(x)$ and slope $N^{(g)}(x)$ . By definition, $ V_\lambda(x) = \inf_{g \in {\mathcal G}} J^{(g)}_\lambda(x). $ Thus, $V_\lambda(x)$ is the lower concave envelope of the family of straight lines $\big\{J^{(g)}_\lambda(x)\big\}_{g \in {\mathcal G}}$ . (See Figure 1 for an illustration.) Thus, we have the following.

Figure 1. An illustration of the plot of $J^{({\cdot})}_\lambda(x)$ versus $\lambda$ for $g \in {\mathcal G} \,:\!=\, \{h_1, h_2, h_3\}$ . Let $\lambda^{ij}_c$ denote the $\lambda$ -value of the intersection of $J^{(h_i)}_\lambda(x)$ and $J^{(h_j)}_\lambda(x)$ . Note that in this plot, for all $\lambda \in \left( -\infty, \lambda^{12}_c \right]$ the policy $h_1$ is optimal; for all $\lambda \in \big[ \lambda^{12}_c, \lambda^{23}_c \big]$ the policy $h_2$ is optimal; and for all $\lambda \in \big[ \lambda^{23}_c, \infty\big)$ the policy $h_3$ is optimal. The lower concave envelope of $J^{(h_i)}_\lambda(x)$ (shown as a thick line) is the value function $V_\lambda(x)$ , which is piecewise linear, concave, increasing, and continuous.

Lemma 2. For any $\lambda^{\prime}, \lambda^{\prime\prime} \in \mathbb{R}$ ,

\begin{align*} (\lambda^{\prime\prime} - \lambda^{\prime}) N^{(g_{\lambda^{\prime\prime}})}(x) \le V_{\lambda^{\prime\prime}}(x) - V_{\lambda^{\prime}}(x) \le (\lambda^{\prime\prime} - \lambda^{\prime}) N^{(g_{\lambda^{\prime}})}(x), \quad \forall x \in \mathcal{X}. \end{align*}

Consequently, $N^{(g_\lambda)}(x)$ is non-increasing in $\lambda$ .

Proof. Recall that $V_\lambda(x) = J^{(g_\lambda)}_\lambda(x) \le J^{(g_{\lambda^{\prime}})}_\lambda(x)$ for any $\lambda^{\prime} \neq \lambda$ . Thus,

(9) \begin{align} V_{\lambda^{\prime\prime}}(x) - V_{\lambda^{\prime}}(x) &= J^{(g_{\lambda^{\prime\prime}})}_{\lambda^{\prime\prime}}(x) - J^{(g_{\lambda^{\prime}})}_{\lambda^{\prime}}(x) \le J^{(g_{\lambda^{\prime}})}_{\lambda^{\prime\prime}}(x) - J^{(g_{\lambda^{\prime}})}_{\lambda^{\prime}}(x) \notag \\ &\stackrel{(a)}= (\lambda^{\prime\prime} - \lambda^{\prime})N^{(g_{\lambda^{\prime}})}(x), \end{align}

where (a) follows from (7). Similarly, we have

(10) \begin{align} V_{\lambda^{\prime\prime}}(x) - V_{\lambda^{\prime}}(x) &= J^{(g_{\lambda^{\prime\prime}})}_{\lambda^{\prime\prime}}(x) - J^{(g_{\lambda^{\prime}})}_{\lambda^{\prime}}(x) \ge J^{(g_{\lambda^{\prime\prime}})}_{\lambda^{\prime\prime}}(x) - J^{(g_{\lambda^{\prime\prime}})}_{\lambda^{\prime}}(x) \notag \\ &\stackrel{(a)}= (\lambda^{\prime\prime} - \lambda^{\prime})N^{(g_{\lambda^{\prime\prime}})}(x), \end{align}

where (a) follows from (7). The result follows from combining the above inequalities.

3.2. Sufficient conditions for indexability

See Appendix B for the proof. The sufficient conditions of Theorem 1 can be difficult to verify. Simpler sufficient conditions are stated below.

See Appendix C for the proof.

Some remarks

1. 1. The sufficient conditions of Theorem 1 and of parts (a), (c), and (d) of Proposition 2 may be viewed as bounds on the discount factor $\beta$ for which an RB is indexable. Numerical experiments to explore such a property are presented in [Reference Niño-Mora25]. A qualitatively similar result was established in [Reference Niño-Mora23, Corollary 5], which showed that a restless bandit process is generalized-conservation-laws (GCL-) indexable for sufficiently small discount factors. GCL-indexability is a sub-class of PCL-indexability, which is a sub-class of Whittle indexability. Thus, Proposition 2 provides a quantitative characterization of the qualitative observation made in [Reference Bertsimas and Niño-Mora9] and generalizes it to a broader class of models.

2. 2. We refer to models that satisfy the sufficient condition of Proposition 2(b) as restless bandits with controlled restarts. Such models arise in various scheduling problems (e.g., machine maintenance, surveillance) where taking the active action resets the state according to a known probability distribution. Specific instances of such models are considered in [Reference Akbarzadeh and Mahajan3, Reference Wang, Ren, Mo and Shi30]. The special case when the active action resets to a specific (pristine) state are considered in [Reference Glazebrook, Mitchell and Ansell18, Reference Glazebrook, Ruiz-Hernandez and Kirkbride19].

3. 3. A different class of restart models have been considered in [Reference Avrachenkov, Ayesta, Doncel and Jacko7, Reference Ayesta, Erausquin and Jacko8, Reference Jacko20], where the passive action resets the state of the arm. Note that such models do not satisfy Proposition 2(b), and additional modeling assumptions are required to establish indexability. See [Reference Avrachenkov, Ayesta, Doncel and Jacko7, Reference Ayesta, Erausquin and Jacko8, Reference Jacko20] for details.

4.1. An efficient implementation using the Sherman–Morrison formula

We now present an efficient implementation of Algorithm 1 using the Sherman–Morrison inverse formula. Suppose $A \in \mathbb {R} ^{n\times n}$ is an invertible square matrix and $u,v\in \mathbb {R} ^{n}$ are column vectors such that $A + u v^{\textsf{T}}$ is invertible. Then the Sherman–Morrison inverse formula is

(15) \begin{equation} \big(A+uv^{\textsf {T}}\big)^{-1}=A^{-1}-{\frac{A^{-1}uv^{\textsf {T}}A^{-1}} {1+v^{\textsf {T}}A^{-1}u}}.\end{equation}

Furthermore, given $b \in \mathbb{R}^n$ , if x is the solution of $Ax = b$ and y is the solution $Ay = u$ , then the solution of $\big(A + u v^{\textsf{T}}\big) \tilde x = b$ is given by the following (see [Reference Egidi and Maponi12, Corollary 2]):

(16) \begin{equation} \tilde x = x - \frac{v^{\mathsf{T}} x}{1 + v^{\mathsf{T}} y} y.\end{equation}

Let $\Phi^{(g)} = \big(I - \beta P^{(g)}\big)^{-1}$ . Note that for any Markov policy g, $\big(I - \beta P^{(g)}\big)$ is invertible, because $\beta P^{(g)}$ is a sub-stochastic matrix and has a spectral radius less than 1. Therefore, the conditions for using the Sherman–Morrison formula are satisfied. Hence, using (16), Equation (8) may be written (in matrix form) as

(17) \begin{equation} D^{(g)} = (1 - \beta) \Phi^{(g)} c^{(g)} \text{ and } N^{(g)} = (1 - \beta) \Phi^{(g)} g.\end{equation}

Now, for any $d \in \{0, \ldots, K_D-1\}$ and a state $y \in {\mathcal X} \backslash {\mathcal P}_d$ , consider the policies $h_d = \bar{g}^{({\mathcal P}_d)}$ and $h_{d,y} = \bar{g}^{({\mathcal P}_{d,y})}$ . Let $e_y$ denote the unit vector with 1 in the yth location, and let $\rho_y$ be a vector given by $ [\rho_y]_x = P_{yx}(1) - P_{yx}(0) $ , for all $x \in {\mathcal X}$ . Then $P^{\left(h_{d,y}\right)} = P^{(h_d)} - e_y \rho^{\mathsf{T}}_y$ . Therefore,

\begin{align*} I - \beta P^{\left(h_{d,y}\right)} = \bigl( I - \beta P^{\left(h_{d}\right)} \bigr) + \beta e_y \rho^{\mathsf{T}}_y. \end{align*}

Let $\Phi^{\left(h_{d}\right)}_{\cdot y}$ denote the yth column of $\Phi^{\left(h_{d}\right)}$ , i.e., $\Phi^{\left(h_{d}\right)}_{\cdot y} = \Phi^{\left(h_{d}\right)} e_y$ . Then, by the Sherman–Morrison inverse formula ((15) and (16)), we have

(18) \begin{align} \Phi^{\left(h_{d,y}\right)} = \Phi^{\left(h_{d}\right)} - \dfrac{\beta\Phi^{\left(h_{d}\right)} e_y \rho^{\mathsf{T}}_y \Phi^{\left(h_{d}\right)} }{1 + \beta \rho^{\mathsf{T}}_y \Phi^{\left(h_{d}\right)}_{\cdot y}}, \end{align}

(19) \begin{align} D^{\left(h_{d,y}\right)} = D^{\left(h_{d}\right)} - \dfrac{\beta\rho^{\mathsf{T}}_y D^{\left(h_{d}\right)} }{1 + \beta\rho^{\mathsf{T}}_y \Phi^{\left(h_{d}\right)}_{\cdot y}} \Phi^{\left(h_{d}\right)}_{\cdot y}, \quad N^{\left(h_{d,y}\right)} = N^{\left(h_{d}\right)} - \dfrac{\beta\rho^{\mathsf{T}}_y N^{\left(h_{d}\right)} }{1 + \beta\rho^{\mathsf{T}}_y \Phi^{\left(h_{d}\right)}_{\cdot y}} \Phi^{\left(h_{d}\right)}_{\cdot y}.\end{align}

Thus, if $\Phi^{\left(h_{d}\right)}$ has been computed, then $\Phi^{\left(h_{d,y}\right)}$ can be computed in $\mathcal{O}\big(K^2\big)$ computations. In addition, if $\Phi^{\left(h_{d}\right)}$ , $D^{\left(h_{d}\right)}$ , and $N^{\left(h_{d}\right)}$ have been computed, then $D^{\left(h_{d,y}\right)}$ and $N^{\left(h_{d,y}\right)}$ can be computed in $\mathcal{O}(K)$ computations.

So we can use (18) and (19) to implement Algorithm 1 in a more efficient manner. However, there is one additional step that needs to be handled, which we explain next.

As $h_{d+1} = \bar{g}^{\big({\mathcal P}_d \cup {\Gamma_{d+1}}\big)}$ , we get

(20) \begin{equation} I - \beta P^{\left(h_{d+1}\right)} = \bigl(I - \beta P^{(h_d)}\bigr) + \beta \sum_{y \in \Gamma_{d+1}} e_y \rho^{\mathsf{T}}_y.\end{equation}

Thus, $I - \beta P^{\left(h_{d+1}\right)}$ is a rank- $|\Gamma_{d+1}|$ update of $I - \beta P^{\left(h_{d}\right)}$ . When $|\Gamma_{d+1}| > 1$ , we can either sequentially apply Equations (18) and (19) for all $y \in \Gamma_{d+1}$ or use the Woodbury formula to compute $\Phi^{\left(h_{d+1}\right)}$ and compute $D^{\left(h_{d+1}\right)}$ and $N^{\left(h_{d+1}\right)}$ using (17). The Woodbury formula is a generalization of the Sherman–Morrison formula. Suppose $A \in \mathbb{R}^{n\times n}$ is an invertible matrix and $U, V \in \mathbb{R}^{n \times m}$ are such that $A + U V^{\mathsf{T}}$ is invertible. Then the Woodbury formula is

\begin{align*} (A + UV^{\mathsf{T}})^{-1} = A^{-1} - A^{-1} U(I + V^{\mathsf{T}} U)^{-1} V^{\mathsf{T}} A^{-1}.\end{align*}

The complexity of sequentially applying the Sherman–Morrison formula is $\mathcal{O}\big(|\Gamma_{d+1}|K^2\big)$ to compute $\Phi^{\left(h_{d+1}\right)}$ and $\mathcal{O}\big(|\Gamma_{d+1}| K\big)$ to compute $D^{\left(h_{d+1}\right)}$ and $N^{\left(h_{d+1}\right)}$ . The complexity of using the Woodbury formula is $\mathcal{O}\big(|\Gamma_{d+1}|^{2.807} + K^2\big)$ to compute $\Phi^{\left(h_{d+1}\right)}$ and $\mathcal{O}\big(K^2\big)$ to compute $D^{\left(h_{d+1}\right)}$ and $N^{\left(h_{d+1}\right)}$ .

We show the complete algorithm to efficiently compute the Whittle index in Algorithm 2, where we use sequential application of Sherman–Morrison formula to compute $\Phi^{\left(h_{d+1}\right)}$ , $D^{\left(h_{d+1}\right)}$ , and $N^{\left(h_{d+1}\right)}$ .

Some remarks

1. 1. The idea of computing the index by iteratively sorting the states according to their index is commonly used in the algorithms to compute the Gittins index; for example, the largest-remaining-index algorithm, the state-elimination algorithm, the triangularization algorithm, and the fast-pivoting algorithm use variations of this idea. See [Reference Chakravorty and Mahajan10] for details.

2. 2. The computational complexity of Algorithm 2 is $\mathcal{O}\!\left(K^3\right)$ , which can be characterized as follows. The algorithm starts with computing $\Phi^{(h_0)}$ , which requires $\mathcal{O}\big(K^{2.807}\big)$ computations (using Strassen’s algorithm), and $D^{(h_0)}$ and $N^{(h_0)}$ , each of which requires $\mathcal{O}\big(K^2\big)$ computations. Then, in the inner for loop, computing each of $D^{\left(h_{d, y}\right)}$ , $N^{\left(h_{d, y}\right)}$ , and $\mu^*_{d,y}$ requires $\mathcal{O}(K)$ computations, and the inner loop is executed $|{\mathcal X} \backslash {\mathcal P}_d|$ times. Afterwards, updating $\Phi^{\left(h_{d+1}\right)}$ , $D^{\left(h_{d+1}\right)}$ , and $N^{\left(h_{d+1}\right)}$ requires $\mathcal{O}\big(|\Gamma_{d+1}|K^2\big)$ , $\mathcal{O}(|\Gamma_{d+1}|K)$ , and $\mathcal{O}(|\Gamma_{d+1}|K)$ computations, respectively. Therefore, the computational complexity of the algorithm is

   \begin{align*} & \mathcal{O}\big(K^{2.807}\big) + \mathcal{O}\big(K^2\big) + \sum_{d = 1}^{K_D} \big( \mathcal{O}(|{\mathcal X} \backslash {\mathcal P}_d| K) + \mathcal{O}\big(|\Gamma_{d+1}| K^2\big) + \mathcal{O}\big(2 |\Gamma_{d+1}| K\big) \big) \\ &\le \mathcal{O}\big(K^{2.807}\big) + \sum_{d = 1}^{K_D} \mathcal{O}\big(K^2\big) + \mathcal{O}\Biggl(\Biggl[ \sum_{d = 1}^{K_D} |\Gamma_{d+1}| \Biggr] K^2\Biggr) \\ &\le \mathcal{O}\big(K^{2.807}\big) + \mathcal{O}\!\left(K^3\right) + \mathcal{O}\!\left(K^3\right) \le \mathcal{O}\!\left(K^3\right), \end{align*}
   where the first inequality uses the fact that $|\mathcal{X}\setminus\mathcal{P}_d| \le K$ and the second inequality uses the fact that $\sum_{d=1}^{K_D}| \Gamma_{d+1}| = K$ .

3. 3. Note that Algorithm 2 computes the Whittle index exactly. In contrast, using binary search [Reference Akbarzadeh and Mahajan2] computes the Whittle index approximately. Let $C_{\max}$ and $C_{\min}$ denote the upper and lower bound on the per-step cost. Then we know that for any state x, $w(x) \in [ C_{\min}, C_{\max} ]$ . Now, suppose we want to compute the Whittle index to an accuracy of $\delta$ . Then the interval $[C_{\min}, C_{\max}]$ needs to be divided into $\log_2( (C_{\max} - C_{\min})/\delta)$ steps. For each step of the binary search, we need to solve the dynamic program (3). There are two main algorithms to solve dynamic programs [Reference Puterman28]: policy iteration, which gives an exact solution, and value iteration, which gives an approximate solution.

   Each iteration of policy iteration has two steps: (i) policy evaluation, which is done by solving a linear system and has complexity $\mathcal{O}\!\left(K^3\right)$ , and (ii) policy improvement, which requires $\mathcal{O}\big(K^2\big)$ computations. Thus, the overall complexity of policy iteration is $\mathcal{O}(N_{\text{PI}} K^3)$ , where $N_{\text{PI}}$ is the number of iterations after which policy iteration converges.

   Each iteration of value iteration requires $\mathcal{O}\big(K^2\big)$ computations. Let $N_{\text{VI}}$ denote the number of iterations for value iteration (typically, the iteration stops when a stopping criterion is met). Then the complexity of approximately solving the dynamic program is $\mathcal{O}\big(N_{\text{VI}}K^2\big)$ .

   Note that the binary search needs to be repeated for each state. Thus, using binary search to compute the Whittle index to an accuracy of $\delta$ has a complexity of $\mathcal{O}\big(\log_2\big( \big(C_{\max} - C_{\min}\big)/\delta\big) N_{\text{PI}} K^4\big)$ if the dynamic program at each step is solved exactly and has a complexity of $\mathcal{O}\big(\log_2\big( \big(C_{\max} - C_{\min}\big)/\delta\big) N_{\text{VI}} K^3\big)$ if the dynamic program at each step is solved approximately.

4.2. Discussion on PCL-indexability

As mentioned earlier, an algorithm very similar to Algorithm 1 was proposed in [Reference Niño-Mora25] for computing the Whittle index for RBs that satisfy a technical condition known as PCL-indexability. The analysis in [Reference Niño-Mora25] is done under the assumption that the system starts from a designated start state distribution $\pi_0$ . For any policy g, define $\mathsf{N}^{(g)} = \sum_{x \in \mathcal{X}} N^{(g)}(x) \pi_0(x)$ , and define $\mathsf{D}^{(g)} = \sum_{x \in \mathcal{X}} D^{(g)}(x) \pi_0(x)$ . Let $x_1,\dots, x_K$ be a permutation of the state space such that the corresponding Whittle indices are non-decreasing: $\lambda_1 \le \dots \le \lambda_K$ . For any $k \in \{1, \dots, K\}$ , let $\bar{\mathcal{P}}_k$ denote the set $\{ x_1, \dots, x_k \}$ .

Now, for any $k \in \{1, \dots, K\}$ and all states $y \in \mathcal{X}\setminus \bar{\mathcal{P}}_k$ , define $\bar h_k = \bar g^{\big(\bar{\mathcal{P}}_k\big)}$ , $\bar h_{k,y} = \bar g^{\big(\bar{\mathcal{P}}_k \cup\{y\}\big)}$ , and define

(21) \begin{equation} \bar \mu_{k,y} = \frac{ \mathsf{D}^{\big(\bar h_{k,y}\big)} - \mathsf{D}^{\big(\bar h_k\big)} } {\mathsf{N}^{\big(\bar h_k\big)} - \mathsf{N}^{\big(\bar h_{k,y}\big)}}.\end{equation}

In [Reference Niño-Mora25] an algorithm called the adaptive greedy algorithm is presented to iteratively identify the sets $\bar{\mathcal{P}}_k$ and compute the corresponding Whittle indices. This algorithm is shown in Algorithm 3.

An RB is said to be PCL-indexable [Reference Niño-Mora25] if it satisfies the following conditions:

1. 1. For any ${\mathcal S} \subseteq {\mathcal X}$ and $y \in {\mathcal X}\backslash{\mathcal S}$ , we have $\mathsf{N}^{\big(\bar{g}^{({\mathcal S})}\big)} - \mathsf{N}^{\big(\bar{g}^{\left({\mathcal S} \cup \{y\}\right)}\big)} > 0$ .

2. 2. The sequence of index values produced by the adaptive greedy algorithm is monotonically non-decreasing.

Finally, the following result is established.

The main differences between our result and that of [Reference Niño-Mora25] are as follows:

1. 1. An implication of the first condition in the definition of PCL indexability is that the denominator in (21) is never zero. In contrast, we do not impose such a restriction and work with the non-empty subset of states for which the denominator in (14) is non-zero.

2. 2. In Algorithm 3, the sets $\{ \bar{\mathcal P}_k \}_{k=1}^{K}$ are constructed by adding states one by one, even when $\bar \mu_{k,y}$ has multiple argmins. In contrast, in Algorithm 1, the sets $\{ \mathcal{P}_d \}_{d=1}^{K_D}$ are constructed by adding all states which have the same Whittle index at once.

3. 3. In Algorithm 3, one has to check that the indices are generated in a non-decreasing order (which is the second condition of PCL-indexability). In contrast, in Algorithm 1, the indices are always generated in an increasing order, and therefore Condition 2 for PCL-indexability is always satisfied.

4. 4. Finally, Theorem 3 only guarantees that Algorithm 3 computes the Whittle index for RBs which satisfy PCL-indexability. Moreover, the second condition for PCL-indexability can only be checked after running Algorithm 3. In contrast, Theorem 2 guarantees that Algorithm 1 computes the Whittle index for all indexable RBs.

We conclude this discussion by revisiting an example from [Reference Niño-Mora25] which is an indexable RB that is not PCL-indexable. For this example, $\mathcal{X} = \{1, 2, 3\}$ ; the transition matrices are

\begin{align*} P(0) = \left[\begin{matrix} 0.3629 & \quad 0.5028 & \quad 0.1343 \\[4pt] 0.0823 & \quad 0.7534 & \quad 0.1643 \\[4pt] 0.2460 & \quad 0.0294 & \quad 0.7246 \end{matrix}\right], \qquad P(1) = \left[\begin{matrix} 0.1719 & \quad 0.1749 & \quad 0.6532 \\[4pt] 0.0547 & \quad 0.9317 & \quad 0.0136 \\[4pt] 0.1547 & \quad 0.6271 & \quad 0.2182\end{matrix}\right]; \end{align*}

the per-step cost is $c(x,0) = 0$ for all $x \in \mathcal{X}$ and $c(1, 1) = -0.44138$ , $c(2, 1) = -0.8033$ , $c(3, 1) = -0.14257$ ; and $\beta = 0.9$ . The corresponding Whittle indices are $[0.18, 0.8, 0.57]$ .

This model is not PCL-indexable, since if $g = [1, 1, 0]$ and $h = [0, 1, 0]$ , then $N^{(g)} = [5.66, 8.24, 4.23]$ and $N^{(h)} = [6.65, 8.59, 4.88]$ . Therefore, for any initial state distribution $\pi_0$ , $\mathsf{N}^{(g)} < \mathsf{N}^{(h)}$ .

However, as the problem is indexable, we can still apply Algorithm 1 to compute Whittle indices without any limitations. The steps are as follows:

1. 1. Initialize $d=0$ and have ${\mathcal P}_0 = \emptyset$ . Thus ${\bar g}^{({\mathcal P}_0)} = [1, 1, 1]$ and we compute $N^{\big({\bar g}^{({\mathcal P}_0)}\big)} = [10, 10, 10]$ and $D^{\big({\bar g}^{({\mathcal P}_0)}\big)} = [{-}6.43, -7.43, -6.51]$ .

2. 2. There are three possibilities for $y \in \mathcal{X} \setminus \mathcal{P}_0 = \{1,2,3\}$ :

• For $y = 1$ , $h_{0,1} = [0, 1, 1]$ . We compute $N^{(h_{0,1})} = [7.88, 9.29, 9.13]$ and $D^{(h_{0,1})} = $ $[{-}6.05, -7.30, -6.35]$ . Therefore, $\Lambda_{0,1} = \big\{x \in \mathcal{X} \,:\, N^{(\bar g(\mathcal{P}_0)}(x) \neq N^{(h_{0,1})}(x) \big\} = \{1, 2, 3\}$ . Now for each $x \in \Lambda_{0,1}$ , we compute $\mu_{0,1}(1) = \mu_{0,1}(2) = \mu_{0,1}(3) = 0.18$ . Therefore, $\mu^*_{0,1} = 0.18$ .

• For $y = 2$ , $h_{0,2} = [1, 0, 1]$ . We compute $N^{\left(h_{0,2}\right)} = [4.58, 2.93, 4.10]$ and $D^{\left(h_{0,2}\right)} = [{-}1.27, -0.7, -0.89]$ . Therefore, $\Lambda_{0,2} = \Big\{x \in \mathcal{X} \,:\, N^{(\bar g^{(\mathcal{P}_0)}}(x) \neq N^{\left(h_{0,2}\right)}(x) \Big\} = \{1, 2, 3\}$ . Now for each $x \in \Lambda_{0,2}$ , we compute $\mu_{0,2}(1) = \mu_{0,2}(2) = \mu_{0,2}(3) = 0.95$ . Therefore, $\mu^*_{0,2} = 0.95$ .

• For $y = 3$ , $h_{0,3} = [1, 1, 0]$ . We compute $N^{(h_{0,3})} = [5.66, 8.24, 4.23]$ and $D^{(h_{0,3})} = [{-}3.64, -6.30, -2.79]$ . Therefore, $\Lambda_{0,3} = \Big\{x \in \mathcal{X} \,:\, N^{(\bar g^{(\mathcal{P}_0)}}(x) \neq N^{(h_{0,3})}(x) \Big\} = \{1, 2, 3\}$ . Now for each $x \in \Lambda_{0,3}$ , we compute $\mu_{0,3}(1) = \mu_{0,3}(2) = \mu_{0,3}(3) = 0.64$ . Therefore, $\mu^*_{0,3} = 0.64$ .

1. Now $\lambda_1 = \min \big\{ \mu^*_{0,1}, \mu^*_{0,2}, \mu^*_{0,3} \big\} = 0.18$ . Therefore, ${\mathcal P}_1 = \{1\}$ , $w(1) = 0.18$ , ${\bar g}^{({\mathcal P}_1)} = [0, 1, 1]$ . We have already computed $N^{\big({\bar g}^{({\mathcal P}_1)}\big)} = [7.88, 9.29, 9.13]$ and $D^{\big({\bar g}^{({\mathcal P}_1)}\big)} = [{-}6.05, -7.30, -6.35]$ .

1. 3. There are two possibilities for $y \in \mathcal{X} \setminus \mathcal{P}_1 = \{2, 3\}$ :

• For $y = 2$ , $h_{1,2} = [0, 0, 1]$ . We compute $N^{(h_{1,2})} = [1.48, 1.52, 2.57]$ and $D^{(h_{1,2})} = [{-}0.21, -0.22, -0.37]$ . Therefore, $\Lambda_{1,2} = \Big\{x \in \mathcal{X} \,:\, N^{(\bar g^{(\mathcal{P}_1)}}(x) \neq N^{(h_{1,2})}(x) \Big\} = \{1, 2, 3\}$ . Now for each $x \in \Lambda_{1,2}$ , we compute $\mu_{1,2}(1) = \mu_{1,2}(2) = \mu_{1,2}(3) = 0.91$ . Therefore, $\mu^*_{1,2} = 0.91$ .

• For $y = 3$ , $h_{1,3} = [0, 1, 0]$ . We compute $N^{(h_{1,3})} = [6.65, 8.59, 4.88]$ and $D^{(h_{1,3})} = [{-}1.22, -0.66, -0.83]$ . Therefore, $\Lambda_{1,3} = \Big\{x \in \mathcal{X} \,:\, N^{(\bar g^{(\mathcal{P}_1)}}(x) \neq N^{(h_{1,3})}(x) \Big\} = \{1, 2, 3\}$ . Now for each $x \in \Lambda_{1,3}$ , we compute $\mu_{1,3}(1) = \mu_{1,3}(2) = \mu_{1,3}(3) = 0.57$ . Therefore, $\mu^*_{1,3} = 0.57$ .

1. Now $\lambda_2 = \min \big\{ \mu^*_{1,2}, \mu^*_{1,3} \big\} = 0.57$ . Therefore, ${\mathcal P}_2 = \{1, 3\}$ , $w(3) = 0.57$ , ${\bar g}^{({\mathcal P}_2)} = [0, 1, 0]$ . We have already computed $N^{\big({\bar g}^{({\mathcal P}_2)}\big)} = [6.65, 8.59, 4.88]$ and $D^{\big({\bar g}^{({\mathcal P}_2)}\big)} = [{-}1.22, -0.66, -0.83]$ .

1. 4. There is only one possibility for $y \in \mathcal{X} \setminus \mathcal{P}_2 = \{2\}$ :

• For $y = 3$ , $h_{3,2} = [0, 0, 0]$ , $N^{(h_{3,2})} = [0, 0, 0]$ , and $D^{(h_{3,2})} = [{-}0.21, -0.22, -0.37]$ . Therefore, $\Lambda_{3,2} = \Big\{x \in \mathcal{X} \,:\, N^{(\bar g^{(\mathcal{P}_2)}}(x) \neq N^{(h_{3,2})}(x) \Big\} = \{1, 2, 3\}$ . Now for each $x \in \Lambda_{3,2}$ , we compute $\mu_{3,2}(1) = \mu_{3,2}(2) = \mu_{3,2}(3) = 0.8$ . Therefore, $\mu^*_{3,2} = 0.8$ .

1. Now $\mu^*_{3,2} = 0.57$ . Therefore, ${\mathcal P}_3 = \{1, 2, 3\}$ and $w(2) = 0.8$ .

Finally, the Whittle indices are $[0.18, 0.8, 0.57]$ .

5.1. Restless bandits with optimal threshold-based policy

Consider an RB $\big(\mathcal{X}, \{0, 1\}, \{P(a)\}_{a \in \{0, 1\}},c, x_0\big)$ where the state space $\mathcal{X}$ is a totally ordered set, i.e., ${\mathcal X} = \{1, \ldots, K\}$ . Let $\mathcal{X}_0 = \{0, \dots, K\}$ ; let $\mathcal{X}_{\ge \ell}$ denote the set of states greater than or equal to the state $\ell$ , and let $\mathcal{X}_{\le \ell}$ denote the set of states less than or equal to the state $\ell$ . We suppose that the model satisfies the following assumption:

1. (P) There exists a non-decreasing family of thresholds $\{\ell_\lambda\}_{\lambda \in \mathbb{R}}$ , $\ell_\lambda \in \mathcal{X}_0$ , such that the threshold policy $g^{(\ell_\lambda)}$ is optimal for Problem 2 with activation cost $\lambda$ .

Several models where (P) holds have been considered in the literature [Reference Akbarzadeh and Mahajan3, Reference Ansell, Glazebrook, Niño-Mora and O’Keeffe5, Reference Avrachenkov, Ayesta, Doncel and Jacko7, Reference Glazebrook, Hodge and Kirkbride15, Reference Glazebrook, Kirkbride and Ouenniche17, Reference Wang, Ren, Mo and Shi30]. A key implication of Property (P) is the following.

As in Section 4, we assume that there are $K_D (\leq K)$ distinct Whittle indices given by $\Lambda^* = \{\lambda_1, \dots, \lambda_{K_D}\}$ , where $\lambda_1 < \lambda_2 < \dots \lambda_{K_D}$ . We also let $\lambda_0 = -\infty$ , and for any $d \in \{0, \dots, K_D\}$ , we let $\mathcal{P}_d = \{ x \in \mathcal X \,:\, w(x) \le \lambda_d \}$ . As stated in the proof of Lemma 4, Property (P) implies that $\mathcal{P}_d = \mathcal{X}_{\le \ell_{\lambda_d}}$ . Therefore, $\Gamma_{d+1} = \big\{\ell_{\lambda_d} + 1, \dots, \ell_{\lambda_{d+1}}\big\}$ . Thus, Theorem 2 simplifies as follows.

Thus, based on Corollary 1, for models that satisfy Property (P), we can simplify Algorithm 2 as shown in Algorithm 4. Instead of computing $\mu^*_{d,y}$ for all $y \in \mathcal{X} \setminus \mathcal{P}_d$ , we can compute it sequentially and break the loop when $\mu^*_{d,y} \neq \lambda_{d+1}$ . Note that this simplification does not change the asymptotic complexity of the algorithm, which is still $\mathcal{O}\!\left(K^3\right)$ .

Remark 1. Note that if the model satisfies additional assumptions such that it is known up front that no two states have the same Whittle index, then we do not need the inner for loop (over y) in Algorithm 4, and can simply compute the Whittle index of the state $\ell$ as

\begin{align*} w(\ell) = \min \frac{D^{\big(\bar g^{\left(\mathcal{X}_{\le \ell + 1}\right)}\big)}(x) - D^{\big(g^{\left(\mathcal{X}_{\le \ell }\right)}\big)}(x)} {N^{\big(\bar g^{\left(\mathcal{X}_{\le \ell}\right)}\big)}(x) - N^{\big(g^{\left(\mathcal{X}_{\le \ell + 1 }\right)}\big)}(x)}, \end{align*}

where the minimum is over all x such that the denominator is non-zero.

In the next section, we present a new model called stochastic monotone bandits, which may be considered as a generalization of monotone bandits [Reference Ansell, Glazebrook, Niño-Mora and O’Keeffe5, Reference Avrachenkov, Ayesta, Doncel and Jacko7, Reference Glazebrook, Hodge and Kirkbride15], and show that these models satisfy Property (P).

5.2. Stochastic monotone bandits

We first start with the definition of a submodular function. Given ordered sets $\mathcal{X}$ and $\mathcal{Y}$ , a function $f\,:\, \mathcal{X} \times \mathcal{Y} \to \mathbb{R}$ is called submodular if for any $x_1, x_2 \in \mathcal{X}$ and $y_1, y_2 \in \mathcal{Y}$ such that $x_2 \geq x_1$ and $y_2 \geq y_1$ , we have $f(x_1, y_2) - f(x_1, y_1) \geq f(x_2, y_2) - f(x_2, y_1)$ .

We say that the RB is stochastic monotone if it satisfies the following conditions:

1. (D1) For any $a \in \{0, 1\}$ , P(a) is stochastically monotone; i.e., for any $x, y \in \mathcal{X}$ such that $x < y$ , we have

   \begin{align*}\sum_{w \in \mathcal X_{\ge z}} P_{xw}(a) \leq \sum_{w \in \mathcal X_{\ge z}} P_{yw}(a)\end{align*}
   for any $z \in \mathcal{X}$ .

2. (D2) For any $z \in \mathcal{X}$ , $S_{zx}(a) \,:\!=\, \sum_{w \in \mathcal X_{\ge z}} P_{xw}(a)$ is submodular in (x, a).

3. (D3) For any $a \in \{ 0, 1 \}$ , c(x, a) is non-decreasing in x.

4. (D4) The cost c(x, a) is submodular in (x, a).

For ease of notation, for any $\ell \in \mathcal{X}_0$ , we let $g^{(\ell)} = \bar g^{(\mathcal X_{\le \ell})}$ denote a policy with threshold $\ell$ (where $\bar{g}^{(\mathcal S)}$ is as defined in (13)).

Proof. For the first part, we note that Conditions (D1)–(D4) are the same as the properties of [Reference Puterman28, Theorem 4.7.4], which implies that there exists a threshold as stated.

For the second part, we first show that for any $\ell \in \mathcal{X}^*$ , $J^{\big(g^{(\ell)}\big)}_\lambda(x)$ is submodular in $(\ell,\lambda)$ for all $x \in \mathcal{X}$ . In particular, for any $k < \ell$ , we have

\begin{align*} J^{\big(g^{(\ell)}\big)}_\lambda(x) - J^{\big(g^{(k)}\big)}_\lambda(x) = D^{\big(g^{(\ell)}\big)}_\lambda(x) - D^{\big(g^{(k)}\big)}_\lambda(x) + \lambda( N^{\big(g^{(\ell)}\big)}_\lambda(x) - N^{\big(g^{(k)}\big)}_\lambda(x)). \end{align*}

Now (D5) implies that the difference $J^{\big(g^{(\ell)}\big)}_\lambda(x) - J^{\big(g^{(k)}\big)}_\lambda(x)$ is non-increasing in $\lambda$ . Therefore, $J^{\big(g^{(\ell)}\big)}_\lambda(x)$ is submodular in $(\ell,\lambda)$ . Consequently, from [Reference Puterman28, Theorem 2.8.2],

\begin{align*}\ell_\lambda = \max \bigg\{ \ell^{\prime} \in \arg\min_{\ell \in \mathcal{X}^*} J^{\big(g^{(\ell)}\big)}_\lambda(x) \bigg\}\end{align*}

is non-decreasing in $\lambda$ .

5.3. Restless bandits with controlled restarts

Consider restless bandits with controlled restarts (i.e., models where $P_{xy}(1)$ does not depend on x). By Proposition 2(c), such models are indexable. In this section, we explain how to simplify the computation of the Whittle index for such models. For ease of notation, we use $P_{xy}$ to denote $P_{xy}(0)$ and $Q_y$ to denote $P_{xy}(1)$ .

Define $ \mathsf{D}^{(g)} = \sum_{x \in \mathcal{X}} Q_x D^{(g)}(x) \quad\text{and}\quad \mathsf{N}^{(g)} = \sum_{x \in \mathcal{X}} Q_x N^{(g)}(x).$ Now, following the discussion of Section 4, we can show that the result of Theorem 2 continues to hold when $\mu_{d,y}$ is replaced by

\begin{align*} \hat{\mu}_{d,y} = \dfrac{\mathsf{D}^{(h_y)} - \mathsf{D}^{\big(\bar{g}^{({\mathcal P}_d)}\big)}} {\mathsf{N}^{\big(\bar{g}^{({\mathcal P}_d)}\big)} - \mathsf{N}^{(h_y)}}. \end{align*}

Therefore, we can replace $\mu_{d,y}(x)$ in Algorithm 1 by $\hat{\mu}_{d,y}$ . Our key result for this section is that $\mathsf{D}^{(g)}$ and $\mathsf{N}^{(g)}$ can be computed efficiently for models with controlled restarts.

For that matter, given any policy g, let $\tau_g$ denote the hitting time of the set $\Pi^{(g)} = \{x \in \mathcal{X} \,:\, g(x) = 1 \}$ . Let

\begin{align*} \mathsf{L}^{(g)} \,:\!=\, \mathbb{E}\Biggl[ \sum_{t=0}^{\tau_g} \beta^t c(X_t, g(X_t)) \Bigm| X_0 \sim Q \Biggr] \quad\text{and}\quad \mathsf{M}^{(g)} \,:\!=\, \mathbb{E}\Biggl[ \sum_{t=0}^{\tau_g} \beta^t \Bigm| X_0 \sim Q \Biggr]\end{align*}

denote the expected discounted cost and expected discounted time for hitting $\Pi^{(g)}$ starting with an initial state distribution of Q. Then, using ideas from renewal theory, we can show the following.

Theorem 4. For any policy g,

\begin{equation*} \mathsf{D}^{(g)} = \dfrac{\mathsf{L}^{(g)}}{\mathsf{M}^{(g)}} \quad\text{and}\quad \mathsf{N}^{(g)} = \dfrac{1}{\beta \mathsf{M}^{(g)}} - \dfrac{1-\beta}{\beta}. \end{equation*}

Proof. The proof follows from standard ideas in renewal theory. By the strong Markov property, we have

(22) \begin{align} \mathsf{D}^{(g)} &= \mathbb{E} \Bigg[ (1-\beta) \sum_{t = 0}^{\tau_g} \beta^{t} c(X_t, g(X_t)) + \beta^{\tau_g+1} \mathsf{D}^{(g)} \Bigm| X_0 \sim Q \Bigg] \nonumber \\ &= (1-\beta) \mathsf{L}^{(g)} + \mathbb{E}\big[ \beta^{\tau_g+1} | X_0 \sim Q \big] \mathsf{D}^{(g)}. \end{align}

Using the definition of $\mathsf{M}^{(g)}$ , we have $\mathbb{E}\big[ \beta^{\tau_g+1} | X_0 \sim Q \big] = 1 - (1-\beta) \mathsf{M}^{(g)}$ . Substituting this in (22) and rearranging the terms we get $\mathsf{D}^{(g)} = \mathsf{L}^{(g)}/\mathsf{M}^{(g)}$ .

For $\mathsf{N}^{(g)}$ , by the strong Markov property we have

\begin{align*} \mathsf{N}^{(g)} &= \mathbb{E}\Big[ (1-\beta) \beta^{\tau_g} + \beta^{\tau_g+1} \mathsf{N}^{(g)} \Bigm| X_0 \sim Q \Big] \nonumber \\ &= \mathbb{E}\big[ \beta^{\tau_g} | X_0 \sim Q \big] \big(1-\beta + \beta \mathsf{N}^{(g)}\big) = \dfrac{ 1-(1-\beta)\mathsf{M}^{(g)} }{\beta}\big(1-\beta + \beta \mathsf{N}^{(g)}\big).\end{align*}

Therefore, we get $\mathsf{N}^{(g)} = \bigl(1 - (1-\beta)\mathsf{M}^{(g)}\bigr) / \beta \mathsf{M}^{(g)}$ .

Given any policy g, we can efficiently compute $\mathsf{L}^{(g)}$ and $\mathsf{M}^{(g)}$ using standard formulas for truncated Markov chains. For any vector v, let $v^{(g)}$ denote the vector with components indexed by the set $\{ x \in \mathcal{X} \,:\, g(x) = 0 \}$ , and let $\tilde v^{(g)}$ denote the remaining components. For example, if $\mathcal{X} = \{1, 2, 3, 4\}$ , $g = (1,0, 1, 0)$ , and $v = [1, 2, 3, 4]$ , then $v^{(g)} = (2, 4)$ and $\tilde v^{(g)}= (1,3)$ . Similarly, for any square matrix Z, let $Z^{[g]}$ denote the square sub-matrix corresponding to the elements $\{ x \in \mathcal{X} \,:\, g(x) = 0\}$ , and let $\tilde Z^{[g]}$ denote the sub-matrix with rows $\{x \in \mathcal{X}\,:\, g(x) = 0 \}$ and columns $\{x \in \mathcal{X} \,:\, g(x) = 1 \}$ . For example, if $g = [1, 0, 1, 0]$ and if

\begin{align*}Z = \left[\begin{matrix} 1 & \quad 2 & \quad 3 & \quad 4 \\[3pt] 5 & \quad 6 & \quad 8 & \quad 8 \\[3pt] 9 & \quad 10 & \quad 11 & \quad 12 \\[3pt] 13 & \quad 14 & \quad 15 & \quad 16 \end{matrix}\right],\end{align*}

then $Z^{[g]} = \left[\begin{matrix} 6 & \quad 8 \\[3pt] 14 & \quad 16 \end{matrix}\right]$ and $Z^{[g]} = \left[\begin{matrix} 5 & \quad 8 \\[3pt] 13 & \quad 15 \end{matrix}\right]$ . Then, from standard formulas for truncated Markov chains, we have the following.

Proposition 3. For any policy g, let $c_0$ and $c_1$ denote column vectors corresponding to $c(\cdot, 0)$ and $c(\cdot, 1)$ . Then

\begin{align*} \mathsf L^{(g)} &= Q^{(g)} \big(I - \beta P^{[g]}\big)^{-1} \Big(c_0^{(g)} + \beta \tilde P^{[g]} \tilde c^{(g)}_1\Big) + \tilde Q^{(g)} \tilde c^{(g)}_1, \\ \mathsf M^{(g)} &= Q^{(g)} \big(I - \beta P^{[g]}\big)^{-1} \Big(\textbf{1}^{(g)} + \beta \tilde P^{[g]} \tilde{\textbf{1}}^{(g)}\Big) + \tilde Q^{(g)} \tilde {\textbf{1}}^{(g)}. \end{align*}

This gives us an efficient method to compute $\mathsf{L}^{(g)}$ and $\mathsf{M}^{(g)}$ , which can in turn be used to compute $\mathsf{D}^{(g)}$ and $\mathsf{N}^{(g)}$ and used in a modified version of Algorithm 1 as explained.

6.1. Experimental setup

In our experiments, we consider restart bandits with $P(1)= [\textbf{1},\textbf{0}, \dots, \textbf{0}]$ . There are two other components of the model: the transition matrix P(0) and the cost function c. We choose these components as follows.

6.2. Experimental details and results

We conduct different experiments to compare the performance of the Whittle index with the optimal policy and the myopic policy for different setups (described in Section 6.1) and for different choices of the size $|\mathcal{X}|$ of the state space, the number n of arms, and the number m of active arms. For all experiments we choose the discount factor $\beta = 0.95$ .

We evaluate the performance of a policy via Monte Carlo simulations over S trajectories, where each trajectory is of length T. In all our experiments, we choose $S= 2500$ and $T = 250$ .

Experiment 3. Comparison of Whittle index with the myopic policy for structured models. We generate the structured models as in Experiment 1 but for $|\mathcal{X}| = 25$ , $n \in \{25, 50, 75\}$ , and $m \in \{1, 2, 5\}$ . In this case, let $\varepsilon_{\small\text{MYP}} = (J({\small\text{MYP}})-J({\small\text{WIP}}))/J({\small\text{MYP}})$ denote the relative improvement of wip compared to myp. The results of $\varepsilon_{\small\text{MYP}}$ for different choices of the parameters are shown in Figure 3.

Figure 2. Relative performance $\alpha_{\small\text{OPT}}$ of wip versus opt for Experiment 2.

Figure 3. Relative improvement $\varepsilon_{\small\text{MYP}}$ of wip versus myp for Experiment 3.

In Figure 3, we observe that wip performs considerably better than myp. In addition to that, the performance of wip is better with respect to myp when $\ell = 4$ , which is more complicated than models where $\ell \in \{1, 2, 3\}$ . However, increasing m does not necessarily contribute to better $\varepsilon_{\small\text{MYP}}$ , as overlap between the choices of the two policies may increase. Note that as ${\mathcal P}_4({\cdot})$ is very different from the rest of the models, the trend of bars in Figure 3d with respect to n varies differently from the rest of the models.

Experiment 4. Comparison of Whittle index with the myopic policy for randomly sampled models. We generate 250 random models as described in Experiment 2 but for $|\mathcal{X}|= 25$ and larger values of n. For each case, $\varepsilon_{\small\text{MYP}}$ is computed. Histograms of $\varepsilon_{\small\text{MYP}}$ for different choices of the parameters are shown in Figure 4.

Figure 4. Relative improvement $\varepsilon_{\small\text{MYP}}$ of wip versus myp for Experiment 4.

The results show that on average, wip performs considerably better than myp, and this improvement is guaranteed as the concentration of data for the sampled models is mostly on positive values of $\varepsilon_{\small\text{MYP}}$ .

Figure 5. Relative improvement $\varepsilon_{\small\text{MYP}}$ of wip versus myp for Experiment 5.