Exploring and Exploiting Model Uncertainty in Bayesian Optimization

We thank Reviewer 44us for the valuable feedback. We have addressed all concerns.

• [W1] Clarify the paper’s positioning in related work for non-BO experts.

Thanks for the suggestion. To make the paper more accessible to non-experts, we will revise Section 2 to clearly highlight the two key components of BO, namely the acquisition policy and the surrogate model, and clarify that our contribution lies in proposing a new, more powerful surrogate. We will also revise the second paragraph of the Introduction and the paragraph of Related Work to better position our $\infty$-GP model as a replacement and enhancement of the standard GP, which serves as a fundamental backbone for BO and many uncertainty quantification tasks.

Here is a brief overview of our contribution. GP is the foundational surrogate model in BO. We propose an enhancement of this backbone: the $\infty$-GP model, which defines an infinite mixture of GPs while remaining tractable in practice. This enables modeling a much richer class of complex reward distributions. The intuition is similar to the universal approximation theorem: one can approximate a complex function using an infinite sum of simple base functions. At the same time, $\infty$-GP inherits the nice properties of standard GPs, such as closed-form updates and extensibility.

• [W2] Concern regarding overlapping confidence bars:

While our method’s confidence intervals partially overlap with those of some baselines, the overlap is minor, and our method consistently achieves significantly better mean performance. BTW, we believe this overlapping phenomenon is fairly common and expected in BO (e.g., as also observed in [1]). One reason for the overlap is that BO algorithms typically encourage exploration, which introduces more randomness and can lead to wider confidence intervals.

[1] Xu, Zhitong, et al. "Standard Gaussian Process is All You Need for High-Dimensional Bayesian Optimization." The Thirteenth International Conference on Learning Representations.

• [Q1] Question regarding the extension to high-dimensional BO:

Since the $\infty$-GP model is defined as a mixture over infinitely many GPs, each layer corresponds to a standard GP. By modifying the base GP in each layer, our framework can seamlessly incorporate any existing GP extensions, including those designed for high-dimensional settings. For example, sparse GPs, Gaussian Markov random fields, additive GPs, GP with low-dimensional subspace embeddings, and high-dimensional BO methods utilizing deep learning.

Indeed, our experiment on the prompt optimization task already provides some empirical evidence on how the $\infty$-GP framework can be extended to high-dimensional settings. In this task, we use UMAP for dimensionality reduction, which can be viewed as a low-dimensional embedding-based BO method, and observe strong empirical performance under this setting.

• [Q2] Clarify how the proposed setting differs from mis-specified BO:

Thanks for the question. Existing mis-specified BO methods typically handle limited types of model mismatch, such as incorrect smoothness assumptions, within a fixed model class like GPs. In contrast, our $\infty$-GP model places a prior over the entire space of reward distributions, allowing us to capture model uncertainty across different function classes and providing much greater flexibility.

• [Q3] Clarify whether Eq. (5) simply expands the model class and how expressive it is:

Thanks for the question. Your understanding of Eq.(5) is essentially correct, but more precisely, we are not simply expanding the model class; instead, we let the reward model (which is a distribution) itself be random and place a prior over the space of reward distributions. Regarding “to what extent is this broader class able to capture arbitrary reward models,” the discussion following Assumption 1 rigorously characterizes the class of reward distributions that can be approximated by our $\infty$-GP model, and shows that it is significantly broader than that of standard GP model. We also note that our paper does not require the reader to have prior knowledge of SDP, as we clearly define the SDP and provide intuitive explanations to aid understanding.

[W3]

The comment in Weakness 3 appears incomplete. We would appreciate it if the reviewer could kindly clarify or complete the sentence, so we can provide a proper and constructive response.

We sincerely thank the reviewers for their thoughtful comments. We have carefully addressed all the questions and concerns raised.

• [W1] Regarding the computational complexity of Gibbs sampling and the choice of truncation number $L$

