Achieving $\tilde{\mathcal{O}}(1/N)$ Optimality Gap in Restless Bandits through Gaussian Approximation

We greatly appreciate the positive comments from the reviewer. We next address the computation issue of SP versus the backward induction suggested by the reviewer.

First, we would like to clarify that our approach applies to a general number of states, time horizon, and number of arms. The stochastic program (SP) (15)-(17) is for this general setting, and a variety of existing algorithms can be used to solve this SP. In our numerical experiments (Section 5), we used the EDDP algorithm Lan (2022) to solve the SP, and our experiment was for a setting with $K=10$ states per arm and a time horizon of $H=5$. In Figure 5 (Appendix G.3 in the supplementary material), we evaluated the time complexity of EDDP with RMAB size as large as $K=5, H=30$ and $K=30, H=5$.

Next, we would like to discuss in detail the computational complexity. Indeed, directly solving the backward induction incurs a complexity of at least $\Theta(N^K)$, which is exponential in the number of states $K$. In contrast, the SP we formulated allows for more efficient algorithms. At an intuitive level, although the SP doesn't exactly decompose the problem in the same way as the LP does, it still takes advantage of state aggregation. More concretely, our Gaussian approximation aggregates the randomness in state transitions into a Gaussian noise ($\mathbf{Z}_h$). One of our key ideas is to make the noises ($\mathbf{Z}_h$'s) independent of actions. This allows the SP to be solved through the Sample Average Approximation (SAA) technique, which uses random samples of the $\mathbf{Z}_h$'s and reduces the problem into a deterministic linear program. The number of random samples required is polynomial in $N$ and $K$. One straightforward approach then is to solve this SAA-induced linear program using the interior point method (IPM) (see, e.g., recent results on the complexity of IPM in Cohen and Song (2021)), which can already lead to a computation complexity that is polynomial in both $N$ and $K$.

That said, directly solving the resulting SAA-induced linear program can lead to an exponential dependence on the time horizon. In our paper, we chose to use the recently developed EDDP algorithm to solve the SAA. In our experiments, EDDP solves the SP efficiently, and in these numerical examples, the computational complexity is linear in both the time horizon and the number of states $K$ (Figure 5 in Appendix G.2 in the supplementary material).

References

• Cohen, M. B., Lee, Y. T., and Song, Z. (2021). "Solving linear programs in the current matrix multiplication time." Journal of the ACM (JACM), 68(1), 1–39.

• Lan, G. (2022). "Complexity of stochastic dual dynamic programming." Mathematical Programming, 191(2), 717–754.

We are delighted to know that we addressed your question. Thanks again for your great questions, and thanks for revising the rating!

We hope our responses to the reviewer's questions below provide clearer intuition and deeper insight into the policy class construction and related issues. We also appreciate the suggestion to enrich the simulation section with additional scenarios. To address this suggestion, we have conducted additional tests on several machine maintenance instances, keeping the problem size and structure the same while varying model parameters for rewards and transition probabilities. All these instances exhibit uniqueness in their LP solutions. The observed total reward gaps in these experiments, displayed in the table below, remain within the same order of magnitude as the example presented in the left subfigure of Figure 3 (Section 5 of the paper) for $N = 100$. We plan to incorporate these additional results, along with the corresponding model parameters, in a future revision of our paper.

Table: Total reward gap across additional machine maintenance instances

• What is the motivation behind considering a Gaussian stochastic system as an approximation of the original system? What other kinds of distribution can be chosen to offer a good approximation? Does a Gaussian stochastic system offer the best approximation?

We consider the Gaussian distribution because by the central limit theorem, the Gaussian distribution is the limiting distribution of the normalized mean of $N$ independent random variables. The central limit theorem also implies that other kinds of distributions can't be such a limiting distribution, so the Gaussian approximation offers the best approximation asymptotically. In particular, the Wasserstein distance bound (Lemma D.1) may not hold under other distributions.

• Why is there a need of a $1/\sqrt{N}$ factor multiplying $\mathbf{Z}_h$ in equation (12)?

The factor $1/\sqrt{N}$ again arises from the central limit theorem (CLT). To see this more clearly, consider $N$ independent and identically distributed random variables $X_1,X_2,\dots,X_N$, each with mean $\mu$ and variance $\sigma^2$. By the CLT, as $N$ becomes large, we have the following approximation

