Nicolas Gast
email: nicolas.gast@inria.fr

Bruno Gaujal
email: bruno.gaujal@inria.fr

Kimang Khun
email: khun.kimang@gmail.com

Testing Indexability and Computing Whittle and Gittins Index in Subcubic Time

Keywords: Whittle Index, Gittins Index, Restless Bandit, Multi-armed Bandit, Sherman-Morrison, Markov Decision Process, Fast Matrix Multiplication

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Markovian bandits form a subclass of multi-armed bandit problems in which each arm has an internal state that evolves over time in a Markovian manner, as a function of the decision maker's actions. In such a problem, at each time step, the decision maker observes the state of all arms and chooses which one to activate. When the state of an arm evolves only when this arm is chosen, one falls into the category of rested Markovian bandits for which an optimal policy (in the discounted case) was found by Gittins [START_REF] Gittins | Bandit processes and dynamic allocation indices[END_REF]. When the state of an arm can also evolve when the arm is not chosen, the problem is called a restless bandit problem, and computing an optimal policy is computationally difficult [START_REF] Papadimitriou | The complexity of optimal queueing network control[END_REF].

In his seminal paper [START_REF] Whittle | Restless bandits: Activity allocation in a changing world[END_REF], Whittle proposed a very efficient heuristic: For each arm, an index function maps each state of the arm to a real number. The Whittle index policy then consists in activating the arms having the highest index first. This heuristic generalizes Gittins index to restless bandits. Contrary to the rested case, the Whittle index policy is in general not optimal. Yet, this policy has been proven to be very efficient over the years: up to a condition called indexability, Whittle index has been shown to be (in the undiscounted case) asymptotically optimal as the number of arms grows to infinity under certain technical assumptions [START_REF] Verloop | Asymptotically optimal priority policies for indexable and nonindexable restless bandits[END_REF][START_REF] Lott | On the optimality of an index rule in multichannel allocation for single-hop mobile networks with multiple service classes[END_REF][START_REF] Weber | On an index policy for restless bandits[END_REF]. Moreover, the heuristic performs extremely well in practice [START_REF] Glazebrook | Some indexable families of restless bandit problems[END_REF][START_REF] Ansell | Whittle's index policy for a multi-class queueing system with convex holding costs[END_REF][START_REF] Glazebrook | An index policy for a stochastic scheduling model with improving/deteriorating jobs[END_REF]. Restless bandits and Whittle index have been applied to many scheduling and resource allocation problems such as wireless communication [START_REF] Aalto | Whittle index approach to opportunistic scheduling with partial channel information[END_REF][START_REF] Liu | Indexability of restless bandit problems and optimality of Whittle index for dynamic multichannel access[END_REF], web crawling [START_REF] Avrachenkov | Whittle index based Q-learning for restless bandits with average reward[END_REF][START_REF] Niño-Mora | A dynamic page-refresh index policy for web crawlers[END_REF], congestion control [START_REF] Avrachenkov | Congestion control of tcp flows in internet routers by means of index policy[END_REF][START_REF] Avrachenkov | Impulsive control for G-AIMD dynamics with relaxed and hard constraints[END_REF], queueing systems [START_REF] Scully | SOAP: One clean analysis of all age-based scheduling policies[END_REF][START_REF] Aalto | Properties of the Gittins index with application to optimal scheduling[END_REF][START_REF] Aalto | On the Gittins index in the M/G/1 queue[END_REF][START_REF] Borkar | Whittle indexability in egalitarian processor sharing systems[END_REF][START_REF] Larrañaga | Asymptotically optimal index policies for an abandonment queue with convex holding cost[END_REF][START_REF] Archibald | Indexability and index heuristics for a simple class of inventory routing problems[END_REF][START_REF] Glazebrook | Index policies for the admission control and routing of impatient customers to heterogeneous service stations[END_REF], and clinical trials [START_REF] Villar | Multi-armed bandit models for the optimal design of clinical trials: benefits and challenges[END_REF].

The above examples show that, when a problem is indexable, computing Whittle index is a very efficient way to construct a nearly-optimal heuristic. This raises a few important questions, that we study in this paper:

• Is testing indexability computationally hard?

• Is there an efficient algorithm to compute Whittle index?

• Is Whittle index harder to compute than Gittins index? Related work. The computation of Gittins index has received a lot of attention in the past, see for instance [START_REF] Chen | Linear programming for finite state multiarmed bandit problems[END_REF][START_REF] Katehakis | The multi-armed bandit problem: decomposition and computation[END_REF][START_REF] Niño-Mora | A (2/3)n 3 fast-pivoting algorithm for the Gittins index and optimal stopping of a markov chain[END_REF][START_REF] Sonin | A generalized Gittins index for a markov chain and its recursive calculation[END_REF] and the recent survey [START_REF] Chakravorty | Multi-armed bandits, Gittins index, and its calculation. Methods and applications of statistics in clinical trials: Planning, analysis, and inferential methods[END_REF]. For a n-state arm, the algorithms having the smallest complexity perform (2/3)n 3 + O(n 2 ) arithmetic operations [START_REF] Chakravorty | Multi-armed bandits, Gittins index, and its calculation. Methods and applications of statistics in clinical trials: Planning, analysis, and inferential methods[END_REF]. Note that in page 4 of [START_REF] Niño-Mora | A fast-pivoting algorithm for Whittle's restless bandit index[END_REF] the author claims that it is unlikely that this complexity can be improved. As we see later, we do improve upon this complexity.

Concerning Whittle index, to the best of our knowledge, there are very few efficient general-purpose algorithms to test indexability, see e.g. [START_REF] Niño-Mora | Characterization and computation of restless bandit marginal productivity indices[END_REF], and most papers studying Whittle index either assume that the studied model is indexable or focus on specific classes of restless bandits for which the structure of arms can be used to show indexability, see e.g. [START_REF] Aalto | Properties of the Gittins index with application to optimal scheduling[END_REF][START_REF] Borkar | Whittle indexability in egalitarian processor sharing systems[END_REF][START_REF] Akbarzadeh | Restless bandits with controlled restarts: Indexability and computation of Whittle index[END_REF][START_REF] Akbarzadeh | Maintenance of a collection of machines under partial observability: Indexability and computation of Whittle index[END_REF]. Assuming indexability, the computation of Whittle index has been considered by a few papers.

The most efficient numerical algorithm to compute Whittle index is recently presented in [START_REF] Niño-Mora | A fast-pivoting algorithm for Whittle's restless bandit index[END_REF]. This algorithm, called fast-pivoting algorithm, performs (2/3)n 3 + O(n 2 ) arithmetic operations 1 if the initialization phase is excluded from the count. This is done by using the parametric simplex method and exploiting the special structure of this linear system to reduce the complexity of simplex pivoting steps. This fast-pivoting algorithm is an efficient implementation of adaptive-greedy algorithm [START_REF] Niño-Mora | Dynamic priority allocation via restless bandit marginal productivity indices[END_REF]. Based on a geometric interpretation of Whittle index, the authors in [START_REF] Akbarzadeh | Conditions for indexability of restless bandits and an O(K 3 ) algorithm to compute Whittle index[END_REF] propose a refinement of the adaptive-greedy algorithm of [START_REF] Niño-Mora | Dynamic priority allocation via restless bandit marginal productivity indices[END_REF] to compute Whittle indices of all indexable restless bandits. For a n-state arm, the refined algorithm achieves a O(n 3 ) complexity by using the Sherman-Morrison formula. The authors also propose a few checkable conditions to test indexability. However, those conditions are not necessary for indexability, which means that if an arm does not verify the conditions, we cannot conclude that the arm is non-indexable and an algorithm to check indexability is still needed. Also, no detailed description is given for adapting those conditions and their algorithm to restless bandit without discount. A thorough comparison between our algorithm and [START_REF] Niño-Mora | A fast-pivoting algorithm for Whittle's restless bandit index[END_REF][START_REF] Akbarzadeh | Conditions for indexability of restless bandits and an O(K 3 ) algorithm to compute Whittle index[END_REF] is given in Appendix E. For continuous-time n-state restless bandits, the work of [START_REF] Ayesta | On the computation of Whittle's index for markovian restless bandits[END_REF] proposes an algorithm to check indexability and compute Whittle index with a complexity exponential in the number of states n of each arm. According to Remark 4.1 of that paper, this complexity can be reduced to O(n 5 ) if the restless bandit is known to be indexable and threshold-based policies are optimal. It is stated that their approach is not applicable for discounted restless bandits.

While computing the Whittle indices of a known arm's model is still a challenge, there is interesting work in trying to learn Whittle index when only the arm's simulator is given and the arm's model is unknown. For instance, [START_REF] Avrachenkov | Whittle index based Q-learning for restless bandits with average reward[END_REF][START_REF] Gibson | A novel implementation of Qlearning for the Whittle index[END_REF][START_REF] Fu | Towards Q-learning the Whittle index for restless bandits[END_REF] use Q-learning algorithm to estimate Whittle index as time evolves in finite-state restless bandits. Moreover, the work of [START_REF] Nakhleh | Neur-WIN: Neural Whittle index network for restless bandits via deep RL[END_REF] uses deep reinforcement learning framework to estimate Whittle indices of the arms with large state space or convoluted transition kernel, assuming a notion of strong indexability. For learning aspect, the work of [START_REF] Gibson | A novel implementation of Qlearning for the Whittle index[END_REF] shows as to learn Whittle index in non-discounted case by maintaining two Q-functions, updating them using Q-learning algorithm, and deducing Whittle index from them when needed. The way that Whittle indices are computed is very close to our work but less efficient than our algorithm since the authors are more interested in learning the index.

Contributions.

In this paper, we investigate Whittle index computation in restless multi-armed bandit problems and present four main contributions.

Our first contribution is to discuss the ambiguities in the classical definitions of indexability. Classical definitions assume that an arm is indexable if the optimal policy is a non-decreasing function of some penalty term λ. While this definition works for most practical cases, it is not always precise enough 1 multiplications and additions of real numbers, regardless of their values because the optimal policy is in general not unique. In our definition, we specify the notion of increasingness that should be used. Our definition guarantees the uniqueness of Whittle indices. Note that our definition is the same as the one used in some recent papers (e.g., [START_REF] Niño-Mora | A fast-pivoting algorithm for Whittle's restless bandit index[END_REF]), but the ambiguity of the classical definition seems rarely mentioned.

Our second contribution is to propose a unified algorithm that computes the Whittle indices for both discounted and non-discounted restless bandits. Our algorithm, which can be viewed as a refinement of the algorithm in [START_REF] Akbarzadeh | Conditions for indexability of restless bandits and an O(K 3 ) algorithm to compute Whittle index[END_REF], tests whether the input arm is indexable or not, and computes Whittle index if the arm is. As a byproduct, our algorithm can compute Gittins index in rested bandits which are a subclass of restless bandits. This algorithm computes the indices in increasing order, and relies on an efficient use of the Sherman-Morrison formula to compute Whittle index in (2/3)n 3 + O(n 2 ) plus subcubic time [START_REF] Strassen | Gaussian elimination is not optimal[END_REF] to solve a linear system of order n. This algorithm can detect on the fly if a computed index violates the indexability condition, which adds an extra (1/3)n 3 + O(n 2 ) arithmetic operations. This later test is optional: the complexity of our algorithm is n 3 + o(n 3 ) when testing indexability and (2/3)n 3 + o(n 3 ) without the test. These two complexities are comparable to the ones excluding the common initialization phase of reduced-pivoting indexability (RPI) and fast-pivoting adaptive greedy (FPAG) algorithms in [START_REF] Niño-Mora | Characterization and computation of restless bandit marginal productivity indices[END_REF]. For discounted problems, our algorithm works for any finite-state arm. For non-discounted problems, our algorithm takes as an input any arm and can output three results: the arm is indexable, non-indexable, or multichain. We show the correctness of the algorithm which proves that for unichain arms, our algorithm have soundness and completeness properties. The possible outputs of our algorithm are summarized in Figure 1. Our third contribution is to show how to reduce the complexity of the above algorithm to obtain the first subcubic algorithm to compute Whittle index. This improvement is made possible by the fact that a linear system can be solved in subcubic time. By carefully reordering the computations, we show that it is possible to reduce the use of the Sherman-Morrison formula at the price of solving more linear systems. The subcubic complexity comes by striking a good balance between having too many or too few linear systems to solve. By using the current fastest matrix multiplication method, our algorithm can test indexability and compute Whittle index in O(n 2.5286 ). Our algorithm is also the first subcubic algorithm to compute Gittins index.

Indexable

Our fourth and last contribution is to provide an open-source implementation of our algorithm in Python with Numba, and to present an empirical evaluation of the performance of our implementation. Our results show that our algorithm is very efficient in computing Whittle index and testing indexability. Moreover, our simulations indicate that the subcubic version of our algorithm not only has an asymptotically small complexity but is also faster in practice than our original (2/3)n 3 algorithm. Testing the indexability and computing indices takes less than one second for n = 1000 states and less than 10 minutes for n = 15000 states. This is 15 to 20 times faster than the original computation times reported in [START_REF] Niño-Mora | A fast-pivoting algorithm for Whittle's restless bandit index[END_REF] (for a Matlab implementation), and about 5 times faster than an optimized implementation of [START_REF] Niño-Mora | A fast-pivoting algorithm for Whittle's restless bandit index[END_REF] (for a Julia implementation).

Road map.

The paper is organized as follows. We introduce the problem and the definition of Whittle index in Section 2. In Section 3, we characterize Whittle index and provide a general idea of how to compute Whittle index. Then, we show, in Section 4, how to use the Sherman-Morrison formula to compute the indices efficiently. We then show how to reduce the complexity of the algorithm by using fast matrix multiplication method in Section 5. We compare the numerical result of different variants of our algorithm in Section 6. We show how to adapt this approach to the discounted case in Section 7. Finally, we conclude in Section 8.

Restless bandits and indexability

Restless bandit arms and multi-armed bandit

In this paper, a restless bandit arm (that we denote later by RB) is a Markov decision process (MDP) with discrete state space [n] := {1, . . . , n} and binary action space {0, 1}, where 0 denotes the action "rest" and 1 denotes the action "activate". The time is discrete and the evolution is Markovian: If the MDP is in state i and action a is chosen, the decision maker earns an instantaneous reward r a i and the arm transitions to a new state j with probability P a ij . We denote this MDP by the pair (r, P ). The name "restless" comes from the fact that an arm put at rest may still transition to a new state.

A restless multi-armed bandit (RMAB) problem is a finite collection of M ∈ N * independent RB arms. At time t, the decision maker observes the state of all arms, and can choose up to m arms to activate, where m ≤ M is a fixed constant. The decision maker then earns a reward that is the sum of the rewards of all arms. The objective of the decision maker is to identify an allocation rule that maximizes the average reward earned over an infinite number of time steps. Such a problem is notoriously difficult to solve, as its complexity grows exponentially with the number of arms [START_REF] Papadimitriou | The complexity of optimal queueing network control[END_REF]. In his seminal paper [START_REF] Whittle | Restless bandits: Activity allocation in a changing world[END_REF], Whittle proposes the following approach: if all arms verify a technical condition known as indexability, then each state i of each arm is associated with a real number λ i , that is now known as the Whittle index of state i. At each time, the decision maker activates the m arms whose Whittle index of their current states are the m greatest indices. As mentioned earlier, this heuristic performs extremely well in practice, see e.g., [START_REF] Glazebrook | Some indexable families of restless bandit problems[END_REF][START_REF] Ansell | Whittle's index policy for a multi-class queueing system with convex holding costs[END_REF][START_REF] Glazebrook | An index policy for a stochastic scheduling model with improving/deteriorating jobs[END_REF]. This shows that if the M arms are all indexable, then one can derive a very efficient allocation rule for RMAB problems by computing the Whittle indices of all arms. The Whittle indices of an arm do not depend on the other arms. This shows that the computational cost of Whittle index policy is linear in the number of arms multiplied by the time to compute the indices for a single arm. Hence, in the remaining of the paper, we focus on a single arm and present a new algorithm to test indexability and compute the index of a given arm.

Indexability and Whittle index

For the remaining of the paper, we consider a single arm (r, P ) that has n states. In this section, we introduce the notion of indexability and discuss some ambiguities that we have found when using the definition of indexability defined in previous works.

Policy and arm structure

An arm is a two-action MDP. Hence, a policy π is a subset of the state space, π ⊆ [n], such that the policy chooses to activate the arm in state i if i ∈ π. We say that π is the set of active states, and we say that state i is passive if i ∈ π. By abuse of notation, we will write π i = 1 if i ∈ π and π i = 0 if i ∈ π, and we denote by P π the transition matrix corresponding to the policy π, i.e., P π ij = P πi ij . Following the classical definitions in the literature [START_REF] Puterman | Markov Decision Processes: Discrete Stochastic Dynamic Programming[END_REF], we say that: • A policy π is unichain if the transition matrix P π induced by π has a unique recurrent class. A policy that is not unichain is called multichain. • An arm is unichain if all policies π ⊆ [n] are unichain. An arm is multichain if it is not unichain, i.e., if there exists a policy π ⊆ [n] that is multichain.

Gain optimality and Bellman optimality

Following a policy π, we denote by g π i the long-run average reward that a decision maker would obtain when starting in state i. In the remainder of the paper, we use the term "gain" to denote the long-run average reward. Let g * i = max π g π i be the maximal gain starting from state i. From [40, Chapter 9], g * is uniquely defined. We say that a policy π is gain optimal if g π i = g * i for all state i.

It is shown in [START_REF] Puterman | Markov Decision Processes: Discrete Stochastic Dynamic Programming[END_REF]Chapter 9] that g * is the optimal gain if and only if there exists a vector h * , called optimal bias vector that satisfies the Bellman optimality equations: for all i ∈ [n],

