Finite-Time Analysis of Whittle Index based Q-Learning for Restless
  Multi-Armed Bandits with Neural Network Function Approximation
Whittle index policy is a heuristic to the intractable restless multi-armed
bandits (RMAB) problem. Although it is provably asymptotically optimal, finding
Whittle indices remains difficult. In this paper, we present Neural-Q-Whittle,
a Whittle index based Q-learning algorithm for RMAB with neural network
function approximation, which is an example of nonlinear two-timescale
stochastic approximation with Q-function values updated on a faster timescale
and Whittle indices on a slower timescale. Despite the empirical success of
deep Q-learning, the non-asymptotic convergence rate of Neural-Q-Whittle, which
couples neural networks with two-timescale Q-learning largely remains unclear.
This paper provides a finite-time analysis of Neural-Q-Whittle, where data are
generated from a Markov chain, and Q-function is approximated by a ReLU neural
network. Our analysis leverages a Lyapunov drift approach to capture the
evolution of two coupled parameters, and the nonlinearity in value function
approximation further requires us to characterize the approximation error.
Combing these provide Neural-Q-Whittle with $\mathcal{O}(1/k^{2/3})$
convergence rate, where $k$ is the number of iterations.
Guojun Xiong, Jian Li
Finite-Time Analysis of Whittle Index based Q-Learning for Restless Multi-Armed Bandits with Neural Network Function Approximation

###### Abstract

Whittle index policy is a heuristic to the intractable restless multi-armed bandits (RMAB) problem. Although it is provably asymptotically optimal, finding Whittle indices remains difficult. In this paper, we present Neural-Q-Whittle, a Whittle index based Q-learning algorithm for RMAB with neural network function approximation, which is an example of nonlinear two-timescale stochastic approximation with Q-function values updated on a faster timescale and Whittle indices on a slower timescale. Despite the empirical success of deep Q-learning, the non-asymptotic convergence rate of Neural-Q-Whittle, which couples neural networks with two-timescale Q-learning largely remains unclear. This paper provides a finite-time analysis of Neural-Q-Whittle, where data are generated from a Markov chain, and Q-function is approximated by a ReLU neural network. Our analysis leverages a Lyapunov drift approach to capture the evolution of two coupled parameters, and the nonlinearity in value function approximation further requires us to characterize the approximation error. Combing these provide Neural-Q-Whittle with \(\mathcal{O}(1/k^{2/3})\) convergence rate, where \(k\) is the number of iterations.

## 1 Introduction

We consider the restless multi-armed bandits (RMAB) problem [56], where the decision maker (DM) repeatedly activates \(K\) out of \(N\) arms at each decision epoch. Each arm is described by a Markov decision process (MDP) [45], and evolves stochastically according to two different transition kernels, depending on whether the arm is activated or not. Rewards are generated with each transition. Although RMAB has been widely used to study constrained sequential decision making problems [5; 38; 12; 64; 29; 36; 35; 25], it is notoriously intractable due to the explosion of state space [44]. A celebrated heuristic is the Whittle index policy [56], which computes the Whittle index for each arm given its current state as the cost to pull the arm. Whittle index policy then activates the \(K\) highest indexed arms at each decision epoch, and is provably asymptotically optimal [54].

However, the computation of Whittle index requires full knowledge of the underlying MDP associated with each arm, which is often unavailable in practice. To this end, many recent efforts have focused on learning Whittle indices for making decisions in an online manner. First, model-free reinforcement learning (RL) solutions have been proposed [10; 23; 53; 8; 28; 57; 59; 61; 58; 3], among which [3] developed a Whittle index based Q-learning algorithm, which we call Q-Whittle for ease of exposition, and provided the first-ever rigorous asymptotic analysis. However, Q-Whittle suffers from slow convergence since it only updates the Whittle index of a specific state when that state is visited. In addition, Q-Whittle needs to store the Q-function values for all state-action pairs, which limits its applicability only to problems with small state space. Second, deep RL methods have been leveraged to predict Whittle indices via training neural networks [40; 41]. Though
these methods are capable of dealing with large state space, there is no asymptotic or finite-time performance guarantee. Furthermore, training neural networks requires to tuning hyper-parameters. This introduces an additional layer of complexity to predict Whittle indices. Third, to address aforementioned deficiencies, [60] proposed Q-Whittle-LFA by coupling Q-Whittle with linear function approximation and provided a finite-time convergence analysis. One key limitation of Q-Whittle-LFA is the unrealistic assumption that all data used in Q-Whittle-LFA are sampled i.i.d. from a fixed stationary distribution.

To tackle the aforementioned limitations and inspired by the empirical success of deep Q-learning in numerous applications, we develop Neural-Q-Whittle, a Whittle index based Q-learning algorithm with _neural network function approximation under Markovian observations_. Like [3; 60], the updates of Q-function values and Whittle indices form a two-timescale stochastic approximation (2TSA) with the former operating on a faster timescale and the later on a slower timescale. Unlike [3; 60], our Neural-Q-Whittle uses a deep neural network with the ReLU activation function to approximate the Q-function. However, Q-learning with neural network function approximation can in general diverge [2], and the theoretical convergence of Q-learning with neural network function approximation has been limited to special cases such as fitted Q-iteration with i.i.d. observations [22], which fails to capture the practical setting of Q-learning with neural network function approximation.

In this paper, we study the non-asymptotic convergence of Neural-Q-Whittle with data generated from a Markov decision process. Compared with recent theoretical works for Q-learning with neural network function approximation [13; 22; 62], our Neural-Q-Whittle involves a two-timescale update between two coupled parameters, i.e., Q-function values and Whittle indices. This renders existing finite-time analysis in [13; 22; 62] not applicable to our Neural-Q-Whittle due to the fact that [13; 22; 62] only contains a single-timescale update on Q-function values. Furthermore, [13; 22; 62] required an additional projection step for the update of parameters of neural network function so as to guarantee the boundedness between the unknown parameter at any time step with the initialization. This in some cases is impractical. Hence, a natural question that arises is

_Is it possible to provide a non-asymptotic convergence rate analysis of Neural-Q-Whittle with two coupled parameters updated in two timescales under Markovian observations without the extra projection step?_

The theoretical convergence guarantee of two-timescale Q-learning with neural network function approximation under Markovian observations remains largely an open problem, and in this paper, we provide an affirmative answer to this question. Our main contributions are summarized as follows:

\(\bullet\) We propose Neural-Q-Whittle, a novel Whittle index based Q-learning algorithm with neural network function approximation for RMAB. Inspired by recent work on TD learning [48] and Q-learning [15] with linear function approximation, our Neural-Q-Whittle removes the additional impractical projection step in the neural network function parameter update.

\(\bullet\) We establish the first finite-time analysis of Neural-Q-Whittle under Markovian observations. Due to the two-timescale nature for the updates of two coupled parameters (i.e., Q-function values and Whittle indices) in Neural-Q-Whittle, we focus on the convergence rate of these parameters rather than the convergence rate of approximated Q-functions as in [13; 22; 62]. Our key technique is to view Neural-Q-Whittle as a 2TSA for finding the solution of suitable nonlinear equations. Different from recent works on finite-time analysis of a general 2TSA [20] or with linear function approximation [60], the nonlinear parameterization of Q-function in Neural-Q-Whittle under Markovian observations imposes significant difficulty in finding the global optimum of the corresponding nonlinear equations. To mitigate this, we first approximate the original neural network function with a collection of local linearization and focus on finding a surrogate Q-function in the neural network function class that well approximates the optimum. Our finite-time analysis then requires us to consider two Lyapunov functions that carefully characterize the coupling between iterates of Q-function values and Whittle indices, with one Lyapunov function defined with respect to the true neural network function, and the other defined with respect to the locally linearized neural network function. We then characterize the errors between these two Lyapunov functions. Putting them together, we prove that Neural-Q-Whittle achieves a convergence in expectation at a rate \(\mathcal{O}(1/k^{2/3})\), where \(k\) is the number of iterations.