$

     \frac{\sum_{i=1}^N X_i-N\mu}{\sigma\sqrt{N}}\approx Z,
 

$

 where $Z$ is a standard Gaussian random variable.
 Rewriting this gives
 

$

   \frac{1}{N}\sum_{i=1}^{N}X_i \approx \mu + \frac{\sigma}{\sqrt{N}}Z,
 

$

 where there is a $1/\sqrt{N}$ factor multiplying $Z$.
 The $1/\sqrt{N}$ factor in our equation (12) comes from the same reason, although for a more complicated setting.

• The meaning of the policy class described in (13) is not clear. What is the physical interpretation behind this class construction? I wonder whether it would include a large number of feasible policies. The main result in the paper (Theorem 4.1) heavily depends on this class of policies. However, the significance of the result can be better understood if we get more insight about the class and its cardinality.

Roughly speaking, the policies in (13) are policies such that when the system state $\mathbf{x}_h$ is in an $\tilde{\mathcal{O}}(1/\sqrt{N})$ neighborhood of the optimal fluid state $\mathbf{x}^{\star}_h$, the action taken is also within an $\tilde{\mathcal{O}}(1/\sqrt{N})$ neighborhood of the optimal fluid action $\mathbf{y}^{\star}_h$. We would like to emphasize that this policy class includes commonly studied fluid-based policies in the literature, e.g., the policies in references Zhang and Frazier (2021), Gast et al. (2023), and Brown and Zhang (2023); and most existing policies for RMAB are fluid-based. Moreover, it also includes the optimal policy for the $N$-armed problem when the uniqueness assumption (Assumption 4.1) is satisfied.

• Please explain what is meant by the round operator in Algorithm 1.

The rounding operation guarantees that the action generated by the algorithm pulls an integer number of arms. It is needed due to the following reason: when computing the first- and second-order approximations, we extend the state and action spaces from their original discrete sets (with granularity $1/N$) to continuous-valued probability simplices. The action computed in this continuous domain must subsequently be converted back to the discrete domain, ensuring that all values multiplied by $N$ are integers and thus applicable to the $N$-armed problem. This rounding procedure is explicitly detailed in Appendix C of the supplementary material.

• Theorem 4.1 holds with rounding. Is it possible to extend this theorem to the case without rounding too?

Without rounding, the policy may output an action that pulls a fraction of an arm, which is typically not allowed in a restless bandit problem. That said, if pulling a fraction of an arm is allowed (say 0.5 arm), e.g. by using a probabilistic method (say pull one arm with probability 0.5), then Theorem 4.1 holds without rounding.

• Theorem 4.3 assumes the existence of a positive constant $\varepsilon$ independent of $N$. How can we guarantee the existence of such an $\varepsilon$?

We expect that such an $\varepsilon$ exists when the fluid optimal action $\mathbf{y}_1^{\star}$ at time step $1$ is not an optimal action for the Gaussian stochastic system. To see this, consider the $N$-independent stochastic program (107), which is derived from problem (14) through re-centering and scaling, as detailed in Appendix G.2 of the supplementary material. If the optimal first-stage solution of (107) differs from zero, which translates to $\mathbf{y}_1^{\star}$ not being an optimal solution, then $\varepsilon$ corresponds to the gap between the optimal objective value of (107) and the objective value (or the $Q$-value) obtained by evaluating the fluid-optimal first-stage action. Equivalently, the value of $\varepsilon$ is the difference between the largest Q-value and the Q-value under the fluid action. If the fluid action is sub-optimal, then the gap $\varepsilon$ is strictly positive and independent of $N$. Otherwise, the fluid optimal action $\mathbf{y}_1^{\star}$ is optimal for the Gaussian stochastic system, and our second-order policy will take the fluid optimal action.

References

• Brown, D. B., and Zhang, J. (2023). "Fluid policies, reoptimization, and performance guarantees in dynamic resource allocation." Operations Research.

• Gast, N., Gaujal, B., and Yan, C. (2023). "Linear program-based policies for restless bandits: Necessary and sufficient conditions for (exponentially fast) asymptotic optimality." Mathematics of Operations Research.

• Zhang, X., and Frazier, P. I. (2021). "Restless bandits with many arms: Beating the central limit theorem." arXiv preprint arXiv:2107.11911.

We greatly appreciate the reviewer's positive comments. We next address the weaknesses and questions raised by the reviewer.