First, we use a truncation technique in Gibbs sampling, i.e. setting the maximum number of mixture component to fixed $L$. By doing so, the complexity of truncated Gibbs sampling is $\mathcal{O}(n^3)$, which is exactly equal to traditional GP (for more details, please refer to our response to W5&Q3).

Second, despite using a finite truncation, we must note that we are not "contradicting the motivation to use an infinite mixture of GP surfaces", and our model is still theoretically an infinite mixture of GPs. This is because, as shown in Eq. (10), when we actually perform posterior prediction using the $\infty$-GP, we do not use infinitely many surfaces, but rather $K_n+1$ surfaces. The truncation level $L$ essentially imposes an upper bound on $K_n$, and in our experiments, we set $L=8$. In practice, this upper bound is rarely reached, and thus the performance is largely insensitive to the choice of $L$. This observation is supported both theoretically and empirically. Theoretically, as shown in Equation (2.6) of [1], the number of active surfaces $K_n$ in a Dirichlet Process grows very slowly with $n$, at a rate of $\log n$. Empirically, the table below reports the best found values (optimal = 0) for the Ackley-HT-NS task under different values of $L$. The performance is not sensitive to $L$. This is because, by the 100th iteration, we observe that only 2–3 surfaces have non-negligible weights in the posterior distribution. The remaining layers are essentially unused, so When $L$ is 6 or greater, changing it has little impact on the result.

[1] Bayesian nonparametrics[M]. Cambridge University Press, 2010.

• [W2 & Q1 & Q2] Prior specification and sensitivity to hyperparameter choices.

Thanks for the questions. The prior specification involves two parts: (i) choosing the prior distribution families, and (ii) setting the hyperparameters for those priors. We explain the two parts one by one below.

(i) We use conjugate priors for $\sigma^2$, $\tau$ and $\nu$ to allow efficient Gibbs sampling with analytical posterior updates. These choices are quite standard and make the inference much easier in practice. For $\phi$, we follow the recommendation from [2] and use a uniform prior, which is shown in their paper to work well and be numerically stable.

(ii) The prior hyperparameters include $a_\nu, b_\nu, a_\tau, b_\tau, a_\sigma, b_\sigma, b_\phi$. We did not include the precision prior hyperparameter $(a_\nu, b_\nu)$ in the appendix, but in our experiments, we simply set it to $(1,1)$. The settings for the other hyperparameters are already listed. We note that, except for $b_\phi$, the results are not sensitive to these prior hyperparameters. Please refer to our response to [W1 & Q2] of Reviewer v6ga for a detailed sensitivity analysis of the prior hyperparameters, where we provide both empirical results and theoretical justifications. We also provide guidance on how to choose the $b_\phi$.

[2] Gelfand A E, Kottas A, MacEachern S N. Bayesian nonparametric spatial modeling with Dirichlet process mixing[J]. Journal of the American Statistical Association, 2005, 100(471): 1021-1035.

• [W3 & W4 & Q5 & Q6] Comparison to additive GP and non-stationary kernel baselines

Thank you for these insightful suggestions. We will include relevant comparisons to BO methods that exploit additive structure and use explicitly non-stationary kernels in the revised version. We would like to clarify that the proposed $\infty$-GP is fundamentally different from additive GP models. Additive BO are designed for high-dimensional problems by partitioning the input dimensions and assigning a single GP surrogate to each subset of dimensions. The resulting model is still a single GP but with an additive kernel. In contrast, our $\infty$-GP framework models the entire input space globally as a mixture of infinitely many GP surfaces. The resulting model is no longer a GP, and can represent reward models that lie far beyond the GP function class. By the way, additive BO can be seamlessly integrated into the $\infty$-GP framework by replacing the standard GP in each surface with an additive GP.

As recommended, we will include the following experimental results comparing our method with TS using additive GP in [4] and non-stationary spectral kernel GPs in [3]. Both alternatives underperform our method. This is because both baselines still operate within the single-GP framework, which is limited in its ability to model complex rewards.

[3] Remes, Sami, Markus Heinonen, and Samuel Kaski. "Non-stationary spectral kernels." Advances in neural information processing systems 30 (2017).