g * i = max a∈{0,1} n j=1 P a ij g * j ( 1
)
g * i + h * i = max a∈{0,1} r a i + n j=1 P a ij h * j . ( 2
)
We say that a policy π is Bellman optimal if there exists2 a bias vector h * that is a solution of (2) and such that π attains the maximum in (2), i.e.: for all i,

n j=1 P πi ij g * j = g * i and π i ∈ arg max a∈{0,1} r a i + n j=1 P a ij h * j . (3)
The notion of Bellman optimality is stronger than the notion of gain optimality: A Bellman optimal policy is gain optimal, but the converse is not true in general. Note that the distinction between gain optimal and Bellman optimal policies is only important for the average reward criterion. This distinction disappears for the discounted case that we discuss in Section 7. The notion of Bellman optimality is equivalent to the notion of canonical optimality, that characterize policies that are optimal for any finite horizon, see [? ].

λ-penalized MDP and definition of indexability

For each λ ∈ R, we define a λ-penalized MDP3 whose transition matrices are the same as in the original MDP and whose reward at time t ≥ 0 when taking action a t in state s t is r at st -λa t . The quantity λ is a penalty for taking action "activate". For λ-penalized MDPs, we define the gain and bias functions as in Section 2.2.2, but these quantities now depend on λ. Hence, we will write them as functions of λ: For instance, the optimal gain is g * (λ), and we will use the notation h * (λ) to denote an optimal bias and π * (λ) to denote an optimal policy.

The classical definition of indexability use in the literature [START_REF] Akbarzadeh | Conditions for indexability of restless bandits and an O(K 3 ) algorithm to compute Whittle index[END_REF][START_REF] Gibson | A novel implementation of Qlearning for the Whittle index[END_REF][START_REF] Nakhleh | Neur-WIN: Neural Whittle index network for restless bandits via deep RL[END_REF] says that an arm is indexable if and only if the optimal policy π * (λ) is non-increasing in λ (for the inclusion order). If an arm is indexable, these papers define the Whittle index of a state i as a real number λ i such that π * (λ) = {i ∈ [n] : λ i > λ}. This definition is ambiguous for two reasons: First, optimal policies are in general not unique. Hence, the notion of π * (λ) being non-increasing is unclear: should all optimal policies be non-increasing or at least one? Second, the notion of optimality for a policy is also unclear: should it mean "gain optimal", "bias optimal" or another notion of optimality?

To solve these ambiguities, in this paper, we use the following definition of indexability.

Definition 1. Given a finite-state arm, let Π * (λ) be the set of Bellman optimal policies for a penalty λ. We say that the arm is indexable if for all λ < λ , and all policies π ∈ Π * (λ) and π ∈ Π * (λ ), then π ⊇ π . This definition says that the function π * (λ) ⊇ π * (λ ) regardless of the choice of Bellman optimal policies. As we show next, it guarantees that the Whittle indices are uniquely defined when they exist. As we detail in Appendix A.1, this is not necessarily the case when we consider other interpretations of the classical definition.

Note that for discounted problems, this definition coincide with the one used in [START_REF] Niño-Mora | A fast-pivoting algorithm for Whittle's restless bandit index[END_REF]. For undiscounted MDPs, we add in addition that the criterion for optimality should be the Bellman optimality.

Definition of Whittle index and characterization of indexability

The proposition below shows that Definition 1 implies that Whittle index is well defined and proposes a characterization of any indexable arm, that we will later use to derive our algorithm.

Lemma 1. In a n-state arm, the following three properties are equivalent: (i) The arm is indexable.

(ii) For all state i ∈ [n], there exists a unique penalty λ i -called the Whittle index of state i -such that if π ∈ Π * (λ) is any Bellman optimal policy for the penalty λ, then

π i = 1 if λ < λ i and π i = 0 if λ > λ i . (iii) There is a non decreasing sequence of penalties µ 0 min := -∞ ≤ µ 1 min ≤ µ 2 min ≤ • • • ≤ µ n min ≤ µ n+1 min := +∞ and a sequence of policies π 1 := [n] π 2 • • • π n+1 := ∅ such that: • If λ ∈ (µ k-1
min , µ k min ), there exists a unique Bellman optimal policy π k . • If k is such that µ k-1 min < µ k min , then all Bellman optimal policies for the penalty µ k-1 min contain π k , and π k contains all Bellman optimal policies for the penalty µ k min .

In the above lemma, we use a subscript "min" in the penalties µ k min in order to be consistent with the same quantities used in Algorithm 1 and 2. The signification of this "min" is because it will be a minimum of values of the form µ k i . We should stress that these quantities (as well as the Whittle index λ i ) can either be finite or infinite. When we say that "a policy π is optimal for the penalty +∞", this means "there exists a penalty λ such that π is optimal for all λ ≥ λ". Also, the last part of the lemma implies that π k is the unique Bellman optimal policy for all penalty λ ∈ (µ k-1 min , µ k min ).

Proof The lemma is a direct consequence of the definition of indexability.

(i) ⇒ (ii) -Assume first that the arm is indexable. Let i ∈ [n] be a state and let λ i = sup{λ : ∃π ∈ Π * (λ) such that π i = 0}. By Definition 1, if π is a Bellman optimal policy for a penalty λ > λ i , then π ⊆ π, which in turn implies that π i = 0. Similarly, if λ < λ i , then π i = 1. This implies (ii).

(ii) ⇒ (iii) -Assume (ii) and let σ k be the state with the kth smallest index (where ties are broken arbitrarily). Let µ k min := λ σ k be the index of the state σ k and let λ ∈ (µ k-1 min , µ k min ). By (ii), any Bellman optimal policy for the penalty λ < λ σ k-1 contains π k := [n] \ {σ 1 , . . . , σ k-1 }. Similarly, π k contains any Bellman optimal policy for the penalty λ > λ σ k-1 . This implies that the policy π k is the unique optimal policy for all λ ∈ (µ k-1 min , µ k min ). (iii) ⇒ (i) -The property (iii) implies that implies that π k is the unique Bellman optimal policy for all λ ∈ (µ k-1 min , µ k min ).

Condition for indexability and basic algorithm

This section aims at providing a basic algorithm to detect whether an arm is indexable or not and if it is the case, to compute the Whittle index of all states. This algorithm tries to construct a sequence of unichain policies π 1 π 2 . . . that satisfy the conditions of Lemma 1. We prove the correctness of our algorithm: if it can construct such a sequence, then the problem is indexable and the computed indices are correct. If the algorithm cannot compute such a sequence of policies, this is either because the problem is not indexable, or because the arm is multichain.

Condition for optimality

In this section, we provide two technical lemmas that we will use in our algorithm. They provide conditions to verify when a given unichain policy is Bellman optimal and if yes, when it is the unique Bellman optimal policy. Let π ⊆ [n] be a unichain policy and fix a penalty λ. By [START_REF] Puterman | Markov Decision Processes: Discrete Stochastic Dynamic Programming[END_REF]Chapter 8], there exists a gain g π and a bias vector h π such that π satisfies Bellman evaluation equations: for all i ∈ [n],

g π (λ) + h π i (λ) = r πi i -λπ i + n j=1 P πi ij h π j (λ). ( 4
) Whittle Index
We denote by α π i the active advantage in state i under policy π, which is the difference between the value in state i of action activate and the one of action rest. It is defined by:

α π i (λ) := r 1 i -r 0 i -λ + n j=1 (P 1 ij -P 0 ij )h π j (λ). ( 5
)
For a unichain policy, Equation (4) uniquely determines the vector h π (λ) up to an additive constant c1 (see [START_REF] Puterman | Markov Decision Processes: Discrete Stochastic Dynamic Programming[END_REF]Chapter 8]). Hence the active advantage vector α π (λ) is uniquely determined for a unichain policy π. As we will see later, the function α π (λ) is affine in λ. Note that despite the name "advantage", α π can be negative.

Our algorithm computes the Whittle index in increasing order, by trying to eliminate states one by one. The following lemma shows that to compute the next Bellman optimal policy, one should look at when the active advantage of a state is equal to 0. In this lemma, π {i} denotes the symmetric difference between π and {i}, i.e., π {i} = π \ {i} if i ∈ π and π {i} = π ∪ {i} if i ∈ π. Also, the active advantage provides necessary and sufficient condition for a unichain policy to be Bellman optimal, and/or to be the unique Bellman optimal policy as shown in the following lemma.

Lemma 2. In a finite-state arm, let π be a unichain policy. Then, for any penalty λ:

(i) π is Bellman optimal if and only if α π i (λ) ≥ 0 for all i ∈ π and α π i (λ) ≤ 0 for all i / ∈ π. (ii) Suppose that π is Bellman optimal and α π i (λ) = 0. Then, π {i} is also Bellman optimal. If, in addition, π {i} is unichain, then α π (λ) = α π {i} (λ). (iii) π is the unique Bellman optimal policy if and only if α π i (λ) > 0 for all i ∈ π and α π i (λ) < 0 for all i / ∈ π.

Proof For the first point (i), one direction of the equivalence is direct: If policy π is Bellman optimal, then α π i (λ) ≥ 0 for all i ∈ π and α π i (λ) ≤ 0 for all i / ∈ π. This is because a bias vector h π that is a solution of Bellman evaluation equations (4) satisfies Bellman optimality equations [START_REF] Papadimitriou | The complexity of optimal queueing network control[END_REF].

We now prove the other direction of Point (i): If α π i (λ) ≥ 0 for all i ∈ π and α π i (λ) ≤ 0 for all i / ∈ π, then policy π is Bellman optimal. Since π is unichain, its gain g π is state independent and satisfies the optimality equations [START_REF] Gittins | Bandit processes and dynamic allocation indices[END_REF]. If α π i (λ) ≥ 0 for all i ∈ π and α π i (λ) ≤ 0 for all i / ∈ π, then any bias vector h π that is a solution of (4) also satisfies [START_REF] Papadimitriou | The complexity of optimal queueing network control[END_REF]. In consequence, g π 1 and h π form a solution of the optimality equations (1) and [START_REF] Papadimitriou | The complexity of optimal queueing network control[END_REF]. From [START_REF] Puterman | Markov Decision Processes: Discrete Stochastic Dynamic Programming[END_REF]Chapter 9], g * is uniquely defined by (1) and (2). Thus, g * i = g π for all i ∈ [n] and the first condition of Bellman policy characterization equation ( 3) is satisfied. Finally, the fact that bias vector h π satisfies (2) fulfills the second condition of (3). That concludes the proof.

