No-Regret Thompson Sampling for Finite-Horizon Markov Decision Processes with Gaussian Processes

Thank you very much for reviewing our paper and for your positive feedback. We are pleased to hear that you found our regret bounds to be a valuable contribution to the community and appreciate the practical recommendations. We address your comments below and hope this will help improve your evaluation of the paper.

Weaknesses

1. Unclear experimental section.
   Thank you for highlighting this. We agree that a more detailed explanation of the link between the experiments and theoretical analysis would improve clarity. Our experiments are specifically designed to validate our theoretical findings by confirming the sublinear growth of cumulative regret, consistent with our $\mathcal{\tilde{O}}(\sqrt{KH\Gamma(KH)})$ bound, and by illustrating how the choice of GP kernel influences learning efficiency in line with the theoretical dependence on $\Gamma(\cdot)$. To emphasize these connections, we will add approximate theoretical reference curves to Figures 1 and 2 in the camera-ready version to visually highlight the sublinear regret growth predicted by our theory relative to the regret growth in practice. We are also open to incorporating additional visualizations if recommended (although we are constrained by this year’s rebuttal format and cannot upload revised figures at this stage).

Questions

1. Are there examples of the algorithm in literature?
   Yes, the algorithm we study is closely related to a recent approach called HOT-GP (Bayrooti et al., 2025), which achieved strong empirical performance on standard benchmarks such as MuJoCo and VMAS. On these tasks, HOT-GP was shown to match or outperform both model-free (e.g., PPO, SAC) and model-based (e.g., MBPO, UCB-based) methods. While this prior work showed the algorithm’s practical utility in simulated robotics settings, it did not provide theoretical guarantees. Our work abstracts the core ideas behind HOT-GP into a form that enables a rigorous regret analysis, serving as the theoretical foundation for this class of algorithms.

Thanks again for your review and we would be happy to further discuss any of these points during the author-reviewer discussion period.

Thank you very much for your time spent reviewing our paper. We are glad that you find the work complete, the paper well-structured, and the proofs rigorous with some results reusable by the community. We address your comments in detail below and hope this will help improve your evaluation of the paper.

1. Significance of contribution and technical novelty.
   Thank you for the opportunity to clarify this. At a high level, both TS and RL are known to be inherently challenging to analyze. For instance, although TS is a natural algorithm dating back to Thompson (1933), its first rigorous analysis in the much simpler setting of K-armed bandits only appeared in the 2010s (Agrawal and Goyal, 2012; Kaufmann et al., 2012a). Extending such analyses to the RL setting, particularly with continuous state and action spaces and minimal structural assumptions, remains a significant open challenge.
   Our work takes a step forward in this direction by providing, to the best of our knowledge, the first sublinear regret bound for TS in continuous state-action MDPs using multi-output GP modeling of the reward and transition functions. We build on tools from the literature while also introducing key theoretical contributions to overcome challenges in this setting through a modular proof strategy. We highlight these innovations by walking through our analysis:

   • Step 1: We decompose the regret into an immediate regret term and a recursive term that captures uncertainty in value propagation through the transition model.
   • Step 2: Bounding the immediate regret term presents two unique challenges as outlined in Section 4. Namely, (i) the target function $Q^\star_h(s, a) = f_R(s, a) + V^\star_{h+1}(f_S(s, a))$ is recursive and compositional, and thus cannot be directly modeled as a GP; and (ii) the TS algorithm samples from the posterior of a proxy $Q_{h,k}$, rather than from the true posterior over $Q^\star_h$. To address these issues, we derive high-probability upper and lower confidence bounds on the proxy and target value functions (Theorem 1). These bounds allow us to relate the sampled proxy values to the true values and quantify the regret incurred. This confidence analysis, which handles recursive composition under GP priors, is novel to our setting.
   • Step 3: We extend the confidence bounds to the multi-step setting to analyze the propagation of uncertainty over the full episode.
   • Step 4: Finally, we bound the uncertainty accumulated across all steps and episodes, which requires controlling the sum of posterior standard deviations over the sequence of inputs. In the multi-output setting, a naive approach would apply the standard elliptical potential lemma independently to each output dimension, resulting in a regret bound that scales with $\sqrt{d \Gamma(T)}$ where $d$ is the number of outputs. Instead, we introduce a new multi-output elliptical potential lemma (Lemma 1) that tracks uncertainty jointly across all outputs using the full output covariance structure. This avoids summing over independent bounds for each output and leads to a strictly tighter regret bound that scales with $\sqrt{\Gamma(T)}$ and is independent of the output dimension $d$. This improvement is essential for scalability and highlights the benefit of using a multi-output GP model rather than independent GPs.

   In summary, while the algorithm we consider follows the general structure of posterior sampling, our analysis addresses a much more complex setting and introduces new analytical tools that are broadly applicable to related problems in RL and Bayesian learning (as noted by the reviewer). We therefore respectfully disagree with the reviewer’s assessment of the significance of our contribution as we believe that our work addresses an important and technically challenging problem of interest to the NeurIPS community.