[4] Gardner, Jacob, Chuan Guo, Kilian Weinberger, Roman Garnett, and Roger Grosse. "Discovering and exploiting additive structure for Bayesian optimisation." In Artificial Intelligence and Statistics, pp. 1311-1319. PMLR, 2017.

• [W5 & Q4] Regarding regret analysis and time complexity analysis

Thanks for your helpful comment. We will include regret and time complexity analysis in the updated version of the paper.

• Regret: The $\infty$-GP-TS indeed maintains sub-linear regret. In fact, we only need one additional step from our established Theorem 1 to derive a regret bound. The instantaneous regret at step $n$ satisfies: $r_n\leq\int |\mu^\ast(x)||\pi_{f_n}(x)-\pi_{f^\ast}(x)|dx\leq2\mu^\ast(x_\ast)\delta(\pi_{f_n},\pi_{f^\ast})\leq2\mu^\ast(x_\ast)\delta(f_n,f^\ast)$, where $\delta$ is total variation distance (TVD), $\pi_{f_n}$ and $\pi_{f^\ast}$ is the TS policy distribution induced by surrogate posterior $f_n$ (defined in Eq.(8)) and true reward $f^\ast$, respectively. The last inequality follows from the data processing inequality. The convergence $\delta(f_n,f^\ast)\to 0$ has already been established in Theorem 1 (noting that L1 norm= 2TVD). Therefore, we obtain $r_n\to 0$, implying sublinear cumulative regret. Moreover, the concrete regret rate can be made explicit by leveraging the posterior contraction result in Theorem 1. Regarding the question about "the key theoretical differences due to operating in the distributional space," we would like to clarify the following: First, our model indeed defines a prior over the distributional space (i.e., $G$ is drawn from an SDP). However, as shown in Eq. (13), we marginalize out the random distribution $G$ over the SDP when computing the posterior predictive distribution. As such, the actual $\infty$-GP-TS algorithm does not explicitly operate over the distributional space. Second, the theoretical difference is as follows: Traditional GP-TS analyzes regret by studying how the posterior predictive mean concentrates around the true function. In contrast, in $\infty$-GP-TS, as shown in Theorem 1 and the above analysis, we bound the regret in terms of the contraction rate of the posterior over the distributional space toward the true reward distribution.

• Time Complexity: By using a truncation technique (see Algorithm 2), the complexity of our method is $\mathcal{O}(n^3)$, which is equal to that of strandard GP-based BO. In short，the complexity is governed by two main components. (i) The computational complexity of Kriging, which is mainly govern by computing the inverse $\Sigma_0(x_{1:n}, x_{1:n})^{-1}$, incurring $\mathcal{O}(n^3)$ cost. However, the $\infty$-GP model enjoys the advantage that all surfaces share the same $\Sigma_0(x_{1:n}, x_{1:n})^{-1}$, allowing this expensive matrix inversion to be computed only once. (ii) The complexity of Gibbs sampling. By truncating the maximum number of surface to $L$, one can sample the $\xi_l(x_{1:n})$ jointly and thus the complexity is $\mathcal{O}(n^3)$. In summary, the total complexity of $\infty$-GP is $\mathcal{O}(n^3)$. The time complexity analysis is supported by the following wall-clock time results (in seconds).

• [Q3] Question about metric that quantifies the reuse of latent surfaces

Thank you for this insightful question. The reuse of latent surfaces can indeed be quantified with the following quantities. First, $n_j$ quantifies how many times the surface $j$ has been reused in past iterations. Second, the posterior of $\nu$ characterizes the decision maker’s model uncertainty and governs the trade-off between introducing new surfaces (exploration) versus reusing existing ones (exploitation). Smaller values of $\nu$ encourage more reuse. Third, the ratio $\frac{{n}_j}{\nu + n}$ quantifies the probability of reusing surface $j$ in the upcoming new iteration.

• [Q7] Minor questions