• Assumes homogeneous arms (i.e., identical transition probabilities and rewards across arms).

Our results can be directly applied to settings with limited heterogeneity. A more detailed response can be found in the answer to the reviewer's question on extending to heterogeneous cases.

• Lacking comments on scalability of the approach with respect to the state-space size and horizon length.

We provide a detailed discussion on scalability later when answering the reviewer's question on scalability, highlighting computational complexity and empirical runtime performance.

• Additional motivation for the degenerate setting would improve significance—are there convincing real-world examples?

Another real-world example is the applicant screening problem studied in Brown and Smith (2020), which explicitly exhibits degeneracy by requiring non-trivial tie-breaking. We will add this example to the revision.

• Since the policy is based on the fluid approximation, it would be helpful to give some idea of how much the SG policy can deviate from the base fluid approximation and still retain optimality or feasibility of computation time – what is the factor at play there?

We thank the reviewer for this insightful question. Our second-order Gaussian policy specifically seeks improvements by considering corrections within a neighborhood of size $\tilde{\mathcal{O}}(1/\sqrt{N})$ around the fluid-optimal state and action. When the LP has a unique solution, we proved in the paper that the optimal policy of the $N$-arm problem is within $\tilde{\mathcal{O}}(1/\sqrt{N})$ neighborhood of the LP (fluid) solution (Lemma D.1 in the supplementary material). Therefore, limiting the deviation within this $\tilde{\mathcal{O}}(1/\sqrt{N})$ neighborhood retains the optimal policy and also guarantees that our second-order approximation has $\tilde{\mathcal{O}}(1/N)$ approximation error (in terms of Wasserstein distance). Expanding the neighborhood would increase the second-order approximation error; while shrinking the neighborhood would exclude the optimal policy. So it turns out that $\tilde{\mathcal{O}}(1/\sqrt{N})$ is the ``exact'' deviation needed for our proposed method. We will add this discussion to the revision.

• Can the gaussian approximation technique be extended to infinite horizon settings?

The algorithm can be directly extended to the infinite horizon settings by considering a second-order approximation around the LP solution to the infinite horizon problem. However, the $\tilde{\mathcal{O}}(1/N)$ optimality gap would require a different proof due to the fundamental difference between finite-horizon MDPs and infinite-horizon average-reward MDPs.

• Can you comment on how you might relax the homogeneity assumption and still maintain the gap improvement?

The heterogeneous setting is more challenging because the arms can no longer be aggregated into a single representative arm as in the homogeneous case. However, if the arms can be divided into different types, and the number of types is independent of $N$, then our result can directly apply. In fact, our machine maintenance problem in Section 5 is such a heterogeneous setting and considers two types of machines (see details in Appendix G.3 of the supplementary materials). We will clarify this in the revision. The ``fully'' heterogeneous setting would likely require new techniques.

• Comment on computational complexity of the approach – how it scales with $|S|$ and $H$?

The computational complexity of our approach stems primarily on solving the Gaussian SP in Line 6 of Algorithm 1. The numerical solution methods are detailed in Appendix G.2 in the supplemented material: After recentering and scaling, we obtain an $N$-independent multi-stage linear stochastic program (problem (107) in the supplementary material). The computational complexity of solving problem (107) depends on the number of states per arm, $|S|$, and the horizon length, $H$. Our numerical evaluations using Explorative Dual Dynamic Programming (EDDP, Lan (2022)) indicate that the runtime scales linearly in both $|S|$ and $H$, as illustrated in Figure 5 in the supplementary material.

• Can this approach be extended to the multi-action/WCMDP setting?

The key is to come up with an equivalent version of the action mapping algorithm in Section D.2.3 (supplementary materials) and then prove it is locally Lipschitz. The proof in this paper utilizes some unique properties of the RMAB. If a similar result can be proven for the more general WCMDP setting, then the result should generalize.

$**References:**$

• Brown, D. B., and Smith, J. E. (2020). "Index policies and performance bounds for dynamic selection problems." Management Science, 66, 3029–3050.

• Lan, G. (2022). "Complexity of stochastic dual dynamic programming." Mathematical Programming, 191(2), 717–754.

We greatly appreciate the reviewer's willingness to reevaluate the score based on our response. We next address the connection to diffusion approximation and whether the $\mathcal{O}(1/N)$ gap can be further reduced.