2. Lack of large-scale experiments.
   We would like to highlight that, as mentioned in the introduction, one of the motivations for our work is the recent empirical success of a closely related algorithm called HOT-GP that is studied by Bayrooti et al. (2025). HOT-GP achieved strong performance on standard benchmarks such as MuJoCo and VMAS and outperformed both model-free (e.g., PPO, SAC) and model-based (e.g., MBPO, UCB-based) methods. While their work demonstrated the algorithm’s practical utility, it did not provide theoretical guarantees. Our goal is complementary: we focus on establishing theoretical bounds for this class of algorithms by deriving the first sublinear regret bounds guaranteeing convergence to the performance of the optimal policy. To support our analysis, we include controlled numerical experiments that align closely with our theoretical framework and are designed to corroborate the analytical findings. Specifically, our experiments confirm sublinear regret and illustrate how the kernel complexity empirically influences the learning efficiency relative to the theoretical dependence on $\Gamma(\cdot)$.

3. Assumptions.
   Below we clarify how Assumptions 1 and 2 relate to theory and practice.

   • Assumption 1: GP models are often used for modeling complex dynamics in RL, particularly in low-data regimes where uncertainty quantification is critical. Notably, GP-based methods such as PILCO (Deisenroth and Rasmussen, 2011) and its successors have demonstrated strong empirical performance in continuous-control robotics. GPs offer flexible representation capacity as the smoothness of functions depends on the GP kernel choice. For example, the Matérn kernel family provides a smoothness parameter $\nu$ that controls function regularity. As shown in Srinivas et al. (2010), functions sampled from Matérn kernels can uniformly approximate any continuous function on compact subsets of $\mathbb{R}^d$, making the GP prior highly expressive. Lower values of $\nu$ yield rougher functions that better capture non-smooth behavior. Thus, this assumption holds for a broad class of reward and transition functions, including those with limited smoothness, so we believe it is not overly restrictive.
   • Assumption 2: This is a mild technical condition required for Theorem 1 to ensure that value functions admit Taylor expansions with controlled error. This holds when the reward and transition functions are differentiable with bounded first and second derivatives.
   • Assumptions 1 and 2 hold in many continuous-control tasks common in RL, such as the Gymnasium Pendulum and MountainCarContinuous tasks. These environments have smooth, deterministic transition dynamics and reward functions, making GP modeling appropriate and ensuring the value functions are sufficiently smooth.

   While these assumptions do not capture all real-world dynamics, our setting still reflects many realistic control problems and provides an important intermediate step toward more general regret bounds. Moreover, our assumptions are relatively relaxed compared to standard assumptions in prior RL theory. For instance, the regret analysis for posterior sampling for RL (Osband et al., 2013; Osband et al., 2017) assumes a finite state-action space and does not extend to continuous MDPs, while RKHS-based regret analyses (e.g., Chowdhury & Gopalan, 2019) rely on predefined feature maps and a fixed kernel class, and linear MDP approaches (e.g., Jin et al., 2020) assume linearity in both rewards and transitions. In contrast, we give a complete analysis under one of the more general settings considered in the literature.

Thanks again for your review and we would be happy to further discuss any of these points during the author-reviewer discussion period.

Thank you very much for reviewing our paper and for your positive and constructive feedback. Below we respond to your comments and questions in detail and we hope this will help improve your evaluation of the work.

Questions

