arXiv CS open access

Sharper Regret Bounds for Time-Varying Gaussian Process Bandits with Constant Exploration

2026 · arxiv_cs

arXiv CS · Papers · 2026

Conceptio Open Archive

Analyze & explore related research in Conceptio

Connect this record directly into your research workflow or AI agent. Explore citations across open repositories, query full text, and export structured proof bundles.

Explore in Search →Connect via MCP / API ↗

arXiv:2608.18863v1 [stat.ML] 19 Aug 2026

Sharper Regret Bounds for Time-Varying Gaussian Process Bandits with Constant Exploration

Matthias Mandl Delft Institute of Applied Mathematics Delft University of Technology Mekelweg 4, 2628 CD, Delft M.J.Mandl-2@tudelft.nl

Hanne Kekkonen Delft Institute of Applied Mathematics Delft University of Technology Mekelweg 4, 2628 CD, Delft H.N.Kekkonen@tudelft.nl

Abstract We study Bayesian optimization in a time-varying environment where the unknown reward function evolves according to a Gaussian process drift model. Existing GP-UCB analyses in this setting typically require the exploration parameter to grow with the horizon to maintain uniform confidence bounds. Using per-round local confidence events, we show that GP-UCB can instead be run with a constant exploration parameter and obtain an expected-regret bound whose coefficient depends on the drift rate. We also derive a sharper time-varying maximum-informatione 1/2 ) and gain bound. For the squared exponential kernel, it yields γ̃T /T = O(ϵ 1/4 e expected average regret O(ϵ ) in the persistent-drift regime. The same constantexploration analysis also yields realized-regret guarantees. Simulations support the predicted logarithmic dependence of the bound-suggested exploration parameter on 1/ϵ.

Introduction

Many sequential decision problems require optimizing an unknown objective through a limited number of costly evaluations. Examples include tuning the hyperparameters of machine learning models [Snoek et al., 2012], adapting control policies in robotics [Calandra et al., 2016], and selecting treatment strategies in clinical trials [Villar et al., 2015]. In these settings, each evaluation reveals only imperfect information about the objective, forcing a fundamental trade-off between exploration and exploitation. Bayesian optimization (BO) provides a framework for sequentially optimizing an unknown objective by combining a probabilistic surrogate model with an acquisition rule. In the standard GP bandit formulation, the objective is a fixed function f : S → R. At each round t, the learner selects xt ∈ S and observes yt = f (xt ) +zt , where zt is observation noise. Gaussian processes (GPs) are commonly used as surrogate models because they provide posterior mean and uncertainty estimates in closed form. When decisions affect performance throughout the learning process, performance is typically meaPT sured by cumulative regret, RT = t=1 (maxx∈S f (x) − f (xt )). This criterion originates in the multi-armed bandit literature [Robbins, 1952, Lai and Robbins, 1985, Auer et al., 2002]. Algorithms such as Gaussian Process Upper Confidence Bound (GP-UCB) [Srinivas et al., 2012] and GP-Thompson sampling [Russo and Van Roy, 2014, Chowdhury and Gopalan, 2017] extend bandit methods to continuous domains and achieve sublinear regret under suitable kernel regularity assumptions, with guarantees typically expressed through the maximum information gain [Srinivas et al., 2012, Vakili et al., 2021, Scarlett et al., 2017]. Preprint.

However, many applications are inherently dynamic. System performance may drift due to environmental changes or evolving user behavior, necessitating a time-dependent reward model ft : S → R. When the objective changes over time, sublinear cumulative regret is generally unachievable. The relevant theoretical goal shifts to bounding the average regret, RT /T , as a function of the environmental drift rate. A principled probabilistic model for this setting was introduced by Bogunovic et al. [2016]. In their time-varying Gaussian process (TV-GP) model, the latent function evolves according to a Markov recursion driven by GP innovations, which induces a separable spatio-temporal kernel. Under this model, the learner updates a temporally discounted GP posterior and selects actions using a modified GP-UCB rule. Their analysis shows that near-optimal regret rates are achievable up to logarithmic factors. Subsequent work has studied related models of nonstationarity, including change-point structure [Brunzema et al., 2022], irregular evaluation times [Imamura et al., 2020], and bounded variation in RKHS norms [Deng et al., 2022]. A separate frequentist line studies nonstationary kernelized bandits under sup-norm variation, including algorithms [Zhou and Shroff, 2021, Hong et al., 2023, Iwazaki and Takeno, 2025] and lower bounds [Cai and Scarlett, 2025]. Other extensions include kernel-agnostic adaptation [Ziomek et al., 2024] and stale-data removal in dynamic Bayesian optimization [Bardou et al., 2024, Bardou and Thiran, 2025]. A recurring theoretical and practical limitation in these analyses is the exploration schedule. Existing regret guarantees typically require the GP-UCB exploration parameter, βt , to grow with time t. This leads to deteriorating performance on long time horizons and introduces schedule tuning parameters. In practice, such schedules are difficult to calibrate and induce unnecessary over-exploration. This paper revisits the TV-GP setting of Bogunovic et al. [2016], where the temporal evolution is governed by a drift parameter ϵ > 0, and shows that GP-UCB can be run with a constant exploration parameter while retaining rigorous performance guarantees. Our analysis uses local confidence events at each round rather than a single event that holds uniformly over the horizon. For any fixed drift level, Theorem 1 gives an expected regret bound with linear dependence on T and a coefficient controlled by γ̃T /T . In the persistent-drift regime, this yields a coefficient that decreases as the drift rate decreases. This differs from standard deterministic uniform-confidence GP-UCB analyses in static environments, where the confidence parameter generally depends on the horizon or time. The analysis also gives a sharper dependence on the drift rate through the time-varying maximum e q ), then in the persistent-drift information gain. If the static information gain satisfies γN = O(N regime 2(1−q) 1−q γ̃T E[RT ] e ϵ 4−q , e ϵ 4−q . =O =O T T e 1/4 ). Our analysis also isolates how the For the squared exponential kernel, the latter becomes O(ϵ bound-suggested constant exploration parameter depends on the drift rate. In the regret bounds of Bogunovic et al. [2016], the prescribed exploration schedule is independent of ϵ. In contrast, balancing the upper-bound scaling suggests that the bound-suggested constant β increases as ϵ decreases. This√does not by itself imply uniformly more exploration, since the UCB rule depends on the product of β and the posterior standard deviation. Rather, when drift is slower, past observations remain informative for longer and the posterior variance is typically smaller. A larger β compensates for this reduced variance by assigning more value to reducing the remaining uncertainty, which is also more consequential because information acquired now remains useful over a longer effective horizon. The dependence of the bound-suggested exploration parameter on the temporal dynamics also gives a simple rule for changes in sampling frequency. If observations are collected every ∆ time units, larger ∆ corresponds to greater effective drift between successive evaluations. Temporal correlations then decay faster, past observations become obsolete more quickly, and the analysis predicts a smaller bound-suggested exploration parameter. The remainder of the paper is organized as follows. Section 2 introduces the time-varying GP model, the domain regularity assumptions, and the temporally discounted GP-UCB policy. Section 3 presents the expected-regret bound and the post-hoc and prescribed-confidence guarantees for realized regret. Section 4 derives the drift-dependent information-gain and regret rates and analyzes the bound-suggested exploration parameter. Section 5 provides a supporting empirical simulation study. Section 6 provides the conclusion. Additional proofs and technical results are deferred to the appendices. 2

Problem setting

We start by introducing the dynamic regret setting. Let S ⊆ Rd be compact. At each round t = 1, 2, . . ., the algorithm selects a query point xt ∈ S and observes zt ∼ N (0, σ 2 ) i.i.d.,

yt = ft (xt ) + zt ,

(1)

where ft : S → R is an unknown reward function that changes with time. Let x⋆t ∈ arg max ft (x) x∈S

denote a maximizer at round t. The instantaneous regret is rt := ft (x⋆t ) − ft (xt ) ≥ 0,

(2)

and the cumulative regret over T rounds is RT :=

T X

rt .

(3)

t=1

In the non-stationary setting, sublinear cumulative regret is not expected in general, so the main objective is to control the linear growth rate and its dependence on the temporal properties of ft . We model the function sequence ft using the time-varying Gaussian process construction of Bogunovic et al. [2016]. Let (gt )t≥1 be i.i.d. draws from GP(0, kS ), where kS is a spatial kernel on S. The latent process is defined by √ √ f1 (x) = g1 (x), ft+1 (x) = 1 − ϵ ft (x) + ϵ gt+1 (x), t ≥ 1, (4) with drift parameter ϵ ∈ (0, 1]. For each fixed t, the marginal law remains GP(0, kS ), while the temporal dependence decays geometrically. Equivalently, {ft (x)}(x,t) is a Gaussian process on S × N with separable kernel √ ′ k (x, t), (x′ , t′ ) = kS (x, x′ ) kT (t, t′ ), kT (t, t′ ) = ( 1 − ϵ)|t−t | . (5) We refer to (5) as the time-decayed kernel, since its temporal factor discounts correlations geometrically with temporal separation. We assume throughout that the prior variance is uniformly bounded, kS (x, x) ≤ 1 for all x ∈ S. Given past data Ht−1 = {(xi , yi )}t−1 y1:t−1 = (y1 , . . . , yt−1 )⊤ , i=1 , Gaussian process regression under the kernel (5) yields a Gaussian posterior for ft (x) with mean and variance −1 µt−1 (x) = k̃t−1 (x)⊤ K̃t−1 + σ 2 I y1:t−1 , (6) −1 2 σt−1 (x) = kS (x, x) − k̃t−1 (x)⊤ K̃t−1 + σ 2 I k̃t−1 (x), (7) where

√ t−1 k̃t−1 (x) := kS (xi , x) ( 1 − ϵ)|t−i| i=1 ,

√ t−1 K̃t−1 := kS (xi , xj ) ( 1 − ϵ)|i−j| i,j=1 .

(8)

These expressions are obtained by evaluating the spatio-temporal kernel (5) at the observed pairs √ (xi , i) and the prediction point (x, t). The factors ( 1 − ϵ)|t−i| account for the fact that older observations are less informative about the current function ft . As in standard GP bandit analysis, we measure cumulative uncertainty by the time-varying maximum information gain [Bogunovic et al., 2016], which controls sums of posterior variances in regret bounds [Srinivas et al., 2012, Lemma 5.4]. For a design set A = {(x1 , 1), . . . , (xT , T )}, let fA = (f1 (x1 ), . . . , fT (xT ))⊤ and yA = (y1 , . . . , yT )⊤ . Since yA = fA + zA with independent Gaussian noise, 1 I(yA ; fA ) = log det I + σ −2 K̃A , (9) 2 3

where K̃A is the covariance matrix of fA induced by (5). We define the time-varying maximum information gain as γ̃T := max I(yA ; fA ), (10) x1:T ∈S T

To control the discretization error incurred on a continuous domain, we impose the following standard GP-UCB regularity condition [Srinivas et al., 2012, Bogunovic et al., 2016]. Assumption 1 (Derivative tail bound). There exist constants a ≥ 1 and b > 0 such that, for each t ≥ 1, each spatial dimension j ∈ {1, . . . , d}, and all L > 0, o n L 2 ∂ft (x) > L ≤ a exp − P sup . (11) (j) b x∈S ∂x This condition holds for many commonly used kernels, including the squared exponential kernel and Matérn kernels with smoothness parameter ν > 2, under appropriate differentiability conditions on kS ; see Ghosal and Roy [2006, Theorem 5]. The TV-GP-UCB selection rule is given in Algorithm 1. At round t, the algorithm selects p (12) xt ∈ arg max µt−1 (x) + βt σt−1 (x) . x∈S

The main result of this paper shows that under the time-varying model (4), one can take βt ≡ β constant and obtain a regret guarantee with improved long-horizon scaling. This contrasts with earlier analyses that rely on time-increasing exploration schedules. Algorithm 1 Constant-β Time-Varying GP-UCB (TV-GP-UCB) Require: Domain S, GP prior (0, kS , σ 2 ), drift parameter ϵ ∈ (0, 1], constant exploration parameter β>0 1: for t = 1, 2, . . . do √ 2: Choose xt ∈ arg maxx∈S µt−1 (x) + βσt−1 (x) 3: Observe yt = ft (xt ) + zt 4: Update posterior mean µt and variance σt2 via (6) and (7) 5: end for

Improved Regret Bound

Static GP-UCB analyses usually prove a single confidence event that holds simultaneously over all rounds and all points in a discretization. To keep this event at fixed probability over a horizon T , the per-round failure probabilities must be summable, which forces the exploration parameter βt to grow with t. Existing analyses of the time-varying GP model (4), in particular Bogunovic et al. [2016], use the same time-increasing exploration schedule. This introduces horizon-dependent tuning, leads to average-regret bounds that grow with T , and obscures the dependence of exploration on the drift rate ϵ. We avoid this global-in-time confidence event. When ϵ > 0, older observations are progressively discounted, so posterior uncertainty need not decrease monotonically. This built-in forgetting prevents the algorithm from becoming permanently overconfident at a fixed point and removes the need for a time-increasing exploration schedule in the analysis. We prove that Algorithm 1 with a constant exploration parameter attains an O(T ) expected cumulative regret bound for every fixed ϵ > 0. The proof combines a per-round UCB regret bound on a local confidence event with an algorithm-agnostic bound on the rare complement event. 3.1

Expected regret with a constant exploration parameter

Theorem 1 (Expected regret for constant-β TV-GP-UCB). Let S ⊆ [0, r]d be compact and convex. Assume the time-varying GP model (4) with drift ϵ ∈ (0, 1], spatial kernel kS satisfying kS (x, x) ≤ 1 for all x ∈ S, and the tail condition (11). Assume observations are generated by zt ∼ N (0, σ 2 )

yt = ft (xt ) + zt , 4

with independent noise. For t ≥ 1, define the spatial span Wt := sup ft (x) − inf ft (x), x∈S

x∈S

2 let µw = E[Wt ], and let σw be any finite sub-Gaussian variance proxy for Wt , so that u2 for all u > 0. P{Wt − µw > u} ≤ exp − 2 2σw

Fix δ ∈ (0, 1) and τ ∈ N+ . If the evaluation points {xt }Tt=1 are selected according to Algorithm 1 with any constant exploration parameter satisfying (13) β ≥ 2 log 6δ + d log(τ ) . then the expected cumulative regret satisfies p E[RT ] ≤ C1 T β γ̃T + T C2 + T δC3 ,

(14)

where γ̃T = γ̃T (S, kS , ϵ, σ) is the size-T maximum information gain (10) under the time-decayed kernel, and r p 8 rdb 3ad C1 = , C2 = log , C3 = µw + 2σw 2 log(e/δ). −2 log(1 + σ ) τ δ For comparison, since rt ≤ Wt for every policy, the policy-independent span bound is RT ≤

T X

E[RT ] ≤ µw . T

Wt ,

t=1

Because each ft has marginal law GP(0, kS ), the quantity µw is independent of ϵ. Thus Theorem 1 improves this policy-independent baseline whenever r γ̃T C1 β + C2 + δC3 < µw . T Under the drift-dependent optimization in Section 4.1, the algorithm-specific coefficient decreases as ϵ decreases in the persistent-drift regime, whereas µw does not. The span constants µw and σw are finite by Borell–TIS, since Wt is the supremum of the centered Gaussian difference process ft (x) − ft (y) over S × S and the sample paths are almost surely bounded on the compact domain. The proof, given in Appendix A, constructs a local confidence event Et at each round rather than a single event over the full horizon. The event combines posterior concentration at the selected point, posterior concentration over a spatial grid of size at most τ d , and the derivative tail bound (11). On Et , the instantaneous regret satisfies p rt ≤ 2 β σt−1 (xt ) + C2 , where C2 is the discretization error. On the complement Etc , we use only the deterministic inequality rt ≤ Wt . The sub-Gaussian tail of Wt gives E Wt 1Etc ≤ δC3 . Therefore, p E[rt ] ≤ 2 β E[σt−1 (xt )] + C2 + δC3 . Summing over t and applying the standard information-gain bound on posterior variances yields (14). We use the notation δ and τ from Bogunovic et al. [2016]: δ controls the confidence level and τ d is the spatial grid size. The distinction is that in their analysis δ is tied to a single confidence event over the full horizon, and therefore determines the time-varying exploration schedule βt . Here δ only controls the local event Et used in the expected-regret decomposition. Thus δ and τ are analysis parameters: each choice gives a corresponding constant β and expected-regret bound. Since γ̃T depends on the drift rate ϵ, the bound-balancing choice of (δ, τ ) generally depends on ϵ as well. This is how the bound-suggested constant exploration parameter inherits a dependence on the e ϵα ), the resulting temporal dynamics. Section 4.2 makes this explicit: under the scaling γ̃T = O(T bound-suggested choice satisfies βopt = O(log(1/ϵ)). 5

3.2

Realized cumulative regret and confidence guarantees

Theorem 1 bounds the expected cumulative regret. For realized regret, we distinguish two uses of a constant exploration parameter. In the post-hoc setting, β is fixed without reference to the reporting confidence level ω, whereas in the prescribed-confidence setting, ω is specified before the run and is used to select a constant β. In both cases, β is independent of t and T . Standard analyses of time-varying GP bandits construct probabilistic guarantees using a timeincreasing exploration schedule, typically √ βt = O(log(t/ω)). This enforces a global union bound over the horizon but introduces an O( log T ) factor in the leading information-gain term. We instead retain a constant exploration parameter and control the accumulated contribution of the local confidence failures. Post-hoc confidence.

Let

p UT := C1 T βγ̃T + T C2 denote the deterministic contribution from the UCB confidence widths and the discretization error. The proof bounds the realized pathwise regret by combining this term with the accumulated span over failed confidence events and applying Markov’s inequality to the latter. Proposition 2 (Realized cumulative regret). Assume the setting and parameter choices of Theorem 1. For any target failure probability ω ∈ (0, 1), the realized cumulative regret RT satisfies the following bound with probability at least 1 − ω: T δC3 RT ≤ UT + . (15) ω The proof is deferred to Appendix C.3. Here δ is the per-round analysis parameter that determines the constant β, whereas ω enters only in the reported bound. Thus the algorithm can be run with β fixed independently of ω, and the bound can subsequently be evaluated at any fixed reporting level ω chosen independently of the realized data. The cost of this post-hoc guarantee is the polynomial 1/ω dependence. Prescribed confidence. If the target failure probability is specified before the run, the same result gives logarithmic dependence on 1/ω while retaining a constant exploration parameter. Corollary 3 (Prescribed-confidence realized regret). Assume the setting of Theorem 1 and suppose that, in the regime under consideration, γ̃T e α ), = O(ϵ α > 0. T Fix ω ∈ (0, 1) and choose l m τ = ϵ−α/2 ,

δ = ωϵα/2 ,

together with the constant exploration parameter l m 6 −α/2 β = 2 log α/2 + d log ϵ = O(log(1/ω) + log(1/ϵ)) . ωϵ Then, with probability at least 1 − ω, p RT e ϵα/2 1 + log(1/ω) + log(1/ϵ) . =O T In particular, β remains independent of t and T . Proposition 2 gives r RT γ̃T δC3 ≤ C1 β + C2 + . T T ω Under the choices above, the three terms are respectively p p e ϵα/2 log(1/ω) + log(1/ϵ) , O O ϵα/2 1 + log(1/ω) + log(1/ϵ) , and

p O ϵα/2 1 + log(1/ω) + log(1/ϵ) ,

which yields (16). 6

(16)

Policy-independent fallback.

Independently of the sampling policy, RT ≤

T X

Wt .

t=1

√

1 − ϵ, this gives with probability at least 1 − ω, r T X 1+ρ RT ≤ Wt ≤ T µw + 8T log(1/ω). 1−ρ t=1

As established in Appendix C.2, with ρ =

(17)

This bound does not use the sampling decisions and is therefore not a replacement for the algorithmspecific guarantees above.

Drift-dependent regret rates and exploration

Theorem 1 gives a constant-exploration regret bound whose leading term is controlled by γ̃T /T . We first combine it with the time-varying information-gain bound in Appendix B to derive explicit drift-dependent rates, and then examine the corresponding bound-suggested choice of β. 4.1

Information-gain and regret rates

The time-varying information-gain bound in Appendix B gives an explicit dependence on the drift rate. Suppose the static maximum information gain satisfies e q ), γN = O(N

q ∈ [0, 1).

At the polynomial scale, the corresponding block length is 2

Nϵ = ϵ− 4−q . In the persistent-drift regime where T exceeds this effective block length, Corollary B.4 gives 2(1−q) γ̃T e ϵ 4−q . =O (18) T Balancing the remaining terms in Theorem 1 as in Section 4.2 then yields 1−q E[RT ] e ϵ 4−q . =O T

(19)

For the squared exponential kernel, γN is polylogarithmic in N [Srinivas et al., 2012, Vakili et al., 2021, Iwazaki, 2025], so q = 0 up to logarithmic factors and γ̃T e 1/2 ), = O(ϵ T

E[RT ] e 1/4 ). = O(ϵ T

For a Matérn kernel with smoothness ν > 2, the static rate q = d/(2ν + d) [Vakili et al., 2021, Iwazaki, 2025] gives 4ν 2ν γ̃T E[RT ] e ϵ 8ν+3d e ϵ 8ν+3d =O , =O . T T These rates are not fixed-T , ϵ → 0 statements. If Nϵ > T , taking Ñ = T in Theorem B.3 gives √ σ −2 ϵ T T 2 − 1 √ √ γ̃T ≤ γT + γT , 2 3 and hence γ̃T → γT for fixed T as ϵ → 0. The constant-exploration result is therefore a positive-drift result and is not claimed to recover the usual static sublinear-regret guarantee at ϵ = 0; the static case requires a separate analysis [Takeno et al., 2025]. 7

Comparison with previous rates. The information-gain improvement is separate from the e q ), the original block perturbation in Bogunovic confidence-schedule improvement. If γN = O(N et al. [2016] gives the first rate below. The O(ϵN 5/2 ) block order from the unit-time specialization of Imamura et al. [2020] gives the second, and Theorem B.3 gives the third: 1−q 2(1−q) 2(1−q) γ̃T e ϵ 3−q , e ϵ 5−2q , e ϵ 4−q . =O O O T p The time-growing confidence schedule contributes a factor of order log(T /ω) to the leading regret term, whereas the prescribed-confidence constant-β choice uses β = O(log(1/ω) + log(1/ϵ)) , so the corresponding confidence factor is independent of T . 4.2

Bound-suggested β dependence on ϵ

Here β denotes any constant exploration parameter satisfying (13), while βopt denotes the boundsuggested value obtained from the scaling choice of the analytical parameters (δ, τ ), rather than a policy-level optimal parameter. The drift parameter ϵ is already required to construct the TV-GP posterior, so using it to guide β introduces no additional model parameter. The scaling βopt = O(log(1/ϵ)) below is a bound-balancing recommendation rather than a validity condition; Theorem 1 holds for any constant β satisfying (13). The analysis of Bogunovic et al. [2016] specifies the time-growing schedule ! r 2 t2 π 2 t2 daπ βt = 2 log + 2d log rdbt2 log , (20) 2δ 2δ where (a, b, r, d) are spatial domain and smoothness constants, and δ is a failure probability. This schedule is determined by the spatial discretization argument and is not optimized as a function of the drift parameter ϵ. Our theorem does not prescribe β as an explicit function of ϵ directly. Instead, ϵ enters through the time-decayed maximum information gain γ̃T . Balancing the upper-bound terms through the analytical parameters (δ, τ ) therefore induces an ϵ-dependent bound-suggested constant βopt (ϵ). To make this dependence explicit, suppose that, for the chosen spatial kernel and in the long-horizon regime where T exceeds the effective block length, the time-decayed information gain satisfies e ϵα ) , γ̃T = O(T α > 0. (21) For fixed T , this scaling should not be interpreted as an ϵ → 0 limit, since γ̃T then approaches the static information gain. Appendix B verifies such a scaling for common kernels through the standard growth rate of the time-invariant information gain. Taking equality in (13), β = 2(log(6/δ)+d log τ ), p while C2 = (rdb/τ ) log(3ad/δ). Thus the three per-round contributions in the expected regret p bound are controlled by β γ̃T /T , C2 , and δC3 . √ Under (21), the information-gain contribution per round scales as β ϵα/2 up to logarithmic factors. For ϵ ∈ (0, 1), the choice below balances the polynomial dependence on ϵ across the information-gain, discretization, and failure-event terms and is sufficient to identify the logarithmic dependence of βopt on 1/ϵ: l m δ = ϵα/2 , τ = ϵ−α/2 . (22) With this choice, the discretization and failure-event contributions match the same polynomial order, up to logarithmic factors: p C2 = O ϵα/2 log(ϵ−1 ) , (23) p δC3 = O ϵα/2 log(ϵ−1 ) . (24) Substituting (22) into (13) at equality gives βopt (ϵ) ≈ 2 log(6) + α(1 + d) log

. (25) ϵ Across environments with the same spatial kernel but different drift rates, the bound suggests that βopt grows logarithmically with 1/ϵ. This scaling should be interpreted together with the posterior variance: smaller ϵ also makes past data more informative and typically reduces σt−1 (x). 8

4.3

Sampling frequency and effective drift

The same calculation gives a simple rule for changing the sampling frequency. Suppose the underlying + process evolves at drift rate ϵ, but observations √ are collected every ∆ ∈ N time units. Then the correlation between successive observations is ( 1 − ϵ)∆ , which corresponds to the effective discretetime drift ϵ∆ = 1 − (1 − ϵ)∆ . (26) Larger ∆ increases ϵ∆ , so past observations become less relevant. Substituting this effective drift into (25) yields the heuristic 1 1 βopt (ϵ∆ ) ∝ log = log . (27) ϵ∆ 1 − (1 − ϵ)∆ This adjustment reflects the fact that less frequent sampling increases the effective drift between observations and reduces the useful lifetime of past data. 4.4

Average-reward perspective

The dynamic setting also changes the interpretation of the benchmark. Let Mt := supx∈S ft (x) denote the round-t optimal reward. In the static case, Mt is the same random variable for all t, so the benchmark remains tied to a single realization of the objective. Regret bounds must therefore account for the joint distribution of the realized function and the algorithm’s decisions. When ϵ > 0, this coupling is weaker in the long run. The sequence (Mt )t≥1 is stationary and has exponentially decaying temporal dependence. As shown in Appendix C.1, T

1X Pr Mt −−−−→ E[M1 ]. T →∞ T t=1 Thus the oracle part of the regret concentrates around a deterministic quantity. Consequently, asymptotic average regret can be viewed as the gap between the deterministic benchmark E[M1 ] and the algorithm’s realized average reward. Cumulative regret remains the relevant finite-horizon criterion, but this perspective separates the intrinsic variability of the moving optimum from the performance loss due to the algorithm’s decisions.

Simulation study

We empirically evaluate constant-β GP-UCB under the TV-GP model on the spatial domain S = [0, 1]2 . The domain is discretized using a 20 × 20 grid, giving a state dimension of N = 400. The latent function uses a radial basis function (RBF, squared exponential) kernel with lengthscale l = 0.1 and prior variance 1.0, while the observation noise variance is fixed at σ 2 = 0.01. We simulate T = 10,000 rounds for each parameter pair (ϵ, β) and record the time-averaged regret RT /T over 100 independent trials. This horizon is used as a long-run stress test rather than as an application-specific benchmark; it is made feasible by the Kalman implementation described below. The simulation setup and compute resources are reported in Appendix E, and the accompanying code archive contains the scripts used to reproduce Figure 1. Exact inference via Kalman filtering. A direct GP implementation would require repeatedly conditioning on the full observation history, making long-horizon simulations computationally expensive. To avoid this bottleneck, we exploit the autoregressive structure of the TV-GP prior. On the finite grid, the latent function can be represented as a state vector ft ∈ R400 , and the TV-GP model becomes a linear Gaussian state-space model. We therefore compute the GP posterior exactly using a Kalman filter. This reduces the per-round posterior update from the O(t3 ) history-dependent cost of naive GP regression to a fixed-grid O(N 2 ) recursion. The full recursion is given in Appendix D. Results. The left panel of Figure 1 displays the average regret RT /T as a function of the constant exploration parameter β for three drift rates, ϵ ∈ {10−1 , 10−2 , 10−3 }. For each drift rate, let β̂opt denote the value of β on the tested grid that minimizes the empirical average regret. The main empirical test is whether β̂opt increases with − log(ϵ) as predicted by Section 4.2. The empirical 9

ε = 10−1

1.5 1.4

β̂opt = 0.84[−log(ε)] + 0.29

1.3 5

0.6

β̂opt

RT /T

0.7

0.5

0.20

0.10

ε = 10−2

ε = 10−3

2 1

constant β heuristic βt

0.15 0

−log(ε)

Figure 1: Simulation study for constant-β GP-UCB under the time-varying GP model. Left: timeaveraged regret RT /T over T = 10,000 rounds as a function of the constant exploration parameter β for three drift rates ϵ ∈ {10−1 , 10−2 , 10−3 }; shaded bands denote 95% confidence intervals over 100 independent trials, and dashed horizontal lines show the logarithmic time-varying-βt heuristic. Right: bootstrapped estimates of the empirical regret-minimizing constant β̂opt versus − log(ϵ) across the full ϵ grid. Points denote mean bootstrap estimates, vertical bars denote 95% percentile intervals, and the dashed line shows a least-squares linear fit used to visualize the log-linear trend. minimizer shifts toward larger values as ϵ decreases: β̂opt ≈ 2.12 for ϵ = 0.1, β̂opt ≈ 4.34 for ϵ = 0.01, and β̂opt ≈ 6.36 for ϵ = 0.001. Because the exact theoretical schedule in Equation (20) is overly conservative in practice [Srinivas et al., 2012, Bogunovic et al., 2016], empirical evaluations often use tuned logarithmic schedules. Following Bogunovic et al. [2016], we compare against βt = c1 log(c2 t) with c1 = 0.8 and c2 = 4. Across all three drift rates, a broad interval of constant β values achieves lower average regret than this time-growing heuristic over the 10,000-step horizon. The right panel of Figure 1 examines this dependence across the full grid of drift rates. The relationship between β̂opt and − log(ϵ) is approximately linear, matching the qualitative prediction βopt (ϵ) ∝ log(1/ϵ) from Equation (25). This agreement is nontrivial: the theory balances terms in an upper bound, and such bounds need not match the exploration scale that performs best in finite simulations. We do not interpret the fitted slope as an estimate of the information-gain exponent α. Such an interpretation would require the analytical discretization choice τ = τ (ϵ) and the global spatial union-bound factor in Equation (25) to be tight. Neither condition holds here: the simulation uses a fixed 20 × 20 action grid, and the adaptive policy concentrates its evaluations near high-reward regions rather than exploring the full domain uniformly. Thus the slope is best viewed as empirical support for the logarithmic dependence of the regret-minimizing constant on the drift rate, not as a quantitative estimate of the information-gain scaling.

Conclusion

Existing theoretical guarantees for Gaussian process bandits in time-varying environments typically require exploration schedules that grow with time. We showed that, under the TV-GP drift model, a constant exploration parameter can instead be analyzed using local confidence events at each round rather than a single horizon-wide confidence event. The resulting expected-regret bound has a drift-dependent leading term controlled by the time-varying information gain, and the same local-confidence analysis yields realized-regret guarantees with β constant in t and T . e q ), We also derived a sharper bound on the time-varying maximum information gain. If γN = O(N e (1−q)/(4−q) ); for the squared exponential the resulting persistent-drift average-regret rate is O(ϵ 10

e 1/4 ). These rates apply when the optimized block length fits within the horizon kernel this gives O(ϵ and are not fixed-T , ϵ → 0 statements. Separately, balancing the remaining terms in the regret bound suggests a bound-suggested constant exploration parameter with βopt (ϵ) ∝ log(1/ϵ). This provides a simple tuning principle when changing the sampling frequency changes the effective drift between observations. The same Markov structure also implies exponential covariance decay for the oracle-value process, so long-run average performance can be compared against a deterministic benchmark. Several extensions remain open. Most importantly, the present analysis relies on the geometric temporal kernel induced by the Markov drift model. Extending constant-exploration guarantees to more general temporal kernels would clarify which forms of nonstationarity are sufficient to prevent posterior overconfidence. Another important direction is to remove the assumption that the drift parameter is known and fixed, allowing β to be adapted online when the rate of temporal variation is unknown or changes over time. Finally, sharper concentration tools for the realized regret of TV-GP-UCB could improve the current high-probability bounds and better reflect the empirical stability observed in the simulations.

Acknowledgments and Disclosure of Funding The first author gratefully acknowledges support from the NWO Spinoza Prize awarded to Aad van der Vaart.

References Robert J Adler and Jonathan E Taylor. Random fields and geometry. Springer Science & Business Media, 2009. Peter Auer, Nicolo Cesa-Bianchi, and Paul Fischer. Finite-time analysis of the multiarmed bandit problem. Machine learning, 47:235–256, 2002. Anthony Bardou and Patrick Thiran. Time-varying Bayesian optimization without a metronome. arXiv preprint arXiv:2501.18963, 2025. Anthony Bardou, Patrick Thiran, and Giovanni Ranieri. This too shall pass: Removing stale observations in dynamic Bayesian optimization. Advances in Neural Information Processing Systems, 37:42696–42737, 2024. Ilija Bogunovic, Jonathan Scarlett, and Volkan Cevher. Time-varying Gaussian process bandit optimization. In Artificial Intelligence and Statistics, pages 314–323. PMLR, 2016. Paul Brunzema, Alexander von Rohr, Friedrich Solowjow, and Sebastian Trimpe. Event-triggered time-varying Bayesian optimization. arXiv preprint arXiv:2208.10790, 2022. Xu Cai and Jonathan Scarlett. Lower bounds for time-varying kernelized bandits. In Proceedings of The 28th International Conference on Artificial Intelligence and Statistics, volume 258 of Proceedings of Machine Learning Research, pages 73–81. PMLR, 2025. Roberto Calandra, André Seyfarth, Jan Peters, and Marc Peter Deisenroth. Bayesian optimization for learning gaits under uncertainty: An experimental comparison on a dynamic bipedal walker. Annals of Mathematics and Artificial Intelligence, 76(1):5–23, 2016. Sayak Ray Chowdhury and Aditya Gopalan. On kernelized multi-armed bandits. In International Conference on Machine Learning, pages 844–853. PMLR, 2017. Yuntian Deng, Xingyu Zhou, Baekjin Kim, Ambuj Tewari, Abhishek Gupta, and Ness Shroff. Weighted Gaussian process bandits for non-stationary environments. In International Conference on Artificial Intelligence and Statistics, pages 6909–6932. PMLR, 2022. Subhashis Ghosal and Anindya Roy. Posterior consistency of Gaussian process prior for nonparametric binary regression. The Annals of Statistics, 34(5):2413–2429, 2006. 11

Kihyuk Hong, Yuhang Li, and Ambuj Tewari. An optimization-based algorithm for non-stationary kernel bandits without prior knowledge. In International Conference on Artificial Intelligence and Statistics, pages 3048–3085. PMLR, 2023. Hideaki Imamura, Nontawat Charoenphakdee, Futoshi Futami, Issei Sato, Junya Honda, and Masashi Sugiyama. Time-varying Gaussian process bandit optimization with non-constant evaluation time. arXiv preprint arXiv:2003.04691, 2020. Shogo Iwazaki. Improved regret bounds for Gaussian process upper confidence bound in Bayesian optimization. Advances in Neural Information Processing Systems, 38:96922–96964, 2025. Shogo Iwazaki and Shion Takeno. Near-optimal algorithm for non-stationary kernelized bandits. In Proceedings of The 28th International Conference on Artificial Intelligence and Statistics, volume 258 of Proceedings of Machine Learning Research, pages 406–414. PMLR, 2025. Tze Leung Lai and Herbert Robbins. Asymptotically efficient adaptive allocation rules. Advances in applied mathematics, 6(1):4–22, 1985. Herbert Robbins. Some aspects of the sequential design of experiments. Bulletin of the American Mathematical Society, 58(5):527–535, 1952. Daniel Russo and Benjamin Van Roy. Learning to optimize via posterior sampling. Mathematics of Operations Research, 39(4):1221–1243, 2014. Jonathan Scarlett, Ilija Bogunovic, and Volkan Cevher. Lower bounds on regret for noisy Gaussian process bandit optimization. In Conference on Learning Theory, pages 1723–1742. PMLR, 2017. Jasper Snoek, Hugo Larochelle, and Ryan P Adams. Practical Bayesian optimization of machine learning algorithms. Advances in neural information processing systems, 25, 2012. Niranjan Srinivas, Andreas Krause, Sham M Kakade, and Matthias W Seeger. Information-theoretic regret bounds for Gaussian process optimization in the bandit setting. IEEE transactions on information theory, 58(5):3250–3265, 2012. Shion Takeno, Yu Inatsu, and Masayuki Karasuyama. Regret analysis for randomized Gaussian process upper confidence bound. Journal of Artificial Intelligence Research, 84(18), 2025. doi: 10.1613/jair.1.19393. Sattar Vakili, Kia Khezeli, and Victor Picheny. On information gain and regret bounds in Gaussian process bandits. In International Conference on Artificial Intelligence and Statistics, pages 82–90. PMLR, 2021. Sofía S Villar, Jack Bowden, and James Wason. Multi-armed bandit models for the optimal design of clinical trials: benefits and challenges. Statistical science: a review journal of the Institute of Mathematical Statistics, 30(2):199, 2015. Xingyu Zhou and Ness Shroff. No-regret algorithms for time-varying Bayesian optimization. In 2021 55th Annual Conference on Information Sciences and Systems (CISS), pages 1–6. IEEE, 2021. Juliusz Ziomek, Masaki Adachi, and Michael A Osborne. Time-varying Gaussian process bandits with unknown prior. arXiv preprint arXiv:2402.01632, 2024.

Improved Expected Regret Bound

This appendix proves Theorem 1. The argument proceeds in three steps: (i) bounding the instantaneous regret rt for Algorithm 1 conditionally on a high-probability event; (ii) bounding the expected regret on the complement event via concentration inequalities for the extrema of Gaussian processes; and (iii) summing the expected regret over t and bounding the sum of predictive variances by the maximum information gain. A.1

Algorithm-specific bound

Lemma A.1 (Instantaneous regret under constant β). Consider the setting of Theorem 1. Fix δ ∈ (0, 1) and τ ∈ N+ . Let S ∗ ⊂ S be a discretization with |S ∗ | ≤ τ d satisfying rd ∥x − [x]∥1 ≤ , for all x ∈ S, (28) τ where [x] ∈ S ∗ is an ℓ1 -closest grid point. Let β and C2 be defined as in Theorem 1. If xt is selected according to Algorithm 1, then for each fixed t, the instantaneous regret satisfies p rt = ft (x∗t ) − ft (xt ) ≤ 2 β σt−1 (xt ) + C2 (29) with probability at least 1 − δ. Proof. We bound the instantaneous regret by conditioning on the intersection of three high-probability events at a fixed time t. Let Et = E1,t ∩ E2,t ∩ E3,t , defined as follows: (i) Pointwise concentration: Given the history Ht−1 , the query point xt is deterministic. Because 2 ft (xt ) ∼ N (µt−1 (xt ), σt−1 (xt )), standard Gaussian tail bounds ensure the event p |ft (xt ) − µt−1 (xt )| ≤ β σt−1 (xt ) (30) holds with probability at least 1 − δ/3 under the specified constant β. (ii) Uniform grid concentration: Applying a union bound over the |S ∗ | ≤ τ d elements of the discretization, the event p |ft (x) − µt−1 (x)| ≤ β σt−1 (x), for all x ∈ S ∗ (31) holds with probability at least 1 − δ/3. (iii) Lipschitz continuity: By the derivative tail assumption (11) and a union bound over the d p spatial dimensions, setting L = b log(3ad/δ) ensures the event |ft (x) − ft (x′ )| ≤ L∥x − x′ ∥1 , for all x, x′ ∈ S holds with an unconditional probability of at least 1 − δ/3 under the prior measure.

(32)

Because events (i) and (ii) hold conditionally with probability 1 − δ/3 for any history Ht−1 , the law of total expectation guarantees they also hold with an unconditional probability of at least 1 − δ/3. Applying a union bound over the unconditional probabilities of these three events yields P(Etc ) ≤ δ, guaranteeing the joint event holds unconditionally with P(Et ) ≥ 1 − δ. Conditioned on the realization of Et , the continuity event (iii) and the spatial discretization bound (28) guarantee that the approximation error for any point in the domain is bounded by L(rd/τ ) = C2 . Evaluating this spatial bound at the instantaneous optimum x∗t yields ft (x∗t ) ≤ ft ([x∗t ]) + C2 . Because the projection [x∗t ] ∈ S ∗ , we apply the uniform concentration event (ii) to bound ft ([x∗t ]), followed by the selection rule of Algorithm 1, which maximizes the upper confidence bound at xt : p ft (x∗t ) ≤ µt−1 ([x∗t ]) + β σt−1 ([x∗t ]) + C2 p ≤ µt−1 (xt ) + β σt−1 (xt ) + C2 . (33) Finally, applying the pointwise concentration event (i) at xt yields p µt−1 (xt ) ≤ ft (xt ) + β σt−1 (xt ). Substituting this inequality into (33) bounds the instantaneous regret: p rt = ft (x∗t ) − ft (xt ) ≤ 2 β σt−1 (xt ) + C2 , (34) which concludes the proof. 13

A.2

Algorithm-agnostic tail control

The span Wt can be written as Wt = sup ft (x) − ft (y) . x,y∈S

Thus Wt is the supremum of a centered Gaussian process indexed by S × S. Since the sample paths are almost surely bounded on the compact domain and kS (x, x) ≤ 1, Borell–TIS implies that 2 Wt − E[Wt ] has a sub-Gaussian upper tail with a finite variance proxy. We denote this proxy by σw . Lemma A.2 (Tail expectation bound). Assume Wt − µw has the sub-Gaussian upper tail u2 for all u > 0. P{Wt − µw > u} ≤ exp − 2 2σw Then, for every event A with P(A) ≤ δ, p E[Wt 1A ] ≤ δ µw + 2σw 2 log(e/δ) . Proof. Write X = (Wt − µw )+ . Then E[Wt 1A ] ≤ µw P(A) + E[X1A ]. Let q = σw

2 log(e/δ). We split E[X1A ] ≤ qP(A) + E[X1{X>q} ].

Using the sub-Gaussian tail and integrating, Z ∞ E[X1{X>q} ] = qP(X > q) +

P(X > u) du. q

Since q = σw

p 2 log(e/δ), we have q2 P(X > q) ≤ exp − 2 2σw

δ . e

Moreover, the standard Gaussian tail integral bound gives Z ∞ u2 q2 δ σ2 σ2 δ exp − 2 du ≤ w exp − 2 = w ≤ q , 2σw q 2σw q e e q where the last inequality uses q ≥ σw . Therefore, E[X1{X>q} ] ≤ 2q

δ ≤ qδ. e

Combining this with qP(A) ≤ qδ yields E[X1A ] ≤ 2qδ.

A.3

Expected cumulative regret

Proof of Theorem 1. Fix δ ∈ (0, 1) and τ ∈ N+ , and define β and C2 as in Theorem 1. By Lemma A.1, for each time step t there exists a confidence√event Et with P(Et ) ≥ 1 − δ such that, conditioned on Et , the instantaneous regret satisfies rt ≤ 2 β σt−1 (xt ) + C2 . Because instantaneous regret is always non-negative, we can bound it pathwise over the entire probability space by introducing the spatial span Wt = maxx ft (x) − minx ft (x) as the worst-case penalty when the confidence event fails: p rt ≤ 2 β σt−1 (xt ) + C2 + Wt 1Etc , (35) 14

where 1Etc is the indicator for the complement event. Taking the unconditional expectation of both sides, and applying Lemma A.2 to bound the expected tail contribution E[Wt 1Etc ] ≤ δC3 , yields the valid deterministic bound: p E[rt ] ≤ 2 β E[σt−1 (xt )] + C2 + δC3 . (36) Summing (36) over t = 1, . . . , T , we obtain: " T # X p σt−1 (xt ) + T C2 + T δC3 . E[RT ] ≤ 2 β E

(37)

t=1

To bound the sum of the posterior standard deviations, we apply the Cauchy-Schwarz inequality strictly inside the expectation: " T # "v # u T u X X t 2 E σt−1 (xt ) ≤ E T σt−1 (xt ) . (38) t=1

t=1

We then bound the sum of predictive variances by the maximum information gain. For every realized query sequence x1:T , the Schur-complement determinant identity gives T 1 1X 2 log 1 + σ −2 σt−1 (xt ) = log det IT + σ −2 K̃T (x1:T ) ≤ γ̃T . 2 t=1 2

(39)

Combining this pathwise identity with the variance-to-logarithm inequality from Srinivas et al. [2012, Lemma 5.4] under the assumption that kS (x, x) ≤ 1 yields almost surely T X

2 σt−1 (xt ) ≤

t=1

2 γ̃T . log(1 + σ −2 )

(40)

PT 2 (xt ) ≤ C1 β γ̃T . Because γ̃T is defined as Defining C1 = 8/ log(1 + σ −2 ), we have t=1 4β σt−1 the maximum information gain over all possible sequences, it is a deterministic upper bound. The expectation of a deterministic bound is simply the bound itself. Substituting this into the cumulative sum yields the final result: p E[RT ] ≤ C1 T β γ̃T + T C2 + T δC3 . (41)

The explicit dependence of γ̃T on T and ϵ for the time-decayed kernel is detailed in Appendix B.

Maximum Information Gain

Theorem B.3 (Time-varying maximum information gain). Let γ̃T := supx1:T ∈S T 12 log det(IT + σ −2 K̃T ) denote the maximum information gain under the separable spatio-temporal kernel √ ′ k̃((x, t), (x′ , t′ )) = kS (x, x′ ) ρ|t−t | with ρ = 1 − ϵ ∈ [0, 1). Assume |kS (x, x′ )| ≤ 1 for all x, x′ ∈ S. Then for any block length Ñ ∈ N+ , ! ! p p l m σ −2 ϵ Ñ Ñ 2 − 1 √ σ −2 ϵ Ñ Ñ 2 − 1 √ T T √ √ γ̃T ≤ γÑ + γÑ ≤ Ñ +1 γÑ + γÑ , Ñ 2 3 2 3 (42) where γÑ := supx ∈S Ñ 12 log det(IÑ + σ −2 KÑ ) is the standard time-invariant maximum infor1:Ñ mation gain for f ∼ GP(0, kS ). Proof. Step 1 (block decomposition). Partition the sequence {1, . . . , T } into m = ⌈T /Ñ ⌉ consecutive blocks, each of length at most Ñ . This partition is used only in the analysis; the GP posterior is 15

neither reset nor truncated at block boundaries. Let (f (i) , y(i) ) denote the latent function values and noisy observations within block i. By the chain rule for mutual information, I(fT ; yT ) =

m X

I(fT ; y(i) | y(<i) ).

i=1 (i)

(i)

Because y = f + z with z drawn independently of (fT , y(<i) ), the variables satisfy the Markov property y(i) ⊥ ⊥ (fT \ f (i) ) | f (i) . Consequently, I(fT ; y(i) | y(<i) ) = I(f (i) ; y(i) | y(<i) ). Furthermore, because the additive Gaussian noise z(i) is independent of the sequence history, the conditional entropy simplifies to H(y(i) | f (i) , y(<i) ) = H(z(i) ). While conditioning reduces entropy (H(y(i) | y(<i) ) ≤ H(y(i) )), it does not universally reduce mutual information. However, due to this independent noise structure, the explicit algebraic cancellation holds: I(f (i) ; y(i) | y(<i) ) = H(y(i) | y(<i) ) − H(z(i) ) ≤ H(y(i) ) − H(z(i) ) = I(f (i) ; y(i) ). It follows that I(fT ; yT ) ≤

m X

(43)

I(f (i) ; y(i) ).

i=1

Pm Taking the supremum over the choice of locations x1:T ∈ S T yields γ̃T ≤ i=1 γ̃ni ≤ m γ̃Ñ , where ni ≤ Ñ is the size of block i and the sequence γ̃n is monotonically non-decreasing. Step 2 (one-block bound). Fix a block of length Ñ . Define its spatial kernel matrix as K := KÑ and its temporal correlation matrix as D := DÑ , where Dij = ρ|i−j| . The spatio-temporal covariance over this block is the Hadamard product K̃ := K ◦ D. By the Schur product theorem, K̃ is positive semi-definite because both K ⪰ 0 and D ⪰ 0. Let 1 denote the Ñ × Ñ all-ones matrix and define the difference matrix A := K̃ − K = K ◦ (D − 1). Let M := IÑ + σ −2 K ≻ 0,

E := σ −2 A,

so that IÑ + σ −2 K̃ = M + E ≻ 0. Applying the congruence factorization M + E = M 1/2 IÑ + B M 1/2 with B := M −1/2 EM −1/2 yields 1 1 1 2 log det(M + E) = 2 log det(M ) + 2 log det(IÑ + B).

The matrix B is symmetric and satisfies IÑ + B ≻ 0. Therefore, every eigenvalue of B is strictly greater than −1. Using the scalar inequality log(1 + u) ≤ u for all u > −1, we bound the log-determinant term by the trace: log det(IÑ + B) =

Ñ X

Ñ X log 1 + λj (B) ≤ λj (B) = tr(B).

j=1

Substituting this trace inequality provides the upper bound 1 1 1 2 log det(M + E) ≤ 2 log det(M ) + 2 tr(B). By the cyclicity of the trace operator, tr(B) = tr(M −1 E). Since Dii = 1, the diagonal entries of

A = K ◦ (D − 1) vanish, and hence tr(E) = 0. Therefore, tr(B) = tr (M −1 − IÑ )E . Using the Frobenius inner product, tr(B) ≤ tr (M −1 − IÑ )E

Write

≤ ∥M −1 − IÑ ∥F ∥E∥F .

1 log det(IÑ + σ −2 K). 2 Let λ1 , . . . , λÑ be the eigenvalues of K and set uj = σ −2 λj ≥ 0. Then 2 X Ñ Ñ X uj ∥M −1 − IÑ ∥2F = ≤ log(1 + uj ) = 2Γ(K). 1 + uj j=1 j=1 Γ(K) :=

Here the inequality follows from

u 1+u

≤ log(1 + u),

u ≥ 0,

since the difference has value zero at u = 0 and derivative (1 + u2 )/(1 + u)3 > 0. Hence p ∥M −1 − IÑ ∥F ≤ 2Γ(K). It remains to bound ∥E∥F . As above, |Dij − 1| = 1 − ρ|i−j| ≤ ϵ|i − j|. Since |Kij | ≤ 1, Ñ X

∥A∥2F ≤ ϵ2

(i − j)2 =

i,j=1

and therefore ∥E∥F = σ

−2

ϵ2 Ñ 2 (Ñ 2 − 1) , 6

p σ −2 ϵ Ñ Ñ 2 − 1 √ . ∥A∥F ≤ 6

Combining the two norm bounds gives p σ −2 ϵ Ñ Ñ 2 − 1 p 1 −2 √ log det(IÑ + σ K̃) ≤ Γ(K) + Γ(K). 2 2 3 √ Since s 7→ s + c s is increasing for s ≥ 0 and Γ(K) ≤ γÑ , p σ −2 ϵ Ñ Ñ 2 − 1 √ √ γ̃Ñ ≤ γÑ + γÑ . 2 3 Summing this bound over the m blocks in Step 1 completes the proof. Remark. The block decomposition follows the approach of Bogunovic et al. [2016]. For unitspaced observations, an O(ϵÑ 5/2 ) block perturbation also follows from the more general analysis of Imamura et al. [2020]. Without the zero-trace step above, the congruence argument recovers this p order directly by combining ∥M −1 ∥F ≤ Ñ with ∥A∥F ≤ ϵÑ 2 , without requiring sign restrictions on the perturbation. The additional identity tr(E) = 0 allows M −1 to be replaced by M −1 − IÑ , yielding the sharper bound in Theorem B.3. B.1

Asymptotic Scaling with Drift

Theorem B.3 gives a two-regime dependence on the drift rate once the block length is chosen according to the growth of the time-invariant maximum information gain. Corollary B.4 (Information gain scaling). Suppose the spatial kernel admits a time-invariant maximum information gain satisfying e Ñ q ) γÑ = O( for some capacity exponent q ∈ [0, 1). Define the polynomial block scale 2

Nϵ := ϵ− 4−q , and choose Ñ = min {T, ⌈Nϵ ⌉} . Then

2(1−q) e T q + T ϵ 4−q . γ̃T = O

In particular, in the persistent-drift regime T ≳ Nϵ , 2(1−q) e T ϵ 4−q , γ̃T = O 17

(44)

(45)

whereas for T ≲ Nϵ one has e q ). γ̃T = O(T For fixed T , choosing Ñ = T in Theorem B.3 gives √ σ −2 ϵ T T 2 − 1 √ √ γT , γ̃T ≤ γT + 2 3

(46)

so the perturbation vanishes as ϵ → 0 and the static information gain is recovered. p Proof. For Ñ ≤ T , Theorem B.3 and Ñ 2 − 1 ≤ Ñ give e T Ñ q−1 + T ϵÑ 1+q/2 . e T Ñ q + ϵÑ 2+q/2 =O γ̃T = O Ñ Balancing the polynomial orders gives Ñ q−1 ≍ ϵÑ 1+q/2 ,

or equivalently

Ñ ≍ ϵ− 4−q = Nϵ .

Thus, when T ≳ Nϵ , choosing Ñ at this scale gives 2(1−q) e T ϵ 4−q . γ̃T = O e q ) give When T ≲ Nϵ , choose Ñ = T . Equation (46) and γT = O(T e T q + ϵT 2+q/2 . γ̃T = O e q ). Combining The condition T ≲ Nϵ implies ϵT 2−q/2 ≲ 1, so the perturbation term is at most O(T the two regimes yields (44). For fixed T , the temporal matrix satisfies Dij → 1 as ϵ → 0. For the continuous kernels considered here on the compact domain S, continuity of the log determinant under maximization therefore gives γ̃T → γT .

Concentration of RT

This appendix analyzes the concentration properties of the time-varying setting. Section C.1 establishes the convergence of the oracle-value process via exponential covariance decay, which justifies the asymptotic average reward perspective discussed in Section 4.4. Section C.2 proves the spanonly concentration fallback, and Section C.3 proves the algorithm-specific realized-regret bound in Proposition 2. C.1

Concentration of the oracle-value process via covariance decay

Define Mt := supx∈S ft (x) for t ≥ 1. Under the TV-GP dynamics, the sequence {ft } is strictly stationary with each ft ∼ GP(0, kS ). Consequently, {Mt } is stationary with a constant expectation E[Mt ]. Assume the functions ft have almost surely bounded sample paths on S and are separable, implying Mt < ∞ almost surely. By the Borell–TIS inequality [Adler and Taylor, 2009, Theorem 2.1.1], Mt is sub-Gaussian around its mean: u2 2 Pr |Mt − E[Mt ]| > u ≤ 2 exp − 2 , σM := sup Var ft (x) ≤ 1. 2σM x∈S This concentration ensures E[Mt2 ] < ∞. The following generic lemma controls the variance of a sum of random variables that exhibit exponentially decaying covariance. Lemma C.5 (Averaging under exponential covariance decay). Let {Xt }t≥1 be random variables with E[Xt2 ] < ∞ and suppose there exist m > 0 and c ∈ [0, 1) such that Cov(Xi , Xj ) ≤ m c|i−j| PT for all i, j. Let ST := t=1 Xt . Then: 18

PT −1 1+c (i) Var(ST ) ≤ m T + 2 h=1 (T − h)ch ≤ m T 1−c . (ii) For any η > 0, ST m 1+c E[ST ] Pr >η ≤ · − −−−−→ 0, T T T η 2 1 − c T →∞ hence STT − E[STT ] → 0 in probability. (iii) For any ω ∈ (0, 1), r p mT 1 + c ST ≤ E[ST ]+ Var(ST )/ω ≤ E[ST ]+ · ω 1−c

with probability at least 1−ω.

Proof. By bilinearity of covariance, grouping terms by lag h = |i − j| yields Var(ST ) =

T X T X

Cov(Xi , Xj ) ≤ m

i=1 j=1

Since (T − h) ≤ T and

T X T X

T −1 X c|i−j| = m T + 2 (T − h)ch .

i=1 j=1

h=1

P∞

h h=1 c = c/(1 − c), we obtain

T +2

T −1 X

∞ X

h=1

(T − h)ch ≤ T + 2T

ch = T

1+c , 1−c

proving (i). Items (ii) and (iii) follow by Chebyshev’s inequality applied to ST /T and ST , respectively.

The next lemma bounds the covariance of suprema under the canonical Gaussian autoregressive coupling that underlies the TV-GP recursion. Lemma C.6 (Covariance contraction for Gaussian suprema).pLet f and g be centered, separable Gaussian processes indexed by a set T , with g(t) = ρf (t) + 1 − ρ2 h(t) for ρ ∈ [0, 1], where h is an independent copy of f . Assume supt∈TVar(f (t)) ≤ σT2 < ∞ and supt∈T |f (t)| < ∞ almost surely. Then Cov supt∈T f (t), supt∈T g(t) ≤ ρ σT2 . Proof. We first prove the claim for a finite set T = {1, . . . , k} and then pass to a general separable T . Let A = (f (1), . . . , f (k)) and C = (g(1), . . . , g(k)). By construction, the joint distribution of (A, C) is a zero-mean Gaussian with cross-covariance Cov(Ai , Cj ) = ρ Σij , where Σij = Cov(Ai , Aj ). Pk nui Let ϕn (u) := n1 log be the smooth approximation of the maximum function. The i=1 e gradient ∇ϕn (u) is the softmax function, ensuring that for any u, ∇ϕn (u) is a probability vector (∥∇ϕn (u)∥1 = 1 with strictly non-negative components). To compute the covariance Cov(ϕn (A), ϕn (C)), we employ the Gaussian interpolation technique standard in the proofs of Gaussian comparison inequalities [see, e.g., Adler and Taylor, 2009, Section 2.2]. Consider the interpolation path parameterized by t ∈ [0, 1] with the block covariance matrix Σ tρΣ Γ(t) = . tρΣ Σ Let (A(t), C(t)) ∼ N (0, Γ(t)). For t = 0, A(0) and C(0) are independent, implying Cov(ϕn (A(0)), ϕn (C(0))) = 0. For t = 1, we recover the joint distribution of (A, C). Differentiating the expectation of ϕn (A(t))ϕn (C(t)) with respect to t and applying multivariate Gaussian integration by parts [Adler and Taylor, 2009, Lemma 2.1.4] yields: k X k X d dΓAC,ij (t) ∂ϕn (A(t)) ∂ϕn (C(t)) E[ϕn (A(t))ϕn (C(t))] = E . dt dt ∂Ai ∂Cj i=1 j=1 19

Since

dΓAC,ij (t) = ρΣij , integrating from t = 0 to t = 1 establishes the exact identity: dt

Z 1 E[⟨∇ϕn (A(t)), Σ∇ϕn (C(t))⟩] dt.

Cov(ϕn (A), ϕn (C)) = ρ 0

We bound the bilinear form inside the expectation via ℓ1 /ℓ∞ duality. Let p = ∇ϕn (A(t)) and q = ∇ϕn (C(t)). Because p and q are probability vectors, X X pi Σij qj ≤ max |Σij | pi qj = max |Σij |. i,j

i,j

Because Σ is positive semi-definite, its maximum absolute entry lies on the diagonal: maxi,j |Σij | ≤ maxi Σii . By assumption, maxi Σii ≤ σT2 . Therefore, the bilinear form is bounded almost surely by σT2 , independent of the dimension k. Substituting this bound into the integral yields: Z 1 Cov(ϕn (A), ϕn (C)) ≤ ρ

E[σT2 ]dt = ρσT2 .

The approximations ϕn (A) and ϕn (C) converge almost surely to ϕ(A) and ϕ(C). The Borell–TIS inequality ensures that {ϕn (A)} and {ϕn (C)} are uniformly square-integrable. Vitali’s convergence theorem permits passing the limit through the covariance, which yields Cov(maxi Ai , maxi Ci ) ≤ ρσT2 . To extend this to a general separable set T , let {Tm }m≥1 be an increasing sequence of finite subsets whose union is dense in a separability set for f and g. This guarantees that supt∈Tm f (t) and supt∈Tm g(t) converge almost surely to their respective suprema over T . Applying the finite-set result yields Cov(supt∈Tm f (t), supt∈Tm g(t)) ≤ ρσT2 for all m. The Borell–TIS inequality provides supm E[(supt∈Tm f (t))2 ] < ∞. Uniform integrability allows taking the limit as m → ∞, which establishes the bound for the general separable set T . p √ In the TV-GP model, the transition rule is ft+h (·) = ρh ft (·) + 1 − ρ2h g̃t,h (·), where ρ := 1 − ϵ and g̃t,h is a centered Gaussian process with the same spatial kernel as ft . The process g̃t,h is independent of ft because it is a linear combination of the independent innovations {gt+1 , . . . , gt+h }. Setting T = S in Lemma C.6 gives σT2 ≤ 1 and Cov(Mt , Mt+h ) ≤ ρh for h ≥ 0. We apply Lemma C.5 by substituting Xt := Mt , m = 1, and c = ρ. For any η > 0, ! ! PT T T E[ t=1 Mt ] 1X 1X Pr Mt − E[M1 ] > η ≤ Pr Mt − > η → 0, T t=1 T t=1 T PT PT since E[ t=1 Mt ]/T = E[M1 ] by stationarity. Equivalently, the time average T1 t=1 ft (x⋆t ) converges in √ probability to E[supx∈S f1 (x)]. The deviation bound in Lemma C.5(iii) provides an explicit O( T ) concentration inequality, with constants depending on ϵ entirely through the factor (1 + ρ)/(1 − ρ). C.2

Span-only concentration fallback

This subsection justifies the span-only fallback bound used in Section 3.2. Recall that Wt = sup ft (x) − ft (y) . x,y∈S

Since rt ≤ Wt for every policy and every round, T X

RT ≤

Wt .

t=1

Define STW :=

T X t=1

Wt .

Equivalently, STW =

T X

sup

ft (xt ) − ft (yt ) .

(x1 ,y1 ),...,(xT ,yT )∈S×S t=1

Thus STW is the supremum of a centered Gaussian process indexed by (S × S)T . √ Let ρ = 1 − ϵ. For any fixed index sequence (x1 , y1 ), . . . , (xT , yT ), the variance of the corresponding Gaussian variable satisfies ! T T h i X X Var ft (xt ) − ft (yt ) = ρ|t−s| kS (xt , xs ) − kS (xt , ys ) − kS (yt , xs ) + kS (yt , ys ) t=1

t,s=1

≤4

T X

ρ|t−s|

t,s=1

1+ρ . 1−ρ Here we used kS (x, x) ≤ 1, which implies Var f (x) − f (y) ≤ 4 and bounds the absolute covariance between two such increments by 4. ≤ 4T

By Borell–TIS, for all u > 0, P STW − E[STW ] > u ≤ exp −

u2 . 8T (1 + ρ)/(1 − ρ) Stationarity gives E[STW ] = T µw . Therefore, with probability at least 1 − ω, r 1+ρ W log(1/ω). ST ≤ T µw + 8T 1−ρ Since RT ≤ STW , we obtain r 1+ρ log(1/ω) RT ≤ T µw + 8T 1−ρ with probability at least 1 − ω. C.3

Concentration of realized regret (Proof of Proposition 2)

Lemma A.1 establishes that the combined confidence and regularity event Et holds with an unconditional probability of P(Et ) ≥ 1 − δ. Because all instantaneous regret terms are non-negative, rt is bounded pathwise by the UCB width on the success event and the spatial span on the failure event: p rt ≤ 2 βσt−1 (xt ) + C2 + Wt IEtc , (47) where Wt = supx∈S ft (x) − inf x∈S ft (x). Summing over T and applying the Cauchy-Schwarz inequality to the standard deviations yields the realized cumulative regret: T X p RT ≤ C1 T βγ̃T + T C2 + Wt IEtc . (48) | {z } t=1 :=UT | {z } :=FT

Since γ̃T deterministically bounds the information gain for any sequence of length T , UT is a valid upper bound almost surely. To bound the stochastic failure penalty FT , we rely on the unconditional properties of the spatial span. By Lemma A.2, the expected product of the sub-Gaussian span and the failure indicator satisfies E[Wt IEtc ] ≤ δC3 . Because instantaneous regret is strictly non-negative, FT ≥ 0. Applying Markov’s inequality directly to the sum yields: T δC3 E[FT ] P FT > ≤ = ω. (49) ω (T δC3 )/ω Substituting this limit into (48) yields the bound in (15). 21

Kalman filtering implementation for the simulation study

This appendix describes the exact posterior updates used in Section 5. On the finite simulation grid {x(1) , . . . , x(N ) }, define ft = (ft (x(1) ), . . . , ft (x(N ) ))⊤ ∈ RN , and let KS ∈ RN ×N be the spatial covariance matrix with entries (KS )ij = kS (x(i) , x(j) ). Restricting the TV-GP recursion (4) to this grid gives the linear Gaussian state-space model

where ρ = xt .

√

ft = ρft−1 + ξt ,

ξt ∼ N (0, ϵKS ),

yt = Ht ft + vt ,

vt ∼ N (0, σ 2 ),

1 − ϵ and Ht ∈ R1×N selects the coordinate corresponding to the sampled grid point

µt|t−1 = ρµt−1|t−1 , St = Ht Pt|t−1 Ht⊤ + σ 2 ,

Lt = Pt|t−1 Ht⊤ St−1 ,

µt|t = µt|t−1 + Lt (yt − Ht µt|t−1 ),

Pt|t = Pt|t−1 − Lt Ht Pt|t−1 .

The acquisition rule is evaluated on the grid using the predictive mean and marginal variances, o n p q it ∈ arg max (µt|t−1 )i + β (Pt|t−1 )ii . i∈{1,...,N }

This recursion is equivalent to GP posterior inference under the spatio-temporal kernel (5) restricted to the grid, but avoids recomputing the posterior from the full observation history. Since N is fixed, each posterior update costs O(N 2 ) rather than growing with t, enabling exact simulations over long horizons.

Simulation reproducibility details

This appendix provides the reproducibility details for the simulation study in Section 5. The experiments use synthetic data generated directly from the TV-GP model in Equation (4). No external datasets, pretrained models, or human-subject data are used. The action space is the fixed 20 × 20 grid on S = [0, 1]2 , giving N = 400 candidate actions. The spatial kernel is the squared exponential kernel with lengthscale l = 0.1 and prior variance 1.0. The observation noise variance is fixed at σ 2 = 0.01. Each reported average regret value is computed over 100 independent TV-GP sample paths of length T = 10,000. For each drift rate ϵ and exploration parameter β, the latent process is simulated according to the recursion in Equation (4). Posterior inference is performed exactly on the finite grid using the Kalman filtering recursion in Appendix D. The constant-β TV-GP-UCB policy selects actions according to n o p q it ∈ arg max (µt|t−1 )i + β (Pt|t−1 )ii . i∈{1,...,N }

The logarithmic baseline uses the schedule βt = c1 log(c2 t),

c1 = 0.8,

c2 = 4,

as described in Section 5. The accompanying code archive contains the scripts used to generate the synthetic sample paths, run constant-β TV-GP-UCB, run the logarithmic time-varying-βt baseline, compute the average regret curves, bootstrap the empirical minimizers β̂opt , and regenerate Figure 1. For the left panel of Figure 1, shaded bands denote approximate 95% confidence intervals for the mean average regret, computed as ±1.96 standard errors over the 100 independent trials. For the right panel, the uncertainty intervals for β̂opt are computed using bootstrap percentile intervals over the independent trials. 22

The computations for Figure 1 were run on a MacBook Pro with an Apple M1 Pro chip, 8 CPU cores, and 16 GB of unified memory. The simulations used CPU multiprocessing and no GPU acceleration. The complete run used to generate the final figure took approximately 48 hours of wall-clock time. Memory usage was modest because the Kalman filter operates on a fixed N = 400 state dimension, but the total runtime is dominated by the large number of independent trials and parameter settings.

Sharper Regret Bounds for Time-Varying Gaussian Process Bandits with Constant Exploration

Analyze & explore related research in Conceptio

Related documents