\(\bullet\) Finally, we conduct experiments to validate the convergence performance of Neural-Q-Whittle, and verify the sufficiency of our proposed condition for the stability of Neural-Q-Whittle.
## 2 Preliminaries

**RMAB.** We consider an infinite-horizon average-reward RMAB with each arm \(n\in\mathcal{N}\) described by a unichain MDP [45]\(\mathcal{M}_{n}:=(\mathcal{S},\mathcal{A},P_{n},r_{n})\), where \(\mathcal{S}\) is the state space with cardinality \(S<\infty\), \(\mathcal{A}\) is the action space with cardinality \(A\), \(P_{n}(s^{\prime}|s,a)\) is the transition probability of reaching state \(s^{\prime}\) by taking action \(a\) in state \(s\), and \(r_{n}(s,a)\) is the reward associated with state-action pair \((s,a)\). At each time slot \(t\), the DM activates \(K\) out of \(N\) arms. Arm \(n\) is "active" at time \(t\) when it is activated, i.e., \(A_{n}(t)=1\); otherwise, arm \(n\) is "passive", i.e., \(A_{n}(t)=0\). Let \(\Pi\) be the set of all possible policies for RMAB, and \(\pi\in\Pi\) is a feasible policy, satisfying \(\pi:\mathcal{F}_{t}\mapsto\mathcal{A}^{N}\), where \(\mathcal{F}_{t}\) is the sigma-algebra generated by random variables \(\{S_{n}(h),A_{n}(h):\forall n\in\mathcal{N},h\leq t\}\). The objective of the DM is to maximize the expected long-term average reward subject to an instantaneous constraint that only \(K\) arms can be activated at each time slot, i.e.,

\[\text{RMAB:}\quad\max_{\pi\in\Pi}\ \liminf_{T\to\infty}\frac{1}{T}\mathbb{E}_{ \pi}\left(\sum_{t=0}^{T}\sum_{n=1}^{N}r_{n}(t)\right),\quad\text{s.t.}\ \sum_{n=1}^{N}A_{n}(t)=K,\quad\forall t. \tag{1}\]

**Whittle Index Policy.** It is well known that RMAB (1) suffers from the curse of dimensionality [44]. To address this challenge, Whittle [56] proposed an index policy through decomposition. Specifically, Whittle relaxed the constraint in (1) to be satisfied on average and obtained a unconstrained problem: \(\max_{\pi\in\Pi}\liminf_{T\to\infty}\frac{1}{T}\mathbb{E}_{\pi}\sum_{t=1}^{T} \sum_{n=1}^{N}\{r_{n}(t)+\lambda(1-A_{n}(t))\}\), where \(\lambda\) is the Lagrangian multiplier associated with the constraint. The key observation of Whittle is that this problem can be decomposed and its solution is obtained by combining solutions of \(N\) independent problems via solving the associated dynamic programming (DP): \(V_{n}(s)=\max_{a\in\{0,1\}}Q_{n}(s,a),\forall n\in\mathcal{N},\) where

\[Q_{n}(s,a)+\beta=a\Big{(}r_{n}(s,a)+\sum_{s^{\prime}}p_{n}(s^{\prime}|s,1)V_{n} (s^{\prime})\Big{)}+(1-a)\Big{(}r_{n}(s,a)+\lambda+\sum_{s^{\prime}}p_{n}(s^{ \prime}|s,0)V_{n}(s^{\prime})\Big{)}, \tag{2}\]

where \(\beta\) is unique and equals to the maximal long-term average reward of the unichain MDP, and \(V_{n}(s)\) is unique up to an additive constant, both of which depend on the Lagrangian multiplier \(\lambda\). The optimal decision \(a^{*}\) in state \(s\) then is the one which maximizes the right hand side of the above DP. The Whittle index associated with state \(s\) is defined as the value \(\lambda_{n}^{*}(s)\in\mathbb{R}\) such that actions \(0\) and \(1\) are equally favorable in state \(s\) for arm \(n\)[3; 23], satisfying

\[\lambda_{n}^{*}(s):=r_{n}(s,1)+\sum_{s^{\prime}}p_{n}(s^{\prime}|s,1)V_{n}(s^{ \prime})-r_{n}(s,0)-\sum_{s^{\prime}}p_{n}(s^{\prime}|s,0)V_{n}(s^{\prime}). \tag{3}\]

Whittle index policy then activates \(K\) arms with the largest Whittle indices at each time slot. Additional discussions are provided in Section B in supplementary materials.

**Q-Learning for Whittle Index.** Since the underlying MDPs are often unknown, [3] proposed Q-Whittle, a tabular Whittle index based Q-learning algorithm, where the updates of Q-function values and Whittle indices form a 2TSA, with the former operating on a faster timescale for a given \(\lambda_{n}\) and the later on a slower timescale. Specifically, the Q-function values for \(\forall n\in\mathcal{N}\) are updated as

\[Q_{n,k+1}(s,a):=Q_{n,k}(s,a)+\alpha_{n,k}\mathds{1}_{\{S_{n,k}=s,A_{n,k}=a\}}\Big{(}r_{n}(s,a)+(1-a)\lambda_{n,k}(s)\\ +\max_{a}Q_{n,k}(S_{n,k+1},a)-I_{n}(Q_{k})-Q_{n,k}(s,a)\Big{)}, \tag{4}\]

where \(I_{n}(Q_{k})=\frac{1}{2S}\sum_{s\in\mathcal{S}}(Q_{n,k}(s,0)+Q_{n,k}(s,1))\) is standard in the relative Q-learning for long-term average MDP setting [1], which differs significantly from the discounted reward setting [45; 1]. \(\{\alpha_{n,k}\}\) is a step-size sequence satisfying \(\sum_{k}\alpha_{n,k}=\infty\) and \(\sum_{k}\alpha_{n,k}^{2}<\infty\).

Accordingly, the Whittle index is updated as

\[\lambda_{n,k+1}(s)=\lambda_{n,k}(s)+\eta_{n,k}(Q_{n,k}(s,1)-Q_{n,k}(s,0)), \tag{5}\]

with the step-size sequence \(\{\eta_{n,k}\}\) satisfying \(\sum_{k}\eta_{n,k}=\infty\), \(\sum_{k}\eta_{n,k}^{2}<\infty\) and \(\eta_{n,k}=o(\alpha_{n,k})\). The coupled iterates (4) and (5) form a 2TSA, and [3] provided an asymptotic convergence analysis.

## 3 Neural Q-Learning for Whittle Index

A closer look at (5) reveals that Q-Whittle only updates the Whittle index of a specific state when that state is visited. This makes Q-Whittle suffers from slow convergence. In addition, Q-Whittle
needs to store the Q-function values for all state-action pairs, which limits its applicability only to problems with small state space. To address this challenge and inspired by the empirical success of deep Q-learning, we develop Neural-Q-Whittle through coupling Q-Whittle with neural network function approximation by using low-dimensional feature mapping and leveraging the strong representation power of neural networks. For ease of presentation, we drop the subscript \(n\) in (4) and (5), and discussions in the rest of the paper apply to any arm \(n\in\mathcal{N}\).