1. How was Assumption 1 used in the proof?
   Assumption 1 is used in the proof to derive the confidence bounds that serve as a key building block in the regret analysis and are specifically used in Step 2 of the proof of Theorem 2 to bound the immediate regret. To elaborate more on this, our analysis is decomposed into four steps and Step 2 bounds the immediate regret, i.e., the loss in the optimal Q-function $Q^{\star}(s_h, a^{\star}_h) - Q^{\star}(s_h, a_h)$, where $a_h$ is the action selected by the Thompson sampling algorithm at step $h$. This step relies on Corollary 1, which in turn is based on Theorem 1. Briefly, Theorem 1 establishes high-probability confidence intervals for composed GP functions and Corollary 1 applies this to recursive value functions in RL. These confidence intervals crucially rely on Assumption 1, which models the joint reward and transition functions as a multi-output GP. To clarify this connection explicitly, we have added this sentence before introducing Theorem 1: “Under Assumption 1, the true reward and transition functions are modeled with a GP prior, allowing us to derive high-probability bounds for compositional functions of the form $v(f)$.” We hope this improves the clarity and would like to note that the full derivation can be found in Appendices A and B and is also outlined in the proof sketch in Section 4.

2. Is it possible to generalize the proof for the case of stochastic rewards and transitions?
   Thank you for this insightful question. Generalizing to stochastic environments is indeed a natural next step and we view our work as a foundation for this direction. In the case of stochastic rewards, the extension is straightforward and does not affect our analysis. The confidence intervals remain valid as long as the GP noise parameter $\lambda^2$ is chosen to upper-bound the variance of the reward noise. This adjustment ensures the posterior remains well-calibrated and the regret analysis proceeds as in the deterministic case. In contrast, generalizing to stochastic transitions introduces new challenges as this setting corresponds to GP regression with input uncertainty, since the state $s$ that is used recursively in the value function would depend on a random variable. Addressing this would require nontrivial extensions beyond the standard GP framework used in our work. Developing such tools is outside the scope of this paper but we believe the techniques that we introduce provide a strong foundation for extending regret analysis to stochastic transition dynamics in future work.

3. Comparing with kernel RL (Chowdhury and Gopalan, 2019; Yang and Wang, 2020; and Vakili and Olkhovskaya, 2023) and their dependence on $H$. All of these works assume that the conditional probability distribution of the next state belongs to the RKHS of a known kernel: $p_h(\cdot \mid s,a) \in \mathcal{H}_k$ for all $(s,a)$ and each step $h$. In contrast, our assumption is that the function $f(\cdot \mid s,a)$, which maps state-action pairs to next states, can be modeled using a GP and this model remains fixed across all steps $h$ within an episode.

   There are also key differences in algorithm design: prior works adopt UCB-based approaches while we analyze Thompson Sampling.

   Our regret bounds are not only tighter by a factor of $H$, but also by a factor of $\sqrt{\Gamma(T)}$ compared to Chowdhury and Gopalan (2019) and Yang and Wang (2020). This is particularly important, as their $O(\Gamma(T)\sqrt{T})$ bounds can become vacuous (i.e., linear in $T$) for many practical kernels such as the Matérn family. A key reason is the confidence intervals used in our analysis, which are tighter due to the Bayesian modeling of transitions via $f$ rather than the non-Bayesian RKHS assumption on $p_h$.

   Roughly speaking, in a non-formal comparison of the regret analyses, the improved confidence intervals contribute a $\sqrt{H}$ factor. In addition, while those works model $p_h(\cdot \mid s,a)$ separately at each step $h$ (i.e., non-homogeneous transitions), our model assumes that $f(\cdot \mid s,a)$ remains fixed across steps within an episode. This shared structure contributes another $\sqrt{H}$ improvement. We will highlight these distinctions in the related work section for better clarity.

   The gap in confidence intervals under RKHS assumptions is also evident in the Bayesian optimization literature, and is formally stated in the following COLT open problem: Open Problem: Tight Online Confidence Intervals for RKHS Elements, Sattar Vakili, Jonathan Scarlett, Tara Javidi, COLT 2021.

Discussing Fan et al. (2018).
Thank you for pointing out the related work by Fan et al. (2018). We will add the reference and include a comparison in the paper. Their analysis builds on a technique from Osband & Van Roy (2017), which was introduced as a general method for establishing Bayesian regret, that is, the expected regret averaged over the distribution of all problem instances. In contrast, our notion of regret is stronger: we establish a high-probability bound that holds uniformly over all instances of the problem. As a result, our analysis requires fundamentally different and novel proof techniques.