The y-axis in Figure 2 denote the original value of the test function. In our experiments, all methods—including both the baseline BO approaches and our $\infty$-GP-TS use sklearn.MinMaxScaler to linearly map observed rewards into the fixed interval $[-5,5]$ before fitting the surrogate. It is a common practice for BO to ensure numerical stability. And all methods use the same affine transform to ensure fair comparisons. After fitting the surrogate, we use the inverse transform to show the final result in the paper. The dimensions for the synthetic function is 20.

I have read the rebuttal and appreciate the authors’ efforts in addressing all my concerns. The additional insights, particularly the discussion on regret bounds and time complexity, strengthen the paper and could enhance its acceptance within the community. While the authors have provided their rationale for choosing the prior distribution families, I remain somewhat unconvinced about the practical feasibility of this approach. Nevertheless, I will retain my score as "Borderline Accept."

We sincerely thank the reviewers for their thoughtful comments. We have carefully addressed all the questions and concerns raised.

• [W1 & Q2] Sensitivity analysis to hyperparameter choices in the prior

Thank you for the helpful comments. We have now included a sensitivity analysis to examine how the hyperparameters in the prior influence the performance. The prior hyperparameters include $a_\nu, b_\nu,a_\tau, b_\tau,a_\sigma,b_\sigma,b_\phi$, and their settings are already provided in the appendix. We note that, except for $b_\phi$, the results are not sensitive to these prior hyperparameters. See the sensitivity analysis table below, where we report the best-found objective values (optimal value=0) under different hyperparameter configurations on Ackley-HT-NS.

There are two main reasons why the performance is insensitive to the choice of $a_\nu,b_\nu,a_\tau, b_\tau,a_\sigma,b_\sigma$: First, as shown in Lines 19, 21, and 22 of Algorithm 2, the posteriors of $\nu$, $\sigma$, and $\tau$ all follow distributions with parameters of the form $(a_\ast + \mathcal{O}(n),\ b_\ast + \mathcal{O}(n))$, where $a_\ast$ and $b_\ast$ refer to $a_\nu, b_\nu, a_\tau, b_\tau, a_\sigma, b_\sigma$. The only exception is the second parameter in the posterior of $\nu$, but it can still be written as $b_\nu$ plus a very large term $-\log(w_L)$. $a_\ast$ and $b_\ast$ are typically set within the range 0–3. Therefore, as $n$ increases, the influence of these prior hyperparameters $(a_\ast,b_\ast)$ becomes negligible. They mainly act as initialization values and quickly fade in importance once $n$ is moderately large. Second, using TS policy, each decision step is based on only one sample from the posterior, which introduces significant randomness. This further reduces the influence of the initial hyperparameter values on the overall performance.

But the choice of $b_\phi$ is critical, as it sets the range over which $\phi$ can vary. We have already included guidance on setting $b_\phi$ in the appendix. In general, this value is set to be quite large (since we use a very small $d$, e.g., 0.01), which ensures a sufficiently wide range for $\phi$ to explore. In practice, we also use an adaptive strategy: we start with a very large $b_\phi$, run a few initialization steps to observe the typical scale of $\phi$, and then adjust $b_\phi$ accordingly. This approach helps set $b_\phi$ in a task-specific way.

• [W2&Q1] Regarding time complexity analysis.

Thanks for your helpful comments. We will include time complexity analysis in the updated version of this paper, and show that, by using a $L$-truncation technique (see Algorithm 2), the complexity of our method is $\mathcal{O}(n^3)$, which is equal to that of strandard GP-based BO.

In short，the complexity is governed by two main components. (i) The computational complexity of Kriging, which is mainly govern by computing the inverse $\Sigma_0(x_{1:n}, x_{1:n})^{-1}$, incurring $\mathcal{O}(n^3)$ cost. However, the $\infty$-GP model enjoys the advantage that all surfaces share the same $\Sigma_0(x_{1:n}, x_{1:n})^{-1}$, allowing this expensive matrix inversion to be computed only once. (ii) The complexity of Gibbs sampling. By truncating the maximum number of surface to $L$, one can sample the $\xi_l(x_{1:n})$ jointly and thus the complexity is $\mathcal{O}(n^3)$. In summary, the total complexity of $\infty$-GP is $\mathcal{O}(n^3)$.