Specifically, given a set of basis functions \(\phi_{\ell}:\mathcal{S}\times\mathcal{A}\mapsto\mathbb{R},\forall\ell=1, \cdots,d\) with \(d\ll SA\), the approximation of Q-function \(Q_{\boldsymbol{\theta}}(s,a)\) parameterized by a unknown weight vector \(\boldsymbol{\theta}\in\mathbb{R}^{md}\), is given by \(Q_{\boldsymbol{\theta}}(s,a)=f(\boldsymbol{\theta};\boldsymbol{\phi}(s,a)),\ \forall s \in\mathcal{S},a\in\mathcal{A}\), where \(f\) is a nonlinear neural network function parameterized by \(\boldsymbol{\theta}\) and \(\boldsymbol{\phi}(s,a)\), with \(\boldsymbol{\phi}(s,a)=(\phi_{1}(s,a),\cdots,\phi_{d}(s,a))^{\intercal}\). The feature vectors are assumed to be linearly independent and are normalized so that \(\|\boldsymbol{\phi}(s,a)\|\leq 1,\forall s\in\mathcal{S},a\in\mathcal{A}\). In particular, we parameterize the Q-function by using a two-layer neural network [13; 62]

\[f(\boldsymbol{\theta};\boldsymbol{\phi}(s,a)):=\frac{1}{\sqrt{m}}\sum_{r=1}^{ m}b_{r}\sigma(\mathbf{w}_{r}^{\intercal}\boldsymbol{\phi}(s,a)), \tag{6}\]

where \(\boldsymbol{\theta}=(b_{1},\ldots,b_{m},\mathbf{w}_{1}^{\intercal},\ldots, \mathbf{w}_{m}^{\intercal})^{\intercal}\) with \(b_{r}\in\mathbb{R}\) and \(\mathbf{w}_{r}\in\mathbb{R}^{d\times 1},\forall r\in[1,m]\). \(b_{r},\forall r\) are uniformly initialized in \(\{-1,1\}\) and \(w_{r},\forall r\) are initialized as a zero mean Gaussian distribution according to \(\mathcal{N}(\mathbf{0},\mathbf{I}_{d}/d)\). During training process, only \(\mathbf{w}_{r},\forall r\) are updated while \(b_{r},\forall r\) are fixed as the random initialization. Hence, we use \(\boldsymbol{\theta}\) and \(\mathbf{w}_{r},\forall r\) interchangeably throughout this paper. \(\sigma(x)=\max(0,x)\) is the rectified linear unit (ReLU) activation function1.

Footnote 1: The finite-time analysis of Deep Q-Networks (DQN) [13; 22; 62] and references therein focuses on the ReLU activation function, as it has certain properties that make the analysis tractable. ReLU is piecewise linear and non-saturating, which can simplify the mathematical analysis. Applying the same analysis to other activation functions like the hyperbolic tangent (tanh) could be more complex, which is out of the scope of this work.

Given (6), we can rewrite the Q-function value updates in (4) as

\[\boldsymbol{\theta}_{k+1}=\boldsymbol{\theta}_{k}+\alpha_{k}\Delta_{k}\nabla _{\boldsymbol{\theta}}f(\boldsymbol{\theta}_{k};\boldsymbol{\phi}(S_{k},A_{k} )), \tag{7}\]

with \(\Delta_{k}\) being the temporal difference (TD) error defined as \(\Delta_{k}:=r(S_{k},A_{k})+(1-A_{k})\lambda_{k}(s)-I(\boldsymbol{\theta}_{k})+ \max_{a}f(\boldsymbol{\theta}_{k};\boldsymbol{\phi}(S_{k+1},a))-f(\boldsymbol {\theta}_{k};\boldsymbol{\phi}(S_{k},A_{k}))\), where \(I(\boldsymbol{\theta}_{k})=\frac{1}{2S}\sum_{s\in\mathcal{S}}[f(\boldsymbol{ \theta}_{k};\boldsymbol{\phi}(s,0))+f(\boldsymbol{\theta}_{k};\boldsymbol{ \phi}(s,1))]\). Similarly, the Whittle index update (5) can be rewritten as

\[\lambda_{k+1}(s)=\lambda_{k}(s)+\eta_{k}(f(\boldsymbol{\theta}_{k};\boldsymbol {\phi}(s,1))-f(\boldsymbol{\theta}_{k};\boldsymbol{\phi}(s,0))). \tag{8}\]

The coupled iterates in (7) and (8) form Neural-Q-Whittle as summarized in Algorithm 1, which aims to learn the coupled parameters \((\boldsymbol{\theta}^{\star},\lambda^{\star}(s))\) such that \(f(\boldsymbol{\theta}^{\star},\boldsymbol{\phi}(s,1))=f(\boldsymbol{\theta}^{ \star},\boldsymbol{\phi}(s,0)),\forall s\in\mathcal{S}\).

```
1:Input:\(\boldsymbol{\phi}(s,a)\) for \(\forall s\in\mathcal{S},a\in\mathcal{A}\), and learning rates \(\{\alpha_{k}\}_{k=1,\ldots,T}\), \(\{\eta_{k}\}_{k=1,\ldots,T}\)
2:Initialization:\(b_{r}\sim\text{Unif}(\{-1,1\}),\mathbf{w}_{r,0}\sim\mathcal{N}(\mathbf{0}, \mathbf{I}_{d}/d),\forall r\in[1,m]\) and \(\lambda(s)=0\), \(\forall s\in\mathcal{S}\)
3:for\(s\in\mathcal{S}\)do
4:for\(k=1,\ldots,T\)do
5: Sample \((S_{k},A_{k},S_{k+1})\) according to the \(\epsilon\)-greedy policy;
6:\(\Delta_{k}\!=\!r(S_{k},A_{k})+(1-A_{k})\lambda_{k}(s)+\max_{a}f(\boldsymbol{ \theta}_{k};\boldsymbol{\phi}(S_{k+1},a))-I(\boldsymbol{\theta}_{k})\!-\!f( \boldsymbol{\theta}_{k};\boldsymbol{\phi}(S_{k},A_{k}))\);
7:\(\boldsymbol{\theta}_{k+1}=\boldsymbol{\theta}_{k}+\alpha_{k}\Delta_{k}\nabla _{\boldsymbol{\theta}}f(\boldsymbol{\theta}_{k};\boldsymbol{\phi}(S_{k},A_{k}))\);
8:\(\lambda_{k+1}(s)=\lambda_{k}(s)+\eta_{k}(f(\boldsymbol{\theta}_{k};\boldsymbol{ \phi}(s,1))-f(\boldsymbol{\theta}_{k};\boldsymbol{\phi}(s,0))\);
9:endfor
10:endfor
11:Return:\(\lambda(s),\forall s\in\mathcal{S}\).
```

**Algorithm 1** Neural-Q-Whittle: Neural Q-Learning for Whittle Index

**Remark 1**.: _Unlike recent works for Q-learning with linear [6; 37; 67] or neural network function approximations [13; 22; 62], we do not assume an additional projection step of the updates of unknown parameters \(\boldsymbol{\theta}_{k}\) in (7) to confine \(\boldsymbol{\theta}_{k},\forall k\) into a bounded set. This projection step is often used to stabilize the iterates related to the unknown stationary distribution of the underlying Markov chain, which in some cases is impractical. More recently, [48] removed the extra projection step and established the finite-time convergence of TD learning, which is treated as a linear stochastic approximation algorithm. [15] extended it to the Q-learning with linear function approximation._
_However, these state-of-the-art works only contained a single-timescale update on Q-function values, i.e., with the only unknown parameter \(\mathbf{\theta}\), while our Neural-Q-Whittle involves a two-timescale update between two coupled unknown parameters \(\mathbf{\theta}\) and \(\lambda\) as in (7) and (8). Our goal in this paper is to expand the frontier by providing a finite-time bound for Neural-Q-Whittle under Markovian noise without requiring an additional projection step. We summarize the differences between our work and existing literature in Table 1._