For the second point (ii), since π is unichain, the optimal gain g * i = g π for all i ∈ [n]. So, π {i} satisfies the first condition of (3). Moreover, α π i (λ) = 0 implies that policy π {i} satisfies evaluation equations (4) for some h π . Since π is Bellman optimal, h π is a solution of (2). So, π {i} satisfies the second condition of (3). We conclude that π {i} is Bellman optimal. Last but not least, if, in addition, π {i} is unichain, then h π is a solution of the evaluation equations (4) for policy π {i}. Consequently, α π (λ) = α π {i} (λ). Points (i) and (ii) also show one direction of the equivalence of (iii): If policy π is the unique Bellman optimal policy, then for all state i, α π i (λ) = 0. The non-trivial property is the other direction of the equivalence. This is a consequence of Lemma 6 that we prove in Appendix B.1.

Note that the main difficulty in proving Lemma 6 is that we do not assume the arm to be unichain: for a unichain arm, the bias of the optimal policy is unique up to a constant vector (see [START_REF] Schweitzer | The functional equations of undiscounted Markov renewal programming[END_REF] and [START_REF] Puterman | Markov Decision Processes: Discrete Stochastic Dynamic Programming[END_REF]Section 8.4]). This implies that if π is an optimal policy and α π i (λ) = 0 for all i, then π is the unique optimal policy. The proof of Lemma 6 that we do in Appendix B.1 does not require the MDP to be unichain, but only the Bellman optimal policy π to be unichain. Lemma 6 shows that α π i (λ) = 0 implies that no other policy can be Bellman optimal (not even multichain policies). In Appendix B.1, we prove that this holds not only for two-action MDPs but also for any MDP with finite state and action spaces.

Fig. 2:

The active advantage as a function of penalty for two unichain examples, one is indexable (Figure 2a) and the other is not (Figure 2b). The red dots mark where the lines change their slope. The parameters of both examples are provided in Appendix A.4.

To illustrate Lemma 1 and Lemma 2, we consider two three-state arms. For each model, we plot in Figure 2 the active advantage α π * (λ) i (λ) as a function of the penalty λ for all states i ∈ {1, 2, 3}. By Lemma 2, we know that an optimal policy should activate all states having a positive advantage and rest all states having a negative advantage. Combined with the characterization of Lemma 1, this shows that:

• The model presented in Figure 2a is indexable: the optimal policy is a non-increasing function of λ and the indices are λ 1 ≈ 0.3, λ 2 ≈ 0.8 and λ 3 ≈ 0.7. • The model presented in Figure 2b is not indexable: the optimal policy π * (0.7) = {3} is not included in π * (0.6) = {1}.

Overview of the algorithm

Our algorithm computes Whittle index in increasing order by navigating through unichain Bellman optimal policies and using the characterization provided by Lemma 1. It follows the graphical construction given in Figure 2a.

It uses the following facts:

• If policy π 1 := [n] is unichain, then α π 1 is decreasing in λ • Similarly, if policy π n+1 := ∅ is unichain, then α π n+1 is decreasing in λ
• For an indexable arm, computing the index can be done by a greedy algorithm that constructs a sequence of penalties

µ 1 min ≤ µ 2 min ≤ • • • ≤ µ n
min and a sequence of unichain policies π 1 • • • π n by looking at where α π k i (λ) intersects horizontal axis for all i ∈ π k . • The arm is indexable if and only if for all k such that µ k-1 min < µ k min , the constructed π k is the largest Bellman optimal policy for the penalty µ k min . In order to compute Whittle index and test indexability, our algorithm needs that all policies π k constructed by the algorithm to be unichain. It does not require the arm to be unichain.

This leads to Algorithm 1, that we write in pseudo-code. This algorithm relies on two subroutines: on Line 7, to compute the next index and on Line 8 to test if a policy is Bellman optimal. We will describe later in the paper how to implement these functions in an efficient manner. Note that in all the paper, we use the superscript k (e.g., π k , µ k , σ k ) to refer to the quantities computed at iteration k. We use the subscripts i or j (e.g., π i , λ i , µ i , π j ) to refer to the quantities related to states i or j.

Note that when µ k min = +∞, the quantity α π k i (µ k min ) defined in Line 7 of Algorithm 1 should be understood as lim λ→∞ α π k i (λ) ∈ R∪{-∞, +∞}. These limits are well defined because the functions αs are affine in λ.

To illustrate how the algorithm works, we plot in Figure 3 the values computed by the algorithm for the two arms represented in Figure 2. For both models (indexable and non-indexable), the algorithm starts with the policy [n] for which the derivative of the active advantage with respect to λ is -1 for all states. It then computes µ 1 min which is the potential index of State σ 1 = 1 for 3a and of State σ 1 = 3 for 3b. The algorithm then moves to iteration 2 and computes µ 2 min > µ 1 min for both models and observes that for both models α π 2 σ 1 (µ 2 min ) < 0 for both models. Then the algorithm moves to iteration 3 and computes µ 3 min > µ 2 min . There are now two cases:

Algorithm 1 Given a n-state arm, test indexability and compute Whittle index (if indexable).

1: Set π 1 := [n], µ 0 min := -∞ 2: if π 1 is multichain then 3: return the arm is multichain 4: end if 5: for k = 1 to n do 6: Compute α π k (λ) 7: Let µ k min := inf{λ ≥ µ k-1 min : ∃i ∈ π k , α π k i (λ) = 0} 8: if µ k-1 min <µ k min and for some i / ∈ π k , α π k i (µ k min ) ≥ 0 then 9:
return the arm is not indexable

10: end if 11: if µ k min = +∞ then 12:
Set

λ i := +∞ for all i ∈ π k 13:
return the arm is indexable and the indices are

{λ i } i∈[n] .
14:

end if 15: Let σ k ∈ π k be such that α π k σ k (µ k min ) = 0 and λ σ k = µ k min 16: Set π k+1 := π k \ {σ k } 17: if π k+1 is multichain then 18:
return the arm is multichain • For 3a, the algorithm verifies that π 4 := ∅ is unichain (α ∅ i is decreasing in λ for all i), terminates, and returns that the arm is indexable.