The computational efficiency is supported by the following wall-clock time results (in seconds). Its runtime is comparable to that of GP-based baselines. This is because, although $\infty$-GP requires expensive Gibbs sampling, the baseline methods also involve hyperparameter inference either through fully Bayesian MCMC or through MLE that relies on gradient-based optimization, both of which incur comparable computational complexity. Notably, $\infty$-GP-TS performs only a single Gibbs sampling run in each optimization iteration, which significantly reduces its overall computational burden.

• [Q3] Theoretical contribution compared with existing works.

We would like to clarify our theoretical contribution. In short, while the convergence of classic GP-based methods typically relies on restrictive and often opaque assumptions, our results substantially relax these conditions. To provide a clearer comparison, we rewrite the observed reward $y(x)$ as a deterministic component $\mu^\ast(x)$ plus a noise term $\epsilon^\ast(x)$.

• In the classic GP literature, $\mu^\ast$ is assumed either to be a sample path drawn from a GP, or to lie in a RKHS associate with certain kernel and have a known RHKS norm bound (e.g. [1,2]), which is difficult to compute in practice. These assumptions are rarely verifiable in practice. For the noise term $\epsilon^\ast(x)$, to utilize concentration results, it is assumed to be independent and (sub-)Gaussian, or even noiseless in some formulations.

• In our result, regarding the noise $\epsilon^\ast(x)$, our framework accommodates a significantly broader class of distributions. For example, Assumption 1 is satisfied by noise of the form $\sum_{l=1}^L w_l(x)\psi_l$, where each $\psi_l$ belongs to distribution families such as Weibull (with shape $k > 1$), Gaussian, Laplace, Gamma, or exponential. Regarding the $\mu^\ast$, our result in Theorem 1 does not require any RKHS assumptions or GP sample path assumption as in the classic GP literature. But we note that, different from existing literature, Theorem 1 assumes a finite decision set. To handle the continuous case, a standard discretization argument (as in [2]) is required. Consequently, we need an additional assumption that $\mu^*(x)$ is Lipschitz continuous, which is still much weaker than the assumptions typically made in the GP literature (those assumptions themselves imply Lipschitz continuity).

[1] Russo D, Van Roy B. Learning to optimize via posterior sampling[J]. Mathematics of Operations Research, 2014, 39(4): 1221-1243.

[2] Srinivas N, Krause A, Kakade S M, et al. Gaussian process optimization in the bandit setting: No regret and experimental design[J]. arXiv preprint arXiv:0912.3995, 2009.

• [W3] Clarification of performance metrics and plot axes

The y-axis in Figure 2 and 3 denote the 'best-found' objective value. All the synthetic benchmarks show the original test function value in the y-axis. The y-axis for the portfolio task denote the portfolio return. The y-axis on both HPOBench and NASBench denote the accuracy on the validation set.

We sincerely thank the reviewers for their thoughtful comments and constructive feedback. We have carefully addressed all the questions and concerns raised, and have made corresponding revisions or clarifications as detailed below.

• [W1]:Missing discussion on existing literature on model uncertainty in hyperparameter learning

Thank you. We will revise the last paragraph of the “Related Work” section to explicitly contrast our $\infty$-GP framework with prior work on model uncertainty in hyperparameter learning and related areas, including Hvarfner et al. (2024) and other relevant studies such as [1,2]. Specifically, we will clarify that those approaches manage model uncertainty within a fixed model class (e.g., GPs with uncertain kernel design), whereas our $\infty$-GP addresses uncertainty over model classes by placing a prior over the distributional space of reward models, offering greater modeling flexibility.

[1] Malkomes G, Garnett R. Automating Bayesian optimization with Bayesian optimization[J]. Advances in Neural Information Processing Systems, 2018, 31.