## 4 Finite-Time Analysis of Neural-Q-Whittle

In this section, we present the finite-time analysis of Neural-Q-Whittle for learning Whittle index \(\lambda(s)\) of any state \(s\in\mathcal{S}\) when data are generated from a MDP. To simplify notation, we abbreviate \(\lambda(s)\) as \(\lambda\) in the rest of the paper. We start by first rewriting the updates of Neural-Q-Whittle in (7) and (8) as a nonlinear two-timescale stochastic approximation (2TSA) in Section 4.1.

### A Nonlinear 2TSA Formulation with Neural Network Function

We first show that Neural-Q-Whittle can be rewritten as a variant of the nonlinear 2TSA. For any fixed policy \(\pi\), since the state of each arm \(\{S_{k}\}\) evolves according to a Markov chain, we can construct a new variable \(X_{k}=(S_{k},A_{k},S_{k+1})\), which also forms a Markov chain with state space \(\mathcal{X}:=\{(s,a,s^{\prime})|s\in\mathcal{S},\pi(a|s)\geq 0,p(s^{\prime}|s,a )>0\}\). Therefore, the coupled updates (7) and (8) of Neural-Q-Whittle can be rewritten in the form of a nonlinear 2TSA [20]:

\[\mathbf{\theta}_{k+1}=\mathbf{\theta}_{k}+\alpha_{k}h(X_{k},\mathbf{\theta}_{k},\lambda_ {k}),\qquad\lambda_{k+1}=\lambda_{k}+\eta_{k}g(X_{k},\mathbf{\theta}_{k},\lambda_{ k}), \tag{9}\]

where \(\mathbf{\theta}_{0}\) and \(\lambda_{0}\) being arbitrarily initialized in \(\mathbb{R}^{md}\) and \(\mathbb{R}\), respectively; and \(h(\cdot)\) and \(g(\cdot)\) satisfy

\[h(X_{k},\mathbf{\theta}_{k},\lambda_{k}) :=\nabla_{\mathbf{\theta}}f(\mathbf{\theta}_{k};\mathbf{\phi}(S_{k},A_{k})) \Delta_{k},\qquad\mathbf{\theta}_{k}\in\mathbb{R}^{md},\lambda_{k}\in\mathbb{R}, \tag{10}\] \[g(X_{k},\mathbf{\theta}_{k},\lambda_{k}) :=f(\mathbf{\theta}_{k};\mathbf{\phi}(s,1))-f(\mathbf{\theta}_{k};\mathbf{\phi}(s,0)),\qquad\mathbf{\theta}_{k}\in\mathbb{R}^{md}. \tag{11}\]

Since \(\eta_{k}\ll\alpha_{k}\), the dynamics of \(\mathbf{\theta}\) evolves much faster than those of \(\lambda\). We aim to establish the finite-time performance of the nonlinear 2TSA in (9), where \(f(\cdot)\) is the neural network function defined in (6). This is equivalent to find the root2\((\mathbf{\theta}^{*},\lambda^{*})\) of a system with _two coupled_ nonlinear equations \(h:\mathcal{X}\times\mathbb{R}^{md}\times\mathbb{R}\to\mathbb{R}^{md}\) and \(g:\mathcal{X}\times\mathbb{R}^{md}\times\mathbb{R}\to\mathbb{R}\) such that

Footnote 2: The root \((\mathbf{\theta}^{*},\lambda^{*})\) of the nonlinear 2TSA (9) can be established by using the ODE method following the solution of suitably defined differential equations [9; 49; 3; 21; 19; 20], i.e., \(\hat{\mathbf{\theta}}=H(\mathbf{\theta},\lambda),\dot{\lambda}=\frac{\eta}{\alpha}G( \mathbf{\theta},\lambda)\), where a fixed stepsize is assumed for ease of expression at this moment.

\[H(\mathbf{\theta},\lambda):=\mathbb{E}_{\mu}[h(X,\mathbf{\theta},\lambda)]=0,\qquad G (\mathbf{\theta},\lambda):=\mathbb{E}_{\mu}[g(X,\mathbf{\theta},\lambda)]=0, \tag{12}\]

where \(X\) is a random variable in finite state space \(\mathcal{X}\) with unknown distribution \(\mu\). For a fixed \(\mathbf{\theta}\), to study the stability of \(\lambda\), we assume the condition on the existence of a mapping such that \(\lambda=y(\mathbf{\theta})\) is the unique solution of \(G(\mathbf{\theta},\lambda)=0\). In particular, \(y(\mathbf{\theta})\) is given as

\[y(\mathbf{\theta})=r(s,1)+\sum_{s^{\prime}}p(s^{\prime}|s,1)\max_{a}f(\mathbf{\theta}; \mathbf{\phi}(s^{\prime},a))-r(s,0)-\sum_{s^{\prime}}p(s^{\prime}|s,0)\max_{a}f(\bm {\theta};\mathbf{\phi}(s^{\prime},a)). \tag{13}\]

\begin{table}
\begin{tabular}{|c|c|c|c|c|} \hline
**Algorithm** & **Noise** & **Approximation** & **Timescale** & **Whittle index** \\ \hline \hline Q-Whittle [3] & _i.i.d._ & \(X\) & _two-timescale_ & ✓ \\ Q-Whittle-LFA [60; 3; 67] & _i.i.d._ & linear & _two-timescale_ & ✓ \\ Q-Learning-LFA [6; 37; 67] & _Markovian_ & linear & _single-timescale_ & ✗ \\ Q-Learning-NFA [13; 15; 22; 62] & _Markovian_ & _neural network_ & _single-timescale_ & ✗ \\ TD-Learning-LFA [48] & _Markovian_ & linear & _single-timescale_ & ✗ \\
2TSA-IID [19; 21] & _i.i.d._ & ✗ & _two-timescale_ & ✗ \\
2TSA-Markovian [20] & _Markovian_ & ✗ & _two-timescale_ & ✗ \\ \hline \hline Neural-Q-Whittle (**this work**) & _Markovian_ & _neural network_ & _two-timescale_ & ✓ \\ \hline \end{tabular}
\end{table}
Table 1: Comparison of settings in related works.
### Main Results

As inspired by [20], the finite-time analysis of such a nonlinear 2TSA boils down to the choice of two step sizes \(\{\alpha_{k},\eta_{k},\forall k\}\) and a Lyapunov function that couples the two iterates in (9). To this end, we first define the following two error terms:

\[\tilde{\mathbf{\theta}}_{k}=\mathbf{\theta}_{k}-\mathbf{\theta}^{*},\qquad\tilde{\lambda}_{ k}=\lambda_{k}-y(\mathbf{\theta}_{k}), \tag{14}\]

which characterize the coupling between \(\mathbf{\theta}_{k}\) and \(\lambda_{k}\). If \(\tilde{\mathbf{\theta}}_{k}\) and \(\tilde{\lambda}_{k}\) go to zero simultaneously, the convergence of \((\mathbf{\theta}_{k},\lambda_{k})\) to \((\mathbf{\theta}^{*},\lambda^{*})\) can be established. Thus, to prove the convergence of \((\mathbf{\theta}_{k},\lambda_{k})\) of the nonlinear 2TSA in (9) to its true value \((\mathbf{\theta}^{*},\lambda^{*})\), we can equivalently study the convergence of \((\tilde{\mathbf{\theta}}_{k},\tilde{\lambda}_{k})\) by providing the finite-time analysis for the mean squared error generated by (9). To couple the fast and slow iterates, we define the following weighted Lyapunov function