• For 3b, the algorithm realizes that α π 3 σ 1 (µ 3 min ) > 0 which shows that this model is not indexable. Note that for the indexable example of Figure 3a, the active advantage function is not a decreasing function of λ. Hence, this example is neither PCL-indexable (defined in [START_REF] Niño-Mora | A fast-pivoting algorithm for Whittle's restless bandit index[END_REF]Definition 3]) nor strongly-indexable (defined in [START_REF] Nakhleh | Neur-WIN: Neural Whittle index network for restless bandits via deep RL[END_REF]). However, this does not prevent our algorithm from working.

Correctness of Algorithm 1

The following result shows that Algorithm 1 is correct.

Theorem 1. Given a n-state arm: (i) if Algorithm 1 outputs "the arm is indexable and the indices are {λ i } i∈[n] ", then the arm is indexable and each λ i is the Whittle index of state i; (ii) if Algorithm 1 outputs "non-indexable", then the arm is non-indexable;

(iii) if Algorithm 1 outputs "multichain", then the arm is multichain. Active Advantage

1 min 1 = 1 1 ={1, 2, 3} 2 min 2 = 3 2 ={2, 3} 3 min 3 = 2 3 ={2} 4 = State 1 State 2 State 3
(a) Indexable arm with 3 states. Note that the algorithm does not compute

α π 4 . It checks if policy π 4 =∅ is unichain or not. If it is, then α π 4
i is decreasing in λ for each i. A direct consequence of Theorem 1 is that for unichain arms, Algorithm 1 provides a full characterization of indexability.

Corollary 2. Given a unichain arm with finite states, Algorithm 1 outputs "the arm is indexable" if and only if it is indexable.

Note that Algorithm 1 does not require the arm to be unichain to work. The required condition is that the Bellman optimal policies {π k } k≥1 that Algorithm 1 uses are all unichain. In particular, there exist examples of arms that are multichain and indexable and for which the algorithm returns "indexable". Similarly, there exist examples of arms that are multichain and non-indexable and for which the algorithm returns "non-indexable". We provide such examples in Appendix A.2.

Our algorithm works by exploring solely unichain policies. It does so because (in general) the bias of multichain policy is not unique. The characterization of Bellman optimal policies is much more difficult for multichain models, and the notion of indexability becomes more elusive (see Example A1 in Appendix A.1). When Algorithm 1 returns "multichain", it means that the algorithm is unable to decide whether the arm is indexable or not (but the algorithm knows that the arm is multichain because it has just found a multichain policy).

Proof of Theorem 1 Proof of (i) -We first prove by induction on k that:

If Algorithm 1 completes iteration k ≥ 1, then π k and π k+1 are unichain and (A) π k is the unique Bellman optimal policy for all λ ∈ (µ k-1 min , µ k min );

(B) π k+1 is Bellman optimal for µ k min and α π k+1 (µ k min ) = α π k (µ k min ).

Base case k = 1: As we prove later in [START_REF] Borkar | Whittle indexability in egalitarian processor sharing systems[END_REF],

∂α π 1 i ∂λ = -1. So, for each i ∈ π 1 , α π 1 i (λ) is decreasing in λ. By definition, µ 1
min is the smallest λ such that one of the α π 1 i (λ) = 0. Hence, for all λ < µ 1 min : α π 1 i (λ) > 0. By Lemma 2, this shows that π 1 is the unique Bellman optimal policy for all λ ∈ (µ 0 min , µ 1 min ), so (A) is true. Moreover, since π 1 is Bellman optimal for the penalty µ 1 min and π 2 is unichain, Lemma 2 implies that π 2 is Bellman optimal for the penalty µ 1 min and α π 2 (µ 1 min ) = α π 1 (µ 1 min ). This shows (B).

Penalty

Active Advantage 4. By the induction hypothesis, π k is Bellman optimal for the penalty µ k-1 min . Moreover, by definition of µ k min , together with

k 1 min k 1 k min k k i ( k 1 min ) > 0 for all i k k i ( k 1 min ) 0 for all i k k is
µ k-1 min < µ k min , α π k i (µ k-1 min ) = 0 for all i ∈ π k . Hence: α π k i (µ k-1 min ) > 0 for i ∈ π k and α π k i (µ k-1 min ) ≤ 0 for i ∈ π k .
Finally, thanks to the test on Line 8 of the algorithm, we have:

α π k i (µ k min ) ≥ 0 for i ∈ π k and α π k i (µ k min ) < 0 for i ∈ π k . ( 6
)
In consequence, for each λ

∈ (µ k-1 min , µ k min ), α π k i (λ) > 0 for i ∈ π k and α π k i (λ) < 0 for i / ∈ π k .
Lemma 2 implies that π k is the unique Bellman optimal policy for each λ ∈ (µ k-1 min , µ k min ). This shows (A). Also, [START_REF] Weber | On an index policy for restless bandits[END_REF] implies that π k is Bellman optimal for the penalty µ k min . Combine this with the fact that π k+1 is unichain, Lemma 2 implies that π k+1 is Bellman optimal for µ k min and α π k+1 (µ k min ) = α π k (µ k min ). This shows (B). So, the induction is also true for iteration k.

This shows that the induction property is true for all k ∈ [n]. In particular, when π n+1 := ∅ is unichain, α π n+1 i is decreasing in λ as we prove later in [START_REF] Larrañaga | Asymptotically optimal index policies for an abandonment queue with convex holding cost[END_REF] that

∂α π n+1 i ∂λ = -1. Combine this with the fact that α π n+1 i (µ n min ) = α π n i (µ n min )
≤ 0 for all i, Lemma 2 implies that π n+1 is the unique Bellman optimal policy for λ > µ n min . We simply set µ n+1 min := +∞. To sum up, there are two cases for which the algorithm outputs that the arm is indexable:

1. if the algorithm goes until the end of iteration n, then the sequence of values {µ k min } k∈[n+1] and of policies {π k } k∈[n+1] satisfies the conditions of Lemma 1(iii) and the arm is indexable. 2. if the algorithm stops at iteration k because µ k min = +∞, then one can set µ k+1 min := • • • := µ n+1 min := +∞ and define a sequence of policies π k+1 • • • π n+1 := ∅ by eliminating all states of π k in an arbitrary order. These sequences satisfies the conditions of Lemma 1(iii) and the arm is indexable.

Proof of (ii) -Our algorithm outputs non-indexable if there exists an iteration

k and a state j / ∈ π k , such that α π k j (µ k min ) ≥ 0 when µ k-1 min < µ k min . We know that π k is Bellman optimal for µ k-1
min , otherwise the algorithm would have stopped before. Assume that:

All Bellman optimal policies for any λ ∈ (µ k-1 min , µ k min ] are included in π k . (7) We will see that this assumption leads to a contradiction. We distinguish two possibilities:

1. π k is Bellman optimal for the penalty µ k min -This implies that α π k i (µ k min ) ≤ 0 for all i ∈ π k which, together with α π k j (µ k min ) ≥ 0, implies that α π k j (µ k min ) = 0. By Lemma 2, this would imply that π k ∪ {j} is Bellman optimal for µ k min . This is in contradiction with [START_REF] Glazebrook | Some indexable families of restless bandit problems[END_REF]. 2. π k is not Bellman optimal for µ k min -In this case, we denote by λ the smallest penalty λ ∈ [µ k-1 min , µ k min ] such that there exists π π k that is Bellman optimal for λ (it exists because π k is not Bellman optimal for µ k min and we assumed ( 7)). By definition of λ, π and π k are both Bellman optimal for the penalty λ.

Let i ∈ π k \ π. By Lemma 2, this implies that α π k i ( λ) = 0.
The problem is that by definition, µ k min is the smallest penalty λ for which there exists i ∈ π k such that α π k i (λ) = 0. This implies that λ = µ k min which in turn implies that π k is optimal for µ k min . This leads to a contradiction. This shows that neither case 1 nor 2 are possible. So, (7) cannot be true. In consequence, the negation of ( 7) is true: there exists λ > µ k-1 min and π ⊆ π k such that π is Bellman optimal for λ. This contradicts Definition 1 and therefore implies that the arm is not indexable.

Proof of (iii) -if our algorithm outputs multichain, then the arm is multichain. This is straightforward based on the definition of multichain MDP.

We should note that by Line 15, it is possible to have µ k min = µ k-1 min . This happens when several states have the same value of Whittle index. This is not problematic because we are sure that σ k = σ k-1 by Line 16.

In the proof of Theorem 1(i), we showed that when policy [n] is unichain, the function α π 1 i (λ) is decreasing in λ which implies that it crosses the line 0 at some finite value µ 1 i . This implies that for an indexable arm, if [n] is unichain then the Whittle index are all strictly larger than -∞. A symmetric argument shows that if policy ∅ is unichain, then all Whittle index are strictly smaller than +∞. This implies the following result.

Corollary 3. Given a unichain arm with n states, if the arm is indexable, then the indices of the n states are finite:

λ i ∈ {-∞, +∞} for all i ∈ [n].
This is not necessarily true for multichain arms (see the discussion in Appendix A.3.)

Naive implementation of Algorithm 1 (in O(n 4 ))

For a given penalty λ, we consider a policy π that is Bellman optimal and unichain. Recall that g * (λ) is the maximal gain, and h π (λ) ∈ R n is a solution of (4). We consider h π (λ) such that h π 1 (λ) = 0. Recall from ( 4) that for all i ∈ [n] :

g * (λ) + h π i (λ) = r πi i -λπ i + n j=1 P πi ij h π j (λ). ( 8
)
The above system is a system of n+1 linear equations with n+1 variables (the additional equation begins with h π 1 (λ) = 0). As π is unichain, the maximal gain g * (λ) and h π (λ) are uniquely determined by the system of linear equations ( 8), together with the condition that h π 1 (λ) = 0. Note that in (8) the sum is for j = 1 to n. Since h π 1 (λ) = 0, it can be transformed into a sum from j = 2 to n. Let us define the vector v π (λ) := [g * (λ) h π 2 (λ) . . . h π n (λ)] which is similar to the vector h π (λ) in which we replaced h π 1 (λ) by g * (λ). We can write Equation ( 8) under a matrix form as:

A π v π (λ) = r π -λπ, ( 9
)
where r π is the reward vector under π:

r π :=[r π1 1 . . . r πn n ] , π:=[π 1 . . . π n ] ,
and A π is the following square matrix:

A π :=      1 1 1 1 . . . 1 1      -    
As we show in Lemma 7, the matrix A π is invertible if and only if policy π is unichain. In consequence, v π is an affine function of λ:

v π (λ) = (A π ) -1 (r π -λπ) = (A π ) -1 r π -λ(A π ) -1 π (11)
For a state i, let δ i := r 1 i -r 0 i , ∆ i1 := 0 and ∆ ij := P 1 ij -P 0 ij for j ∈ {2, . . . , n}. By definition of the advantage function of (5), we have:

α π (λ) = δ -λ1 + ∆v π (λ).
For each active state i, we want to find the smallest penalty

µ k i ≥ µ k-1 min such that α π k i (µ k i ) = 0. Suppose that π k-1 and π k are unichain and π k-1 is Bellman optimal for µ k-1 min . By Lemma 2, α π k (µ k-1 min ) = α π k-1 (µ k-1 min ). Let d π k := -(A π k ) -1 π k . By (11), α π k (λ) is a linear function of λ whose derivative is -(1 -∆d π k ). In particular, α π k (λ) = α π k-1 (µ k-1 min ) -(λ -µ k-1 min )(1 -∆d π k ). Thus, α π k i (λ) = 0 if and only if α π k-1 i (µ k-1 min ) = (λ -µ k-1 min )(1 - n j=2 ∆ ij d π k j ). (12)
Recall that

α π k-1 i (µ k-1 min ) is non-negative for active state i ∈ π k . The value µ k i is the smallest λ ≥ µ k-1
min that satisfies Equation [START_REF] Avrachenkov | Whittle index based Q-learning for restless bandits with average reward[END_REF]. There are three cases:

1. if α π k-1 i (µ k-1 min ) = 0, then µ k i := µ k-1 min ; 2. if α π k-1 i (µ k-1 min ) > 0 and (a) if 1 - n j=2 ∆ ij d π k j > 0, then µ k i := µ k-1 min + α π k-1 i (µ k-1 min ) 1 - n j=2 ∆ ij d π k j ; ( 13
) (b) if 1 - n j=2 ∆ ij d π k j ≤ 0, then µ k i := +∞.
This shows that, for a given k, computing µ k min of Line 7 can be done in O(n 3 ): A first part in O(n 3 ) to compute the inverse of matrix A π and to compute d π , plus some smaller order terms to compute the solutions of [START_REF] Avrachenkov | Whittle index based Q-learning for restless bandits with average reward[END_REF]. Similarly, the test in Line 8 of Algorithm 1 can also be implemented in

O(n 3 ) by using α π k (µ k min ) = α π k-1 (µ k-1 min ) -(µ k min -µ k-1 min )(1 -∆d π k
) with the convention that when µ k min = +∞, +∞ × 0 = 0 and +∞ × x = sign(x)∞ for any x = 0. This leads to an overall complexity of O(n 4 ) for Algorithm 1 that contains n loops each having a O(n 3 ) complexity. If at some iteration k the matrix A π k is not invertible, then Lemma 7 implies that π k is multichain. In consequence, the algorithm outputs multichain and stops. We integrate this in the newer version of our algorithm below.

The (2/3)n 3 + o(n 3 ) algorithm

This section describes a way to implement Algorithm 1 efficiently using O(n 3 ) operations. The main idea is to use the Sherman-Morrison formula to compute in O(n 2 ) the active advantage vector α π k (λ) associated to π k from the one associated to π k-1 . This leads to a O(n 3 ) algorithm. Once this main idea is in place, we show how to avoid unnecessary computations to obtain an algorithm that performs (2/3)n 3 + o(n 3 ) arithmetic operations.

Additional notations

In order to obtain a more efficient and compact algorithm, for an iteration k and a state i, we define

y k i := n j=2 ∆ ij d π k j and z k i := α π k i (µ k min )
, where d π k j , ∆ ij and α π k i are as in [START_REF] Niño-Mora | A dynamic page-refresh index policy for web crawlers[END_REF]. Equation ( 13) can be rewritten as

µ k i = µ k-1 min + z k-1 i 1 -y k i . ( 14
)
The above equation can be used to compute µ k i and µ k min easily from y k i and z k-1 i . Indeed, from the previous section, we have α

π k i (µ k min ) = α π k-1 i (µ k-1 min ) - (µ k min -µ k-1 min )(1 - n j=2 ∆ ij d π k j ) which translates into z k i = z k-1 i -(µ k min -µ k-1 min )(1 -y k i ). (15)
This shows that the critical values to compute are the variables y k i . In the remainder of this section, we show that the quantity y k i can be computed efficiently by a recursive formula.

Application of the Sherman-Morrison formula

To compute

y k+1 i := n j=2 ∆ ij d π k+1 j
, we need to compute the quantities d π k+1 := -(A π k+1 ) -1 π k+1 . This requires the inverse of A π k+1 . By definition of π k , two policies π k and π k+1 differ by exactly one state: π k+1 = π ke σ k where e j denotes the column vector with a 1 in jth coordinate and 0's elsewhere. Also by definition of A π in [START_REF] Aalto | Whittle index approach to opportunistic scheduling with partial channel information[END_REF], the two matrices A π k and A π k+1 differ only at the row σ k :

A π k+1 = A π k + e σ k ∆ σ k ( 16
)
where ∆ σ k is a row vector defined as in the previous section. One can efficiently compute the inverse of matrix A π k+1 from the one of A π k by using the Sherman-Morrison formula, which says that if A ∈ R n×n is an invertible square matrix and p, q ∈ R n are two column vectors, then the matrix A + pq is invertible if and only if 1 + q A -1 p = 0 and if A + pq is invertible, then:

(A + pq ) -1 = A -1 - A -1 pq A -1 1 + q A -1 p . Let X k ij := ∆ i (A π k ) -1 e j .
Following [START_REF] Scully | SOAP: One clean analysis of all age-based scheduling policies[END_REF], we can apply the Sherman-Morrison formula with matrix A π k , and vectors p = e σ k and q = ∆ σ k . After some simplification, we get:

X k+1 ij := ∆ i (A π k+1 ) -1 e j = ∆ i (A π k + e σ k ∆ σ k ) -1 e j = X k ij - X k iσ k 1 + X k σ k σ k X k σ k j ( 17
)
In particular,

X k+1 iσ k = X k iσ k 1 + X k σ k σ k .
Before computing X k+1 ij , we need to verify that π k+1 is unichain. With the help of Lemma 7 and the Sherman-Morrison formula, this can be done easily:

π k+1 is unichain if and only if 1 + X k σ k σ k = 0. For y k+1 i := ∆ i d π k+1 = -∆ i (A π k+1 ) -1 π k+1 , we use π k+1 = π k -e σ k
and apply the Sherman-Morrison formula to get:

y k+1 i = -∆ i (A π k + e σ k ∆ σ k ) -1 (π k -e σ k ) = -∆ i (A π k ) -1 (π k -e σ k ) + ∆ i (A π k ) -1 e σ k ∆ σ k (A π k ) -1 1 + ∆ σ k (A π k ) -1 e σ k (π k -e σ k ) = y k i + X k iσ k + X k iσ k (-y k σ k -X k σ k σ k ) 1 + X k σ k σ k = y k i + X k iσ k (1 -y k σ k ) 1 + X k σ k σ k = y k i + (1 -y k σ k )X k+1 iσ k (18)
The above formula indicate how to compute y k+1 from y k . To complete this analysis, let us show that y 1 = 0. For a given policy π, the vector d π satisfies the same equation as Equation (9) but replacing r πi i -λπ i by -π i . This implies that for π = π

1 = [1 . . . 1] , one has d π = [-1 0 . . . 0] as d π 1
is the long-run average reward of a Markov reward process whose reward is negative one in all states and d π 2 , . . . , d π n is the bias of this process. This shows that for all i, one has

y 1 i := n j=2 ∆ ij d π 1 j = 0.
Moreover, by [START_REF] Liu | Indexability of restless bandit problems and optimality of Whittle index for dynamic multichannel access[END_REF], one has

α π 1 (λ) = δ -λ1 + ∆(A π 1 ) -1 r π 1 -λ ∆(A π 1 ) -1 π 1 =:y 1 = δ -λ1 + X 1 r π 1 . ( 19
)
Finally, for π n+1 = [0 . . . 0] , one has

α π n+1 (λ) = δ -λ1 + ∆(A π n+1 ) -1 r π n+1 . ( 20
)

Detailed algorithm

Equation ( 14) shows how to compute µ k i from the values of y k i and z k-1 i while ( 18), ( 17) and [START_REF] Avrachenkov | Impulsive control for G-AIMD dynamics with relaxed and hard constraints[END_REF] show how to compute the values of y, z and X recursively in k. In order to compute µ k min and σ k , one only needs to compute the values

µ k i for i ∈ π k . Once µ k min = min i∈π k µ k i is determined, if µ k-1
min < µ k min , then Line 8 of Algorithm 1 can be performed, based on [START_REF] Avrachenkov | Impulsive control for G-AIMD dynamics with relaxed and hard constraints[END_REF], by checking if

z k i ≥ 0 for some i ∈ [n] \ π k .
This leads to Algorithm 2 that can be decomposed as follows: 1. In Lines 1 to 8, we initialize the various variables. The main complexity of this part is to compute the matrix X 1 , which is equivalent to solving the linear system XA π 1 = ∆. It can be done by inverting the matrix A π 1 and multiplying this by the matrix ∆. This can be done in a subcubic complexity by using for instance Strassen's algorithm [START_REF] Strassen | Gaussian elimination is not optimal[END_REF]. 2. We then enter the main loop:

• We update the vectors µ, z by using Equations ( 14) and ( 15), and test indexability. This costs O(n) operations per iteration, thus O(n 2 ) in total. • We update the vector X according to [START_REF] Aalto | Properties of the Gittins index with application to optimal scheduling[END_REF]. The "naive" way to do so is to use Subroutine 3. At iteration k this costs 2kn arithmetic operations if we test indexability, and 2 k l=1 (n -l) if we do not test indexability. The total complexity of computing X is n 3 + O(n 2 ) arithmetic operations if we test indexability and (2/3)n 3 +O(n 2 ) if we do not. A detailed study of the arithmetic complexity is provided in Appendix C, where we also provide details on how to efficiently implement the algorithm, including how to optimize the cost of memory access.

• In Line 11, we update the vector y by using Equation [START_REF] Aalto | On the Gittins index in the M/G/1 queue[END_REF], which costs O(n) per iteration. 3. Testing if π n+1 is unichain can be done in O(n 2 ) by Tarjan's strongly connected component algorithm. Hence, the total complexity of this algorithm is n 3 + o(n 3 ) if we test indexability and (2/3)n 3 + o(n 3 ) if we do not test indexability. Without testing the indexability, our algorithm has the same main complexity as [START_REF] Niño-Mora | A fast-pivoting algorithm for Whittle's restless bandit index[END_REF]. However, the algorithm of [START_REF] Niño-Mora | A fast-pivoting algorithm for Whittle's restless bandit index[END_REF] computes Whittle index only for an arm that is PCLindexable. Hence we can claim that our algorithm is the first algorithm that computes Whittle index with cubic complexity for general indexable restless bandits. It is also the first algorithm with cubic complexity that tests nonindexability of restless bandits composed of unichain arms. Note that the ties are broken arbitrarily for Lines 7 and 15 of Algorithm 2.

Algorithm 2 Given a n-state arm, test indexability and compute Whittle indices (if indexable).

1: Set π 1 := [n], k 0 = 1, ∆ i1 := 0, ∆ ij := P 1 ij -P 0 ij , ∀i ∈ [n]
, j ∈ {2, . . . , n}, y 1 = 0 0 0 and A π 1 is defined by [START_REF] Aalto | Whittle index approach to opportunistic scheduling with partial channel information[END_REF]. 2: if π 1 is multichain then 3:

return the arm is multichain 4: end if 5: Set X 1 := ∆(A π 1 ) -1 6: Set µ 1 := r 1 -r 0 + X 1 r 1 7: Let σ 1 := arg min i∈π 1 µ 1 i , λ σ 1 = µ 1 min := µ 1 σ 1 , and π 2 := π 1 \ {σ 1 } 8: Set z 1 = µ 1 -µ 1 min 1 9: for k = 2 to n do 10: Update X(k -1)
Here we call Subroutine 3 or Subroutine 4

11:

Set

y k = y k-1 + (1 -y k-1 σ k-1 )X k :σ k-1 12: for i ∈ π k do 13: Set µ k i :=          µ k-1 min , if z k-1 i = 0 µ k-1 min + z k-1 i 1-y k i , if z k-1 i >0 and 1-y k i >0 +∞, otherwise 14:
end for 15:

Let σ k := arg min i∈π k µ k i and µ k min := µ k σ k 16: Set z k = z k-1 -(µ k min -µ k-1 min )(1 -y k ) 17: if µ k-1 min < µ k min and z k i ≥ 0 for some i ∈ [n] \ π k then 18:
return the arm is not indexable Remark: Equivalently, instead of calling Subroutine 3 at Line 10, one could do the following update for all j ∈ π k+1 and i ∈ [n] (or i ∈ π k+1 if we do not test indexability):

X k+1 ij = X k ij - X k iσ k 1 + X k σ k σ k X k σ k j . ( 21
)
Subroutine 3 Update X(k)

1: for = 1 to k -1 do 2:
for i ∈ [n] do or i ∈ π +1 if we do not test indexability.

3:

X +1 iσ k = X iσ k -X +1 iσ X σ σ k 4:
end for 5: end for

6: if 1 + X k σ k σ k = 0 then 7:
return the arm is multichain 8: end if 9: for i ∈ [n] do or i ∈ π k+1 if we do not test indexability.

10:

X k+1 iσ k = X k iσ k 1 + X k σ k σ k 11: end for
This iterative update is very close to the one used in [START_REF] Niño-Mora | A fast-pivoting algorithm for Whittle's restless bandit index[END_REF][START_REF] Akbarzadeh | Conditions for indexability of restless bandits and an O(K 3 ) algorithm to compute Whittle index[END_REF] for discounted restless bandit. This results in an algorithm that has the same total complexity as Subroutine 3 (of n 3 + O(n 2 ) or (2/3)n 3 + O(n 2 ) with or without the indexability test) because both algorithms will have computed the same values of X k ij . The reason to use Subroutine 3 is that, as we will see in the next section, not all values of X k ij are needed at iteration k: in particular, for < k, the computation of the value X iσ k-1 has no interest per say and is only useful because it allows to recursively compute X k iσ k-1 . In the section below, we show how to reduce the cost by avoiding the computation of X iσ k-1 when is much smaller than k. We comment more on the differences with [START_REF] Niño-Mora | A fast-pivoting algorithm for Whittle's restless bandit index[END_REF][START_REF] Akbarzadeh | Conditions for indexability of restless bandits and an O(K 3 ) algorithm to compute Whittle index[END_REF] in Appendix E and in particular we explain why our approach can be tuned into a subcubic algorithm while (21) cannot.

The subcubic algorithm 5.1 Main idea: recomputing X k from A π k periodically

The main computational burden of Algorithm 2 is concentrated on two lines: on Line 5 where we compute X 1 by solving a linear system, and on Line 10 where we compute the column vector

X k :σ k-1 from X 1 :σ k-1 .
The remainder of the code runs in O(n 2 ) operations and is therefore negligible for large matrices. In fact, the computation of X 1 is done by solving a linear system, which can be computed by using a subcubic algorithm (like Strassen [START_REF] Strassen | Gaussian elimination is not optimal[END_REF]). In this section, we show how to optimize our algorithm by reducing the complexity of the update X() function, at the price of recomputing the full matrix X k0 from A π k 0 periodically.

At iteration k of Algorithm 2, the quantities X k iσ k-1 are used at Line 11 to obtain the values y k i . The matrix X k is defined as X k := ∆(A π k ) -1 . It also Whittle Index satisfies Equation [START_REF] Aalto | Properties of the Gittins index with application to optimal scheduling[END_REF], that is, for an iteration and states i and j, we have:

X +1 ij = X ij -X +1 iσ X σ j . ( 22
)
The way Subroutine 3 is implemented is to initialize X 1 := ∆(A π 1 ) -1 and then use [START_REF] Glazebrook | Index policies for the admission control and routing of impatient customers to heterogeneous service stations[END_REF] recursively to compute the column vector

X k :σ k-1 from X 1 :σ k-1
at each iteration.

Here, we propose an alternative formulation which consists in recomputing the whole matrix X k := ∆(A π k ) -1 every K iterations. In the meantime, we use [START_REF] Glazebrook | Index policies for the admission control and routing of impatient customers to heterogeneous service stations[END_REF] to compute the values X iσ k-1 for ∈ {k 0 + 1, . . . , k} where k 0 is the iteration at which we recomputed the whole matrix X k0 := ∆(A π k 0 ) -1 . This can be implemented by replacing the call to Subroutine 3 at Line 10 with a call to Subroutine 4.

Subroutine 4 Update X FMM(k)

1: if k is an iteration at which we recompute the whole X k+1 then 2:

Set k 0 := k + 1 3: if A π k 0 is not invertible then 4:
return the arm is multichain

5: end if 6: Set X k0 := ∆(A π k 0 ) -1 7: else 8: for = k 0 to k -1 do 9: for i ∈ [n] do
or i ∈ π +1 if we do not test indexability.

10:

X +1 iσ k = X iσ k -X +1 iσ X σ σ k 11:
end for 12:

end for 13:

if 1 + X k σ k σ k = 0 then 14:
return the arm is multichain

X k+1 iσ k = X k iσ k 1 + X k σ k σ k 18:

end for 19: end if

To see why this can be more efficient, we illustrate in Figure 5 the pairs ( , σ k-1 ) for which we compute the vector X :σ k-1 , either for Subroutine 3 or Subroutine 4. In each case, a vertical blue line indicates that we recompute the whole matrix X k by solving a linear system. The gray zone corresponds to the values ( , σ k-1 ) for which we compute X :σ k-1 using Equation ( 22) and the red squares represent the vector X k :σ k-1 used at Line 11 of Algorithm 2.

For Subroutine 3, we do one matrix inversion at the beginning and then compute for all ( , σ k-1 ) with ≤ k because the red square at value (k, σ k-1 ) is computed starting from the vertical blue line at value (1, σ k-1 ). For Subroutine 4, we do n/K (here n/K = 4) full recomputation of X k , which correspond to the vertical blue lines. We gain in terms of operations because the surface of the gray zone to compute is divided by n/K .

σ k σ k K K σ K σ 2K
. . . with Subroutine 4 when we do not test the indexability. For Subroutine 4, we solve more linear systems (each vertical blue line corresponds to solving a linear system) but we reduce the gray zone to compute. A linear arrow corresponds to the internal loop of Line 8 of Subroutine 4.

Note that the y-axis of Figure 5 is ordered by increasing value of σ k (and not by increasing value of k). The value of σ k is computed at iteration k but unknown before iteration k. This explains why in Subroutine 4, when we recompute the matrix (X k ij ) at an iteration k (vertical blue lines in Figure 5(b)), we recompute it for all i, j and not just i, j ∈ {σ k , . . . , σ K } (which are the only values that we will use): Indeed, σ k , . . . , σ K are unknown at iteration k.

A subcubic algorithm for Whittle index

We now assume to have access to a subcubic matrix multiplication algorithm that satisfies the following property: (FMM) There exists an algorithm to multiply a matrix of size n × n by a matrix of size 4 n × n γ that runs in O(n ω(γ) ), where ω : [0, 1] → [2, 3] is a nondecreasing function. Going back to Algorithm 2 where Line 10 is Subroutine 4, we now assume that we recompute the whole matrix X k every O(n γ ) iterations. The new algorithm has a subcubic complexity: Theorem 4. Given a n-state arm, Algorithm 2 with Subroutine 4 checks indexability and computes Whittle (and Gittins) index in time at least Ω(n 2.5 ) and at most O(n 2.5286 ) when choosing γ = 0.5286. We believe that Theorem 4 is the first theoretical result that shows that Whittle index can be computed in subcubic time. As we show in Section 7, this algorithm can be directly extended to discounted index. As a byproduct, we also obtain the first subcubic algorithm to compute Gittins index.

Proof The algorithm starts by computing X 1 which can be done in O(n ω (1) ). Then, there are n 1-γ times that we do: 1. We fill the "gray" mini matrices by using [START_REF] Glazebrook | Index policies for the admission control and routing of impatient customers to heterogeneous service stations[END_REF]. This amounts to three for loops of size n (for i), n γ (for k) and n γ (for ). Hence, each small gray matrix costs O(n 1+2γ ).

At the end of a cycle, we recompute the full inverse by updating (A

π k ) -1 from (A π k-n γ ) -1 .
As we show in Lemma 3 (stated below), this can be done in O(n ω(γ) ).

This implies that the algorithm has a complexity:

O(n ω(γ) ) + n 1-γ O(n 1+2γ ) + O(n ω(γ) ) = O(n max{2+γ,1-γ+ω(γ)} ).
To compute the optimal value of γ minimizing this expression requires the knowledge of the function ω(γ) which is not known. The current state of the art only gives a lower bound (ω(γ) ≥ 2 ) and an upper bound described in [START_REF] Gall | Improved rectangular matrix multiplication using powers of the Coppersmith-Winograd tensor[END_REF].

It is shown in [START_REF] Gall | Improved rectangular matrix multiplication using powers of the Coppersmith-Winograd tensor[END_REF]Table 3] that γ = 0.5286 is the smallest currently known value of γ for which ω(γ) < 1 + 2γ. This implies that the complexity is at most O(n 2.5286 ).

