Gaussian Process Thompson Sampling

• Gaussian Process Thompson Sampling is a Bayesian optimization method that uses a GP surrogate and random posterior sampling to efficiently explore global optima.
• It leverages advanced strategies like rootfinding and local optimizers to overcome challenges in high-dimensional, multimodal spaces.
• The method offers strong theoretical regret bounds and scalable implementations, including sparse and distributed variations, for practical applications.

Gaussian Process Thompson Sampling (GP-TS) is a Bayesian decision-making policy for global optimization, sequential design, and bandit problems in continuous or structured action spaces. It combines Gaussian process (GP) regression—as a nonparametric Bayesian surrogate for an unknown, expensive-to-evaluate objective function—with the Thompson sampling paradigm, which stochastically explores the decision space by sampling functions from the posterior and acting greedily with respect to these samples. GP-TS is of significant interest due to its strong theoretical guarantees, conceptual simplicity, and empirical performance, especially in high-dimensional and multi-modal settings.

1. The GP-TS Framework: Prior, Posterior, and Decision Rule

The GP-TS framework begins by positing a GP prior on the unknown objective:

$f(x) \sim \mathcal{GP}(m(x), k(x,x'))$

where $m(x)$ is the mean function (typically $0$) and $k(x,x')$ is a positive-definite kernel. After each sequential observation $y^i = f(x^i) + \epsilon_i$ ($\epsilon_i \sim \mathcal{N}(0, \sigma_n^2)$), the posterior remains a GP with mean and covariance functions: $$m_{t-1}(x) = k(x,X) [K + \sigma_n^2 I]^{-1} y, \quad k_{t-1}(x,x') = k(x,x') - k(x,X)[K + \sigma_n^2 I]^{-1} k(X,x')$$ At each iteration $t$, a sample function $f_t(\cdot)$ is drawn from the GP posterior. The next evaluation point $x_t$ is selected as: $x_t = \arg\max_{x \in \mathcal{X}} f_t(x)$ This maximization is critical: since $f_t$ is a typical sample from the current posterior, globally maximizing $f_t$ translates the stochastic exploration of Thompson sampling into a functionally rich, spatially aware action selection strategy (Adebiyi et al., 2024).

2. Challenges and Strategies for Global Maximization

Posterior samples from a GP are generally highly non-convex and exhibit many local extrema; in high dimensions, $f_t$ can have a complex modal structure with strong spatial correlations. Theoretical regret guarantees for GP-TS rely on the ability to perform exact (or sufficiently accurate) global maximization of each $f_t$. In practice, this step is arduous and usually dominates the computational cost.

Recent advances have addressed this challenge by exploiting the structure of GP priors, particularly for separable kernels. By representing the sampled function as a decoupled sum $f_t(x) = f_0(x) + b(x)$, where $f_0$ is a separable prior sample and $b(x)$ is a smooth data-dependent "correction" (via Matheron's rule), all critical points of $f_0$ can be cheaply identified via univariate rootfinding in each input dimension. The Cartesian product of these critical points forms an exploration set $S_e$ that covers the local optima of the prior sample. One then initializes a local optimizer (e.g., L-BFGS) from each point in $S_e$ as well as from data points (exploitation set $S_p$), refining toward the global maximizer of the perturbed function $f_t$ (Adebiyi et al., 2024).

The core algorithm is as follows:

1. Fit GP posterior $(m_{t-1}, k_{t-1})$.
2. Draw $f_0$ by pathwise or Mercer/Bochner expansion.
3. Solve $f_i'(x_i) = 0$ in every coordinate to obtain all roots per dimension.
4. Form $S_e$ from Cartesian products corresponding to prior maxima/minima.
5. Union $S = S_e \cup S_p$.
6. From each $s \in S$, run a gradient ascent on $f_t(x)$ to a local maximizer.
7. Select $x_t = \arg\max_s f_t(x_s)$.
8. Sample $y_t = f(x_t) + \epsilon$, update data.

Benchmarks show this method accelerates and improves global maximization compared to random multi-start or evolutionary optimizers, especially as dimension grows (Adebiyi et al., 2024).

3. Theoretical Regret, Concentration, and Rates

Regret bounds for GP-TS are typically of the form:

$$\mathrm{BayesRegret}(T) = \tilde O(\sqrt{T \gamma_T})$$

where $\gamma_T$ is the maximum information gain of the GP kernel. Recent frequentist analyses have refined this picture by using fractional (tempered) posteriors (variance inflation by $\alpha^{-1}$, $\alpha \in (0,1)$), yielding kernel-agnostic regret bounds in terms of $\gamma_T$ and a posterior contraction rate $\epsilon_t$: $$\mathbb{E}[\text{Regret}(T)] = \tilde O\left(\sqrt{ \gamma_T \cdot \sum_{t=1}^T t \,\epsilon_t^2 } \right)$$ For squared-exponential kernels, the rate is $O(\sqrt{T}\,\mathrm{polylog}\,T)$; for Matérn-$\nu$ kernels, $O(T^{(2\nu+3d)/(2(2\nu+d))})$; both match minimax lower bounds up to polylogarithmic factors (Roy et al., 16 Feb 2026). Strict exponential concentration of the iterates $x^t$ around the global maximizer $x^*$ can also be established under regularity and sharpness conditions: $$\Pr(\|x^t - x^*\| > \epsilon) \leq C t^{d/2} \delta_\epsilon^{-d} \exp(-c \delta_\epsilon^2 t)$$ Here, $\delta_\epsilon$ captures the steepness of the function near the optimum; adaptive rates follow from its magnitude (Basu et al., 2017).

4. Extensions: Adaptive Priors, Sparse GPs, and Distributed GP-TS

Practical deployment often requires adaptation to unknown or misspecified priors. Recent work introduces prior selection strategies for GP-TS: Prior-Elimination GP-TS (PE-GP-TS) eliminates bad priors using predictive error tests, while HyperPrior GP-TS (HP-GP-TS) maintains a hyperposterior distribution over the prior set, sampling from both prior and function posteriors. These methods achieve cumulative regret matching UCB-style elimination in $O(\sqrt{T \beta_T \hat \gamma_T})$ (Sandberg et al., 3 Feb 2025).

GP-TS’s cubic scaling in traditional GP inference and sampling can be circumvented by using sparse approximations such as SVGP or finite-feature decompositions. Scalable TS (STS) leverages inducing points and finite spectral expansions to drive per-iteration costs to $O(m^2 t + (M+m)BN)$, maintaining regret guarantees equivalent to the standard GP-TS for appropriate spectral decay and error budgets (Vakili et al., 2020).

In distributed and decentralized optimization, multi-agent GP-TS over a communication network exploits local GPs and partial data fusion. Regret bounds accelerate with the graph's clique number (fully-connected yields a $\sqrt{M}$ speed-up), and performance is sensitive to connectivity (Zerefa et al., 2024).

5. Algorithmic Variants and Alternative Inference Techniques

Alternative techniques for implementing GP-TS arise when exact global maximization is infeasible. Sequential Monte Carlo (SMC) methods approximate the distribution of the maximizer directly, using particles and importance weighting rather than challenging inner-loop optimization of function samples. The resulting SMC-TS algorithm samples maximizer candidates, challenges them via posterior draws, and uses systematic resampling to focus computation. This paradigm avoids acquisition function optimization entirely, yielding competitive regret with low overhead, especially in moderate dimensions (Bijl et al., 2016).

Further, adaptive grid-based TS algorithms efficiently target parameter–trajectory pairs in high-noise or stochastic simulation settings by coupling GP surrogates over $(x, s)$ (input, seed) and refining candidate sets using surrogate-informed filtering and MCMC densification. These techniques substantively accelerate trajectory discovery versus summary-statistics-based BO methods (Fadikar et al., 20 Oct 2025).

6. Applications and Empirical Insights

Empirical evaluations across high-dimensional deterministic and stochastic test functions demonstrate that appropriately initialized and globally optimized GP-TS (notably, rootfinding-based variants) achieve or surpass the performance of widely used competitors such as GP-UCB and Expected Improvement (EI). On problems like 2D Schwefel and 10D Lévy, rootfinding GP-TS finds near-global optima rapidly, reduces simple regret faster, and balances exploration/exploitation more effectively than random Fourier features or genetic algorithm (GA) starts. In adaptive trajectory-based simulation contexts, grid-based GP-TS with common random numbers identifies high-likelihood single-realization trajectories more efficiently (Adebiyi et al., 2024, Fadikar et al., 20 Oct 2025).

7. Open Problems and Ongoing Directions

Despite impressive advances, practical challenges persist: (i) ensuring truly global maximization for high-dimensional non-separable GPs; (ii) balancing computational resources between surrogate modeling and sample optimization; (iii) advancing theoretical regret guarantees for approximate or sparse GP-TS; and (iv) understanding optimal tradeoffs between exploration and exploitation under severe model misspecification or complex, structured priors (Adebiyi et al., 2024, Vakili et al., 2020, Sandberg et al., 3 Feb 2025). Recent progress in regret analysis without discretization (Roy et al., 16 Feb 2026) and in data-driven, structure-exploiting initializations (Adebiyi et al., 2024) points toward more robust and scalable future variants of Gaussian Process Thompson Sampling.