\[M(\mathbf{\theta}_{k},\lambda_{k}):=\frac{\eta_{k}}{\alpha_{k}}\|\tilde{\mathbf{\theta }}_{k}\|^{2}+\|\tilde{\lambda}_{k}\|^{2}=\frac{\eta_{k}}{\alpha_{k}}\|\mathbf{ \theta}_{k}-\mathbf{\theta}^{*}\|^{2}+\|\lambda_{k}-y(\mathbf{\theta}_{k})\|^{2}, \tag{15}\]

where \(\|\cdot\|\) stands for the the Euclidean norm for vectors throughout the paper. It is clear that the Lyapunov function \(M(\mathbf{\theta}_{k},\lambda_{k})\) combines the updates of \(\mathbf{\theta}\) and \(\lambda\) with respect to the true neural network function \(f(\mathbf{\theta};\mathbf{\phi}(s,a))\) in (6).

To this end, our goal turns to characterize finite-time convergence of \(\mathbb{E}[M(\mathbf{\theta}_{k},\lambda_{k})]\). However, it is challenging to directly finding the global optimum of the corresponding nonlinear equations due to the nonlinear parameterization of Q-function in Neural-Q-Whittle. In addition, the operators \(h(\cdot),g(\cdot)\) and \(y(\cdot)\) in (10), (11) and (13) directly relate with the convoluted neural network function \(f(\mathbf{\theta};\mathbf{\phi}(s,a))\) in (6), which hinders us to characterize the smoothness properties of theses operators. Such properties are often required for the analysis of stochastic approximation [15; 19; 21].

To mitigate this, **(Step 1)** we instead approximate the true neural network function \(f(\mathbf{\theta},\mathbf{\phi}(s,a))\) with a collection of local linearization \(f_{0}(\mathbf{\theta};\mathbf{\phi}(s,a))\) at the initial point \(\mathbf{\theta}_{0}\). Based on the surrogate stationary point \(\mathbf{\theta}_{0}^{*}\) of \(f_{0}(\mathbf{\theta};\mathbf{\phi}(s,a))\), we correspondingly define a modified Lyapunov function \(\hat{M}(\mathbf{\theta}_{k},\lambda_{k})\) combining updates of \(\mathbf{\theta}\) and \(\lambda\) with respect to such local linearization. Specifically, we have

\[\hat{M}(\mathbf{\theta}_{k},\lambda_{k}):=\frac{\eta_{k}}{\alpha_{k}}\|\mathbf{\theta }_{k}-\mathbf{\theta}_{0}^{*}\|^{2}+\|\lambda_{k}-y_{0}(\mathbf{\theta}_{k})\|^{2}, \tag{16}\]

where \(y_{0}(\cdot)\) is in the same expression as \(y(\cdot)\) in (13) by replacing \(f(\cdot)\) with \(f_{0}(\cdot)\), and we will describe this in details below. **(Step 2)** We then study the convergence rate of the nonlinear 2TSA using this modified Lyapunov function under general conditions. **(Step 3)** Finally, since the two coupled parameters \(\mathbf{\theta}\) and \(\lambda\) in (9) are updated with respect to the true neural network function \(f(\mathbf{\theta};\mathbf{\phi}(s,a))\) in (6) in Neural-Q-Whittle, while we characterize their convergence using the approximated neural network function in Step 2. Hence, this further requires us to characterize the approximation errors. We visualize the above three steps in Figure 1 and provide a proof sketch in Section 4.3. Combing them together gives rise to our main theoretical results on the finite-time performance of Neural-Q-Whittle, which is formally stated in the following theorem.

**Theorem 1**.: _Consider iterates \(\{\mathbf{\theta}_{k}\}\) and \(\{\lambda_{k}\}\) generated by Neural-Q-Whittle in (7) and (8). Given \(\alpha_{k}=\frac{\alpha_{0}}{(k+1)},\eta_{k}=\frac{\eta_{0}}{(k+1)^{4/3}}\), we have for \(\forall k\geq\tau\)_

\[\mathbb{E}[M(\mathbf{\theta}_{k+1},\lambda_{k+1})|\mathcal{F}_{k-\tau }]\leq\frac{2\tau^{2}\mathbb{E}[\hat{M}(\mathbf{\theta}_{\tau},\lambda_{\tau})]}{ (k+1)^{2}}+\frac{1200\alpha_{0}^{3}}{\eta_{0}}\frac{(C_{1}+\|\hat{\mathbf{\theta}} _{0}\|)^{2}+(2C_{1}+\|\tilde{\lambda}_{0}\|)^{2}}{(k+1)^{2/3}}\] \[+\frac{2\eta_{0}c_{0}^{2}}{\alpha_{0}(1-\kappa)^{2}}\|span(\Pi_{ \mathcal{F}}f(\mathbf{\theta}^{*})-f(\mathbf{\theta}^{*}))\|^{2}+\bigg{(}\frac{2}{(k+ 1)^{2/3}}+2\bigg{)}\,\mathcal{O}\Big{(}\frac{c_{1}^{3}(\|\mathbf{\theta}_{0}\|+| \lambda_{0}|+1)^{3}}{m^{1/2}}\Big{)}, \tag{17}\]

_where \(C_{1}:=c_{1}(\|\mathbf{\theta}_{0}\|+\|\lambda_{0}\|+1)\) with \(c_{1}\) being a proper chosen constant, \(c_{0}\) is a constant defined in Assumption 3, \(\tau\) is the mixing time defined in (22), \(span\) denotes for the span semi-norm [47], and \(\Pi_{\mathcal{F}}\) represents the projection to the set of \(\mathcal{F}\) containing all possible \(f_{0}(\mathbf{\theta};\mathbf{\phi}(s,a))\) in (18)._

The first term on the right hand side (17) corresponds to the bias compared to the Lyapunov function at the mixing time \(\tau\), which goes to zero at a rate of \(\mathcal{O}(1/k^{2})\). The second term corresponds to the accumulated estimation error of the nonlinear 2TSA due to Markovian noise, which vanishes at the rate \(\mathcal{O}(1/k^{2/3})\). Hence it dominates the overall convergence rate in (17). The third term captures the distance between the optimal solution \((\mathbf{\theta}^{*},\lambda^{*})\) to the true neural network function \(f(\mathbf{\theta}_{k};\mathbf{\phi}(s,a))\) in
(6) and the optimal one \((\mathbf{\theta}^{*}_{0},y_{0}(\mathbf{\theta}^{*}_{0}))\) with local linearization \(f_{0}(\mathbf{\theta}_{k};\mathbf{\phi}(s,a))\) in (18), which quantifies the error when \(f(\mathbf{\theta}^{*})\) does not fall into the function class \(\mathcal{F}\). The last term characterizes the distance between \(f(\mathbf{\theta}_{k};\mathbf{\phi}(s,a))\) and \(f_{0}(\mathbf{\theta}_{k};\mathbf{\phi}(s,a))\) with any \(\mathbf{\theta}_{k}\). Both terms diminish as \(m\to\infty\). Theorem 1 implies the convergence to the optimal value \((\mathbf{\theta}^{*},\lambda^{*})\) is bounded by the approximation error, which will diminish to zero as representation power of \(f_{0}(\mathbf{\theta}_{k};\mathbf{\phi}(s,a))\) increases when \(m\to\infty\). Finally, we note that the right hand side (17) ends up in \(\mathcal{O}(1/k^{2})+\mathcal{O}(1/k^{2/3})+c\), where \(c\) is a constant and its value goes to \(0\) as \(m\to\infty\). This indicates the error bounds of linearization with the original neural network functions are controlled by the overparameterization value of \(m\). Need to mention that a constant step size will result in a non-vanishing accumulated error as in [15].