Connection to diffusion approximation:

In the current version of the paper, we only briefly mention the connection to diffusion approximation in Appendix B of the supplementary material. We agree that this connection warrants further elaboration, and we plan to expand this discussion following the outline below.

First, there is a rich literature on diffusion approximations for continuous-time queueing systems, particularly in the heavy-traffic regime (see, e.g., seminal work [3], [4], [12], [13] in the reference list of the supplementary material). This body of work typically analyzes queueing systems under a fixed policy, and the focus is on characterizing the convergence or the rate of convergence of the centered and scaled queueing system to the diffusion limit. In contrast, our work addresses a different goal for a different application: we focus on finding a near-optimal policy for a discrete-time restless bandit problem.

A more closely related line of work, which builds on the diffusion approximation literature, is on approximate diffusion control for queueing systems or stochastic control in general (see, e.g., seminal work and textbook such as Beneš et al. (1980), Harrison (1985), and Harrison and Wein (1989)). In this line of work, the problem is approximated by a diffusion control problem, which is then solved via the associated Hamilton-Jacobi-Bellman (HJB) equation. Recent work by Braverman et al (2020) shows that the HJB equation naturally arises under the second-order approximation of value functions of general Markov decision processes (MDPs).

While we can apply the same diffusion control to our problem, solving the HJB equation is known to be notoriously difficult in general. This is because the diffusion coefficients often depend on the control, which makes the HJB a highly nonlinear PDE. This actually motivated us to consider a different second-order approximation.

A key innovation in our work is to replace (state, action)-dependent diffusion terms with (state, action)-independent diffusion terms. We construct our Gaussian approximation system such that the covariances of the Gaussian noises (analogous to the diffusion coefficients) are independent of the actions, by leveraging the optimal fluid solution. This proposed Gaussian approximation allows us to tap into stochastic programming techniques, where (state,action)-independent randomness is required for efficient algorithms like the EDDP algorithm we used in this paper. In other words, our Gaussian approximation can be viewed as a further ``approximation/simplification'' of the diffusion approximation, but it has the same order of approximation error in our problem and can be solved efficiently using existing stochastic programming methods.

We note that our results are possible because in the restless bandit problem, the noise terms are on a smaller scale compared to the deterministic terms as the number of arms scales, due to the central limit theorem. This is reflected in the $1/\sqrt{N}$ factor for $\mathbf{Z}_h$ in the SP in equations (14)-(17). This scale separation makes it possible to use the (state,action)-independent covariances generated by the fluid optimal solution as a surrogate for the action-dependent covariances, and leads to a critical result (Lemma D.1) that shows that the optimal policy of the $N$-arm problem is within $\tilde{\mathcal{O}}(1/\sqrt{N})$ neighborhood of the fluid solution under a uniqueness assumption. We do not usually have such a property in traditional diffusion control.

In summary, our approach may be viewed as a further approximation of the traditional ``diffusion approximation''. Based on the proximity result (Lemma D.1) for restless bandits, our approximation reduces the problem to simpler stochastic programming (instead of solving nonlinear PDEs) and guarantees $\tilde{\mathcal{O}}(1/N)$ optimality gap.

We hope these expanded discussion addresses the reviewer's concern, and we will incorporate a detailed version of this discussion into our paper.

The $\mathcal{O}(1/N)$ gap.

We thank the reviewer for this question. Our conjecture is that any asymptotic methods that reduce the problem to continuous state and action spaces will have an $\Omega(1/N)$ gap in a general setting. This is because rounding procedure (like the one detailed in Appendix C), which is necessary to map a continuous action to a discrete action in the N-arm system, will result in an $\mathcal{O}(1/N)$ error in general, except for some special cases, such as when the optimal policy is a priority policy.

References:

• Beneš, V. E., Shepp, L. A., and Witsenhausen, H. S. (1980). "Some solvable stochastic control problems." Stochastics, 4(1), 39–83.

• Braverman, A., Gurvich, I., and Huang, J. (2020). "On the Taylor expansion of value functions." Operations Research, 68(2), 631–654.

• Harrison, J. M. (1985). "Brownian Motion and Stochastic Flow Systems." Wiley, New York.

• Harrison, J. M., and Wein, L. M. (1989). "Scheduling networks of queues: Heavy traffic analysis of a simple open network." Queueing Systems, 5(4), 265–279.