As for the lower bound, ω(γ) ≥ 2 implies that the complexity of the algorithm is at least Ω(n 2.5 ).

In the next lemma, B plays the role of A π k and A the role of A π k-n γ . Note that as required in the lemma, exactly n γ rows and columns are changed between the two.

Lemma 3. Assume (FMM).

Let A be a square matrix whose inverse A -1 has already been computed, and let B be an invertible square matrix such that A -B is of rank smaller than n γ . Then, it is possible to compute the inverse of B in O(n ω(γ) ).

Proof The matrix B can be written as

B = A + U CV where U is a n × n γ matrix, C is n γ × n γ and V is n γ × n.
The Sherman-Morrison-Woodbury formula [START_REF] Woodbury | Inverting modified matrices[END_REF] states that

B -1 = (A + U CV ) -1 = A -1 -A -1 U C -1 + V A -1 U -1 V A -1 .
This shows that B -1 can be computed by:

• Computing D := A -1 U and E := V A -1 : this takes O(n ω(γ) ). • Computing F := C -1 + V A -1 U
-1 : as this is the inversion of a n γ × n γ matrix, it can be done in O(n γω (1) ) where γω(1) ≤ ω(γ).

• Computing G := DF and then GE: this again takes O(n ω(γ) ).

Hence, computing B -1 can be done in O(n ω(γ) ) operations for the inversion and all multiplications plus an additional O(n 2 ) term for the subtraction and the addition.

As ω(γ) ≥ 2, this concludes the proof of the lemma.

The subcubic algorithm in practice

The complexity of O(n 2.5286 ) given in Theorem 4 is mainly of theoretical interest. The value γ = 0.5286 is obtained by using the best upper bound on ω(γ) known today which is based on the Coppersmith-Winograd algorithm and its variants. The Coppersmith-Winograd algorithm (or its variants) are, however, known as a galactic algorithm: the hidden constant in the O() is so large that their runtime is prohibitive for any reasonable value of n. Hence, the existence of these algorithms is of theoretical interest but has limited applicability. This does not discard the practical improvement provided by Subroutine 4 which is based on the mere fact that multiplying two matrices (or inverting a matrix) is faster than three nested loops even for matrices of moderate size. To verify this, we launched a detailed profiling of the code of Algorithm 2 with the non-optimized Subroutine 3. It shows that for a problem of dimensions 5000, the update of Line 10 takes more than 90% of the computation time, the initialization of X 1 on Line 5 takes about 5% of the time and the rest of the code takes less than 1% of the running time. Now, if inverting the full matrix takes about 5% of the execution time, and updating the gray zone takes 95%, then by doing 5 updates, one can hope to obtain an algorithm whose running time is roughly 5 × 5 + 95/5 ≈ 43% the one of the original implementation. As we observe in Section 6, this is close to the gain that we obtain in practice. A general way to choose the best number of updates is used in the numerical section. It is based on the following reasoning. For large matrices (say n ≥ 10 3 ), the fastest implementations of matrix multiplication and inversion are based on Strassen's algorithm [START_REF] Huang | Strassen's algorithm reloaded[END_REF][START_REF] Huang | Practical fast matrix multiplication algorithms[END_REF]. As we report in Appendix D, the time to solve a linear system of size n by using the default installation of scipy seems to run in O(n 2.8 ). By replacing the function ω(γ) used in Theorem 4 by a more practical bound (ω(γ) = 2.8), the best value for γ becomes γ = 0.9. This indicates that our algorithm can be implemented in O(n 2.9 ) by doing O(n 0.1 ) recomputation of X k from A π k . Note that even for very large values of n (like n = 15000), n 0.1 remains quite small, e.g., 15000 0.1 ≈ 2.6. In practice, we observe that updating int(2n 0.1 ) times (the notation int(x) indicates that it is rounded to the closest integer) gives the best performance among all algorithms, as reported in the next section.

Numerical experiments

In complement to our theoretical analysis, we developed a python package that implements Algorithm 2 and gives the choice of using the variant of Subroutine 3 or of Subroutine 4 to do the "update X()" function. This package relies on three python libraries: scipy and numpy for matrix operations, and numba to compile the python code. To facilitate its usage, this package can be installed by using pip install markovianbandit-pkg.

All experiments were conducted on a laptop (Macbook Pro 2020) with an Intel Core i9 CPU at 2.3 GHz with 16GB of Memory using Python 3.6.9 :: Anaconda custom (64-bit) under macOS Big Sur version 11.6.2. The version of the packages are scipy version 1.5.4, numpy version 1.19.5 and numba version 0.53.1. The code of all experiments is available at https://gitlab.inria.fr/ markovianbandit/efficient-whittle-index-computation.

In all of our experiments, the way we generated arms guarantees that they are almost surely unichain because all the elements on their diagonal as well as on the upper and lower diagonals are positive.

Time to compute Whittle indices