**Remark 2**.: _A finite-time analysis of nonlinear 2TSA was presented in [39]. However, [39] required a stability condition that \(\lim_{k\to\infty}(\mathbf{\theta}_{k},\lambda_{k})=(\mathbf{\theta}^{*},\lambda^{*})\), and both \(h\) and \(g\) are locally approximated as linear functions. [19; 60] relaxed these conditions and provided a finite-time analysis under i.i.d. noise. These results were later extended to Markovian noise [20] under the assumption that \(H\) function is strongly monotone in \(\mathbf{\theta}\) and \(G\) function is strongly monotone in \(\lambda\). Since [20] leveraged the techniques in [19], it needed to explicitly characterize the covariance between the error caused by Markovian noise and the parameters' residual error in (14), leading to the convergence analysis much more intrinsic. [15] exploited the mixing time to avoid the covariance between the error caused by Markovian noise and the parameters' residual error, however, it only considered the single timescale \(Q\)-learning with linear function approximation. Though our Neural-Q-Whittle can be rewritten as a nonlinear 2TSA, the nonlinear parameterization of \(Q\)-function caused by the neural network function approximation makes the aforementioned analysis not directly applicable to ours and requires additional characterization as highlighted in Figure 1. The explicit characterization of approximation errors further distinguish our work._

### Proof Sketch

In this section, we sketch the proofs of the three steps shown in Figure 1 as required for Theorem 1.

#### 4.3.1 Step 1: Approximated Solution of Neural-Q-Whittle

We first approximate the optimal solution by projecting the Q-function in (6) to some function classes parameterized by \(\mathbf{\theta}\). The common choice of the projected function classes is the local linearization of \(f(\mathbf{\theta};\mathbf{\phi}(s,a))\) at the initial point \(\mathbf{\theta}_{0}\)[13; 62], i.e., \(\mathcal{F}:=\{f_{0}(\mathbf{\theta};\mathbf{\phi}(s,a)),\forall\mathbf{\theta}\in\Theta\}\), where

\[f_{0}(\mathbf{\theta};\mathbf{\phi}(s,a))=\frac{1}{\sqrt{m}}\sum_{r=1}^{m}b_{r}\mathds {1}\{\mathbf{w}_{r,0}^{\mathsf{T}}\mathbf{\phi}(s,a)>0\}\mathbf{w}_{r}^{\mathsf{T }}\mathbf{\phi}(s,a). \tag{18}\]

Then, we define the approximate stationary point \(\mathbf{\theta}^{*}_{0}\) with respect to \(f_{0}(\mathbf{\theta};\mathbf{\phi}(s,a))\) as follows.

**Definition 1**.: _[[13; 62]] A point \(\mathbf{\theta}^{*}_{0}\in\Theta\) is said to be the approximate stationary point of Algorithm 1 if for all feasible \(\mathbf{\theta}\in\Theta\) it holds that \(\mathbb{E}_{n,\pi,\mathcal{P}}[(\Delta_{0}\cdot\nabla_{\mathbf{\theta}}f_{0}(\mathbf{ \theta};\mathbf{\phi}(s,a)))^{\mathsf{T}}(\mathbf{\theta}-\mathbf{\theta}^{*}_{0})]\geq 0, \forall\mathbf{\theta}\in\Theta,\) with \(\Delta_{0}:=[r(s,a)+(1-a)\lambda^{*}-I_{0}(\mathbf{\theta})+\max_{a^{\prime}}f_{0} (\mathbf{\theta};\mathbf{\phi}(s^{\prime},a))-f_{0}(\mathbf{\theta};\mathbf{\phi}(s,a))]\), where \(I_{0}(\mathbf{\theta})=\frac{1}{2S}\sum_{s\in\mathcal{S}}[f_{0}(\mathbf{\theta};\mathbf{ \phi}(s,0))+f_{0}(\mathbf{\theta};\mathbf{\phi}(s,1))]\)._

Figure 1: Neural-Q-Whittle operates w.r.t. true neural function \(f(\cdot)\) with its finite-time performance given in Theorem 1 (indicated in dashed lines). Our proofs operate in three steps: (i) Step 1: Obtain local linearization \(f_{0}(\cdot)\) and define Lyapunov function \(\hat{M}(\cdot)\) w.r.t. \(f_{0}(\cdot)\). (ii) Step 2: Characterize the finite-time performance w.r.t. \(\hat{M}(\cdot)\) using Lyapunov drift method. Since Neural-Q-Whittle is updated w.r.t. \(f(\cdot)\), we need to characterize the gap between \(f(\cdot)\) and \(f_{0}(\cdot)\). (iii) Step 3: Similarly, we characterize the approximation errors between \(M(\cdot)\) and \(\hat{M}(\cdot)\).
Though there is a gap between the true neural function (6) and the approximated local linearized function (18), the gap diminishes as the width of neural network i.e., \(m\), becomes large [13; 62].

With the approximated stationary point \(\boldsymbol{\theta}_{0}^{*}\), we can redefine the two error terms in (14) as

\[\hat{\boldsymbol{\theta}}_{k}=\boldsymbol{\theta}_{k}-\boldsymbol{\theta}_{0}^ {*},\qquad\hat{\lambda}_{k}=\lambda_{k}-y_{0}(\boldsymbol{\theta}_{k}), \tag{19}\]

using which we correspondingly define a modified Lyapunov function \(\hat{M}(\boldsymbol{\theta}_{k},\lambda_{k})\) in (16), where

\[y_{0}(\boldsymbol{\theta})\!=\!r(s,\!1)\!+\!\sum_{s^{\prime}}\!p(s^{\prime}|s, 1)\max_{a}\!f_{0}(\boldsymbol{\theta};\boldsymbol{\phi}(s^{\prime},a))\!-\!r (s,\!0)\!\!-\!\!\sum_{s^{\prime}}p(s^{\prime}|s,0)\max_{a}\!f_{0}(\boldsymbol {\theta};\boldsymbol{\phi}(s^{\prime},a)). \tag{20}\]

#### 4.3.2 Step 2: Convergence Rate of \(\hat{M}(\boldsymbol{\theta}_{k},\lambda_{k})\) in (16)

Since we approximate the true neural network function \(f(\boldsymbol{\theta};\boldsymbol{\phi}(s,a))\) in (6) with the local linearized function \(f_{0}(\boldsymbol{\theta};\boldsymbol{\phi}(s,a))\) in (18), the operators \(h(\cdot)\) and \(g(\cdot)\) in (10)-(11) turn correspondingly to be

\[h_{0}(X_{k},\boldsymbol{\theta}_{k},\lambda_{k})=\nabla_{\boldsymbol{\theta} }f_{0}(\boldsymbol{\theta}_{k};\boldsymbol{\phi}(S_{k},A_{k}))\Delta_{k,0},\ g_{0}( \boldsymbol{\theta}_{k}):=f_{0}(\boldsymbol{\theta}_{k};\boldsymbol{\phi}(s,1 ))-f_{0}(\boldsymbol{\theta}_{k};\boldsymbol{\phi}(s,0)), \tag{21}\]

with \(\Delta_{k,0}:=r(S_{k},A_{k})+(1-A_{k})\lambda_{k}-I_{0}(\boldsymbol{\theta}_{ k})+\max_{a}f_{0}(\boldsymbol{\theta}_{k};\boldsymbol{\phi}(S_{k+1},a))-f_{0}( \boldsymbol{\theta}_{k};\boldsymbol{\phi}(S_{k},A_{k})).\)

Before we present the finite-time error bound of the nonlinear 2TSA (9) under Markovian noise, we first discuss the mixing time of the Markov chain \(\{X_{k}\}\) and our assumptions.

**Definition 2** (Mixing time [15]).: _For any \(\delta>0\), define \(\tau_{\delta}\) as_