Lack of large scale experiments.
We would like to highlight that, as mentioned in the introduction, one of the motivations for our work is the recent empirical success of a closely related algorithm called HOT-GP that is studied by Bayrooti et al. (2025). HOT-GP achieved strong performance on standard benchmarks such as MuJoCo and VMAS and outperformed both model-free (e.g., PPO, SAC) and model-based (e.g., MBPO, UCB-based) methods. While their work demonstrated the algorithm’s practical utility, it did not provide theoretical guarantees. Our goal is complementary: we focus on establishing theoretical bounds for this class of algorithms by deriving the first sublinear regret bounds guaranteeing convergence to the performance of the optimal policy. To support our analysis, we include controlled numerical experiments that align closely with our theoretical framework and are designed to corroborate the analytical findings. Specifically, our experiments confirm sublinear regret and illustrate how the kernel complexity empirically influences the learning efficiency relative to the theoretical dependence on $\Gamma(\cdot)$.

Thanks again for your review and we would be happy to further discuss any of these points during the author-reviewer discussion period.

Thank you very much for taking the time to review our paper. We appreciate your comments and are glad that you find the paper to be clear and well-structured and the theoretical analysis to be rigorous. We address the weaknesses highlighted in detail below and hope our responses will help improve your evaluation of the paper.

Weaknesses

1. Relying on a GP to model the environment dynamics restricts the application to a narrow class of smooth, well-behaved true functions that may not be practical.
   We agree that GP modeling introduces an assumption about the true underlying environment that does not always hold (as explicitly stated in the Limitations section), however we believe it is not overly restrictive, especially when compared to existing regret analyses. To clarify this, we respond in several parts:

   • GP models are often used for modeling complex dynamics in RL, particularly in low-data regimes where uncertainty quantification is critical. Notably, GP-based methods such as PILCO (Deisenroth and Rasmussen, 2011) and its successors have demonstrated strong empirical performance in continuous-control robotics.
   • GPs offer flexible representation capacity as the smoothness of functions depends on the GP kernel choice. For example, the Matérn kernel family provides a smoothness parameter $\nu$ that controls function regularity. As shown in Srinivas et al. (2010), functions sampled from Matérn kernels can uniformly approximate any continuous function on compact subsets of $\mathbb{R}^d$, making the GP prior highly expressive. Lower values of $\nu$ yield rougher functions that better capture non-smooth behavior. Thus, our setting can model a broad class of reward and transition functions, including those with limited smoothness.
   • Our assumptions are milder than many related works on regret bounds for RL. For instance, the regret analysis for posterior sampling for RL (Osband et al., 2013; Osband et al., 2017) assumes a finite state-action space and does not extend to continuous MDPs, while RKHS-based regret analyses (e.g., Chowdhury & Gopalan, 2019) rely on predefined feature maps and a fixed kernel class, and linear MDP approaches (e.g., Jin et al., 2020) assume linearity in both rewards and transitions. In contrast, our multi-output GP model flexibly models correlations between reward and transitions without relying on discretization, linearity, or fixed feature representations.
   • Our assumption that the true joint reward and transition function is drawn from a multi-output GP enables a tractable regret analysis in continuous state-action spaces. While GPs may not capture all real-world dynamics, such as sharp discontinuities, this setting still reflects many realistic control problems (e.g., MuJoCo environments) and provides an important intermediate step toward more general regret bounds.

2. Empirical evaluation is insufficient to demonstrate the practical value of the algorithm.
   We would like to highlight that, as mentioned in the introduction, one of the motivations for our work is the recent empirical success of a closely related algorithm called HOT-GP that is studied by Bayrooti et al. (2025). HOT-GP achieved strong performance on standard benchmarks such as MuJoCo and VMAS and outperformed both model-free (e.g., PPO, SAC) and model-based (e.g., MBPO, UCB-based) methods. While their work demonstrated the algorithm’s practical utility, it did not provide theoretical guarantees. Our goal is complementary: we focus on establishing theoretical bounds for this class of algorithms by deriving the first sublinear regret bounds guaranteeing convergence to the performance of the optimal policy. To support our analysis, we include controlled numerical experiments that align closely with our theoretical framework and are designed to corroborate the analytical findings. Specifically, our experiments confirm sublinear regret and illustrate how the kernel complexity empirically influences the learning efficiency relative to the theoretical dependence on $\Gamma(\cdot)$.

Thanks again for your review and we would be happy to further discuss any of these points during the author-reviewer discussion period.