To test the implementation of our algorithm, we randomly generate restless bandit arms with n states where n ∈ {100, 1000, . . . , 15000}. In each case, both transition matrices are uniform probabilistic matrices: for each row of each matrix, we generate n i.i.d. entries following the exponential distribution and divide the row by its sum. This means that all matrices are dense. We use dense matrix since it is the worst case for computational interest. By running our algorithms on sparse matrix, we would expect to have faster running time. Note that all tested matrices are indexable. This is coherent with [START_REF] Niño-Mora | Dynamic priority allocation via restless bandit marginal productivity indices[END_REF][START_REF] Gast | Exponential convergence rate for the asymptotic optimality of Whittle index policy[END_REF] that report that for uniform matrices, the probability of finding a non-indexable example decreases very rapidly with the dimension n. Finally, reward vectors were generated from random Uniform[0,1) entries.

We record the runtime of the different variants of our algorithm and report the results in Table 1. Note that these results present the whole execution time of the algorithm, including the initialization phase in which X 1 is computed. For each value of n, we run Algorithm 2 with four variants:

• The first two columns correspond to the O(n 3 ) algorithm (that uses Subroutine 3 for Line 10), either (a) with the indexability test, or (b) without the indexability test. • The last two columns correspond to the subcubic algorithm (that uses Subroutine 4 for Line 10 with int(2n 0.1 ) updates), either (c) with the indexability test, or (d) without the indexability test. Our numbers show that our algorithm can compute the Whittle index in less than one second for n = 1000 states and slightly less than 7 minutes for n = 15000 states with variant (d). As expected, not doing the indexability test does improve the performance compared to doing the indexability test (here by a factor approximately 1/3 for Subroutine 3 and 1/5 for Subroutine 4). More importantly, this table shows that the time when using the subcubic variant, Subroutine 4, diminishes the computation time by about 40% to 50% compared to when using Subroutine 3. Note that for n = 100, using the Subroutine 3 is slightly faster than using the Subroutine 4 (while both takes only a few milliseconds). This indicates that the subcubic algorithm becomes interesting when n is large enough (say n ≥ 2000).

To give a visual idea of how the various variants of the algorithms compare, we plot in Figure 6a the runtime of the four variants along with two variants of the algorithms of [START_REF] Niño-Mora | A fast-pivoting algorithm for Whittle's restless bandit index[END_REF]: FPA-Matlab (the original matlab implementation), and FPA-Julia: a Julia's implementation of the algorithm provided by the authors. We choose to compare to this algorithm as it was the one with the smallest complexity up to now. The numbers for FPA-Matlab are the ones reported in [START_REF] Niño-Mora | A fast-pivoting algorithm for Whittle's restless bandit index[END_REF] and they are comparable to the ones that we obtained on our machine with the same algorithm. FPA-Julia is significantly faster. Hence, we plot in Figure 6b the runtime of each variant divided by the runtime of FPA-Julia. For large n, our best implementation is about 4 to 6 times faster than the best one of [START_REF] Niño-Mora | A fast-pivoting algorithm for Whittle's restless bandit index[END_REF]. For instance, for n = 15000, our implementation takes about 7 minutes to compute the index (or 9 minutes when checking indexability) whereas FPA-Julia takes 33 minutes (and does not check indexability on the fly). In our implementation, not testing indexability reduces the computation time of 15 -20% for Subroutine 4 or 25 -30% for Subroutine 3 compared to the version that tests indexability.

It should be clear that the comparison of the computation times of our implementation versus the ones of FPA-Matlab or FPA-Julia has its limits, because we do not use the same programming language. The influence of the choice of programming language is clear when comparing FPA-Matlab and FPA-Julia: while both codes are similar, the compiled Julia code is about 5 times faster. To obtain a fairer comparison with our algorithm, we tried to Whittle Index rewrite the algorithm of [START_REF] Niño-Mora | A fast-pivoting algorithm for Whittle's restless bandit index[END_REF] in Python with Numba but our implementation was significantly slower than the one of FPA-Julia. We do not know if this is by lack of optimization of our code or because Julia is indeed faster. Regardless of the choice programming langages, our implementation has nevertheless a few advantages that can explain why it is faster: our algorithm has some technical advantage (a simpler internal loop, and the use of Subroutine 4), and our implementation is optimized to improve data access pattern (which consists mostly in sorting on the fly the array according to the permutation σ, see Appendix C). Fig. 6: Numerical result over 7 simulations: in each simulation, we run the algorithm over randomly generated RBs with the state size ranging over {1000, . . . , 15000}. We plot the average runtime over 7 simulations. The solid lines represent the result of Subroutine 3 and the dashed-dot lines represent the one of Subroutine 4. The marker "+" indicates that algorithms test indexability and the triangles indicates that algorithms do not test indexability.

Statistics of indexable problems

To the best of our knowledge, our algorithm provides the first indexability test that scales well with the dimension n. We used this to answer a very natural question: given a randomly generated arm, how likely is it to be indexable? This question was partially answered in [START_REF] Niño-Mora | Dynamic priority allocation via restless bandit marginal productivity indices[END_REF] that shows that when generating dense arms, the probability of generating a non-indexable arm is close to 10 -n for n ∈ {3, . . . , 7}. This suggests that most arms are indexable. Below, we answer two questions: what happens for larger values of n, and more importantly, what happens when the state transition matrices are not dense?

To answer these questions, we consider randomly generated arms where the matrices P 0 and P 1 are b-diagonal matrices with b non-null diagonals. In particular, b = 3 corresponds to tridiagonal matrices, b = 5 corresponds to pentadiagonal matrices and b = 7 corresponds to septadiagonal matrices. We also compare with the classical case of dense matrices (which corresponds to b = 2n -1). For each model, the entries are generated from the exponential distribution for each row and we divide the row by the sum of generated entries for row. We vary n from 3 to 50 and for each case, we generate 100000 arms. We report in Table 2 the number of indexable arms for each case. Note that pentadiagonal matrices are dense matrices for n = 3 and septadiagonal matrices are dense matrices for n = 4 and do not make sense for n = 3, which is why no numbers are reported. Based on these results, we can assert that dense models are essentially always indexable which conforms with the data reported in [START_REF] Niño-Mora | Dynamic priority allocation via restless bandit marginal productivity indices[END_REF]. The situation is, however, radically different for sparse models: the number of indexable problems decreases quickly with the number of states. For instance, there are only 1 823 indexable 50-state problems among 100 000 generated tridiagonal models (i.e. around 1.8% are indexable). Note that a tridiagonal model is a birth-death Markov chain which is frequently used for queueing systems. Hence, it is very important to check the indexability of the problem because it is not a prevalent property for sparse models. This also calls for new efficient policies in restless multi-arm bandit problems that are not based on Whittle indices.

Extension to the discounted case

The model described in Section 2 corresponds to the definition of Whittle index for a time-average criterion, for which Whittle index is known to be asymptotically optimal [START_REF] Weber | On an index policy for restless bandits[END_REF] for restless bandits. Yet, Whittle index can also be defined for the discounted case [START_REF] Niño-Mora | A fast-pivoting algorithm for Whittle's restless bandit index[END_REF][START_REF] Akbarzadeh | Conditions for indexability of restless bandits and an O(K 3 ) algorithm to compute Whittle index[END_REF]. Notably, the discounted Whittle index simplifies into Gittins index when the bandit is rested (i.e., when P 0 = I and r 0 = 0). In this section, we show how to adapt our algorithm to the discounted case. As a by product, we obtain the first subcubic algorithm to compute Gittins index rested bandit.

Discounted Whittle index

We now consider a λ-penalized MDP in which the instantaneous reward received at time t ≥ 0 is discounted by a factor β t , where β ∈ (0, 1) is called the discount factor: when executing action a in state i at time t ≥ 0, the decision maker earns a reward β t (r a i -λa). For a given policy π, we denote by u π i (λ) the expected sum of discounted rewards earned by the decision maker when the MDP starts in state i at time 0. The vector u π (λ) = [u π 1 (λ) . . . u π n (λ)] is called the value function of the policy π. From [START_REF] Puterman | Markov Decision Processes: Discrete Stochastic Dynamic Programming[END_REF], it satisfies Bellman's equation, that is, for all state i we have:

u π i (λ) = r πi i -λπ i + β n j=1 P πi ij u π j (λ). ( 23
)
The above equation is a linear equation, whose solution is unique because β < 1. It is given by:

u π (λ) = (I -βP π ) -1 (r π -λπ). ( 24
)
For a given penalty λ and a state i, we denote by u * i (λ) := max π u π i (λ) be the optimal value of state i. A policy π is optimal for the penalty λ, i.e., π

∈ Π * (λ), if for all state i ∈ [n], u * i (λ) = u π i (λ)
. By [START_REF] Puterman | Markov Decision Processes: Discrete Stochastic Dynamic Programming[END_REF] such a policy exists, |Π * (λ)| > 0. As mentioned in Section 2.2.3, the distinction between Bellman optimal and gain optimal disappears in discounted MDP in which we are concerned with maximizing the value function. Similarly to the time-average criterion studied before, a β-discounted RB is called indexable if for all penalty λ < λ , all π ∈ Π * (λ) and π ∈ Π * (λ ), one has π ⊇ π .

Analogy between the time-average and the discounted versions

Let π be a policy and i ∈ π be an active state. Similarly to average reward model studied before, the advantage of action activate over action rest in state i right before following policy π is given by,

α π i (λ) := r 1 i -r 0 i -λ+β n j=1 (P 1 ij - P 1 ij )u π j (λ). Then α π i (λ) = 0 if and only if λ = δ i + n j=1 ∆ij u π j (λ).
where δ i is defined as for the average reward model and ∆ is such that for all5 states i, j ∈ [n]: ∆ij := β(P 1 ij -P 0 ij ). To finish the derivation of the algorithm, one should note that the value function u(λ) plays the same role as the vector v(λ) defined for the average reward model. In particular, the definition of u in Equation ( 24) is the analogue of the definition of v in [START_REF] Liu | Indexability of restless bandit problems and optimality of Whittle index for dynamic multichannel access[END_REF] up to the replacement of the matrix A π in (11) by the matrix I -βP π . This means that similarly to v, the value function u(λ) is affine in λ.

Hence, following the same development in Section 4, we can modify Algorithm 2 to compute the discounted Whittle index by modifying only the initialization phase: ∆ := β(P 1 -P 0 ) and X 1 := ∆(I -βP π 1 ) -1 .