\[\tau_{\delta}=\!\min\{k\geq 1:\|\mathbb{E}[h_{0}(X_{k},\boldsymbol{\theta}, \lambda)|X_{0}=x]-H_{0}(\boldsymbol{\theta},\lambda)\|\leq\delta(\|\boldsymbol {\theta}-\boldsymbol{\theta}_{0}^{*}\|+\|\lambda-y_{0}(\boldsymbol{\theta}_{0 }^{*}\|\|)\}. \tag{22}\]

**Assumption 1**.: _The Markov chain \(\{X_{k}\}\) is irreducible and aperiodic. Hence, there exists a unique stationary distribution \(\mu\)[32], and constants \(C>0\) and \(\rho\in(0,1)\) such that \(d_{TV}(P(X_{k}|X_{0}=x),\mu)\leq C\rho^{k},\forall k\geq 0,x\in\mathcal{X},\) where \(d_{TV}(\cdot,\cdot)\) is the total-variation (TV) distance [32]._

**Remark 3**.: _Assumption 1 is often assumed to study the asymptotic convergence of stochastic approximation under Markovian noise [4; 9; 15]._

**Lemma 1**.: _The function \(h_{0}(X,\boldsymbol{\theta},\lambda)\) defined in (21) is globally Lipschitz continuous w.r.t \(\boldsymbol{\theta}\) and \(\lambda\) uniformly in \(X\), i.e., \(\|h_{0}(X,\boldsymbol{\theta}_{1},\lambda_{1})\!-\!h_{0}(X,\boldsymbol{ \theta}_{2},\lambda_{2})\|\leq L_{h,1}\|\boldsymbol{\theta}_{1}-\boldsymbol{ \theta}_{2}\|+L_{h,2}\|\lambda_{1}-\lambda_{2}\|,\forall X\in\mathcal{X}\), and \(L_{h,1}=3,h_{h,2}=1\) are valid Lipschitz constants._

**Lemma 2**.: _The function \(g_{0}(\boldsymbol{\theta})\) defined in (21) is linear and thus Lipschitz continuous in \(\boldsymbol{\theta}\), i.e., \(\|g_{0}(\boldsymbol{\theta}_{1})-g_{0}(\boldsymbol{\theta}_{2})\|\leq L_{g} \|\boldsymbol{\theta}_{1}-\boldsymbol{\theta}_{2}\|\), and \(L_{g}=2\) is a valid Lipschitz constant._

**Lemma 3**.: _The function \(y_{0}(\boldsymbol{\theta})\) defined in (20) is linear and thus Lipschitz continuous in \(\boldsymbol{\theta}\), i.e., \(\|y_{0}(\boldsymbol{\theta}_{1})-y_{0}(\boldsymbol{\theta}_{2})\|\leq L_{y} \|\boldsymbol{\theta}_{1}-\boldsymbol{\theta}_{2}\|\), and \(L_{y}=2\) is a valid Lipschitz constant._

**Remark 4**.: _The Lipschitz continuity of \(h_{0}\) guarantees the existence of a solution \(\boldsymbol{\theta}\) to the ODE \(\dot{\boldsymbol{\theta}}\) for a fixed \(\lambda\), while the Lipschitz continuity of \(g_{0}\) and \(y_{0}\) ensures the existence of a solution \(\lambda\) to the ODE \(\dot{\lambda}\) when \(\boldsymbol{\theta}\) is fixed. These lemmas often serve as assumptions when proving the convergence rate for both linear and nonlinear 2TSA [31; 39; 18; 24; 19; 17; 27]._

**Lemma 4**.: _For a fixed \(\lambda\), there exists a constant \(\mu_{1}>0\) such that \(h_{0}(X,\boldsymbol{\theta},\lambda)\) defined in (10) satisfies_

\[\mathbb{E}[\dot{\boldsymbol{\theta}}^{\intercal}h_{0}(X,\boldsymbol{\theta}, \lambda)]\leq-\mu_{1}\|\hat{\boldsymbol{\theta}}\|^{2}.\]

_For fixed \(\boldsymbol{\theta}\), there exists a constant \(\mu_{2}>0\) such that \(g_{0}(X,\boldsymbol{\theta},\lambda)\) defined in (11) satisfies_

\[\mathbb{E}[\dot{\lambda}g_{0}(X,\boldsymbol{\theta},\lambda)]\leq-\mu_{2}\| \hat{\lambda}\|^{2}.\]

**Remark 5**.: _Lemma 4 guarantees the stability and uniqueness of the solution \(\boldsymbol{\theta}\) to the ODE \(\dot{\boldsymbol{\theta}}\) for a fixed \(\lambda\), and the uniqueness of the solution \(\lambda\) to the ODE \(\dot{\lambda}\) for a fixed \(\boldsymbol{\theta}\). This assumption can be viewed as a relaxation of the stronger monotone property of nonlinear mappings [19; 15], since it is automatically satisfied if \(h\) and \(g\) are strong monotone as assumed in [19]._

**Lemma 5**.: _Under Assumption 1 and Lemma 1, there exist constants \(C>0\), \(\rho\in(0,1)\) and \(L=\max(3,\max_{X}h_{0}(X,\boldsymbol{\theta}_{0}^{*}),y_{0}(\boldsymbol{\theta}_{0 }^{*}))\) such that_

\[\tau_{\delta}\leq\frac{\log(1/\delta)+\log(2LCmd)}{\log(1/\rho)}.\]
**Remark 6**.: \(\tau_{\delta}\) _is equivalent to the mixing time of the underlying Markov chain satisfying \(\lim_{\delta\to 0}\delta\tau_{\delta}=0\)[15]. For simplicity, we remove the subscript and denote it as \(\tau\)._

We now present the finite-time error bound for the Lyapunov function \(\hat{M}(\boldsymbol{\theta}_{k},\lambda_{k})\) in (16).

**Theorem 2**.: _Consider iterates \(\{\boldsymbol{\theta}_{k}\}\) and \(\{\lambda_{k}\}\) generated by_ Neural-Q-Whittle _in (7) and (8). Given Lemma 1-4, \(\alpha_{k}=\frac{\alpha_{0}}{(k+1)},\eta_{k}=\frac{\eta_{0}}{(k+1)^{4/3}}\), \(\hat{C}_{1}:=c_{1}(\|\boldsymbol{\theta}_{0}\|+\|\lambda_{0}\|+1)\) with a constant \(c_{1}\),_

\[\mathbb{E}[\hat{M}(\boldsymbol{\theta}_{k+1},\lambda_{k+1})| \mathcal{F}_{k-\tau}]\leq\frac{\tau^{2}\mathbb{E}[\hat{M}(\boldsymbol{\theta }_{\tau},\lambda_{\tau})]}{(k+1)^{2}}+\frac{600\alpha_{0}^{3}}{\eta_{0}}\frac {(C_{1}+\|\hat{\boldsymbol{\theta}}_{0}\|)^{2}+(2C_{1}+\|\hat{\lambda}_{0}\|) ^{2}}{(k+1)^{2/3}}\\ +\frac{\mathcal{O}\Big{(}c_{1}^{3}(\|\boldsymbol{\theta}_{0}\| \!+\!|\lambda_{0}|\!+\!1)^{3}m^{-1/2}\Big{)}}{(k+1)^{2/3}},\quad\forall k\geq\tau. \tag{23}\]

3.3 Step 3: Approximation Error between \(M(\boldsymbol{\theta}_{k},\lambda_{k})\) and \(\hat{M}(\boldsymbol{\theta}_{k},\lambda_{k})\)

Finally, we characterize the approximation error between Lyapunov functions \(M(\boldsymbol{\theta}_{k},\lambda_{k})\) and \(\hat{M}(\boldsymbol{\theta}_{k},\lambda_{k})\). Since we are dealing with long-term average MDP, we assume that the total variation of the MDP is bounded [47].

**Assumption 2**.: _There exists \(0<\kappa<1\) such that \(\sup_{(s,a),(s^{\prime},a^{\prime})}\|p(\cdot|s,a)-p(\cdot|s^{\prime},a^{ \prime})\|_{TV}=2\kappa\)._

Hence, the Bellman operator is a span-contraction operator [47], i.e.,

\[span(\mathcal{T}f_{0}(\boldsymbol{\theta}_{0}^{*})-\mathcal{T}f(\boldsymbol{ \theta}^{*}))\leq\kappa\;span(f_{0}(\boldsymbol{\theta}_{0}^{*})-f(\boldsymbol {\theta}^{*})). \tag{24}\]

**Assumption 3**.: \(\|\boldsymbol{\theta}_{0}^{*}-\boldsymbol{\theta}^{*}\|\leq c_{0}\|span(f_{0 }(\boldsymbol{\theta}_{0}^{*})-f(\boldsymbol{\theta}^{*}))\|\), _with \(c_{0}\) being a positive constant._

**Lemma 6**.: _For \(M(\boldsymbol{\theta}_{k},\lambda_{k})\) in (15) and \(\hat{M}(\boldsymbol{\theta}_{k},\lambda_{k})\) in (16), with constants \(c_{1}\) and \(c_{0}\) (Assumption 3),_

\[M(\boldsymbol{\theta}_{k},\lambda_{k})\leq 2\hat{M}(\boldsymbol{\theta}_{k}, \lambda_{k})+\frac{2\eta_{k}c_{0}^{2}}{\alpha_{k}(1-\kappa)}\|span(\Pi_{ \mathcal{F}}f(\boldsymbol{\theta}^{*})-f(\boldsymbol{\theta}^{*}))\|+2 \mathcal{O}\Big{(}\frac{c_{1}^{3}(\|\boldsymbol{\theta}_{0}\|\!+\!|\lambda_{0} |\!+\!1)^{3}}{m^{1/2}}\Big{)}.\]

## 5 Numerical Experiments

We numerically evaluate the performance of Neural-Q-Whittle using an example of circulant dynamics [23; 3; 8]. The state space is \(\mathcal{S}=\{1,2,3,4\}\). Rewards are \(r(1,a)=-1,r(2,a)=r(3,a)=0,\) and \(r(4,a)=1\) for \(a\in\{0,1\}\). The dynamics of states are circulant and defined as

\[P^{1}=\begin{bmatrix}0.5&0.5&0&0\\ 0&0.5&0.5&0\\ 0&0&0.5&0.5\\ 0.5&0&0&0.5\end{bmatrix}\text{ and }P^{0}=\begin{bmatrix}0.5&0&0&0.5\\ 0.5&0.5&0&0\\ 0&0.5&0.5&0\\ 0&0&0.5&0.5\end{bmatrix}.\]

This indicates that the process either remains in its current state or increments if it is active (i.e., \(a=1\)), or it either remains the current state or decrements if it is passive (i.e., \(a=0\)). The exact value of Whittle indices [23] are \(\lambda(1)=-0.5,\lambda(2)=0.5,\lambda(3)=1,\) and \(\lambda(4)=-1\).

In our experiments, we set the learning rates as \(\alpha_{k}=0.5/(k+1)\) and \(\eta_{k}=0.1/(k+1)^{4/3}\). We use \(\epsilon\)-greedy for the exploration and exploitation tradeoff with \(\epsilon=0.5\). We consider a two-layer neural network with the number of neurons in the hidden layer as \(m=200\). As described in Algorithm 1, \(b_{r},\forall r\) are uniformly initialized in \(\{-1,1\}\) and \(w_{r},\forall r\) are initialized as a zero mean Gaussian distribution according to \(\mathcal{N}(\textbf{0},\textbf{I}_{d}/d)\). These results are carried out by Monte Carlo simulations with 100 independent trials.

Figure 2: Convergence of Neural-Q-Whittle.
**Convergence to true Whittle index.** First, we verify that Neural-Q-Whittle convergences to true Whittle indices, and compare to Q-Whittle, the first Whittle index based Q-learning algorithm. As illustrated in Figure 1(a), Neural-Q-Whittle guarantees the convergence to true Whittle indices and outperforms Q-Whittle [3] in the convergence speed. This is due to the fact that Neural-Q-Whittle updates the Whittle index of a specific state even when the current visited state is not that state.

Second, we further compare with other White index learning algorithms, i.e., Q-Whittle-LFA [60], WIQL [8] and QWIC [23]in Figure 3. As we observe from Figure 3, only Neural-Q-Whittle and Q-Whittle-LFA in [60] can converge to the true Whittle indices for each state, while the other two benchmarks algorithms do not guarantee the convergence of true Whittle indices. Interestingly, the learning Whittle indices converge and maintain a correct relative order of magnitude, which is still be able to be used in real world problems [60]. Moreover, we observe that Neural-Q-Whittle achieves similar convergence performance as Q-Whittle-LFA in the considered example, whereas the latter has been shown to achieve good performance in real world applications in [60]. Though this work focuses on the theoretical convergence analysis of Q-learning based whittle index under the neural network function approximation, it might be promising to implement it in real-world applications to fully leverage the strong representation ability of neural network functions, which serves as future investigation of this work.

**Convergence of the Lyapunov function defined in (15).** We also evaluate the convergence of the proposed Lyapunov function defined in (15), which is presented in Figure 1(b). It depicts \(\mathbb{E}[M(\mathbf{\theta}_{k},\lambda_{k})]\) vs. the number of iterations in logarithmic scale. For ease of presentation, we only take state \(s=4\) as an illustrative example. It is clear that \(M(\mathbf{\theta}_{k},\lambda_{k})\) converges to zero as the number of iterations increases, which is in alignment with our theoretical results in Theorem 1.

**Verification of Assumption 3.** We now verify Assumption 3 that the gap between \(\mathbf{\theta}_{0}^{*}\) and \(\mathbf{\theta}^{*}\) can be bounded by the span of \(f_{0}(\mathbf{\theta}_{0}^{*})\) and \(f(\mathbf{\theta}^{*})\) with a constant \(c_{0}\). In Figure 4, we show \(c_{0}\) as a function of the number of neurons in the hidden layer \(m\). It clearly indicates that constant \(c_{0}\) exists and decreases as the number of neurons grows larger.

## 6 Conclusion

We presented Neural-Q-Whittle, a Whittle index based Q-learning algorithm for RMAB with neural network function approximation. We proved that Neural-Q-Whittle achieves an \(\mathcal{O}(1/k^{2/3})\) convergence rate, where \(k\) is the number of iterations when data are generated from a Markov chain and Q-function is approximated by a ReLU neural network. By viewing Neural-Q-Whittle as 2TSA and leveraging the Lyapunov drift method, we removed the projection step on parameter update of Q-learning with neural network function approximation. Extending the current framework to two-timescale Q-learning (i.e., the coupled iterates between Q-function values and Whittle indices) with general deep neural network approximation is our future work.

Figure 4: Verification of Assumption 3 w.r.t the constant \(c_{0}\).

Figure 3: Convergence comparison between Neural-Q-Whittle and benchmark algorithms.
## Acknowledgements

This work was supported in part by the National Science Foundation (NSF) grant RINGS-2148309, and was supported in part by funds from OUSD R&E, NIST, and industry partners as specified in the Resilient & Intelligent NextG Systems (RINGS) program. This work was also supported in part by the U.S. Army Research Office (ARO) grant W911NF-23-1-0072, and the U.S. Department of Energy (DOE) grant DE-EE0009341. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the funding agencies.