[2] Malkomes G, Schaff C, Garnett R. Bayesian optimization for automated model selection[J]. Advances in neural information processing systems, 2016, 29.

• [W2] Missing comparison with state-of-the-art industrially relevant methods such as HEBO.

Thanks for your helpful comments. In response, we have included HEBO as an additional baseline in our updated experiments (see the table below, which reports the best-found performance after 100 iterations). Our method consistently outperforms HEBO in Portfolio task, HPOBench and NASBench.

First, while HEBO incorporates input and output warping to address non-stationarity (in fact, in our experiments, all baselines for nonstationary tasks already adopt input warping techniques by default), it still relies on a single GP model as the surrogate. As such, the model class that a single GP can represent remains limited, and it cannot handle heavy-tailed noises. Our results show that warping alone is insufficient. Second, HEBO ensembles multiple base acquisition functions (e.g., EI, UCB) to improve search efficiency. However, as our results demonstrate, in highly complex reward landscapes involving model misspecification, the choice of a more expressive surrogate model (such as $\infty$-GP) is often more critical than acquisition efficiency. Importantly, the complex structure of our $\infty$-GP surrogate is learned directly from data, rather than hand-crafted through warping techniques.

• [W3] Lack of Regret Analysis

Thanks for your helpful comment. We will include regret analysis in the updated version of the paper. In fact, we only need one additional step from our established Theorem 1 to derive a regret bound. The instantaneous regret at step $n$ satisfies: $r_n\leq\int |\mu^\ast(x)||\pi_{f_n}(x)-\pi_{f^\ast}(x)|dx\leq2\mu^\ast(x_\ast)\delta(\pi_{f_n},\pi_{f^\ast})\leq2\mu^\ast(x_\ast)\delta(f_n,f^\ast),$ where $\delta$ is total variation distance (TVD), $\pi_{f_n}$ and $\pi_{f^\ast}$ is the TS policy distribution induced by surrogate posterior $f_n=f(y(x_n)|\mathcal{H}_n)$ (defined in Eq.(8)) and true reward $f^\ast$, respectively. The last inequality follows from the data processing inequality. The convergence $\delta(f_n,f^\ast)\to 0$ has already been established in Theorem 1 (noting that L1 norm= 2TVD). Therefore, we obtain $r_n\to 0$, implying so-called no-regret or sublinear cumulative regret.

Moreover, the concrete regret rate can be made explicit by leveraging the posterior contraction result in Theorem 1, requiring only a standard discretization argument (e.g., from [3]) to handle the continuous decision space.

[3] Srinivas N, Krause A, Kakade S M, et al. Gaussian process optimization in the bandit setting: No regret and experimental design[J]. arXiv preprint arXiv:0912.3995, 2009.

• [W4&Q3] Regarding Experimental Details:

(1)HPObench details: We would like to clarify that we used the "kc1" dataset for the MLP task, which is widely adopted in the AutoML and HPO benchmark community. The OpenML task ID for kc1 is 3917. This dataset consists of 2,109 training instances and 22 features (after preprocessing), and it originates from a real-world software defect prediction task.

(2)Prompt optimization details: We apply UMAP to compress high-dimensional prompt embeddings into a 10-dimensional manifold-preserving latent space. The specific UMAP configuration is as follows: we set the number of components to 10 (i.e., n_components=10), and the number of neighbors to 15 (n_neighbors=15). We also used the 'cosine' distance metric, as it is more suitable for semantic embeddings than the default Euclidean metric. The minimum distance parameter was set to 0.0 to encourage tighter clusters in the projected space. Other dimensionality reduction techniques such as PCA and t-SNE were considered, but we chose UMAP because it preserves both local and global structure, is computationally efficient for large text collections, and allows flexible adjustment of clustering granularity via the n_neighbors parameter.