Note that we still have y 1 = 0 0 0 because (I -βP π 1 ) -1 π 1 = 1 1-β 1 1 1 (value function in a β-discounted Markov reward process with reward equals to 1 in all states) and ∆1 1 1 = 0 0 0. Also, if Subroutine 4 is used, Line 6 should be changed to X k0 := ∆(I -βP π k 0 ) -1 . Last but not least, in the discounted case, we no longer need to check if the optimal policies are unichain because the matrix (I -βP π ) is invertible for any policy π as long as β < 1 (from Perron-Frobenius' Theorem).

Gittins index

The notion of "restless" bandit comes from the fact that even when the action "rest" is taken, the Markov chain can still change state and generate rewards. When this is not the case (i.e., when P 0 = I and r 0 = 0), an arm is no longer restless and is simply called a Markovian bandit (or a rested Markovian bandit if one wants to emphasize that it is not restless).

In a discounted rested bandit, the notion of Whittle index coincides with the notion of Gittins index (In fact, Whittle index was first introduced as a generalization of Gittins index to restless bandit in [START_REF] Whittle | Restless bandits: Activity allocation in a changing world[END_REF]). In such a case, there is no notion of indexability: a discounted rested bandit is always indexable. Its index can be computed by Algorithm 2 without testing indexability. The best known algorithms to compute Gittins index runs in (2/3)n 3 +O(n 2 ) [START_REF] Chakravorty | Multi-armed bandits, Gittins index, and its calculation. Methods and applications of statistics in clinical trials: Planning, analysis, and inferential methods[END_REF]. When using fast multiplication, our algorithm computes Gittins index in O(n 2.5286 ) which makes it the first algorithm to compute Gittins index in subcubic time.

Note that when P 1 = I, it is possible to compute X 1 = ∆/(1-β) without having to solve a linear system which means that, when P 1 = I, our algorithm with the variant Subroutine 3 has complexity (2/3)n 3 +O(n 2 ). This shows that our cubic algorithm can compute the Gittins index also in (2/3)n 3 + O(n 2 ) instead of (2/3)n 3 + o(n 3 ).

Conclusion

In this paper, we propose a univocal definition of indexability and present an algorithm that is efficient for detecting the non-indexability and computing the Whittle index of all indexable finite-state restless bandits whose arms are all unichain. With no assumptions on the structure of arms, this algorithm can still test the indexability and compute the Whittle index of some arms that are multichain and remains efficient if it is able to do so. Our algorithm is based on the efficient application of the Sherman-Morrison formula. This is a unified algorithm that works for both discounted and non-discounted restless bandits, and can be used for Gittins index computation. We present a first version of our algorithm that runs in n 3 + o(n 3 ) arithmetic operations (or in (2/3)n 3 +o(n 3 ) if we do not test the indexability). So, we conclude that Whittle index is not harder to compute than Gittins index. The second version of our algorithm uses the fastest matrix multiplication method and has a complexity of O(n 2.5286 ). This makes it the first subcubic algorithm to compute Whittle index or Gittins index. We provide numerical simulations that show that our algorithm is very efficient in practice: it can test indexability and compute the index of a n-state restless bandit arm in less than one second for n = 1000, and in a few minutes for n = 15000. These numbers are provided for dense matrices. One might expect to have more efficient algorithms if the arm has a sparse structure. We leave this question for future work.

(GO) "optimal" means gain optimal. (BO) "optimal" means Bellman optimal. The most problematic choice is the notion of increasingness: Interpretation (∃) is more permissive: For instance, consider an arm with two states and assume that the optimal policy is {1, 2} for λ < 0 and is either {1} or ∅ for λ > 0. Interpretation (∃) says that the arm is indexable while interpretation (∀) says that this arm is not indexable. If the arm is indexable, what should the index of state 1 be? Any choice λ 1 ∈ [0, +∞] seems reasonable. Saying that the arm is not indexable clarifies the situation. This is why we choose interpretation (∀) in our paper.

In our paper, we choose the combination (∀-BO) because we believe that, for a problem that has transient state, the notion of Bellman optimality is more meaningful than the notion of gain optimality. Also, our combination (∀-BO) allows for more problems to be indexable compared to (∀-GO) and is easier to characterize. We illustrate these different definitions in Figure A1. For example (a): • The gain optimal policies are {1, 2} and {2} for λ < 0, and {1} and ∅ for λ > 0: According to the interpretation (∃), the problem should be indexable but the index for state 1 is unclear. According to the interpretation (∀), the problem should not be indexable. • The Bellman optimal policy is {1, 2} for λ < 0, and ∅ for λ > 0. According to our definition, (∀-BO), the problem is indexable and the indices are λ 1 = λ 2 = 0. For example (b), the Bellman optimal and gain optimal policies are identical and equal to the gain optimal policies of example (a). Hence, example (b) is not indexable according to our definition. The output of our algorithm for this problem is "multichain".

Note that if the distinction between (BO) and (GO) disappears for discounted problems, the distinction between (∀) and (∃) remains.

A.2 Possible outputs for multichain arms

When running our algorithm on a multichain arm, it outputs "multichain" if one of the policies π k is multichain. This suggests that our algorithm will not necessarily output "multichain" for all multichain arms because it explores only a small subset of the policies. In Figure A2, we provide two examples that illustrate this case. The two are multichain arms for which our algorithm is able to identify the indexability. A2: Two examples of multichain arms for which our algorithm does not return "multichain" but returns "indexable" (a) or "non-indexable" (b).

In the first example shown in Figure A2(a), the arm is multichain because the policy {3} has two recurrent classes: {1, 2} and {3}. Yet, this arm is indexable and the indices are {11, 8, -10}. Our algorithm will output that this arm is indexable because it will explore the sequence of policies π 1 , π 2 , π 3 , π 4 , where

• π 1 = {1, 2, 3} is the unique Bellman optimal policy for λ < -10;

• π 2 = {1, 2} is the unique Bellman optimal policy for λ ∈ (-10, 8); • π 3 = {1} is the unique Bellman optimal policy for λ ∈ (8, 11);

• π 4 = ∅ is the unique Bellman optimal policy for λ > 11. All these policies are unichain, and the policy {3} will never be explored. Hence, our algorithm will output "indexable" for this case and will compute the indices.

In the second example, shown in Figure A2(b), we construct a nonindexable arm by taking the non-indexable 3-state example shown in Figure 3b (parameters are given in (A1)) to which we add an extra state "4". For this state, the active action has a very high reward (1000) and leads to the nonindexable recurrent class. The passive action has a very low reward and stays in state 4. Any policy that does not activate 4 is multichain. The algorithm will start by exploring policies that activate the state 4. As for the original example presented in Figure 3b, our algorithm will realize that the arm is non-indexable when exploring values around λ ≈ 0.70. The algorithm will stop and answer "not indexable" before trying the passive action for state 4 because the active avantage for state 4 is larger than 2000 -λ. The output of the algorithm is thus "non-indexable".

A.3 Multichain arms and infinite indices

Consider the two examples of Figure A3. The examples are multichain because policy ∅ has two irreducible classes for example (a) and policy {1, 2} has two irreducible classes for example (b). These two problems are indexable:

• For (a), the Bellman optimal policy for λ < 0 is {1, 2} and {1} for λ > 0.

The indices are λ 2 = 0 and λ 1 = +∞. • For (b), the Bellman optimal policy is {2} for λ < 0 and {1, 2} for λ > 0.

The indices are λ 1 = 0 and λ 2 = -∞. For the first example, our algorithm returns the correct indices because the constructed policies are π 1 := {1, 2} π 2 := {1} and they are both unichain. For the second example, our algorithm will start with the policy {1, 2} and will stop by saying that this example is multichain.

A.4 Parameters for the example of Figure 2

The numerical data of the indexable arm presented in Figure 2a

Definition and notation

We are given a MDP S, A, r, P with finite state and action spaces. As shown in [START_REF] Puterman | Markov Decision Processes: Discrete Stochastic Dynamic Programming[END_REF]Chapter 9], for such a MDP, the optimal gain g * is a vector that satisfies the multichain optimality equations: for each i ∈ S:

max a∈A j∈S P a ij g * j -g * i = 0 (B2) max a∈A r a i -g * i + j∈S P a ij h j -h i = 0. (B3)
This system uniquely determines the optimal gain that we denote by g * . However, vector h is not uniquely determined by the system. In the following, we denote by H the set of bias vector h such that (g * , h) is a solution of the optimality equations (B2)-(B3). A policy π is Bellman optimal if there exists a bias h ∈ H such that policy π attains the maximum (B3), i.e., for all i ∈ S:

π i ∈ arg max a∈A r a i -g i + j∈S P a ij h j -h i = arg max a∈A r a i + j∈S P a ij h j .
For a given policy π : S → A, we denote by r π and P π the reward vector and state transition matrix under policy π: r π i = r πi i and P π ij = P πi ij . Let g π , h π ∈ R |S| be a solution of the following system:

g π -P π g π = 0 (B4) r π -g π + P π h π -h π = 0. (B5)
The vector g π is uniquely determined by this system of equations and is called the long-run average reward or gain of policy π. The vector h π is unique up to an element of the null space of (I -P π ). Such a vector h π is called a bias of policy π. If (g π , h π ) is a solution of (B4) and (B5), then the advantage of action a over the action π i when the MDP is in state i is given by:

B a i (h π ) := r a i + j∈S P a ij h π j -g π i -h π i = r a i -r πi i + j∈S (P a ij -P πi ij )h π j
We recall the two notions of optimality:

• A policy π is gain optimal if its gain g π equals g * .

• A policy π is Bellman optimal if g * = P π g * and if there exists a bias h π that is a solution of (B5) and also (B3). Recall that a Bellman optimal policy is also gain optimal.

Note that by the definition the advantage B(•), a policy is Bellman optimal if g * = P π g * and if there exists h ∈ H such that B πi i (h) = 0.

Useful notations

Let π : S → A be a policy. We say that a state is recurrent for π if it is recurrent for the Markov chain whose transition matrix is P π . In other words, a state i ∈ S is recurrent if when the chain starts in i at time 0, it almost surely visits the state i at some time t ≥ 1. We denote by R π the set of recurrent states of policy π.

We also define P π as the Cesaro limit of the sequence {(P π ) t } t=1 :

P π := lim T →∞ 1 T T t=1 (P π ) t-1 .
From [40, Section A.4 of Appendix A], the matrix P π exists and has the following properties:

• P π is a stochastic matrix, and satisfies P π P π = P π P π = P π . • For all state i, j, if j ∈ R π , then P π ij = 0. • If π is unichain, then the rows of P π are identical.

B.1.2 Characterization of gain optimal policies

The following lemma characterizes gain optimal policy by showing that the policy must satisfies (B3) on their recurrent states.

Lemma 4. Let π : S → A be a policy and recall that R π is the set of recurrent states of policy π. The three properties below are equivalent. (i) P π g * = g * and for all h ∈ H, B πi i (h) = 0 for all i ∈ R π (ii) P π g * = g * and for some h ∈ H, B πi i (h) = 0 for all i ∈ R π (iii) π is gain optimal. = 0.

Proof (i) ⇒ (ii) is trivial. (ii) ⇒ (iii): By definition of B a i (h), we have r π i -g * i = h i -j∈S P π ij h j for
By Theorem 8.2.6 of [START_REF] Puterman | Markov Decision Processes: Discrete Stochastic Dynamic Programming[END_REF], the average reward of π is P π r π . The above equation shows that P π r π = P π g * . Moreover, the assumption P π g * = g * implies that P π g * = g * which in turn implies that P π r π = g * . This shows that the average reward of π is g * and therefore π is gain optimal.

(iii) ⇒ (i): If π is gain optimal, then P π g * = g * and P π (r πg * ) = 0. The latter rewrites as i∈S P π ki (r π i -g * i ) = 0 for all state k. Let h ∈ H be an optimal bias. For all state k, we have

i∈S P π ki B πi i (h) = i∈S P π ki (r π i -g * i + j∈S P π ij h j -h i ) = 0. (B6)
As h satisfied (B3), for all action a, we have B a i (h) ≤ 0 for all states i ∈ S and in particular B πi i (h) ≤ 0. This shows that for any state i such that P π ki > 0, one must have B πi i (h) = 0. Such state i are the recurrent states of π. This shows that B πi i (h) = 0 for all i ∈ R π .

B.1.3 Characterization of Bellman optimal policies

The previous lemma shows that a policy is gain optimal if and only if the actions for the recurrent states of the policy satisfy (B3). The following lemma shows the relationship between two policies that are unichain and satisfy (B3) on all states.

Lemma 5. Suppose that two policies π and θ are Bellman optimal, unichain and have at least one common recurrent state: R π ∩ R θ = ∅.

Then for any h π and h θ solutions of (B5) for π and θ respectively, there exists a constant c such that for all state i: h

π i -h θ i = c. Moreover, in this case, B θi i (h π ) = B πi i (h θ ) = 0 for all i.
Proof Since π and θ are Bellman optimal, h π , h θ ∈ H. In consequence, we have

h π ≥ r θ -g * + P θ h π .
By Lemma 4 (i), the above inequality is an equality for all i ∈ R θ because θ is gain optimal.

As h θ satisfies (B5), we have

h θ -h π ≤ r θ -g * + P θ h θ -(r θ -g * + P θ h π ) = P θ (h θ -h π ),
with equality for all state i ∈ R θ . This shows that for all t, h θ -h

π ≤ (P θ ) t (h θ -h π ) which implies that h θ -h π ≤ P θ (h θ -h π ) with equality for all states i ∈ R θ . Similarly, h π -h θ ≤ P π (h π -h θ ) with equality for any state i ∈ R π . Let c π i = j∈S P π ij (h π j -h θ j ) and c θ i = j∈S P θ ij (h π j -h θ j )
. By what we have just shown, for all state i, we have

c θ i ≤ equality if i ∈ R θ h π i -h θ i ≤ equality if i ∈ R π c π
As both policies are unichain, c π i and c θ i do not depend on i. Moreover, if there exists i ∈ R θ ∩ R π , we have c π i = c θ i =: c. In consequence, h π i -h θ i = c for all state i.

B.1.4 Unicity of Bellman optimal policy

Lemma 6. Let π be a Bellman optimal policy that is unichain. If π is not the unique Bellman optimal policy, then there exists a state i and an action a = π i such that B a i (h π ) = 0.

C.1 Arithmetic complexity of Subroutine 3 and memory usage

Recall that in Subroutine 3, we compute the values X ij by doing the update (for all iteration k, for all = 1 to k and for all i ∈ [n] or all i ∈ π +1 if we do not test indexability):

X +1 iσ k = X iσ k - X iσ 1 + X σ σ X σ σ k (C7)
If we test indexability, there are Let W i := X iσ /(1 + X σ σ ) and V i := X iσ k . Using this, Equation (C7) can be rewritten as:

V i = V i -W i V σ . (C8)
This results in the following loop at iteration k:

• Initialize V i from X 1 iσ k . • For all ∈ {1, . . . , k -1}, and all i ∈ [n] (or i ∈ π +1 ), apply (C8).

• Compute W ik = V i /(1 + V σ k ) Note that the value of V is not necessary for iteration k (only the values of W i are needed). This shows that the algorithm can be implemented with a memory O(n 2 ).

C.2 Speedup when not checking the indexability: First found go last

When the indexability is not tested, the update (C8) is computed for all i ∈ π +1 . This creates inefficiencies (due to inefficient cache usage) because the elements V i are not accessed sequentially.

To speedup the memory accesses, our solution is to sort the items during the execution of the algorithm. At iteration k, the algorithm computes σ k . When this is done, our implementation switches all quantities in positions σ k and n -k + 1. These quantities are δ, y, z, W and X. For instance, once σ 1 is found, we know that the state at position n is state n and we do the following switches:

δ σ 1 , δ n → δ n , δ σ 1 y 1 σ 1 , y 1 n → y 1 n , y 1 σ 1 z 1 σ 1 , z 1 n → z 1 n , z 1 σ 1
W σ 1 : , W n: → W n: , W σ 1 : and X σ 1 : , X n: → X n: , X σ 1 : . To do so, we need an array to store all states such that at iteration k, the first n -k states of the array are the active states. We will need to track the position of each state in such array.

E.1 Comparison with [34]

The paper presents an algorithm that computes Whittle indices in O(n 3 ) (no explicit constant before n 3 is given) for all indexable problems. Despite following a different approach, our algorithm for computing Whittle index can be viewed as a refinement of this work. Let us recall once again that our approach also allows one to check the indexability of general restless bandits.

In the following, we show how we can refine the work of [START_REF] Akbarzadeh | Conditions for indexability of restless bandits and an O(K 3 ) algorithm to compute Whittle index[END_REF] to obtain an algorithm that is exactly the same as ours. Let D π and N π be two vectors defined as in [START_REF] Akbarzadeh | Conditions for indexability of restless bandits and an O(K 3 ) algorithm to compute Whittle index[END_REF] (we use the same notation, D π and N π , as the cited paper), D π = (1 -β)(I -βP π ) -1 r π , and N π = (1 -β)(I -βP π ) -1 π.

Then, we have D π -λN π = (1 -β)u π (λ) where u π (λ) is defined as in [START_REF] Chen | Linear programming for finite state multiarmed bandit problems[END_REF]. In our proposition, at each iteration k, we compute µ k i by Line 13. Instead, it is defined in [START_REF] Akbarzadeh | Conditions for indexability of restless bandits and an O(K 3 ) algorithm to compute Whittle index[END_REF] by two steps:

1. for all state j ∈ [n] such that N π k \{i} j =N π k j , one needs to compute

µ k ij = D π k \{i} j -D π k j N π k \{i} j -N π k j 2. compute µ k i = arg min j∈[n]:N π k \{i} j =N π k j µ k ij .
From [START_REF] Akbarzadeh | Conditions for indexability of restless bandits and an O(K 3 ) algorithm to compute Whittle index[END_REF]Theorem 2], in an indexable problem, for state σ k , there exists a state j ∈ [n] such that N π k \{σ k } j = N π k j . Now, suppose that for any active state i ∈ π k , there exists j ∈ [n] such that N π k \{i} j = N π k j . Using the Sherman-Morrison formula, we have6

D π k \{i} -D π k = - (1 -β)δ i + ∆i D π k 1 + ∆i [(I -βP π k ) -1 ] :i [(I -βP π k ) -1 ] :i ,
and

N π k \{i} -N π k = - (1 -β) + ∆i N π k 1 + ∆i [(I -βP π k ) -1 ] :i [(I -βP π k ) -1 ] :i .
Then, for any j ∈ [n] such that N

π k \{i} j = N π k j , µ k ij = (1 -β)δ i + ∆i D π k (1 -β) + ∆i N π k
which does not depend on j. Then, we simply have

µ k i = (1 -β)δ i + ∆i D π k (1 -β) + ∆i N π k . Also, we have ∆i N π k = (1 -β) ∆i (I -βP π k ) -1 π k = -(1 -β)y k i and ∆i D π k = (1 -β) ∆i u π k (µ k-1 min ) + µ k-1 min ∆i N π k = (1 -β)z k-1 i -(1 -β)µ k-1 min y k i .
So, replacing these terms in µ k i , we get the formula in Equation ( 14) of our work.

Note that the algorithm of [START_REF] Akbarzadeh | Conditions for indexability of restless bandits and an O(K 3 ) algorithm to compute Whittle index[END_REF] was only developed for the discounted case. Our approach for the time-average reward case is different because we use the active advantage function defined in (5) instead of working with the expected discounted total reward D π k and total number of activations N π k under policy π k . Note that the counterpart of D π k in undiscounted MDP is the average reward g π k and as we have seen in Appendix A.1, utilizing average reward optimality is not rich enough for undiscounted MDPs with transient states. In addition, our code is also optimized to avoid unnecessary computation and to reduce memory usage. Finally, the way we do the update of our matrix X makes it possible to obtain a subcubic algorithm whereas their approach does not (see also below).

E.2 Comparison with the algorithm of [29]

The algorithm [START_REF] Niño-Mora | A fast-pivoting algorithm for Whittle's restless bandit index[END_REF] has the best complexity up to date for discounted restless bandit. There is a square matrix A that plays a similar role as the square matrix X in our proposed algorithm. The most costly operations in the algorithm of [START_REF] Niño-Mora | A fast-pivoting algorithm for Whittle's restless bandit index[END_REF] is to update their matrix A at each iteration and it is done by Equation ( 21) that we recall here (using the same notation A as the cited paper):

for i, j ∈ π k , A k+1 ij = A k ij - A k iσ k A k σ k σ k A k σ k j . (E9)
This incurs a total complexity of (2/3)n 3 + O(n 2 ) arithmetic operations. As mentioned in Section 4.3, if we updated X k+1 as given by ( 21), our algorithm would also have a (2/3)n 3 +O(n 2 ) complexity but this version of update cannot be optimized by using fast matrix multiplication.

E.3 Their approach cannot be directly transformed into a subcubic algorithm

In addition to all the previously cited differences, one of the major contribution of our algorithm with respect to [START_REF] Niño-Mora | A fast-pivoting algorithm for Whittle's restless bandit index[END_REF][START_REF] Akbarzadeh | Conditions for indexability of restless bandits and an O(K 3 ) algorithm to compute Whittle index[END_REF] is that the most advanced version of our algorithm runs in a subcubic time. The approach7 used in [START_REF] Niño-Mora | A fast-pivoting algorithm for Whittle's restless bandit index[END_REF][START_REF] Akbarzadeh | Conditions for indexability of restless bandits and an O(K 3 ) algorithm to compute Whittle index[END_REF] is to update the full matrix X k+1 at iteration k, by using (E9). This idea is represented in Figure E5(a): for a given , their algorithm compute X :σ k for all k (i.e., the full vertical lines represented by arrows). Our first Subroutine 3 uses an horizontal approach based on [START_REF] Glazebrook | Index policies for the admission control and routing of impatient customers to heterogeneous service stations[END_REF], which we recall here:

X +1 iσ k = X iσ k - X iσ 1 + X σ σ X σ σ k .
At iteration k + 1, we use X 1 :σ k to compute all values of X :σ k up to = k + 1. This is represented in Figure E5(b). Our approach can be used to obtain the subcubic algorithm illustrated in Figure E5(c) by using subcubic algorithms for multiplication.

σ k σ k σ k K K σ K σ 2K
. . . This leads to the next fundamental question: why should the computation of Whittle index be harder than matrix inversion (or multiplication)? To us, the main difference is that when computing Whittle indices, the permutation σ is not known a priori but discovered as the algorithm progresses: σ k is only known at iteration k. Hence, while all terms of the matrices X k ij are not needed, it is difficult to know a priori which ones are needed and which ones are not. Hence, a simple divide and conquer algorithm cannot be used. This is why when recomputing X k+1 in Subroutine 4, we recompute the whole matrix (the vertical blue line) and not just the part that will be used to compute the gray zone: we do not know a priori what part of X k+1 ij will be useful or not.

21 :

 21 return the arm is indexable and the indices are {λ i } i∈[n] .

  Non-indexable arm with 3 states. The algorithm stops at iteration 3 because α π 3 3 (µ 3 min ) > 0 (the green line at the zone circled with red ellipse).

Fig. 3 :

 3 Fig. 3: The active advantage α π k i (λ) computed by the algorithm, for the two examples of Figure 2.

Fig. 5 :

 5 (a) Computation load of Subroutine 3. (b) Computation load of Subroutine 4 Illustration of the improvement proposed by replacing Subroutine 3

8

  , numbers from [29] FPA-Julia (a) Alg.3 with index. test (b) Alg.3 without test (c) Alg.4 with index. test (d) Alg.4 without test (a) Runtime of our implementations as a function of the state size n. Speed-up over FPA-Julia [29] (a) Alg.3 with index. test (b) Alg.3 without test (c) Alg.4 with index. test (d) Alg.4 without test (b) Speedup of each implementation over the data from [29] as a function of the state size n.

  Ambiguous examples: All transitions are deterministic and labels on transitions indicate rewards. Solid black arrows correspond to the action "activate" and dashed red arrows to the action "rest".

  non-indexable arm. Fig.

  Example of an indexable multichain problem with infinite Whittle index. Transitions are deterministic and labels on edges indicate rewards (for the λ-penalized arm). Solid black transitions correspond to action "activate" and dashed red to the action "rest".

  of the non-indexable arm presented in Figure2b is

n k=1 kn = n 3 / 2 +

 32 O(n 2 ) such updates. If we do not test indexability, there aren k=1 k =1 (n -) = n 3 /3 + O(n 2 )such updates. Below, we show each update of Equation (C7) can be done in two arithmetic operations (one addition and one multiplication), which leads to the complexity of n 3 + O(n 2 ) (or (2/3)n 3 + O(n 2 )) arithmetic operations for the computation of all the needed X k ij . We also show how to reduce the memory size to O(n 2 ).

  Comparison of the computation load of (E9) used in[START_REF] Niño-Mora | A fast-pivoting algorithm for Whittle's restless bandit index[END_REF][START_REF] Akbarzadeh | Conditions for indexability of restless bandits and an O(K 3 ) algorithm to compute Whittle index[END_REF] with the one of Subroutine 3 and Subroutine 4.

  Bellman optimal for k 1

	min
	Active State
	Passive State
	Fig. 4: Illustration of what happens when µ k-1 min < µ k min and when the test of Line 8 is successful (recall that the function α π k
	Suppose that the induction is true until iteration k -1 and that the algorithm
	completes iteration k. If µ k-1 min = µ k min , (A) is trivial and (B) is a direct consequence of the definition of µ k min and of Lemma 2. Consider now that µ k-1 min < µ k min and observe
	what happens in Figure

i (λ) is affine in λ). The black lines are the advantage functions α π k i (λ) of active state i ∈ π k . The dashed red lines are the advantage functions α π k i (λ) of passive state i ∈ π k .

  return the arm is indexable and the indices are{λ i } i∈[n] .

	19:	end if
	20:	if µ k min = +∞ then
	25: end for
	26: if π n+1 is multichain then
	27:	return the arm is multichain
	28: end if
	29:

21:

Set λ i := +∞ for all i ∈ π k 22: return the arm is indexable and the indices are {λ i } i∈[n] . 23: end if 24:

Set λ σ k = µ k min and π k+1 := π k \ {σ k }

Table 1 :

 1 Running time (in seconds) of the variants of Algorithm 2

		O(n 3 ) algorithm (Subroutine 3)	Subcubic algorithm (Subroutine 4)
	n	(a) With index. test (b) Without test	(c) With index. test (d) Without test
	100	0.006	0.004	0.007	0.005
	1000	0.2	0.2	0.2	0.2
	2000	1.9	1.2	1.1	1.0
	3000	7.2	5.0	3.2	2.6
	4000	17	12	8	6
	5000	34	25	16	12
	6000	60	43	27	20
	7000	95	68	42	33
	8000	142	99	63	49
	9000	199	141	92	70
	10000	275	191	122	95
	11000	361	257	164	133
	12000	471	339	225	190
	13000	620	428	286	243
	14000	790	553	408	317
	15000	965	685	501	403

Table 2 :

 2 Number of indexable problems among 100 000 randomly generated problems.

	Problem size n Tridiagonal 5-diagonal 7-diagonal	Dense
	3	98 731	-	-	99 883
	4	95 067	99 655	-	99 931
	5	89 198	99 309	99 902	99 969
	10	54 129	90 377	98 914 100 000
	30	7 094	29 699	66 143 100 000
	50	1 823	9 332	32 069 100 000

  any recurrent state i of π. Multiply this with P π ki and sum over i ∈ S (if i is not recurrent,

	then P π ki = 0) gives
	P π ki (r π i -g * i ) =	P π ki h i -		P π ki		P π ij h j
	i∈S	i∈S	i∈S	j∈S
		=	j∈S	P π kj hj since	P π P π =	P π .

If the MDP is unichain, then the bias vector h * is unique up to an additive constant. This is in general not the case for multichain MDPs.

not to be confused with β-discounted MDPs, where the discount is on rewards and not on actions.

We write n γ which is possibly non-integer. For the sake of simplicity, we write n γ and n1-γ but they should be understood as n γ and n 1-γ respectively.

Note that the definition of ∆ is identical to ∆ except when j = 1, for which ∆i1 := 0 but ∆i1 := β(P 1 i1 -P 0 i1 ).

the expression of D π\{i} and N π\{i} given by Equation (18) in[START_REF] Akbarzadeh | Conditions for indexability of restless bandits and an O(K 3 ) algorithm to compute Whittle index[END_REF] are erroneous.

Equation (18) of[START_REF] Akbarzadeh | Conditions for indexability of restless bandits and an O(K 3 ) algorithm to compute Whittle index[END_REF], which is central to their algorithm is the same as the above equation (E9).

Acknowledgments. This work is supported by the French National Research Agency (ANR) through REFINO Project under Grant ANR-19-CE23-0015.

Declarations

Appendix A Examples and counterexamples

In this section, we provide a few examples to illustrate the ambiguities in the classical definition of indexability, and to illustrate what can happen for some multichain arms. We also provide the parameters of arms presented in Figure 2.

A.1 Discussion on the definition of indexability

The classical notion of indexability used in the literature is to say that the optimal policy π * (λ) should be non-increasing in λ. Yet, we argue that this definition has two problems:

1. What does "increasing" mean when π * (λ) is not unique? Two possibilities are: for all penalties λ < λ : (∃) there exist policies π, π with π optimal for λ and π optimal for λ such that π ⊇ π ; (∀) for all policies π, π such that π is optimal for λ and π is optimal for λ , we have π ⊇ π . 2. What notion of "optimality" should be used? Two possibilities are:

Proof Let θ = π be another Bellman optimal policy. Since θ is gain optimal and h π ∈ H, Lemma 4 implies that B θi i (h π ) = 0 for all i ∈ R θ . If there exists i ∈ R θ such that θ i = π i , then the proof is concluded. Otherwise, θ i = π i for all i ∈ R θ . This show that π and θ coincide for all recurrent states of θ and that R θ = R π . Moreover, as π is unichain, θ is also unichain. Hence, Lemma 5 implies that B θi i (h π ) = 0 for all i. Since θ = π, there exists at least one state i ∈ S such that θ i = π i .

B.2 Unichain property

Lemma 7. Given a two-action MDP [n], {0, 1}, r, P , let P π be the transition matrix under a Bellman optimal policy π. Policy π is unichain if and only if the matrix

Proof A π is not invertible if there exists a column vector u = 0 such that u A π = 0.

We prove that such u does not exist when policy π is unichain. Let u ∈ R n be an arbitrary vector such that u A π = 0. Then, we have

Combining the above equation with j P π ji = 1, we get:

where we used that n i=1 u i = 0 to obtain the last equality. This shows that

The set of vector u such that u P π is a vector space. It is of dimension 1 if and only if π is unichain, in which case the vector u verifying u P π = u are multiples of a stationary distribution under policy π [START_REF] Puterman | Markov Decision Processes: Discrete Stochastic Dynamic Programming[END_REF]. Thus, if the policy π induces a unichain Markov chain, then n i=1 u i = 0 implies u = 0. If policy π is not unichain, there exists u = 0 such that n i=1 u i = 0.

Appendix C Implementations

Appendix D Analysis of the experimental time to solve a linear system

In this section, we report in Figure D4 the time taken by the default implementation to solve a linear system of the form AX = B where A and B are two square matrices. To obtain this figure, we generated random (full) matrices where each entry is between 0 and 1 and use the function scipy.linalg.solve from the library scipy. The reported numbers suggest that the complexity of the solver is closer to O(n 2.8 ) than to O(n 3 ), although we agree that the difference between the O(n 2.8 ) and the O(n 3 ) curves is small. Note that this is in accordance with the papers [START_REF] Huang | Strassen's algorithm reloaded[END_REF][START_REF] Huang | Practical fast matrix multiplication algorithms[END_REF]

Appendix E Detailed comparison with [34] and [29]

In this section, we compare our algorithm with two main related works for finite-state restless bandits problem.