Regarding the comment on “missing comparisons to domain-specific BO baselines,” to the best of our knowledge, there is no established class of BO algorithms specifically tailored for prompt optimization that fundamentally differs from general-purpose BO methods. Notable works, including [4] and [5], all rely on standard BO algorithms such as GP-UCB or GP-EI. Compared with our implementation, their main difference lies in the way they embed or reduce the dimensionality of the prompt space. To ensure a fair comparison focused on the BO algorithm itself rather than the embedding method, we adopt UMAP across all methods.

(3) Hyperparameter setting: In our experiments, we set the number of MCMC iterations to B = 500 and the truncation number of mixture components to L = 8.

[4] Chen L, Chen J, Goldstein T, et al. Instructzero: Efficient instruction optimization for black-box large language models[J]. arXiv preprint arXiv:2306.03082, 2023.

[5] Sabbatella A, Ponti A, Giordani I, et al. Prompt optimization in large language models[J]. Mathematics, 2024, 12(6): 929.

• [Q1] Clarify difference between $\infty$-GP and spectral mixture kernel:

The spectral mixture kernel (SMK) enhances the expressiveness of a single GP by constructing flexible kernels through spectral mixtures. However, these methods still assume the reward process lies within the GP family, thus limiting their modeling flexibility. In contrast, $\infty$-GP places a prior over the space of reward distributions, allowing us to model rewards that are not themselves GPs. This enables model uncertainty over function classes, rather than within a fixed class. Compared to using a highly flexible kernel within a standard GP (e.g., SMK), the $\infty$-GP offers several distinct advantages.

• First, it provides greater modeling flexibility and accommodates not only non-stationary rewards, but also heavy-tailed, heteroskedastic, and other non-Gaussian noise structures, as supported by our theoretical guarantee in Theorem 1.
• Second, SMK can be directly incorporated into $\infty$-GP framework by replacing the base kernel. Since $\infty$-GP is a mixture over infinitely many GPs, it is compatible with any GP extensions, including custom kernel designs and other architectural improvements.
• Third, unlike standard GPs, $\infty$-GP quantifies not only value uncertainty but also model uncertainty.

• [Q2] question regarding the role of the noise term

Let us further explain the role of the noise term $\epsilon$. As defined in Eq. (6), $\xi$ is discretely distributed (though with infinite support), so it cannot directly model a continuous reward distribution. By mixing $\xi$ with a Gaussian $\epsilon$, the resulting distribution of $y(x)=\xi(x)+\epsilon(x)$ becomes continuous. Here, the term $\epsilon$ does not mean that we impose a Gaussian assumption on the true measurement noise, as we directly model the overall distribution of the observation $y$ via Eq. (4), rather than modeling the expectation and measurement noise separately as in classic GP-based BO. Therefore, we say the term $\epsilon$ is introduced solely for “modeling purposes”, without real physical meaning.

Dear Reviewer BoMR

Thank you once again for your valuable comments on our submission. As the discussion phase is approaching its end, we would like to kindly confirm whether we have sufficiently addressed all of your concerns (or at least part of them). Should there be any remaining questions or areas requiring further clarification, please do not hesitate to let us know. If you are satisfied with our responses, we would greatly appreciate your consideration in adjusting the evaluation scores accordingly.

We sincerely look forward to your feedback.

Dear Authors, Thank you for your response and for providing the additional experiments and theoretical insights. I appreciate the effort. I have two remaining points that, if addressed, would make me happy to raise my score.

• Regret Bound Derivation: I find the statement from the instantaneous regret to the final sublinear regret bound unclear. Could you please provide a more detailed derivation or explanation showing how the instantaneous regret convergence rate is sufficient to construct the overall sublinear cumulative regret bound? Strictly speaking, instantaneous regret converging to zero is a necessary condition for the sublinear cumulative regret. For example $\lim\frac{T}{2^t}\rightarrow 0$, while $\lim \frac{\sum \frac{T}{2^t}}{T} > 0$

• Statistical Significance: For the empirical comparison involving HEBO and $\infty$-GP-TS, could you please report the standard error over the multiple runs? Alternatively, the results of a formal significance test would also be sufficient.