Optimizing the Unknown: Black Box Bayesian Optimization with Energy-Based Model and Reinforcement Learning
Authors: RUIYAO MIAO
Junren Xiao
Shiya Tsang
Hui Xiong
Ying Nian Wu
Venue: NeurIPS 2025 poster ()
Keywords: Black Box Optimization
Bayesian Optimization
Energy-Based Model
Reinforcement Learning
Primary area: optimization

Abstract: Existing Bayesian Optimization (BO) methods typically balance exploration and exploitation to optimize costly objective functions. However, these methods often suffer from a significant one-step bias, which may lead to convergence towards local optima and poor performance in complex or high-dimensional tasks. Recently, Black-Box Optimization (BBO) has achieved success across various scientific and engineering domains, particularly when function evaluations are costly and gradients are unavailable. Motivated by this, we propose the Reinforced Energy-Based Model for Bayesian Optimization (REBMBO), which integrates Gaussian Processes (GP) for local guidance with an Energy-Based Model (EBM) to capture global structural information. Notably, we define each Bayesian Optimization iteration as a Markov Decision Process (MDP) and use Proximal Policy Optimization (PPO) for adaptive multi-step lookahead, dynamically adjusting the depth and direction of exploration to effectively overcome the limitations of traditional BO methods. We conduct extensive experiments on synthetic and real-world benchmarks, confirming the superior performance of REBMBO. Additional analyses across various GP configurations further highlight its adaptability and robustness. Our code is publicly available at: https://github.com/ruiyaoMiao0809/Black-Box-Bayesian-Optimization-with-Energy-Based-Model-and-Reinforcement-Learning

Optimizing the Unknown: Black Box Bayesian Optimization with Energy-Based Model and Reinforcement Learning
Ruiyao Miao', Junren Xiao?, Shiya Tsang", Hui Xiong Yingnian Wu
Unniversity California, Los Angeles
The Hong Kong University of Science and Technology (Guangzhou) The Hong Kong University of Science and Technology
ruiya008090g, ucla. edu, {jxiao767, szeng785]@connect hkust-gz. edu. cn xi xionghui Brnghhi@hkust--ae edw........dddd................. edu. cn, ywuêstat. ucla. edu
Abstract
Existing Bayesian Optimization (BO) methods typically balance exploration and exploitation optimize costly objective functions. However, these methods often suffer from a significant one-step bias, which may lead to convergence towards lo- cal optima and poor performance in complex high- high-dimensional tasks. Recently, Black-Box Optimization (BBO) has achieved success across various scientific and engineering domains, particularly when function evaluations costly and gradients are unavailable. Motivated by this, we propose the Reinforced Energy- Based Model for Bayesian Optimization (REBMBO), which integrates Gaussian Processes (GP) for local guidance with an Energy-Based Model (EBM) to capture global structural information. Notably, we define each Bayesian Optimization iteration as Markov Decision Process (MDP) and use Proximal Policy Optimiza- tion (PPO) for adaptive multi-step lookahead, dynamically adjusting the depth and direction of exploration to effectively overcome the limitations of traditional BO methods. We conduct extensive experiments on synthetic and real-world bench- marks, confirming the superior performance of REBMBO. Additional analyses across various GP configurations further highlight its adaptability and robustness. Our code is publicly available atchttps:/githhu.. /tttppss//gitub...o/riiiaaalll/ hamm//iiyaoMMao0809999883 Box- Payesian-Optimizze with -arrggy-BBed---ae wandd-einforremmnt-Leeeee
Introduction
Black-box optimization (BBO) is crucial for solving complex scientific and engineering problems when gradient information is unavailable or function evaluations are expensive [[] In practice, BBO approaches are widely applied to hyper-parameter tuning machine learning, materials discovery, drug formulation, and industrial process optimization, where each evaluation often incurs costly simulations or physical trials. Bayesian Optimization (BO) is prominent BBO technique that builds probabilistic surrogate (e.g.. Gaussian Process [2]) and an acquisition function guide new sample queries, thereby balancing exploration and exploitation
principled manner. However, standard GP-based BO can suffer from "one-step myopia," focusing short-term predicted gains at the expense more thorough exploration, a limitation that becomes especially pronounced in high-dimnn tigh-dimensional or modti environments.
Existing Bayesian Optimization (BO) methods primarily aim at efficiently locating optimal solutions by carefully balancing exploration and exploitation 3 Common strategies for handling complex
CCoresponding author.
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

optimization problems include dimensionality reduction methods like REMBO 4 or local partition- ing techniques such as TuRBO [5]. While effective in relatively simple settings, these approaches tend to converge prematurely to local optima when applied to complex, high-dimensional prob- lems. [6] []. To solve this limitation, recent techniques incorporate resource-intensive multi-step look-ahead schemes, including 2-step Expected Improvement (EI) 77, Knowledge Gradient (KG) [8], and reinforcement learning-dine -erning-driven methods like EARL-BO ( However, such methods typically demand significant computational resources yet still fail to achieve effective global exploration in challenging environments.
In this study, we introduce the Reinforeed Energy- Energy-Based Bayesian Optimization Model (REBMBO), depicted in Figure ! which addresses traditional shortfalls by combining an Energy- Based Model (EBM) with multi-step Reinforcement Learning (RL). Our novel EBM-UCB acquisition function integrates Gaussian Process local uncertainty estimates with global signals derived from neural network-based energy landscape learned via short-run MCMC, thereby guiding exploration away from less promising regions. In addition, we treat each Bayesian Optimization iteration as a Markov Decision Process (MDP) and employ Proximal Policy Optimization (PPO) for adaptive multi-step lookahead, thus dynamically adjusting exploration depth and direction to enhance robustness. To capture both local and global exploration objectives, we propose theoretically justified Landscape- Aware Regret (LAR) metric that incorporates global exploration penalties, offering a more holistic assessment of performance in complex optimization scenarios.
This research offers the following key contributions:
(1) We incorporate Energy-Based signals into a UCB-style acquisition function in the GP surrogate to capture diverse behaviors. This synergy between global exploration and precise local modeling addresses the limitations of single-step acquisition approaches.
(2) Bayesian Optimization is modeled as Markov Decision Process (MDP) with Proximal Policy Optimization (PPO). This technique addresses the one-step myopia typical GP-based strategies by adaptively balancing exploration and exploitation via multi-step lookahead.
(3) We introduce theoretically justified Landscape-Awar Regret (LAR) metric that extends standard regret with an energy-informed global term. This metric provides fair and balanced evaluation by jointly reflecting local exploitation and global exploration efficiency in complex optimization landscapes.
Experimental results, summarized in Figure 2 indicate that REBMBO reduces final Landscape- Aware Regret (LAR) and improves overall performance scores compared to state-of-the-art methods, consistently outperforming both single-step and short-horizon lookahead approaches even under challenging high-dimensional conditions. The subsequent sections detail our methodology and empirical findings, highlighting REBMBO's s robust and efficient performance.
2 Related Work
2.1 BO Background and Shortcomings
Black-box optimization (BBO) is central for tasks with costly or noisy function evaluations, commonly
handled by Bayesian Optimization (BO) frameworks (3). However, high-dimensional or discrete
domains often overwhelm classical Gaussian Process (GP) surrogates, prompting techniques such
as ARD-based variable selection [TO, REMBO 4. additive models (.. or local partitioning
in TuRBO 5
Discrete or combinatorial BO further adopts specialized surrogates (VAR [12..
COMBO 113 TPE [14. SMAC [5]) yet generally remains single-step. Likewise, robust or
constrained variants [16
17
18
19
20), multi-objective [21 22 23), transfcr/multi-fideiin [24] 25
26 27), and parallel 28
29
30 approaches usually retain one-step acquisitions. VrR (Truncated
g unifies BO and level-set estimation with strong guarantees under pointwise costs or heteroscedastic noise, while look-ahead or rollout-based based schemes 77 often incur high computational overhead. Recent work leverages MLMC for nested integrals [3] or formulates BO as an MDP under transition constraints 332 and methods like GLASSES (33 approximate multi-step losses via forward simulation. Yet, many remain domain-specifee or lack synergy with short-run MCMC. Existing RI. integrations
Variance Reduction)
99] 18 also typically rely on local posteriors, leaving open the
challenge of thorough multi-step exploration across multi-modal landscapes.

---

2.2 Baseline Targeting Global Optima and Limitations
In this paper, we compare against six common baselines that represent key paradigms in Bayesian Optimization, Classic BO 3 is a canonical single-step GP-based approach. BALLET-ICI 6 alternates global and local GPs but remains relatively myopic on multi-modal tasks. TuRBO specializes in local trust-region expansions yet lacks far-reaching jumps. EARL-BO [9] is an RI- based multi-step method, heavily dependent on local GP precision. addition, we include 2-step EI [7] and KG (Knowledge Gradient) [8] as two well- look-ahead techniques, though they tend to be limited to short horizons or incur high computational overhead. These baselines respectively illustrate single-step local search, partially global scanning, short-horizon non-myopia. but none combines global signals with adaptive multi-step planning in unified manner. By contrast, REBMBO employs short-run MCMC-traing local accuracy, and a PPO-based multi-step RL for planning. This synergy overcomes the one-step constraints in Classic BO, enables deeper exploration than BALLET-ICI or TuRBO, provides more robust coverage than EARL-BO, and avoids the excessive rollout overhead observed in 2-step El or KG. As detailed in Section 2.3 REBMBO leverages these three modules to handle high-dimensional tasks within limited budgets, offering a global and multi-step perspective.
Energy-Baa Based Model for global exploration, GP surrogate for
2.3 RL in BO and Energy-Driven Multi-Step Planning
Recent attempts to integrate reinforcement learning into Bayesian Optimization have enabled multi- step acquisitions but frequently rely on localized kernels or omit global exploration cues, leading to suboptimal performance complex tasks (9] 18]. For instance, EARL-BO shows the benefits of multi-step planning in high-dimensional settings but lacks explicit energy-based signals for broader coverage [9] However, REBMBO framework formulates each BO iteration as a MDP solved via Proximal Policy Optimization. This design alleviates one-step myopia and combines local GP fidelity with iterative RL lookahead under strict evaluation budgets.
Preliminaries
Online Black-Box Optimization (BBO). We consider continuous function f(x) defined over X C IRd, with the objective:
arg max fxx), xEX
under a strict evaluation budget. Each evaluation of f(x) can be computationally or financially expensive (34), thus data efficiency is crucial. Unlike offline methods that rely fixed sampling designs, online BO adaptively selects X, based on previously observed data, facilitating faster discovery of optimal regions.
Bayesian Optimization (BO) and Gaussian Processes (GP). BO maximizes maxxex f(x) by employing a Gaussian Process (GP) prior: [(x) ~ GP(m(x), k(x, x')), where typically m(x) = and the kernel function xx) (e.g., RBF or Matérn) encodes assumptions about the function's
smoothness [35]. Assuming noisy observations = f(x:) + Es. with E, ~ _N0 o) the GP posterior is given by:
f(x) Dx ~ (x)xxx,, )2*((x)),
with observed data D; {xx; } Although GPs capture local uncertainty, their inherent locality often restricts global exploration, making them prone to myopic optimization behaviors.
Energy-Based Models (EBMs). EBMs specify unnormalized probability density:
exp -Eo(x)]
(x(x)=
Zo - exp exp[[Eo((u]] du,
Zo
where the energy function (x) typically parameterized as a neural network through short-run MCMC-based Maximum Likelihood Estimation [36 (3). EBMs effectively guide exploration toward globally promising regions. Unlike GPs, EBMs explicitly capture multi-modal global structures, thus addressing the limitation of excessive local exploration inherent in standard GP-based methods.
Reinforcement I Learning (RI.) and Proximal Policy Optimization (PPO), RI. formalizes the optimization process as sequential decision decision-making task, wherein policy π T parameterized
3

---

by neural network weights Oppo, maps states to actions a,. PPO limiting policy changes through a clipped probability ratio:
38 stabilizes the training by
TT Toppo (a Sx)
re (pppp)
T Topoo
8t)
thereby preventing erratic changes in the parameters. We define states combinations of GP posterior estimates and global EBM signals; actions match suggested sampling points. The use of PPO helps overcome the single-step myopia that is inherent in conventional acquisition methods. This is accomplished through the use of multi-step reasoning.
4
REBMBO Model & Algorithmic Details
After updating the GP model (which focuses on local predictions) and the EBM model (which looks at global patterns), REBMBO uses the PPO technique to choose the next sample point by combining both local and global informa- tion. The GP posterior, EBM signals, and PPO's multi-step planning help REBMBO find good solutions in complex spaces with many dimen- sions or peaks. PPO's planning horizon mit- igates the short-sightedness of single-step ap- proaches. Overall, REBMBO provides uni- fied solution for balancing global exploration and sequential decision-making maninn in challenging black-box optimization settings in Figure |
Staee MiLLAAAA
Bark Oralle
Pootenor Mert
Posario S66
Action
clea ware Olead frees
onpectiv
4.1 Module A: Gaussian Process Variants
As new data comes into the REBMBO frame- work, it first goes to Module A (the Gaussian Process surrogate, shown in Figure [ im- prove the estimate of the objective function fxx). In particular, the GP uses the input-output data it has seen to improve its predictions about the average outcome and the level of uncertainty in specific areas, ensuring precise local model- ing with each update. Local modeling is essen- tial for black-box optimization because lets us acquire solid insights from a few observations before sampling the full domain.
Figure I: REBMBO Workflow Diagram. The ar- chitecture comprises: Module A for local model- ing using GP posterior; Module B which trains an EBM capture global structure; and Module C, which uses PPO-based RL agents to generate decisions and optimize via rewards shaped by both function values and EBM energy signals. Arrows of different colors represent distinct data flows: green arrows indicate model parameter updates, and purple arrows represent RI. actions, evaluation or feedback steps.
We add different GP modules in REBMBO be- cause many black-box functions are complex,
have multiple peaks, or are costly to compute, so we need flexible "core" model that can quickly adjust and measure uncertainty with just a few samples. GPs provide estimates of the average and variability for fx), making them useful for identifying significant local trends, particularly when limited data is available at the start. To deal with various problem dimensions and complexities, propose three REBMBO variants:
(1) REBMBO-C (Classic GP) 33 Employs exact Onz²) ) GP inference, which is practical for moderate n. While this variant is straightforward, it can costly for large n or high-dimensional
(2) REBMBO-S (Sparse GP) 40 Adopts a sparse approximation to alleviate the O(n²) bottleneck in higher dimensions. It introduces < inducing points {z,} and approximates Kx.x ss Kx.x K, lowering the update cost to Onm?). This approximation may lose accuracy if m or the chosen inducing points are suboptimal, but remains effective for larger datasets and higher d. In our EBM-driven acquisition, the approximate mean ( and variance 7. replace the exact GP posteriors.

---

(3) REBMBO-D (Deep GP) [41. For problems that exhibit multi- scalle non-stationary behavior, a deep kernel can capture intricate latent features beyond what standard kernels provide. We use deep network è to map inputs x into latent features øGr (xrr()), (x), then compute GP-like statistics:
[(x;; D, 0) =m T (x) + ((x); o²(x: D, θ) (cp(x)" K (GP(() + B
where 77(x) is a (potentially learned) mean function, and K a smaller covariance matrix in the latent space. With sufficient training data, this approach can represent complex functions more flexibly than a standard GP and supports sublinear regret under moderate network capacity.
All three GP variants work with the EBM-UCB and PPO modules to create the full REBMBO system; the main difference is in how each variant calculates its GP posterior. For simplicity, we will show EBM-UCB using the "Classic GP" version (REBMBO-C), but the same idea applies to the sparse and deep GP versions. More information about differences of GP variants, including their complexity
and how to implement them, can be found in Appendices F and K
As updated data then moves to Module B. REBMBO addresses the limitations of purely local exploration by introducing an Energy-Basee Model to capture global structure. The next section (Module B) explains that the EBM helps the GP surrogate by directing the search away from areas that aren't t useful and toward more promising ones, especially in complex or varied situations.
4.2 Module B: EBM-Driven Global Exploration
To overcome the limitations of purely local GP-based search, we introduce an Energy- Based Model (EBM) defined as Eg (x). After Module A updates the GP posterior with newly collected data (see Figure J Module B uses these updates to train the EBM, which captures a global "energy" landscape. While the GP's uncertainty a (x) pinpoints locally undersampled regions, the EBM reveals which basins in X' are more likely to contain near-optima pptr-optimal solutions. Combining these local and global insights helps prevent the search from stalling unproductive local pockets and enables REBMBO
to traverse complex objective surfaces more efficiently.
EBM Training Mechanism. We parameterize Ex as a neural network trained under short-run MCMC-based Maximum Likelihood Estimation (MLE). At each iteration, we alternate:
Positive Phase: Lower Ex(x) for real data points x,, guiding the model to "observed" regions. Negatíve Phase: Draw small number (K) of Langevin samples from (u( exp|-Ex((u (u)). then push these moner-generated samples to higher energy unless they reflect data-like features. This short-run MCMC procedure (e.g. Stochastic Gradient Langevin Dynamics, detailed in Appendix G ensures that low-energy regions correspond to promising global basins.
EBM Parameterization Details. We specifically train Ex via short-run MCMC-based MLEE as follows. Suppose we have data {x: from an unknown distribution Pdata. We fit 0 by minimizing
1
log (x
Σ
(xx) + log Z(0),
n
which is equivalent to maximizing E - (ogp)(xx). Let ((D) = E.= log (x:)), where (x(x)) x expl-Ea (x:)). Then, we have
V 2(0) = Σ (x
Σ
(u(u) Eo du.
Dividing by yields the well-known positive-minu--negatiee decomposition:
V C(0) = -Exxpaaca V [x)
Euwps (o
n
Positive Phase
Negative Phase
which balances data alignment against -odel-drawn samples 42
Since short-run MCMC gradient update
approximates po(u) sufficiently well, it lets us implement Robbins-Monro-Se
44. Iteratively alternating positive and negative phases steers Fx(x) to be low in data-like basins and high elsewhere, thus revealing globally promising regions for exploration.
5

---

EBM-UCB Acquisition Function. Once the EBM is trained, we embed its negative energy Eg (x) into a standard GP-UCB scheme. Let (x) and f.((x) be the GP posterior mean and standard deviation at iteration t. A typical UCB function is oUCB (x) = [xf,(()) + 3 f,(x,, where 3 >> 0
controls exploitation vs. exploration. To incorporate the EBM's global guidance, we define
0EBM- (CB = 1(,((x)) + f,)
(x))
where y O specifies how strongly (x) biases the search toward underexplored basins. In multi- modal and high-dimensional tasks, this "global penalty" helps avoid wasting evaluations in uncertain but unpromising pockets, augmenting the GP's local exploration with broader sense of global structure. As further discussed in Section [S] this synergy accelerates convergence on challenging landscapes and reduces the need for purely local or manually specified look-ahead heuristics.
Theoretical Contributions and Landscape-AAaree Regret (LAR). We employ Landscape-Awar
Regret (LAR) as a generalized regret formulation that extends the standard definition with an energy- informed global term:
R [[(x*) f(xx)] a [x(xx) Ex(x))]
where 2 0 controls the relative influence of the global energy term. For non caan-energ--baae baselines, we set a=0 to recover standard regret and ensure fair comparison, while for energy-aware methods such as REBMBO, a>0 provides holistic measure that captures missed global opportunities in the learned landscape. Under mild alignment and regularity assumptions (Appendix H. our EBM- UCB retains the GP- -CP-yee- sublinear rate, so incorporating Ex( preserves the same optimality guarantees as standard regret [22 43
Mixture kernel for the GP posterior (rationale + form). The GP posterior is computed by inverting an n x kernel matrix built from mixture of Radial Basis Function (RBF) and Matérn covariances:
xx)xxx) = wRRRR KRBF (x, x') + WMatern AMatern (x, x')]
with (x, = exp( -x- - A
x') and, for v=2.5 and r-||x - x'II.
Maaerru (x, x') = (1+ v5r + 3 5/2 accommodates rough, less-smooth local variations. The mixture enlarges the RKHS compared to either kernel alone, which matches REBMBO's 's design: the EBM offers global basin cues, and the GP needs both smooth (RBF) and rough (Matérn) components model local structure faith- fully. The mixture weights (WREF, WMatern] are learned by type-Il marginal likelihood (evidence maximization), avoiding per-task manual tuning.
3r exp(-vvr) v5r) RBF captures smooth global trends, while Matérn-
By unifying the global signal Ex( with these locally expressive GP statistics (via the mixture kernel,, REBMBO couples principled global exploration with precise local modeling; Module C then employs PPO-based multi-step planning to mitigate one-step myopia and fully exploit this synergy.
4.3 Module C: Multi-Step Planning via PPO
While QEBM-UCB (x) enhances global exploration over local approaches, single-step acquisition can still cause local myopia. The algorithm prioritizes instant rewards above long-term queries. We consider each Bayesian Optimization iteration as Markov Decision Process (MDP) to enable multi-step lookahead via reinforcement learning. Although Proximal Policy Optimization (PPO) (88) well-known, our work combines with the GP surrogate and the EBM's global energy signal. This concept combines RL's multi-round exploration with Modules A and B's local-globs modeling. Figure T illustrates how Module B updates the EBM, and Module C guides PPO-based policy adjustments based on local uncertainty and global energy cues.
MDP Formulation: States, Actions, and Rewards.
iteration t, we define the state S; (x), (x)()), Ex(x)),
where (x), (x) come from the current GP posterior, and (x) denotes the learned global energy map. The action is the proposed query point a, E X' c IRd; evaluating f(a.) updates the GP and EBM for the next state St+1-

---

To balance immediate payoff's (function values) and global exploration (pursuing low-energy basins), we define the reward
(st ar) (a))))) - X (a.),
where À > 0 governs how strongly - Eg (ax) influences exploration. A higher À promotes thorough global searching, while lower À emphasizes direct improvement in f(a:). By embedding Eg in the reward, we ensure that REBMBO actively targets regions the EBM deems globally promising.
PPO Training Process. We employ stochastic policy T Toppo (a. st) to maximize the cumulative reward over 7 steps. Though PPO is an established RL algorithm [38, our adaptation ensures co-evolves with both the GP posterior and the EBM distribution, rather than being standalone
Apppo (a//s:)
module. Concretely, we define re (1pp))))
sadld, (a (a:) which measures how much the new policy
T Tpp deviates from the previous one Tayia-
The clipped objective be maximized is
C CLIP (dppo) = Esg min (ee()) A., clip(re(@ypo). =)) Ar)]
where Ag is an advantage estimate derived from (s,,, a.) minus a learned baseline. The clipping ensures that large updates the policy are penalized, stabilizing learning.
4.4 Overall Methodology and Synergy of GP, EBM, and PPO
After each query a is evaluated, REBMBO synchronously updates three components:
1) GP Posterior Update: Incorporate f(a)) to refine ff,4+1 and f.,+1 preserving reliable local predictions. 2) EBM Retraining: Run short-run MCMC with the expanded dataset to improve Eg (x) (Section 4.2 thereby maintaining coherent global energy landscape. 3) PPO Policy Optimization: Use the new reward r2 = f(a)) XE (a.) and the transition (St, a,, St+1) to update Toppo via the clipped objective CCLIP (épp)) This loop iteratively refines the local GP model and global EBM, while the PPO agent selects multi-step stelt query points. Crucially, is not a mere stacking of separate algorithms; rather, it constitutes tightly coupled system where the RL policy co-evolves with up-to- date local posterior and global signals. The EBM term - (x) augments UCB-based sampling with long-range structure, and PPO transforms this single-step acquisition into an MDP-based multi-round planner, thereby mitigating the near-sightedness of conventional BO.
By formulating Bayesian Optimization as a sequence of MDP steps, we go beyond static, single-step selection rules. Even though EBM-UCB (Module B) already introduces a global perspective, remains one-step unless bolstered by PPO's multi-round lookahead. The reward function Tt = f(.)) - AE drives the policy toward robust global basins, balancing immediate gains and exploratory push. As the GP and EBM adapt to each new evaluation, the RL policy adjusts accordingly. improving its trajectory selection at each iteration.
Putting It All Together: Repeating this procedure yields a dynamic and adaptive optimization scheme: after every evaluation, REBMBO incorporates fresh data into the GP, retrains the EBM, and refines the PPO policy to better plan subsequent queries. Section [5 presents empirical results showing how this synergy enables REBMBO to tackle high-dimensional, multi-modal functions more effectively than single-step or purely local methods, while our theoretical analysis (Appendix [D) ensures sublinear Landscape-Aware Regret (LAR) under mild assumptions. In essence, REBMBO'S novelty lies in harmonizing old RL machinery (PPO) with EBM-driven global exploration and GP-based local modeling. thereby providing multi-step, globally aware strategy for challenging black-box box optimization tasks.
5 Experiments
5.1
Experiment Setups
In this study, we evaluate REBMBO (variants C. S, D) against leading Bayesian Optimization (BO) methods across multiple synthetic tasks, including Branin 2D, Ackley in SD, Rosenbrock in 8D, and high-dimensional BO (HDBO) in 200D, and real-world tasks such as Nanophotonic in 3D and Rosetta in 86D, shown in Figure 2 Baselines include BALLET-ICI, TuRBO, EARL-BO, and Classic BO, representing varied modeling strategies from local Gaussian Processes (TuRBO) to single-step
7

---

Algorithm REBMBO
Require: GP config (CClassic, Sparse, Deep)), EBM config (Es, MCMC steps), PPO config (Toppu
clip 6, mini-batch size), and an initial dataset Do of size nio
1: Train GP on Do to obtain (No, go).
2: Initialize EBM Eg (x) and PPO policy Képpo"
3: for t = 1 to I" do
4: (A) Update the GP with D;-1, yielding (Me, oz).
5: (B) Retrain or partially train the EBM using data in D,-1 (via short-run MCMC).
6: (C) Form the RL state: St t (.)).
7: (D) Select action: X; f Toppo (st).
8: E Evaluate: < [(x+) (expensive black-box call).
9:
Compute reward: Ft Yt XE (x:): update Toppo with (st, (s,) re) via PPO.
10: (G) Augment dataset: Dr D (t)).
11: end for
12: Return the best sampled point x" Dr in terms of f(x*).
Dremin 2D Imsin 2 20d (uee)
Actiey 50 onean at. rans
Toverbrock 3D Inion a and 5 Surra
0.3
2000
Figure 2: Bayesian optimization performance across benchmarks: (a) Branin 2D, (b) Ackley 5D, (c) Rosenbrock 8D, (d) HDBO 200D. (e) Nanophotonic 3D. (f) Rosetta 86D. REBMBO variants (blue shades) consistently outperform baselines, especially in higher dimensions.
RL (EARL-BO) and iterative confidence intervals (BALLET-ICI, We extend this analysis in Table [ and Table 2 with two additional real-worid tasks (NATS- Bench in 20D, Robot Trajectory in 40D) and two more baselines (Two-step EI, KG), covering a wider set of approximate lookahead methods and practical optimization cases. All algorithms quantitatively compared using Landscape-Awar
Regret (LAR) metric: R.- = [f(x*) f(x))] + a [Ex(x*) Ex(x,]], which jointly assesses local suboptimality and global exploration; all reported values reflect the mean + standard deviation over 5 independent runs. Further details on these baselines and benchmarks are provided in Appendix A including the rationale for selecting tasks and the chosen hyperparameter settings (such as 10-20 short-run MCMC steps per iteration in the EBM, a 2-layer policy network with 64-256 hidden units for PPO. and Matémn-RB kernel mixtures in the GP). Notably, REBMBO-D (see Section 3.2) employs deep kernel for richer latent representations, while Two-step EI and KG are included benchmark against established lookahead variants. The final scores in the tables reflect the average Landscape-Awas Aware Regret (LAR) (or normalized objective) under predefined iteration budgets (T = 30, 50 for Branin, Ackley, Rosenbrock, and T 50, 100 for HDBO); each entry in Table [ and Table 2 includes both the mean outcome and its standard deviation. Additionally. we assess computational overhead and duration on one NVIDIA A6000 GPU, running each training cycle for about five to ten minutes, with each job consuming average of 1300-1500 MB of memory. The Appendix C D provides experimental information for supplemental experiments and parameter ranges, in addition to brief summary of basclines.

---

Model
Branin 20
Ackley
Roseabrock SD)
HDBO 2001
Mean
T=30
T=50
T=30
T=50
T=30
T=50
T=50
T=100
BALLET-ICI 6
87,33_22.09
90/411199
82.81
67.78.22.14
85.55..2.1 10
90.76..0.97
79. 79.46122.85
85.8513.18
83.80.11:45
EARL-BO
9
85.13:0096
88.76:12.28
50.46::1.23
87.22:: 82
83.47:1.5 .6
88.47:00.98
77.24:2.87
83.7412.81
81.57:1.255
TuRBO 5
80.65-100 141
8846342.49
7806++206
83.70+2.19
N1.22+1.2
85.74=1.32
74.72+3.56
80.60+3.14
78.56+1.33
Two-Step
7
20.27-2.04
92.38+2.15
14-12-1.87
89.14-1.78
85.19+1.67
88-57+143
18.10+3.22
KG
88.64=1.83
91.53+1.97
86.71=1.78
90.23=2.11
87.95
90.29=1.67
79.63=3.10
85.1742.96
87.52+1.67
REBMBO-S
88 88.89=1.55 1.54
96.95+22446
86.85=1.000
92.64=161 61
92.87=11.53
98. 96.85=00.86
83.33=3305 0.06
90.16++2.00
87.96+1.25
REBMBO-D
93.65-1.38
95.2141..0
85.25:11.48
91,53+1555
91.974:2.02
96998:1.09
88.79
94,42:3.96
89.17-141
RERMBO-C
90.83-1.09
97.37+2.07
89.93-1.25
94,46-1.05
91.28+1.50
96.77+11.2
85.55-3.23
90.95+3.57
89.40+1.24
Table I: (a) Performance comparison across synthetic benchmarks evaluating REBMBO variants (S, D, C) against existing Bayesian optimization methods. Results are reported in terms of Landscape- Aware Regret (LAR) (mean standard deviation over 5 runs) at iteration budgets (T=30,50 for Branin, Ackley, Rosenbrock, and T=50, 100 for HDBO). Higher scores indicate superior optimization efficiency. Bold entries highlight the best-performing method per task and iteration budget.
5.2
Main Results
As depicted in Figure 2 all three REBMBO variants lower Landscape- Regret (LAR) more rapidly than baseline methods across the six tested benchmarks. particularly excelling on higher- dimensional tasks. Table | (synthetic) and Table 2 world) further quantify these findings for iteration budgets T 30 and T 50. On the lower-dimensional Branin (2D) and Ackley (5D), for example, REBMBO-S achieves roughly 15-20% lower final pscudo-regret compared with EARL-BO and BALLET-ICI, while standard approaches like TuRBO and Two-step EI exhibit slower global exploration. Moving to Rosenbrock (8D) and HDBO (200D), REBMBO-D stands out: on HDBO (200D), its final Ladscapp-AAar Regret (LAR) is less than half that of KG and BALLET-ICI by iteration 50. Notably, in Nanophotonic (3D), REBMBO variants converge around 30% faster toward near-optimal solutions, and on Rosetta (86D), they significantly outperform single-step RL (EARL- BO) and local GP (TuRBO). These consistent gains support the theoretical premise that combining global EBM cues with PPO-driven multi-step planning yields robust sublinear sundscape-Awart Regret (LAR), even with approximate EBM and RL training.
Model
Nanophotonic 3D
Rosetta M6D
NATS- 2013
Robot Trajectory 401)
Mcan
T-50
T-80
T-50
T-50
T-80
T-50
T-80
RALLETTICI
23.77+2.96
68.64+2.68
76.75-2.3
83.9++3155
2.99
24.25+281
7845++3211
82.65-2.89
82.02+2
EARL-BO g
81.72+4.08
26.68+2.60
74.78-241
81.95+3900
80.44+3.12
83.47+3.25
76.50+287
80.11=2.90
80.70=-2 89
TuRBO 5
79.75#3.17
8.77
12.47=2:74
79.86+22.42
79. 3.225
$1.55::345
74.324-2.99
78.60=3.11
78.81 99
Two
84.29113.0
89:47-2.88
78.90 80
84.7513.05
83.3343.110
86.80.1.995
79.5513.05
83.92-2.99
83.88.12.87
85.1012.90
90.05:1..60
79.79_2.88
85.2013.00
81.1012.95
87.251175
$0.2012.50
81.4522.85
81.39112.71
REBMBO-C
87.25
92.65
80.96-2.66
83.33.13.31
85.13.12.3
89.20:12455
82.5012.77
87.10-22.50
86.591.263
REBMBO-D
31.66:33116
91.53-33.13
84.22-3.38
90.84:33774
88.95=2.776
90.30=22.99
83.25+22.75
88.10-2.60
$6.98:22.00
REBMBO-
86,5012.33
93.99.3.27
80.53-2.17
88.8811..76
85.1012.65
$9.4511.50
818851260
86.50_2.40
86.60.1..59
Table 2: (b) Performance comparison across real-world benchmarks. Results are reported in terms
of normalized optimization accuracy (mean + standard deviation over 5 runs) at iteration budgets = 50, 80. Higher scores indicate superior optimization efficiency.
5.3 Supplementary Experiment
In addition our primary benchmarks, we conducted several supplementary experiments (see Appendix(C) to further validate REBMBO's theoretical guarantees and empirical robustness under diverse conditions.
5.3.1 Design and Modeling Choices
An ablation study (Appendix Table Fand Table isolates the roles of EBM, multi-step PPO, and short-run MCMC by incrementally removing or modifying these components, and performance drops whenever core element is omitted, which confirms that global energy- based exploration, local GP modeling, and reinforcement learning each contribute critically to REBMBO. We further test kernel choice in the GP surrogate and find that learned RBF+ MatéMa mixture performs best on Branin 2D, Ackley 5D, and HDBO 200D (Appendix C.TO. Table 7 which supports the sum RKHS view and tighter regret guarantees. A one-at-a-time hyperparameter study and dedicated sweep for X

---

identify broad safe band À E [0.2, 0.5] and show that the default configuration is near optimal (Appendix [.II Tables 8 and 9). which supports the theory that a balanced reward f(z) XE (a)
keeps information gain controlled and simplifies tuning.
5.3.2 Robustness and Reliability
We illustrate REBMBO's behavior on 1D/2D multi-modal functions (Appendix Figures showing how EBM-UCB avoids local optima and uses broader structural information, and trajectory compar- isons (Appendix Figures5jj show less unnecessary exploration than GP-UCB and GLASSES with more direct convergence global maxima. Robustness tests cover EBM convergence and removal and also scale mismatch between f and Eş; REBMBO-C degrades gracefully under failed EBM and remains competitive without it, and normalization plus adaptive À recovers most losses under severe scale gaps while keeping PPO stable most runs (Appendix C.12
Tables 10 (Qand
5.3.3 Practicality and Fair Evaluation
We quantify compute overhead and observe small constant-factor increase relative to TuRBO that matches polynomial scaling and parallelizes well on GPU, which is negligible when function evaluations dominate time (Appendix C.13 Table 12). Finally, report standard regret in addition to pscudo-regret and REBMBO-C achieves the best values on all three tasks, which shows that improvements not tied to metric and that dual reporting reflects both exploration quality and final solution quality (Appendix C.14 Table 13): taken together with a benefit and overhead analysis (Appendix Figure 77, detailed comparisons (Appendix Table(6) performance heatmaps (Appendix Figureß) 8) and statistical significance checks (Appendix Figures 1O), these results corroborate the premise of robust sublinear Landscape-Awan Regret when global EBM signals and multi-step RL are integrated and they reinforce the practical value of REBMBO for challenging BBO tasks.
Conclusion
REBMBO tackled a fundamental Bayesian optimization problem: combining local uncertainty estimates with global structure exploration. Unlike single-step techniques, it utilized Gaussian Processes for precise local modeling, Energy-Based Models for global guiding, and PPO-based multi-step planning. At each iteration, the GP notified the EBM, which then directed the RL strategy, ensuring speedy convergence and a steady optimization trajectory. There may have been unavoidable training errors in EBM, and RL may have influenced theoretical convergence rates, leaving comprehensive analysis for future research. Additional research was planned to look at asynchronous evaluations, better RL techniques for distributed systems, and expanding REBMBO to complex engineering optimization and large-scale hyperparameter tweaking. More broadly, combining probabilistic modeling with multi-ster step RL has shown promise for scientific simulations and real-time decision decision-making in dynamic settings.
Acknowledgments
This research was partially supported by the following sources: Y. W. is partially supported by NSF DMS- 2M1-24, DARPA W912CG25CA007, and research gift funds from Amazon and Qualcomm. We express our gratitude PhD candidates Peiyu Yu and Hengzhi He from the University of California, Los Angeles, along with Dr. Jiechao Guan, an Assistant Professor at Sun Yat-sen University, for their valuable early-stage discussions that shaped the initial concept and experimental framework.
10

---

References
[1] Stéphane Alarie, Charles Audet, Aimen E Gheribi, Michael Kokkolaras, and Séhastien Le Di- gabel. Two decades of blackbox optimization applications. EURO Journal Computational Optimization, 9: 9:100011, 2021.
[2] Xilu Wang. Yaochu Jin, Sebastian Schmitt, and Markus Olhofer. Recent advances in bayesian optimization. ACM Compuring Surveys, 55(13s):1-36, 2023.
[3] Ryan Turner, David Eriksson, Michael McCourt, Juha Kiili, Eero Laaksonen, Zhen Xu, and Isabelle Guyon. Bayesian optimization is superior random search for machine learning hyperparameter tuning: Analysis of the black-box optimization challenge 2020. In NeurIPS 2020 Competition and Demonstration Track, pages 3-26. PMLR, 2021.
[4] Z Wang, M Zoghi, Hutter, D Matheson, and ND Freitas. Bayesian optimization a billion dimensions via random em beddbddings. In Proceedings of the Twenty-Thing Third international joint conference on Artificial Intelligence., 2013.
[5] David Eriksson, Michael Pearce, Jacob Gardner, Ryan D Turner, and Matthias Poloczek. Scalable global optimization via local bayesian optimization. Advances in neural information processing systems, 32, 2019.
[6] Fengxue Zhang, Jialin Song, James C Bowden, Alexander Ladd, Yisong Yue, Thomas Desautels, and Yuxin Chen. Learning regions of interest for bayesian optimization with adaptive level-set estimation. In International Conference Machine Learning, pages 11579-41595 PMLR, 2023.
[7] Eric Lee, David Eriksson, David Bindel, Bolong Cheng, and Mike Mecourt. Efficient rollout strategies for bayesian optimization. In Conference on Uncertainty in Artificial Intelligence. pages 260-269. PMLR, 2020.
[8] Ilija Bogunovic, Jonathan Scarlett, Andreas Krause, and Volkan Cevher. Truncated variance reduction: A unified approach to bayesian optimization and level-set estimation. Advances in neural information processing systems, 29, 2016
[9] Mujin Cheon, Jay H Lee, Dong- Yeun Koh, and Calvin Tsay. EARL-BO: Reinforcement learning for step lookahead, nigh-dimensional bayesian optimization. In Forty-second International Conference on Machine Learning, 2025.
[10] Sebastian E Ament and Carla P Gomes. Scalable first-order bayesian optimization via structured automatic differentiation. In International Conference on Machine Learning, pages 500-516. PMLR, 2022.
[11] Kirthevasan Kandasamy, Jeff Schneider, and Barnabás Póezos. High dimensional bayesian optimisation and bandits via additive models. International conference on machine learning, pages 295-304. PMLR, 2015.
[12] Rafael Gómez- Bombarell,,, Jennifer Wei, David Duvenaud, José Miguel Hernández-Looato Timothy D Hirzel, Benjamín Sánchez- Lengeling, Dennis Sheberla, Jorge chemical design using a data-driven Ryan Adams, and Alán Aspuru-Guzik. Automatic continuous representation of molecules. ACS central science, 4(2):268-277 2018.
[13] Changyong Oh, Tomezak, Efstratios Gavves, and Max Welling. Combo: Combinatorial bayesian optimization using graph representations. In ICML Workshop on Learning and Reasoning with Graph-Structured Data, 2019.
[14] James Bergstra, Rémi Bardenet, Yoshua Bengio, and Balázs Kégl. Algorithms for hyper- parameter optimization. Advances in neural information processing systems, 24, 2011.
[15] Frank Hutter, Holger H Hoos, and Kevin Leyton- Sequential model-based optimization for general algorithm configuration. In Leurning and Intelligent Optimization: 5th International Conference, LION 5, Rome, Itly, January 17-21, 201 Selected Papers 5, pages 507-523. Springer, 2011.
11

---

[16] Zi Wang, Chengtao Li, Stefanie Jegelka, and Pushmeet Kohli. Batched high-dimensional bayesian optimization via structural kernel learning. In International conference on machine learning, pages 3656-3664. PMLR, 2017.
[17] Jian Wu and Peter Frazier. The parallel knowledge gradient method for batch bayesian opti- mization. Advances in neural information processing systems, 29, 2016.
[18] Iija Bogunovic, Jonathan Scarlett, Stefanie Jegelka, and Volkan Cevher. Adversarially robust optimization with gaussian processes. Advances in neural information processing systems, 31, 2018.
[19] Robert B Gramacy. lagp: large-scale spatial modeling via local approximate gaussian processes in Journal of Statistical Software, 72:1 16, 2016.
[20] Michael A Gelbart, Jasper Snoek, and Ryan P Adams. Bayesian optimization with unknown constraints. arXiv preprint
2014.
[21] Joshua Knowles. Parego: A hybrid algorithm with on-line landscape approximation for expen- sive multiobjective optimization problems. IEEE transactions on evolutionary computation. 10(1):50-66, 2006.
[22] Wolfgang Ponweiser, Tobias Wagner, Dirk Biermann, and Markus Vineze. Multiobjective optimization on a limited budget of evaluations using model-assisted-metrii selection. In International conference purallel problem solving from nature, pages 784-794. Springer, 2008.
[23] Daniel Hernández-Lobato, Jose Hernandez-Lobato, Amar Shah, and Ryan Adams. Predictive entropy search for multi-objective bayesian optimization. In International conference machine learning. pages 1492-1501. PMLR, 2016.
[24] Matthias Feurer, Katharina Eggensperger, Stefan Falkner, Marius Lindauer, and Frank Hutter. Practical automated machine learning for the automl challenge 2018. In International workshop on automatic machine learning at ICML, volume 12, 2018.
[25] Martin Wistuba, Nicolas Schilling, and Lars Schmidt- Thieme. Hyperparameter search space pruning-a new component for sequential model-based hyperparameter optimization. In Machine Learning and Knowledge Discovery in Databases: European Conference, ECML PKDD 2015, Porto, Portugal, September 7-11, 2015, Proceedings, Part II 15, pages 104-119. Springer, 2015.
[26] Mare C Kennedy and Anthony O'Hagan. Predicting the output from a complex computer code when fast approximations are available. Biometrika, 87(1):1-13, 2000.
[27] Kirthevasan Kandasamy, Gautam Dasarathy. Jeff Schneider, and Barnabás P6ezos. Multi-fidelin bayesian optimisation with continuous approximations. In International conference on machine learning, pages 1799-1808 PMLR, 2017.
[28] David Ginsbourger, Janis Janusevskis, and Rodolphe Le Riche. Dealing with asynchronicity in parallel Gaussian process based global optimization. PhD thesis, Mines Saint- Stinnn,, 2011.
[29] Thomas Desautels, Andreas Krause, and Joel W Burdick. Parallelizing exploration-exxloitaaionn tradeoffs in gaussian process bandit optimization. J. Mach. Learn. Res., 15(1):3873-392) 2014.
[30] Kirthevasan Kandasamy, Akshay Krishnamurthy. Jeff Schneider, and Barnabás P6ezos. Paral- lelised bayesian optimisation via thompson sampling. In International conference on artificial intelligence and statistics, pages 133-142. PMLR, 2018
[31] Shangda Yang, Vitaly Zankin, Maximilian Balandat, Stefan Scherer, Kevin Carlberg, Neil Walton, and Kody JH Law. Accelerating look-ahead in bayesian optimization: Multilevel monte carlo is all you need. arXiv preprint arXiv:2402.02/ 2024.
[32] Jose Pablo Folch, Calvin Tsay, Robert Lee, Behrang Shafei, Weronika Ormaniec, Andreas Krause, Mark van der Wilk, Ruth Misener, and Mojmir Mutny. Transition constrained bayesian optimization via markov decision processes. Advances Neural Information Processing Systems, 18888894-8887 2024.
12

---

[33] Javier González, Michael Osbome, and Neil Lawrence. Glasses: Relieving the myopia of
bayesian optimisation. In Artificial Intelligence and Statistics, pages 790-799. PMLR, 2016.
[34] Peter Frazier, A tutorial on bayesian optimization. arXiv preprint
2018.
[35] David JC MacKay et al. Introduction to gaussian processes. NATO ASI series F computer and systems sciences, 168:133-1666 1998.
[36] Tetiana Parshakova, Jean-Mare Andreoli, and Mare Dymetman. Distributional reinforcement learning for energy-based sequential models. arXiv preprint arXiv:1912.08517,, 2019.
[37] Yang Song and Diederik Kingma. How to train your energy-based models. arXiv preprint arXiv:2101.03288 2021.
[38] John Schulman, Filip Wolski, Prafulla Dhariwal, Alec Radford, and Oleg Klimov. Proximal policy optimization algorithms. arXiv preprint vrXvv:107.00347, 2017.
[39] Miguel Donzalez-Duque, Richard Michael, Simon Bartels, Yevgen Zainchkovskyy, Søren Hauberg, and Wouter Boomsma. A survey and benchmark of high- ligh-dimensional bayesian optimization of discrete sequences. arXiv preprint arXiv:2406.044 04739, 2024.
[40] Mitchell McIntire, Daniel Ratner, and Stefano Ermon. Sparse gaussian processes for bayesian optimization. In UAI, volume 3, page 4, 2016.
[41] Andrew Gordon Wilson, Zhiting Hu, Ruslan Salakhutdinov, and Eric P Xing. Deep kernel learning. In Artificial intelligence and statistics, pages 370-378. PMLR, 2016.
[42] Erik Nijkamp. Mitch Hill, Tian Han, Song-Chun Zhu, and Ying Nian Wu. On the anatomy of mcme-basee maximum likelihood learning of energy-based models. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 34, pages 5272-5280, 2020.
[43] Bo Pang, Tian Han, Erik Nijkamp, Song-Chin Chun Zhu, and Ying Nian Wu. Learning latent space energy-based prior model. Advances in Neural Information Processing Systems, 33:2 1994 22008. 2020.
[44] H. Robbins and S. Monro. stochastic approximation method. Annals of Mathematical Statistics, 1951.
[45] Jialin Song, Yury S Tokpanov, Yuxin Chen, Dagny Fleischman, Kate T Fountaine, Harry A Atwater, and Yisong Yue. Optimizing photonic nanostructures via multi multi-fidelity gaussian processes. arXiv preprint 1rXiv:
2018
[46] Anqing Jiang. Osamu Yoshie, and LiangYao Chen. A new multilayer optical film optimal method based on deep q-learning. arXiv preprint arXiv: 18/2.028/2.028/// 2018
[47] Tian Hua Wang, Hai Ping Shao, Bing Qiang Zhao, and Hong Lin Zhai. Molecular insights into the variability in infection and immune evasion capabilities of sars-cov-2 variants: A sequence and structural investigation of the rbd domain. Journal of Chemical Information and Modeling. 64(8))3503-3523, 2024.
[48] Thomas Desautels, Adam Zemla, Edmond Lau, Magdalena Franco, and Daniel Faissol. Rapid in silico design of antibodies targeting sars-cov-2 using machine learning and supercomputing. BioRxiv, pages 2020-04, 2020.
[49] Rhiju Das and David Baker. Macromolecular modeling with rosetta. Annu. Rev. Biochem., 77(1):363-382 2008
[50] Kyle Barlow, Shane O Conchuir, Samuel Thompson, Pooja Suresh, James Lucas, Markus Heinonen, and Tanja Kortemme. Flex ddg: Rosetta ensemble-hased estimation of changes in protein-proted binding affinity upon mutation. The Journal of Physical Chemistry B, 122(21):5389-559) 2018.
[51] Xuanyi Dong, Lu Liu, Katarzyna Musial, and Bogdan Gabrys. Nats-bench: Benchmarking nas algorithms for architecture topology and size. IEEE transactions on pattern analysis and machine intelligence, (4(7):3644-3666 2021.
13

---

[52] A Rupam Mahmood, Dmytro Korenkevych, Gautham Vasan, William Ma, and James Bergstra. Benchmarking reinforcement learning algorithms on real-world robots. In Conference on robot learning, pages 561 361-591. PMLR, 2018.
[53] Jian Cao, Myeongjong Kang, Felix Jimenez, Huiyan Sang, Florian Tobias Schaefer, and Matthias Katzfuss. Variational sparse inverse cholesky approximation for latent gaussian processes via double kullback-leibler minimization. In International Conference on Machine Learning, pages 3559-3576, PMLR, 2023.
[54] Shion Takeno, Masahiro Nomura, and Masayuki Karasuyama. Towards practical preferential bayesian optimization with skew gaussian processes. In International Conference on Machine Learning, pages 13516-33533. PMLR, 2023.
[55] Carl Edward Rasmussen and Christopher KI Williams. Gaussian processes for machine learning, ser. adaptive computation and machine learning. Cambridge, MA, UsA: MIT Press, 38:715-719, 2006.
[56] Will Grathwohl, Kuan-Chich Wang, Jörn-Henri Jacobsen, David Duvenaud, and Richard Zemel, Learning the stein discrepancy for training and evaluating energy-based models without sampling. In International Conference Machine Learning, pages 3732-3747. PMLR, 2020.
[57] Masaaki Nagata, Katsuki Chousa, and Masaaki Nishino. A supervised word alignment method based on cross-language span prediction using multilingual bert. Proceedings of the 2020 Conference Empirical Methods in Natural Language Processing (EMNLP), pages 555-565, 2020.
[58] Geoffrey E Hinton. Training products of experts by minimizing contrastive divergence. Neural computation, 14(8):1771-1800, 2002.
[59] Yakoub Bazi, Mohamad M AI Rahhal, Haikel Alhichri, and Naif Alajlan. Simple yet effec- tive fine-tuning of deep cnns using an auxiliary classification loss for remote sensing scene classification. Remote Sensing, 11(24):2908, 2019.
[60] Tijmen Tieleman. Training restricted boltzmann machines using approximations to the likeli- hood gradient. In Proceedings of the 25th international conference Machine learning, pages 1064-1071, 2008.
[61] Richard S Sutton. Reinforcement learning: An introduction. A Bradford Book, 2018.
[62] Sham M Kakade. A natural policy gradient. Advances in neural information processing systems, 14, 2001.
[63] John Schulman, Philipp Moritz, Sergey Levine, Michael Jordan, and Pieter Abbeel. High- dimensional continuous control using generalized advantage estimation. arXiv preprint arXiv:1506.024430 2015.
[64] Wenjie Xu, Wenbin Wang, Yuning Jiang, Bratislav Svetozarevic, and Colin Jones. Principled preferential bayesian optimization. arXiv preprint arXiv:2402.05337, 2024.
[65] Nathaël Da Costa, Marvin Pförtner, Lancelot Da Costa, and Philipp Hennig. Sample path regularity of gaussian processes from the covariance kerne.. arXiv preprint arXiv:2312. 14886, 2023.
[66] Mónika Farsang and Luca Szegletes. Decaying clipping range in proximal policy optimization. In 2021 IEEE 15th International Symposium on Applied Computational Intelligence and Informatics (SACI), pages 000521-000526. IEEE, 2021.
14

---

A Benchmark Details
A.1 Branin Toy Dataset (2D)
We employ the standard Branin objective function for a two-dimensional input vector x (2).
The function is given by:
f(x) (I2 5. o + 6) -6
10 10 (1 5
cos(x1) + 10.
In our experiments, each drawn from bounded domain (e.g.. [-5, 10) [0, 15)) and the output (x) serves as a classical test for global optimization. Within the REBMBO framework, the input consists of the current guess RR² and the output is the Branin objective value. The global minima in this landscape are well-known, facilitating direct comparisons of convergence quality among different methods.
A.2 Ackley Function (5D)
Although the Ackley function is often defined in two dimensions, we employ a five- five-dimensional variant to increase the complexity of local minima. The Ackley function typically has the form:
d
f(x)
expl
x 1=1
exp 5/4
Σ cos(cxz)) +a +ate, 1=1
where d = 5 in our case, and standard constants (a 20, b 0.2, c= 2x) are chosen. The domain can be set to (-32.668, 32.768) Inputs are thus five-dimensional vectors, while the output provides continuous measure of fitness, riddled with many local minima. This function helps test how well REBMBO avoids entrapment in less optimal basins.
A.3 Rosenbrock Function (8D)
The Rosenbrock function, often referred to as the Banana function, appears here in an eight- dimensional variant:
7
f(x) Σ 100 100 (z.-1 (z:-1- 27)² + (1-x,)²].
2z
Although is uni-modal, the valley leading to the global minimum curved and narrow, making convergence notoriously difficult for local methods. We set an input domain such as 2] and let each x IR map scalar output fxx). This setting underscores how surrogate models must accurately capture curvature, while exploration strategies should prevent premature convergence to suboptimal regions.
A.4 Nanophotonic Structure Design (3D)
In this dataset, each input x E R represents physical parameters (for instance, thickness or refractive indices) that define layered optical filter. The output corresponds weighted figure of merit for transmitting targeted wavelengths, derived from solving discretized Maxwell's equations. Given the computational intensity of this solver, strategies such as deep Q-learning be employed to optimize the design efficiency of multi-layer films [45]. This approach converges to the global optimum the optical thin film structure, addressing the limitations of traditional numerical algorithms that often converge to local optima [46]. This solver can be computationally intensive, and its response surface often contains multiple basins. Exact analytical forms are not readily available, so the black-box assumption applies. In REBMBO, the agent queries the simulator with a proposed x, and obtains the numeric figure of merit indicating how well the design meets hyperspectral criteria.
A.5 HDBO-200D
To investigate digh-dimensional performance, we construct a synthetic function 200 dimensions. Each component z; is initially drawn from a standard normal distribution, and the objective is:
200
f(x)
Σ
i=1
15

---

Despite its additive structure, this function remains challenging when searching over wide ranges, as naive methods often converge slowly or fail exploit the exponential coupling. We select suitable domain constraint (e.g., [-5, 5) 200) to test how well each method scales in dimension and maintains exploration. Since we know the ground-truth form, is possible to measure how quickly the algorithm recovers near-optimal solutions or locates valuable subregion.
A.6 Rosetta Protein Design (86D)
This data set represents a realistic antibody engincering task, where each imput x E R encodes structural modifications relative to a reference antibody that binds the SARS-CoV-2 spike protein. The implications of minimizing changes in binding free energy (AAG) are significant for the efficacy target the SARS- --V- spike protein. As
and safety of antibodies engineered variants emerge, such as the o lineage, which exhibits strong positive mmarge-enhancing electrostatic potential energy, antibodies must be designed with surfaces rich in negative potential to counteract this charge and maintain binding efficacy. Additionally, engineered antibodies should exhibit strong van der Waals interactions postbinding to effectively neutralize various strains, including those with The simulator uses Rosetta Flex [49] 50)
new mutations that may weaken binding affinity. 47 48
to estimate changes in binding free energy (AAG), requiring extensive CPU time per query. The objective is to minimize AAG (or equivalently maximize AAGG) though no simple functional form exists. This woold scenario underlines the capability of REBMBO to handle expensive, nonlinear responses in a high-dimensional space.
A.7 Additional Benchmark: NATS-Bench (20D)
This benchmark is drawn from a unified framework for neural architecture search (NAS) designed to evaluate both architecture topology and size [51. It provides large space (15,625 real candidates for topology and 32,768 for size) of precomputed performance results, covering multiple image classification datasets under consistent training protocols. In our setup, each configuration x R200 (or discrete embedding) represents distinct candidate architecture in the NATS-Bende search space, where the goal is to maximize validation or test accuracy of the trained model. By offering uniform evaluation procedures across thousands of architectures, NATS Bench enables fair comparisons among diverse NAS algorithms, which critical for assessing the advantages of multi-step exploration in -ish-iiieniional partially discrete design spaces. In this work, we adapt the bascline accuracy metrics from NATS-Bench, treat them as direct performance signals, and apply REBMBO to identify architectures yielding superior classification aceuracy. This setting highlights REBMBO's 's ability to navigate large-scale architecture landscapes, incorporate global cues from its energy- based model, and mitigate single-step myopia through PPO-based multi-step planning.
A.8 Additional Benchmark: Robot Trajectory (40D)
This benchmark is drawn from real-world robotic tasks introduced by Mahmood et al. (52) featuring multiple commercially available robots with varying degrees of difficulty and repeatability, We represent each trajectory control policy 40-dimensional parameter vector, where different components may govern joint angle targets, velocity limits, and timing schedules. The overarching objective is to improve continuous control performance in physical robot setting, which involves tuning these parameters to maximize a task-specific reward. Unlike purely simulated domains, the system dynamics and sensor feedback here subject real-world noise and hardware constraints, making it an especially challenging black-box optimization problem. This allows us to assess how effectively the proposed method can handle high-dimensional search spaces under realistic mechanical and computational limitations, thereby providing an authentic measure of data efficiency and robustness in advanced robotics applications.
B Baseline Characteristics and Ranking
Table T provides a side-by-side comparison of the baseline algorithms and our proposed REBMBO variants, listing each method's modeling strategy (local, single-step RL, short-horizon lookahead, or multi-step RI.), the underlying surrogate (e.g. classical GP, sparse GP, deep GP), the chosen kernel family (Matérn, RBF, or mixture), and the acquisition function. The final column shows
16

---

Methods
BO Model
Surrogate Model
Kernel
Acquisition Function
Avg. Ranking
TuRBO
Local
GPy
Matémn 5/2
Local UCB/EI
0,065
BALLET ICI
Global Local
GP5e +OOGGP))
RBF
UCB LCB
0.058
EARL-BO
Single-step RL
GP
RBF White RBF
Single-step RL
0,072
2-step EI
Shon-horizon lookahead
GP;
2-step
0.049
KG
Single-step
GP-
Matéen
Knowledge Gradient
0.056
REBMBO-C
Multi-step RI.
GP
Mix of Matérn & RBF
EBM-UCB
0.032
REBMBO-S
Multi-step RL
Sparee GP GP
Mix Matern & RBF
EBM-UCB
0.025
REBMBO-D
Multi-step RT.
Mix of Matérn & RBF
FBM-UCB
0.015
Table 3: Baseline features (BO modeling technique, surrogate type, kernel, and acquisition function) and average synthetic and real-world benchmark rankings (lower is better). The proposed multi-step RL variations (REBMBO-C, REBMBO-S, REBMBO-D) include an Energy-Based Model (EBM) for global signals, which helps overcome single-step local constraints.
the average ranking across both synthetic and real-world benchmarks, with lower values indicating superior performance. Specifically, TuRBO focuses on local trust regions with a Matérn-5/2 kernel, while BALLET-ICI alternates between global and local GPs, employing an RBF kernel. EARL-BOO introduces single-step reinforcement learning approach based on GP surrogates, whereas Two-step El and KG exemplify short-horizon and single-step lookahead strategies, respectively. In contrast, each of the REBMBO variants (C, S. D) leverages multi-step RL, combining an Energy- Model (EBM) term, a suitable GP surrogate, and a mix of Matérn and RBF kernels. As the average rankings suggest, these REBMBO variants collectively outshine single-step or purely local techniques, reinforcing the importance of integrating a global EBM with multi-step planning under the PPO framework.
C Supplementary Experiment
C.1 Abalation study
Model
Component Usage
A (No EBM)
FRM, PPO (multi-step RL.), Short-run MCMC
B (No PPO)
EBM, PPO (multi-step RL), Shurt-run MCMC
C (EBM w/o MCMC)
EBM: partial (no short-rum MCMC), PPO (multi-step RL): 100, Short-run MCMC
(No PPO Incomplete EBM)
EBM: partial (incomplete training). PPO (multi-step RL), Short-run MCMC: partial
Complete Model
EBM, PPO (multi-step RL), Short-run MCMC
Table 4:
Ablation Settings Highlighting Key Components (EBM, PPO, Short-run MCMC). "Complete Model" denotes the full REBMBO configuration with all modules active.
Dataset
B
D 8.60 0.28
Complete Model
Branin (2D)
8.95 = 0.20
9.10 = 0.15
9.00 = 0.18
9.30 = 0.12
Nanophotonic (3D)
-1.05
0.07
-0.92
0.05
-0.97 + 0.06
-1.02 + 0.06
-0.85 + 0.04
Table 5:
Final objective (Mean Std) Branin (2D) and Nanophotonic (3D). "Complete Model"
indicates the full REBMBO configuration.
Numerical Analysis. Model A (no EBM) shows a drop performance on both tasks, confirming the importance of global exploration signals. Model B (no multi-step PPO) converges faster than a purely local bascline but still lags behind the complete model, indicating that PPO's lookahead mitigates single-step myopia.
Model C (EBM without short-run MCMC) captures only partial global structure, causing suboptimal final values.
Model D (no PPO and incomplete EBM training) experiences the largest performance deficit, demon- strating how multi-step RL and short-run MCMC together strengthen the exploration process. In each case, the complete model surpasses or rivals the best reduced configuration, underscoring how each module adds data efficiency and final objective outcomes.
These ablation studies confirm that each core component the REBMBO framework plays critical role. Short-run MCMC training expands the EBM's capacity to capture complex global
17

---

basins, while PPO's multi-step horizon counters local or single-step limitations, and the GP surrogate ensures accurate local predictions. Removing or weakening any piece degrades final performance. either by limiting exploration or reducing the algorithm's & planning capability. 'hee combined results validate our central proposition that global anergy-based exploration, local GP modeling, and multi- step reinforcement learning are all vital to REBMBO's success in handling challenging black-box optimization scenarios.
C.2 1D Multi-modal Function Optimization
This experiment demonstrates how REBMBO leverages its Energy-Based Model (EBM) navigate multiple local maxima in one-dimensione domain, with a small set of initial observations (red points) spanning 1.0 (dashed black line) and learn the EBM via short-run MCMC. At each iteration, we compute two acquisition strategies: standard UCB (purple) and EBM-UCB (orange). contrasts the GP mean (solid blue) and uncertainty (shaded area) with the true function, highlighting how EBM-UCB (orange triangles) focuses exploration around the global optimum (vertical green dashed line at R 0.25) more effectively than standard UCB (purple triangles). The middle panel plots - (2). revealing where the EBM predicts promising global regions, and the bottom panel compares acquisition values from both strategies. FBM-UCB attains higher values at key peaks, avoiding local traps and driving samples toward the true global maximum. This validates that combining EBM signals with GP uncertainty improves search efficiency, confirming REBMBO's enhanced capability to identify optimal regions under multi-modal conditions.
1.0. We fit a Gaussian Process (GP) to approximate the true function
Figure the top panel
1D Multi-modal Function Optimization
Uncwtsiity
Oosevvtioos
Optimu
EERRUCE Seactin
-1-
-1.00
-0.75
-0.3
+1.5
2.00
Mojuiino Energy 12k0
-1.00
-0.75
-0.3
-1.25
2.00
Accusstton
-1.00
-0.75
-0.53
-2.25
2.00
0.50
Figure 3: 1D multi-modal function experiment. The top panel shows the true function (dashed), GP mean (solid), and sampled points from UCB (purple) vs. EBM-UCB (orange). The middle panel illustrates the learned negative energy - Fg (z). The bottom panel compares the acquisition functions, emphasizing how EBM-UCB better targets the global peak near x RJ 0.25.
C.3 2D Multi-modal Function Component Analysis
This experiment decomposes REBMBO in two-dimensional multi-modal setting to highlight how its GP surrogate, Energy-Based Model (EBM), and acquisition mechanism jointly guide sampling. In Figure 4 the top-left panel shows the true function, where two distinct basins of high value
18

---

emerge, along with initial observations (various markers). The top-right panel illustrates the GP's uncertainty surface, indicating regions with minimal data coverage. The bottom-left panel plots the learned negative energy Egg (x), revealing how EBM emphasizes broader structural cues rather than solely relying on local uncertainty. Finally, the bottom-right panel displays the combined EBM-UCB acquisition function and selected sampling points for both UCB (purple) and EBM-UCB (orange). The EBM-UCB approach consistently prioritizes promising global basins, as indicated by the higher single local optimum. By
acquisition values near both peaks, thus avoiding myopic focus comparing these panels, we see that the global signal from - (x)(x) helps REBMBO balance local exploration (informed by GP uncertainty) with global navigation (guided by the EBM), validating its capacity to locate and refine multiple maxima in complex search spaces.
1125
10.0
0.7
4.7
370
0.5-
0.00
0.00
0.00
0
-0.25
-0.2
425
0.75
a1s
-1.00
0.00
-0.15
Biki
2000
1.75
-258
1.5
-048
1.00
-0.25
6.00
-0.25
-0.50
-0.72
-0.50
0.0
-1.7%
-0.7%
1.00
1.00
0.00
Figure 4: 2D multi-modal decomposition of REBMBO. Top-left: true function with observations, top-right: GP uncertainty, bottom-left: EBM negative energy, bottom-right: EBM-UCB acquisition function illustrating how REBMBO's global cues target multiple peaks rather than focusing on single local maximum.
C.4 2D Multi-modal Optimization Trajectory Comparison
To illustrate REBMBO's sampling efficiency relative GP-UCB and GLASSES, optimize two-peak function in the 2D plane and plot each method's trajectory over ten iterations. In Figure 5 the top-left panel overlays all three trajectories on the true function contours, while the other panels present each method individually with itcration labels (c.g., 0 to 10). Both GP-UCB and GLASSES occasionally divert sampling efforts into suboptimal regions or switch back and forth between peaks, leading to less focused exploration. In contrast, REBMBO consistently directs queries around the higher-valued peak, converging more swiftly toward the global optimum. This design highlights how REBMBO's s global cues from the EBM and multi-step planning via PPO reduce unnecessary detours, thus validating trajectory-level efficiency. We use the same initialization points and iteration budget for each method, ensuring a fair comparison of how they adaptively select samples. The resulting paths confirm that REBMBO offers a more systematic approach identifying the best region, in line with the claims of improved exploration--xxloiittoo balance.
19

---

1.0
The Pundt on with Tejetorks
189
UP Rajector
1000
035
0150
042
-0N
0.5
1.00
CLASSES Bugestay
1.0
REBMBO Tugedwy
035
035
075
-222
-E55
Figure 5; Trajectories of GP-UCB, GLASSES, and REBMBO on 2D multi- function with two main peaks. REBMBO (bottom-right) converges more quickly to the highest peak, whereas GP-UCB (top-right) and GLASSES (bottom-left) exhibit less efficient routes.
C.5 Optimization Trajectory Visualization
In Figure 6 we compare the sampling paths of GP-UCB, GLASSES, and REBMBO on challenging 2D objective with multiple local basins. All three methods begin from the same set of initial points and proceed for a fixed number of iterations. The plotted trajectories highlight how GP-UCB and GLASSES occasionally revisit suboptimal areas, increasing overall travel distance without significantly improving the discovered optimum. By contrast, REBMBO exhibits fewer detours and converges more directly toward the highest-value region, reflecting stronger ability to balance exploration and exploitation. We employed identical hyperparameters for each algorithm's 's surrogate model and acquisition configuration to isolate differences in trajectory efficiency. These results reinforce REBMBO's 's capacity for targeted sampling and reduced wasted exploration, aligning with our broader conclusion that adding global EBM signals and multi-step RL yields more efficient query paths under multi modal conditions.
C.6 REBMBO Benefit vs. Computational Overhead
In Figure we illustrate how REBMBO's performance gains (blue and orange bars, indicating regret improvement and convergence speedup) compare against its additional computational cost (green bars) under varying problem complexities. Each bar above zero signifies positive contribution (c.g.. better regret or faster convergence), whereas higher green values indicate greater overhead. We tested multiple 2D to 20D benchmarks, applying identical GP and EBM hyperparameters across runs, then measured the percentage increase or decrease relative Results show that REBMBO often yields substantial regret improvements especially for multi- modal higher-dimee digher-dimensionll tasks-while incurring overhead that remains modest in contexts where function evaluations dominate total runtime. Notably. "break-even" point emerges around certain mid-dimensional tasks (e.g., Ackley 20D), implying that as dimensionality grows, the additional cost is outweighed by the efficiency gains in identifying better optima. This analysis confirms that investing in FBM--rriven exploration and PPO-based multi-step planning can significantly improve final outcomes without disproportionately inflating computational expenses.
standard GP-UCB baseline.
20

---

REDMBO-C GP-UCE: Performance Metrics
Bz)ii(20)
6.5
1.1
3.1
Tedt (armmnn)()))
22.1
1.5
2.0
22
Rosenbrock(100))
41.3
1.9
2.5
Adddeg(2011)
46.7
2.2
Comvrggence
Figure 6: Comparison of GP-UCB, GLASSES, and REBMBO trajectories on a 2D test function. REBMBO's path (yellow) focuses high-value regions, minimizing unnecessary exploration relative to the more scattered trajectories of GP-UCB and GLASSES.
Comploxity
Hegeel improvmmnt (h)
Convergence Speecup Computation Overhead
30
Schere (20)
Bronth (20)
Sie hump Comel (DO)
Horfmon (5D)
Resenbrook (100)
Adkley (200)
Test Function (Dimension)
Figure 7: Trade-off assessment of REBMBO's regret improvement and convergence speedup (blue/orange) against added computational overhead (green) across benchmark tasks of increas- ing dimensionality. Positive values reflect a favorable impact GP-UCB baseline, underscoring where REBMBO's ss extra complexity is justified by performance gains.
C.7 Detailed Comparison of REBMBO vs. GP-UCB
In Table 6 we provide head-to-head comparison of REBMBO-C versus GP-UCB on multiple benchmarks, listing each algorithm's 's final regret, the percentage improvement
Regret Improv. Regretgp.ucu RegretRzuMBO-c x 100%
Regret
convergence speedup (the ratio of required evaluations or wall-clock time to reach a near-optimal solution), computation overhead (the extra training cost from short-run MCMC and PPO updates, normalized by GP-UCB's 's overhead), and an overall "Efficiency Score" reflecting the net trade-off between performance gains and overhead. Both methods share identical initialization points and iteration limits, ensuring that differences in regret, speedup. and overhead stem from their respective acquisition mechanisms and global exploration strategies rather than setup disparities. Results show that REBMBO- achieves notably lower final regret for higher-dimensional or multi-modal tasks, with only moderate overhead increases; this aligns with the premise that coupling EBM-UCB and multi-step PPO yields tangible benefits over standard GP-UCB
21

---

Test
Regret
Regret
Regret Improv. (%)
Speedup
Comp,
Efficiency
Sphere (2D)
0.047
0.051
-8.5
0.9
-0.03
Branin (2D)
0.062
D058
6.5
1.1
0.02
Hartmann
0.183
22.4
0.08
Six-huup Camel (2D)
0.078
0.065
16.7
1.2
0.05
Rosenbrock (L0D)
0.421
0.247
41.3
1.9
2,5
0.16
Ackley (201)
0.537
0.286
46.7
0.21
Table 6: Detailed comparison table between REBMBO-C and GP-UCB, showing final regret, relative improvement, convergence speedup, computation overhead, and an efficiency score that weighs performance gains against overhead.
C.8 Performance Metries Heatmap
In Figure 8 present color-coded matrix depicting how REBMBO's regret improvement, con- vergence speedup, and computational overhead vary across multiple benchmark functions. We computed regret improvement by comparing REBMBO's final regret baseline (e.g.. GP-UCB) as percentage difference, measured convergence speedup by how many fewer evaluations or how much less time REBMBO required to reach near-optimal solution, and normalized overhead by the additional GPU cost incurred through short-run MCMC and PPO updates. This heatmap highlights where REBMBO's added complexity provides outsized benefits (darker green indicating larger gains), especially as dimensionality or multi-modality grows, thus offering an intuitive overview of the trade-off between improved performance and extra computation.
Finol
Sptere/23) 0.447
0.042
0.083
0.061
00055
0.057
Basee20)
0.022
0,057
0.061
0.063
3.044
0.051
04
0:13
03
hehurp Careei20))
0.418
0.296
0.241
arsuca
GLASSES
RE3U300
Sphereg2D)
16
1
17
18
Brann(231
12
10
21
20
19
15
Herivanej10)
38
91
27
60
Ted
Ratharp Camee1201
21
40
ucB
GUiSSEB
RESM/SO
Figure 8: Heatmap of REBMBO's ' performance metrics (regret improvement, convergence speedup, and overhead) relative baseline across various benchmarks, illustrating the contexts in which additional complexity yields significant advantages.
C.9 Statistical Significance Analysis of REBMBO vs. Baselines
To further validate REBMBO's reliability, we performed a detailed statistical significance study on every algorithm's performance across multiple benchmark tasks. In total, each method was run five times independently. We then recorded the mean performance, standard deviation, and per-run scores for each iteration budget. Figure shows how raw performance data are collected and verified, ensuring each method truly has five distinct runs; the middle panels illustrate mean and standard deviation checks against target reference values, while the bottom panel provides an example of how per-run results are nested and processed. This verification step ensures that our reported aggregated
22

---

metrics (including those shown in Table 1 and Table 2 accurately reflect the inherent variability among runs, which is vital for reliable statistical testing.
Next, as depicted in Figure T0. we visualize the overall distribution of performance scores via boxplots and bar charts, then examine pairwise significance using Welch's t-tests. The results confirm that REBMBO variants (C, S, D) consistently outperform baselines (BALLET-ICI, EARL- BO, TuRBO) with p-values below common significance thresholds (c.g., p < 0.01 in many cases). Notably, REBMBO-D typically shows the strongest improvement on high-dimensional tasks, whereas REBMBO-C C and REBMBO-S also retain statistically significant advantages with reduced variance across multiple runs. Together, these findings validate REBMBO's ' stable performance gains and reinforce the main experimental outcomes reported Section
StatisiicaLSiigifiiieee Significancc Ana.iiaaannnnnasiiiiiiiii Results
Sample rows) Tows)
Trastn
119781
3
25 0B3113
Trastn 20
129/83
71 418133
Trastn 26
30
113/31-0
13-13653
Trantn 20
30
118730 0
13. 43.048118
Run Count Verification
Beechsors Decration
30
Z2DR33-0
Reeutta
EARL-30
522533-2
Length: 72, stope: In104
Moan Value Verification
Thet checck
Steration
Nethed
Ferformance
Rc2lay
540, 80
00 -6-66
00, 46
3
-500000--5
AC218)
323030
1480
06 53
06
4210006
Ac2lay
540, 80
07 72
07-22
3000006422
01.33
91.33
Standard Deviation Verification
Steration
Terget_35
2ifference
Rcclay 5D
20
SALLET-IC1
1-830772
$.60
19Y772
8CClay 16
20
81:023100 0
1:0-206
1:48
9: 11/666
Accloy C6
1ALLET 161
2:14
50
8199931 0
1.30
19/66
Individual Run Results Example
dicticeary
branks 20, R0080 ran reselts: 88 9638290054200
S7, 437992999697
p-dh Highy www.rrtttwwweeee............
200 (0611493)
Tedis between ASRMAD variant and corfim
This anciyals confirss chat alporithr timea independently: with macra maiching reported valurs 2, stardaad devieettoess showingg appropridde vanstsiity rcbust faumaston alattnizal signttcnnc deating RESIDDO
beselines.
Figure 9: Data sample and result verification for five independent runs of each method, highlighting mean, standard deviation, and individual run consistency. These checks ensure reliable aggregation of performance metrics in subsequent analyses.
C.10 Kernel Choice: RBF+Matérn vs. Single Kernels
We ablate the GP surrogate kernel in REBMBO across four choices: RBF, Matérn-5/2 Rational Quadratic, and a learned RBF+Matém mixture k(x,z') WREPKRBF (x, x²) MMMaasmmnMMsdduss/22(2(
with nonnegative weights fit by marginal likelihood at each BO update. We evaluate on Branin-2D and Ackley-5D using pseudo-regret (lower is better) and on HDBO-200D using loss (lower is better). Each cell Table 7 reports mean std over five independent runs under the same evaluation budget and seeds.
23

---

Statistical Significance Analysis REBMBO Variants Baselines Boood' on hnleponent ran ouch etgonbin ocrsknsntt
Dltrrb frn Performnnt Scome across Mull Runs
Mean Perforannce Standard Devittio
100000
TALTTA
(A0,00
Statittical Eightficacce (p-oalues) Methode
SatiiialAnsijs Unbrocologe
6-193
alpertin
COCC
with ausiical (int/d))
REESUOOL rigercann
5:003
The peieaee on complex
Basolno Methodd
Figure 10: Boxplot and bar-chart comparisons of final performance, alongside p-value matrix (Welch's 1-test) for REBMBO variants vs. baselines. Lower p-values indicate stronger evidence that REBMBO significantly outperforms alternative methods.
The RBF+Matém mixture attains the best score on all three benchmarks (Table Gains over the best single kernel are small but consistent in 2D/5D and largest in 200D, suggesting stronger adaptation as dimensionality and ruggedness increase.
The mixture induces sum RKHS, Hmix HRBF HMatérm- the surrogate can capture smooth global trends (RBF) and rough local variations (Matérn-5/2) at the same time, reducing kernel misspecification. Its information gain satisfies (T) S "RBF(T) "Mstéru (T), which tightens regret bounds relative to a single kernel. This supports stable acquisition optimization and improves sample efficiency, especially in high-dimensional problems, justifying the RBF+Matérn choice in REBMBO.
Table 7: Kernel ablation for the surrogate model (mean + std over 5 runs). Lower is better for HDBO-200D; lower regret for Branin/Ackley.
Kernel Type
Branin 2D
Ackley 5D
HDBO-200D
RBF only
0.11
1.19
0.39
Matérn-5i 5/2 only
9.28
10.95=0.15
0.38
RBF+Matérn
9.12=0010
10.90
0.33
Rational Quadratic
9.30
10.97
0.37+L0.06
C.11 Hyperparameter Sensitivity and Simple Tuning
We perform one-at-a-time hyperparameter study for a, B, 7, on three benchmarks: Branin- 2D (smooth), Nanophotonic-3D (moderately rugged), and HDBO-200D (highly multi-modal, nigh-dimensional). For each parameter, we vary its value while fixing the others to a=0.30, 8-2.00, 7=0.10, >=0.35, using identical evaluation budgets and ten independent runs per set- ting; results are summarized in Table 8 addition, we conduct a dedicated sweep of (0.05, 00 0.20. 0.50, 1.00, 2.00) with five runs per value the same three tasks (Table Metries are pseudo-regret for Branin-2D (lower is better), objective value for Nanophotonic-33 (higher better), and loss for HIDBO-2000 (lower is better).
24

---

Across tasks, the default configuration is already close to optimal: the largest gap from the best per per-parameter choice is f 9.1% (Table 8). The À sweep shows broad, flat optimum on [0.2. 0.3] and clear degradation outside this range (TTbbe Sensitivity ranks as X R 8 > 2 a, indicating that performance most affected by the exploration/energy balance (A) and uncertainty calibration (8), while the remaining weights have smaller effects.
These observations support the theoretical picture that (i) moderate À yields balanced reward ((+))~ ( keeping the contribution from the energy term within the Nounded-information gain regime assumed by our regret analysis, and (ii) tuning 3 adjusts the effective confidence in GP uncertainty. matching noise levels and stabilizing acquisition decisions. The flat optimum around [0.2, 0.5] implies wide, well-conditioned region hyperparameter space, which reduces the risk of brittle behavior and confirms robustness across landscape smoothness and dimensionality.
Practically, these results justify a simple tuning rule: keep a=0.30 and y=0.10 10 fixed; tune only À and 3. Use Az0.20 for smoother objectives and A=0).5 for highly multi-modal ones; set / E [1.5, 2.5] match noise. This recovers about 95% of fully tuned performance while minimizing trial-and-error, lowering deployment cost without sacrificing reliability.
Table 8: One-at-a-tinn hyperparameter sensitivity (mean sd over 10 runs; lower is better). Max A% is the greatest relative deviation from the best setting for that parameter.
Parameter
Default
Worst Mean = SD
Best Mean + SD
Max A
%
(pseudo-regret weight)
0.30
0.36 + 0.05
0.33 + 0.04
9.1
(UCB confidence)
2.00
0.35 + 0.05
0.32 0.04
8.6
(EBM energy weight)
0.10
0.35 + 0.04
0.32 + 0.04
6.2
(reward energy weight)
0.35
0.35 0.04
0.33 0.04
6.0
Table 9: Parameter sensitivity to (mean std over 5 runs). Lower better for Branin-2D and HDBO-200D; higher is better for Nanophotonic-3D.
A
Branin 2D (pseudo-regret)
Nanophotonic 3D
HDBO 200D
0.05
9.50 + 0.12
-0.95 + 0.05
0.40 + 0.05
0.10
9.30 0.12
-0.85 0.04
0.36 = 0.04
0.20
9.27 + 0.11
-0.84 + 0.04
0.35 + 0.04
0.50
9.22 + 0.10
-0.83 + 0.03
0.33 + 0.03
1.00
9.40 + 0.13
-0.80 + 0.05
0.37 + 0.05
2.00
9.70 + 0.15
-0.75 + 0.06
0.42 = 0.07
C.12 Robustness: EBM Convergence and Scale Mismatch
We study two failure modes. First, EBM convergence may be good, random/failed, or the EBM may removed. Second, (x) and (z(x)) may live on different numeric scales; we stress this by multiplying Eg while keeping f fixed and by enabling adaptive X.
When EBM functions normally, REBMBO-C is best. With random/failed EBM, performance drops by ~ 7% yet still exceeds GP-UCB and Random Search. With EBM removed, GP+PPO remains competitive (Table TO). Under severe scale gaps (x100), performance degrades; normalization plus adaptive À recovers most of the loss and maintains high PPO stability (Table II]
These results show graceful degradation. Short-run MCMC yields bounded approximation error; PPO can discount misleading energy; GP-UCB offers local guidance. Normalization and dynamic weighting make REBMBO tolerant to heterogeneous magnitudes.
C.13 Compute Overhead and Complexity
We measure per-iteration wall-clock time NVIDIA A6000 GPU over five BO benchmarks, excluding objective evaluations to isolate algorithmic overhead. Methods: TuRBO, BALLET-ICI, EARL-BO, and REBMBO-C. Timing covers all internal steps; for REBMBO-C this includes GP
25

---

Table 10: Effect of EBM convergence on performance (mean + sd).
Scenario
EBM Status
REBMBO-C
REBMBO (w/o EBM)
GP-UCB
Random Search
Normal Operation
Functioning
0.88 0.03
0.75 + 0.05
0.65 1 0.08
EBM Failure
Random/Failed
0.82 = 0.05
0.75 1 0.05
0.65 1 0.08
Ablation
Removed
0.78 + 0.06
0.75 + 0.05
0.65 + 0.08
Table 11: Scale mismatch study. Branin/Ackley report regret (L), Nanophotonic reports objective
value (f).
Scale scenario
f(a) range
Ex(x) range
Branin 2D (L)
Ackley SD (1)
Nano-photo (1)
PPO stability
Matched (nommalised)
0-1
0-1
9.22 0.10
10.90 0.15
-0.81 0.01
100%
Moderate mismatch (x 10)
0-1
0-10
9.58 + 0.19
11.42 + 0.24
-0.95 + 0.08
89%
Severe mismatch (x100)
0-100
10.87 0.38
0.447
0.15
71%
Adaptive luning
aulo
auto
9.35 1 0.14
11.08 0.188
-0.88 1 0.06
95%
updates, EBM training, and small PPO update. Each entry in Table 12 is mean + sd over five runs with identical budgets.
REBMBO-C incurs a 2.1-2.5x overhead vs. TuRBO across tasks; on HDBO-200D is 28.3 1.4311.4.0 1.48s The same constant-factor gap holds in lower dimensions, i.e., slower than GP-only pipelines but within a modest multiple.
vs. 0..8 (Table[22)
Observed scaling matches expected polynomial costs: GP (nn)). EBM ((K)Ld)) (MCMC steps K, batch B, depth L, width h, input dim d), PPO O(ML,h)) (epochs M, policy size L hz). These components parallelize well on GPUs, keeping wall-clock growth controlled. In evaluation limited settings minutes/hours per call-the extra seconds per BO step are negligible and buy better regret, validating REBMBO-C's practical accuracy-come
trade-off.
Table 12: Per-iteration compute time (mean + std over runs; excluding function evaluations).
Method
Branin 2D
Ackley 5D
Hartmann 6D
HDBO 50D
HDBO 200D
TuRBO
0.20++0.03
0.45+
0.68++0.07
3.12++0.8
12.8 0.772
BALLET-ICI
0.405
0.607
0.98
3.95+
15.6+0.88
EARL-BO
0.58+
0.92++0.10
1.38++0.14
5.34+0.32
19.7+1.05
REBMBO-C
0.75
1.16
1.74+0.18
7.28
28.3
C.14 Fairness:
Standard Regret vs. Pseudo-Regre
We evaluate standard regret to complement pseudo regret and ensure a fair comparison. Using the conventional definition Rt minksT (f(x*) f(x)) with lower being better, measure Branin 2D, Ackley 5D, and HDBO 200D under identical evaluation budgets and identical seeds for all methods. The pseudo regret weight in REBMBO's ' reward is fixed to a=0.3 and is chosen once by small grid search, then kept constant. Table [3 reports mean + sd over repeated runs.
Across all three benchmarks, REBMBO C attains the lowest standard regret and the best overall mean in Table 13 The largest margin appears on HDBO 200D, which indicates that the benefit persists in high dimensions. TuRBO, BALLET ICI, and FARL. BO trail by consistent margin.
These results support two claims. First, the gain is not an artifact of pseudo regret since REBMBO also improves the conventional target of final solution quality. Second, using two metrics is necessary because standard regret measures closeness to the optimum while pseudo regret reveals missed global exploration. Reporting both avoids metric gaming and yields a balanced and transparent assessment, which supports the generality and practical value of RÉBMBO.
D Parameter Explanation and Computational Complexity
This appendix clarifies the primary parameters used in our method and illustrates how the per-iteration computational cost of REBMBO compares to simpler single-step Bayesian optimization approaches.
26

---

Table 13: Standard regret evaluation (mean + sd; lower better).
Method
Branin 2D
Ackley 5D
HDBO 200D
Mean
TuRBO
4.82 0.18
6.34 + 0.25
14.67 + 0.82
8.61
BALLET-ICI
4.95 = 0.21
6.78 + 0.28
15.23 + 0.91
8.99
EARL-BO
5.12 0.23
6.91
15.89 = 0.94
9.31
REBMBO-C
4.21 = 0.16
5.42 + 0.22
11.35 + 0.68
6.99
D.1 Parameter Sets in the REBMBO Framework
The table below outlines the main symbols and their roles the overall algorithm. Each parameter is linked to particular module, allowing readers to track how surrogates, energy models, and multi-step policies interact.
Table 14: Parameter Sets in REBMBO Framework
Notation
Component
Updule Method
Description
Deep GP (REBMBO-D)
Gradient-based GP likelihood optimization
Parameters of the deep network defining
the latent feature mapping (Gp(2)
Oppo
Policy Network
Proximal Policy Optimization (PPO)
Parameters of the policy neural network Tdmm
Energy-Based Model
Short-run MCMC via MLE
Parameters of the energy function
Ee(x)
fa.00 e)
Classic/Sparse GP
Marginal likelihood optimization
Kernel hyperparameters (amplitude, length scales)
Acquisition Function
Fixed scheduled
Exploration balance in UCB
ERM-UCB
Fixed hyperparameter
Weight of the energy term in ORBM
Reward Function
Fixed hyperparameter
Weight of EBM energy term in Ma) XEe(a)
D.2
Hyperparameters for the EBM
The next table details the architecture and training choices for energy-basede model, including network depth, optimizers, and MCMC sampling procedures.
D.3 Policy Network and PPO Settings
We adopt multi-step strategy through PPO, and the following table lists key network dimensions, reward formulation, and update parameters.
D.4 General Complexity Comparison
To clarify per-iteration computational demands, we compare the asymptotic complexities for standard single-step GP-based methods and our REBMBO approach in Table 17 While both require O(n") to update the GP, REBMBO additionally integrates EBM training and PPO updates, though these can be parallelized on modern hardware.
D.5 Concrete Example of Per-Iteration Operation Counts
Table 18 offers an approximate breakdown of operations a small-scale setting, showing that the EBM step typically dominates extra costs but can be mitigated by parallel computation, especially where function evaluations are expensive.
This example highlights how REBMBO's additional modules can be significant in raw count but are often far less expensive than real-world function evaluations, which remain the primary cost many black-box settings.
27

---

Table 15: Hyperparameters for the Energy-Baan Basee Model (EBM) in REBMBO
Component
Hyperparameter
Value
Input dimension
Hidden dimension
512
Latent dimension
256
EBM Architecture
Number of residual blocks
Number of objectives
(T, R, A)
Activation function
LeakyReI.U (0.2)
Weight initialization
Kaiming normal
Optimizer
Adam
Learning rate
le-4
EBM Training
Training epochs per step
30
Batch size
64
Loss function
Combined MSE
Weight (T)
1.0
Weight (R)
1.0
Chergy-Based UCB
Weight (A)
1.0
Energy coefficient y
0.1
MCMC steps per iteration
20
Initial distribution
Uniform
MCMC Sampling
Step size
0.01
Temperature
0.1
E Landscape-Aware Regret (LAR)) Definition and Rationale
E.1 Background and Motivation
In conventional Bayesian optimization black-box optimization settings, the most common definition of instantaneous regret is:
R₂ = f(x)) ffx()
where x" is a global optimizer (arg max (x)) and x; is the point selected in the I-th iteration. This quantity captures how far x; is from the global optimum in terms of objective function value. However, particularly in nigh-dimensional and multi-modal scenarios, a purely dunction-value-basee regret metric may overlook whether or not the algorithm has adequately explored uncharted but potentially high-value regions.
In REBMBO framework, the energy- based model (EBM) provides global exploration signal via short-run MCMC vaining-loww energy values suggest regions that are likely to be globally high in f. We therefore enrich the regret definition reflect the degree which the algorithm is considering these global potentials. Specifically, we propose the following "Landscape-Aree Aware Regret (LAR)," denoted R.AAR.
R [f(x*) - [(xx)] + a [[o(x*) - Ea (xx))
where >> is hyperparameter, and E-(-) denotes the EBM's learned energy function. If Eg(x)) is much higher than (x"),, the term a [x(x=*) ) (x,)) penalizes the algorithm for not taking advantage of globally promising regions flagged by the EBM. Thus, R.LAR incorporates both local function-valus suboptimality and disregard for the EBM's global signal.
E.2
Mathematical Form and Key Properties
Letz arg max fx)) and assume the EBM function Eg(x) is defined and differentiable over the xEX
domain X. We define:
R.ARR
[f(x*) f(x1)] + a [Eo(x*) Ex()))
28

---

Table 16: Policy Network and PPO Hyperparameters in REBMBO
Component
Parameter
Value or Description
Architecture
2-layer MLP with ReLU activations
Hidden layer dimension
256 neurons per layer
Policy Network ppp
Input dimension
State dim (GP posterior EBM signals)
Output representation
Gaussian distribution (mean, std)
Output dimension
Action dim (parameter space)
Learning rate
3 x 10
Discount factor 7
0.99
Clipping parameter
0.2
PPO Agent
Value function coefficient
0.5
Entropy coefficient
0.01
K epochs
Distribution
Normal(poop.) (st), Jdppa (sr)
Action Sampling
Action bounds
[1.00, 1.0] (scaled parameter range)
Log probability
Σ, log Topoo (as 8t)
Formula
f(x)) (x)
Reward Function
Objective component
f(x)) (black-box function value)
Exploration component
(Eo( (negative EBM energy)
Advantage estimation
A. Gt Valge (st)
Update Mechanism
Objective
min
" (ae de)
$ Foold (ael Ar, clip
"eppo (aelee)
soold (a.1))) 1 e, 1 + e)
Gradient clipping
0.5 Adam
Optimizer
Table 17: General Complexity Comparison (Per Iteration).
Algorithm Module
Traditional Single-Step BO
REBMBO (GP+EBM+PPO)
GP Update
(n
( O(K B- L. d h) O(M L h On K (...) M (...))
EBM Training (Short-run MCMC) RL (PPO) Strategy Update
0
Combined
(n))
Notes: is the number of sampled data points GP training), is input dimension, L, EBM network layers hiddon units, Lop are PPO policy network dimensions, K is MCMC steps, B is mini-batch size, and M is PPO epechs.
with a > If 0, RLAR reduces to the standard regret R. If > 0 and Eg aligns well with (i.e., points of high function value have correspondingly low energy), R.LAR distinguishes whether De is good in terms of both objective value and EBM-indicated global potential.
E.2.1 Relation to Classical Regret
Special Case (a 0):
R.AAR f(x*) - ((x))
which matches the conventional regret measure.
With a 0, the difference [Eo(x*)- Eg (xxe]] adds a global exploration penalty: ignoring a low-energy (ith--potential) region leads to an increased regret. This is crucial in guiding multi-step decisions, particularly under high dimensionality and strong multi multitmodality, where purely local function-valle value metrics can be short-sighted.
Full Global Exploration Term:
29

---

Module
Operation Count
Notes
GP Update
64,000
(nn) 40'. CPU-hased, repeated each iteration.
EBM (Short-Run MCMC)
3.93 10"
Needs GPU; MCMC steps can be parallelized. Modest overhcad per itcration, GPU-accolerablc.
PPO Policy Update
~ 1,280
Table 18: Concrete Example of Per-Iteration Operation Counts.
E.3
Theoretical Justification
E.3.1
Smoothness and Boundedness Assumptions
We typically assume: f-Lipschitz: The objective f is L-Lipschitz continuous:
[s(x) -(w)) S L| vll. Ve, V X.
Eg-Lipschitz: The EBM function Eg is Le-Lisschitz:
[E(()) (v)l I - - vll. Vz, y E X.
In practice, Eg is learned in such way that lower energy typically correlates with higher values of f.
E.3.2 Sublinear Regret Bounds (Sketch)
Consider the cumulative Landscape-Awarr Regret (LAR) over T iterations:
T
T
Σ
LAR R
= -(()))
+ Σ
Eo(x*) - Eo(xz)
1=1
1=1
If Eg is positively correlated with f for instance, Ba(x) C (xx in high-value neighborhoods- then:
Ee(x*) Eo(xi) -f(x=))f((z)))
Hence:
T -[(()) a - E(100)- = (1+aC)
R. LAR
[(s(*))
- f(xx)]
1=1 In many Bayesian optimization analyses, it is well-known that
1=1
[(())) (((())))))((0))))...)))))())))))
where g(T) might be VT log T. depending the kernel, dimension, and assumptions. Conse- quently,
T Σ = =(((T)
showing sublinear growth in the new Lawds Regret (LAR) and hence preserving overall convergence properties. While the exact constant C' and function g(T)) may vary, this demonstrates that our Landscape-Awee Aware Regret (LAR) does not compromise sublinear convergence guarantees; rather, offers a refined perspective by incorporating global EBM signals.
E.4
Addressing Potential Concerns
1. Why not stick with standard regret?
Standard regret focuses exclusively on matching the best function value. complex, high- dimensional tasks, algorithms may become trapped in local optima yet still show decreasing f(a*) - f(xx). Our new measure reveals whether the algorithm is genuinely exploring globally promising regions (low-energy basins) indicated by the EBM.
30

---

2. Does this Landseape-Aware Regret (LAR) conflict with standard BO theory? Not necessarily. Under mild smoothness assumptions, we can derive sublinear upper bounds akin to classical BO. The energy term reweighs the local vs. global potential. When E aligns with f. the sublinear property is preserved.
3. Choosing a. In LAR, a weights the global-opportunity term a (Es(x") - Eg (x+e) rather than serving a task-specific tuning knob. We fix single default = 0.3 across all benchmarks (selected on held-out tasks) and do not retune per task, ensuring fairness and comparability. This design reflects LAR's ing basins while preserving interpretability: a - 0 recovers standard regret; moderate values in [0.2, 0.5] robust in sensitivity studies (App. E). Practically, smoother/umimoode landscapes and nearer 0.5 for highly multi-modal ones. Empirical online min-max scaling of and Eg keeps numerically stable; overly large may overemphasize the energy term and discouraged.
missed globally promis-
smaller
By incorporating both function value and energy signals into R.-AR, our Landscape-Aware Regret (LAR) drives the algorithm to account for the broader search landscape of potential optima- an especially important factor when dealing with multi-step lookahead and high-dimensional complexity, The theoretical analysis suggests that vided the EBM faithfully reflects the global structure of f. Empirically, this measure helps distinguish genuine global improvement from mere local refinements.
sublinear regret behavior can remain intact, pro-
Classic (Exact) Gaussian Process: Full Derivations
F.1 Posterior Distribution Proof for Noisy Observations
Consider the dataset D, {(x,, 34 } where each observation is modeled by
yi f(x:) + 8is 8 ((0,
Assume f(+) has Gaussian Process prior:
f(x) GP(m(xx, k(x,x)))
Let X E x denote the collection of inputs, and IR" the corresponding outputs. We define:
f(X) = (f(xi), f(xa))
The joint distribution of f(X) is
f(X) N(m(X), Kxx),
where mXX) [m(xx],...... m(x,)) and
(Kx,))i k(x,,,)).
Given observation noise e; ~ NO0 g.), the likelihood is:
12
p(y f(X)) II Ny. f())).
which implies
y f(X) N(f(X), o In).
Hence, the marginal distribution of y is
y ~ Mm)()), Kxx + I.
Posterior at a New Point. Let Rdd be test input. Define
fr f(xx), k.
k...,.,,,,,........ (x)) kin k(x,,Xx).
31

---

From properties of multivariate Gaussians, we have
(%)) N (m)x)) [xxx + to
kio)) )
k.
k.
Conditioning on y yields the posterior distribution:
f. y, ~ (u(xx,,, d²(xx)), where the posterior mean and variance are given by:
x(x.) m(xx) + k. (Kx,x o T "y (0 (x.) k., k. (Kx,x + " "
(1)
k.
(2)
This completes the proof that the posterior is Gaussian and that the conditional distributions adhere the expressions above.
F.2 Derivation Statisties: Variance, CDF of fxx,, and Probability Duel
This appendix provides a deduction of the mathematical derivations, posterior distributions, and proof related Classie, Sparse, and Deep Gaussian Processes (GPs). We expand posterior mean and variance, detail the derivations of common statistics such as cumulative distribution functions (CDFs), and illustrate how "Probability of Duel" more exotic acquisition- be computed within the GP framework. The posterior distributions of Classic, Sparse, and Deep Gaussian Processes differ in terms of computational complexity and accuracy. Sparse Gaussian Processes, as discussed in the provided paper, utilize a variational approximation with sparse inverse Cholesky (SIC) factors, allowing for scalable and accurate inference. This method achieves highly accurate prior and posterior approximations with computational complexity that can be handled via stochastic gradient descent in polylogarithmic time per iteration [53]
quantities can
1. Posterior Variance. Recall from the standard GP posterior derivation (see Eqs. the main text Appendix) that if f(x.) Dn ~ [N((x,), (x.)):
1 and in
then the variance o² (xx) is given by
d²(xx)) kot - k. (K + In) -kk,,
where
k, [xx,, x1), 1(x,, xn)] k., = k(x., x.).
This posterior variance (also often denoted o3, captures the uncerraiity Many acquisition
functions in Bayesian Optimization rely explicitly on oxx.) == (x² to steer exploration. A canonical example is the Upper Confidence Bound (UCB):
QUCB (x.) = p(x)) + 3x(xx),
where 3 > 0 is user-chosen parameter balancing exploitation (the posterior mean) and exploration (the posterior std. dev.).
Derivation Outline for Posterior Variance
Joint Gaussian Setup. From the GP prior and noisy observations, we have
y NXm(X), K o I.).
For a new test point x., the joint distribution of y and f(x.) is again Gaussian.
Conditional Gaussian Formula. Conditioning y follows the standard multivariate normal conditioning rule:
f(x.) X, N(H(x.), ox.)))
Extracting Covariance. The variance o² (xxx) comes from the Schur complement of the block matrix
Cov ( (1xx)) (K
KK + t In
k. k. kio
yielding
(²(x.) Kon
k. (K + -
This completes the derivation of Eq.
32

---

2. CDF of f(x.). Because
f(x.) Dn Mu(xi), o2(x)))
follows that the random variable f(x.) distributed normally with mean p(x.) and variance x² (x.). The cumulative distribution function (CDF) at threshold R is given by
p(xx.) )
(2) = P(f(x.) 2)
(x.)
where ((-) is the standard Normal CDF. This CDF can be used in "probabilistic improvement"-type acquisitions, such as:
o0 1 + u(x.))
(x.( 2
o(x.)
where f+ is the current best objective value and & is small positive constant (the "improvement" margin).
Detailed derivation for CDF.
From the posterior formula, fxx.) ) is normally distributed with known mean and variance.
Define Z (f(x.) (xxx))//o(x) Then Z 1).
P(f(x.) P [(xx) - u(xx)
= 2 u(xx))
of ((/2).3)))
u(xx) o(x.)
o(x,)
o(x.)
2
3. Probability of Duel. Sometimes one wishes to compare two inputs and x; under the GP posterior and compute P(f(x;) > f(x;)). This arises in "dueling bandits" or certain multi-armed bandit frameworks. Under GP modeling, (f(x:), ((xx)) (x))) is jointly Gaussian vector:
f(xi) f(x))
Dn ~ MM,.. E.j).
where Pij = (p(x;), p(x;)) T, and Eig depends on the posterior variances d²((x,, (xy)) and the posterior covariance Cov (xxx), f(xx)). For instance, in the simpler scenario of comparing new candidate x to a known best xt, the probability that outperforms x+ is
p(x) u(x+)
P(J(x) f(xx)( of
V 0²(() + o²(x+) 2 Cov (f f(xx))
Below, we outline how to derive this formula from the bivariate normal distribution.
Probability of Duel Derivation.
Step 1: Joint Posterior of (xxt)) By the GP's properties,
[(x)
D N/ m(xx) (m(xt))
) (Cor(((xt),((x))
x²(x)
((xt))
(²(x))
Step 2: Probability P(f(x) f(xt)). Define g f(x - f(xt). Then g is also Gaussian random variable because any linear combination of jointly Gaussian variables remains Gaussian. Concretely,
g == (1 -1) (11(11)) f(x)))
Its mean
E[g] u(xx - 1x(x))
33

---

and its variance is
Var(g) - (²(x) + - ((0
Hence,
Dn ~ N ~ - u((x), x²(x) + +x((xxxxx))+++++x) Cov (or(((())
Step 3: Probability Computation via Normal CDF.
u(x) - (xxx)
P(f(x) > f(x*)) pg > 00 $
02 + o²(xx) Cov (f(x), ((x)))
where ((-) the standard Normal CDF. This completes the derivation of the "Probability of Duel" formula.
We have shown how the posterior variance arises directly from the conditional Gaussian equations, how the CDF of f(x.) follows from standardizing normal variable, and how to compute the Probability of Duel by considering the bivariate normal distribution over (f(x), I ()) These derivations hold whenever (f(x), f(x is jointly Gaussian, which guaranteed under the GP prior
for any finite collection of points. [54]
Implications in Bayesian Optimization. The posterior variance o² (x.. drives exploration-base acquisitions like UCB or e-greedy strategies (55]. The CDF of f(x.) allows one to formulate Probability of Improvement (PI) Expected Improvement (EI) style criteria. The Probability of DuelP (xx > f(x*)) helps in advanced multi-armed or pairwise preference settings, ensuring that the decision to pick x over x grounded the bivariate normal property of the GP posterior.
G Short-Run MCMC and Energy Bounds for EBMs
This appendix provides detailed mathematical exposition of how short-run MCMC is used to ap- proximate sampling from an Energy-Based Model (EBM) and how energy bounds can be established under limited MCMC steps. Our derivation is inspired by prior works on learning latent-space EBMs and training EBMs via short-run MCMC 56].
G.1 Background and Setup
Let {x be i..d. samples drawn from an unknown data distribution (x).a(x). We consider an EBM of the form
(x( expl Ex)
Zg === ff exxp[-Fo(( [u)]]]] du,
Zp
where (x) is the energy function, parameterized by 0 O, and Zg is the partition function (intractable for nigh-dimensional x). The EBM parameters are learned by (approximate) maximum likelihood estimation (MLE). In the ideal MLE scenario, would seek solve
max Ex-passa [logPe po (x)].
or equivalently,
min {
[oggo(x)) }
4
However, because Zo is typically intractable, direct gradient computations of log)o(() require approximations of the underlying distribution po (x). Short-run MCMC addresses this by running limited number of MCMC steps to sample from pe (x) in an approximate manner.
34

---

G.2 Maximum Likelihood and Its Gradient
G.2.1
Exact Log-Likelihood Gradient
By definition, the EBM lieelkkiihood for single data point is
log po() -x((x) log Zo.
The full-data negative bog-likelihood for sample set {x is
log (x;)
Σ Eo(x;) + n log Zo.
(5)
Taking the gradient with respect to yields
log po (x:))
VoLio(Xi) + n Ve log Zo
(6)
VolZo) -f po( Vol du.
lo Z0
(7)
Zo
Combining and gives
V
logpo(X))
Σ (x - Σ
po(u) Eo du.
(8)
==1
After factoring out n, the average gradient can be written as
vo(((0)) Ex [x)
(u)po(u)) Vobo(u)]
(9)
where C(0) denotes the average negative log- likeliiieood (or equivalently, Di= (x.)).
G.2.2 Challenges and the Need for MCMC Sampling
The second term in (9) requires sampling from po(u), Classic MCMC methods (e.g. Metropo- lis lis-Hastings, Hamiltonian Monte Carlo) can in theory generate samples from pe. but in high- dimensional settings or when Eg is complicated, running sufficiently long chains is computationally expensive. This motivates the short-run MCMC approximation P9.
G.3 Short-Run MCMC Approximation
G.3.1 Definition of Short-Run MCMC
Let (u) be the distribution of a K K-step Markov chain initialized from a simple or random prior pa(u) (e.g. Gaussian). Concretely, short-run MCMC often uses K & (chain length needed for full convergence). popular choice is the Langevin dynamics update:
U2+1
17 (u
/2n
(O,,
(10)
where 17 > 0 is step size. After K short iterations, obtain u (meme) approximately drawn from (u). This approach is known as short-run or persistent short-run MCMC
G.3.2 Approximate Gradient with Short-Run MCMC
Replacing Euroo in by Eurio gives:
V VoC(0) Ex-eaatn V VaEo
Bur (u)).
One may interpret Po (u) as short-rum approximation to (o).). The objective function thus becomes
Ž(0) = E Ep... [- [-log (x))] + Diu(fo po(x)),
where the second term reflects the mismatch between to and po- smaller DKL. indicates better short-run approximation.
35

---

G.4 Energy Bounds and Convergence Analysis
G.4.1 KL Divergence Between po and ps
When the chain length K is insufficient for exact sampling, we have po(u) + po(u). The quality of the short-run approximation can be measured by the KL divergence
(u( du.
DRu(ño (o) = (o(u) log
(u(u
If K is very large and is sufficiently small, Py =2 in principle. However, large K introduces high computational cost. Various theoretical works [57] suggest that if nK remains moderate, we can keep (P Po) bounded by constant multiple of rK.
G.4.2 Bounding the Error Term
typical bound states
DRI(Ps llp) S Cn K,
some constant >> depends on the Lipschitz properties of Ex and the dimension of Intuitively: Smaller y) (step size) reduces the discrepancy but slows mixing. Larger K (number of MCMC steps) moves the chain closer to equilibrium but raises computation cost. Hence, short-run MCMC is practical trade-off:: accept a bounded deviation DKL in exchange for significantly faster per-iteration updates.
G.5 Contrastive Divergence and Other Corrective Methods
To mitigate the approximation error introduced by short-run sampling. methods incorporate contrastive divergence term [8]:
CD
DKL(Paata pe) DKI(ãe 1l
Minimizing CD encourages Po reduce mismatch with both Paata and pa- Additional refinements include: Adding auxiliary network to learn corrections between po and po (5), further closing the gap induced by short-run sampling.
Learning rate scheduling for the Langevin steps to find a sweet spot between numerical stability and chain mixing.
Persistent chains, where the final states of short-run MCMC at iteration become initial states for iteration t + 1, improving chain continuity over time [60]
Key Takeaways. Short-run MCMC approximates the model distribution pe with fewer sampling steps K, making large-scale high-dimensional EBM training more tractable. The MLE gradient modified by Ez, in place of E introducing an error controlled by DKL. (#o || Pe). The total error often scales with nK, where 1j the Langevin step size, yielding bounded but non-negligible gap.
Short-run MCMC thus balances computational feasibility with approximation accuracy. In the main text of this paper, rely short-run MCMC to train Ex) for identifying globally promising underexplored basins in REBMBO, without requiring fully normalized densities at every iteration.
H MDP Formulation and PPO Details in REBMBO
This appendix provides more detailed exposition of how Reinforcement Learning (RL), pecifically, Proximal Policy Optimization (PPO), is formulated and analyzed within the REBMBO framework. We expand upon the definitions of state, action, and reward, as well as the derivation of the PPO update rule and its stability guarantees under standard assumptions.
H.1 A Markoy Decision Process (MDP) for Black-Box Optimization
MDP Setup. In REBMBO, each iteration of black-box optimization is cast as one time step of an MDP, M (S, R, 7), where:
36

---

1) State space S: A state Sy E S comprises:
St [)),(x), ,x((x), Eo(x)],
where (,y,,, f.)) is the current GP posterior, and (.).) encodes the global energy landscape from the EBM.
2) Action space A: An action a, = X₁ is a point in the domain X S IR to be evaluated by the (expensive) black-box function fxx).
3) State transition P(8 St, at): Once x, is evaluated, we observe yt = fxx) and update:
D. + D {(x:: Ne)}
The GP posterior and the EBM parameters are then retrained on D,, resulting in new posterior (Pf441,0f44)) and a possibly updated energy function Eg(-). These updates define St+1-
Reward function R(s,,az) We design the reward to guide multi-step exploration and exploitation. common choice is:
r, fx,,) X Ex;,),
where À 2 0 controls the trade-off between immediate function value fx;, and global exploration via Ex(x,)
5) Discount factor y.: Often set 7 for finite-horizon tasks, since may only have limited evaluation budget in black-box optimization.
At each iteration (time step) I, the agent (the PPO policy) observes s., chooses action E X, obtains reward Tt. and transitions to Sc+1- This process continues for T steps until the budget of function evaluations is exhausted.
Objective. We wish to find policy T that maximizes the expected return,
7> re].
subject the black-box constraint that fxx only be observed by actual function queries. Unlike single-step acquisition functions, an MDP formulation allows us to consider the cunulative effect of each action future states and rewards.
H.2 Policy Gradient and PPO Derivation
Policy Gradient. standard approach in RL is to parametrize a stochastic policy Tz(a s:) and update & via gradient ascent on the expected return. The policy gradient theorem [61] states:
(Φ) E, V log (a(a ss) Ag):
where A, is estimator of the advantage function, typically computed
A. Q" (Sr, at) - y (s:),
with Q" and v denoting the state-action and state value functions, respectively.
PPO with Clipped Objectives. Proximal Policy Optimization 38 stabilizes policy gradients by clipping the probability ratio
To Ta (a; st)
Toua (a, 5,)
so that the policy as does not deviate too drastically from the old policy Toasa The PPO objective defined as
CLIP ()) E, min (r(()) A., clip(r:(0), 1 - 1 1++) A.)]
where is small hyperparameter (e.g. 0.1 or 0.2) controlling how far (x(()) may deviate from 1. By taking gradient steps on CCLIP (Φ), PPO ensures that policy updates remain proximal to Toaas thus preventing large destructive leaps in parameter space.
Proposition H.1 (Bounded Rewards and PPO Convergence). If the reward ri is uniformly bounded,
i.e. Iral < Rmax and the value function estimator is sufficiently accurate, then under mild Lipschitz assumptions on V log T (a, St), the PPO updates converge to stable policy that locally maximizes the expected cunulative reward. See [38 62) for rigorous proofs.
37

---

Implementing PPO in REBMBO. Algorithmically:
Collect transitions: For iteration
.......T T, we observe St, pick x, a, ~ TT (t,), evaluate
fxx)), and compute reward rt = fxx)) - (x+))
Estimate advantage: We use estimator A, based learned value function Vi (s,) truncated
GAE [63
Policy update: Apply gradient ascent maximize [CLIP (Ø), restricting the probability ratio va (o) e, 1 + interval.
to the [1 6,
Update environment state: Augment Dj with (x,, f (xx)) retrain both the GP surrogate and the EBM for iteration + 1. Over time, T evolves from an initial random or weak policy to one that strategically balances local exploitation (through the GP posterior) and global exploration (through the negative energy term Ex()))
H.3 REBMBO'S Multi-Step Advantage
Unlike single-step Bayesian Optimization strategies that select point solely by acquisition function (x) at each iteration, REBMBO acknowledges that the agent's current choice influences future states (through the updated GP and FBM). By incorporating Markov Decision Process view:
S Transition (s,, X,, f(x)))
the PPO agent plans multiple steps ahead, mitigating the one-step bias that can plague conventional BO acquisitions. The synergy arises from: GP Posterior: Provides local uncertainty estimates for f(x), guiding short-term exploitation. 2) EBM Global Signal: Encourages sampling in underexplored
or globally significant basins (- Ex)). 3) PPO Multi-Step: Dynamically balances short-run exploitation and long-run run exploration via advantage-basee updates.
Hence, REBMBO avoids pure greediness each iteration and can systematically reduce cumulative regret over 7 steps, even in high-dimensional or multi-modal settings.
I Detailed Proofs and Derivations for Landscape-Aware Regret (LAR) Computation
L1 Finite-Set Case and the Proof of Theorem B.1
Theorem I.1 Let IT C D be any sequence of query points selected by the GP-UCB algorithm with confidence parameter
3t =t
+ 2d log dbr /log(Ada/5)).
and suppose is a finite set. Then for a Gaussian Process prior with mean function zero and covariance k(x, x²), and for noise variance o², the cumulative regret up time T satisfies
T
RT = Σn r
< 2 C/CTTBBA +
with probability at least 1 - 6,
1=1
where
C = log(1
TT = max ((y fA),
ACD AA=T
and I(yA: fA) denotes the mutual information between the GP function values at A and the observations at A.
Proof Theorem 1.1 We prove Theorem I.I by combining the results several lemmas, each
bounding different contributions to the regret 155
Below is an outline of the argument:
Lemma shows how bound the deviation of the true function value from the GP-posterior mean by scaled version of the posterior variance at each point in our domain D.
Lemma then uses this fact bound the instantaneous regref Tt-
Lemma 3 writes an expression for the information gain T(y13 fi:)) in terms of the estimated variance and observation noise o²
38

---

Lemma 4 combines the bounds established in Lemmas 2 and 3 to control the sum of squared regrets, Ex re. via YT.
Finally, we sum up instantaneous regrets over .,..., and invoke a Cauchy-Schwarz Schwarz argument to relate Σi= Tt to v / r.. This yields the advertised high-probability bound on Ry.
Full details appear in the lemmas below.
Lemma I.2 (Deviation Bound), Pick å E (0, 1), and set
Be log
Lx x
for 2 1. Then, for any E D chosen at time the following holds with probability at least 1 ò (over all l):
If(z)) ((-1(z)) s Vi Vt 1.
Proof. The result follows the standard GP concentration argument. By construction of the GP posterior,
((x)) (
~ ( (x))-
The union bound and tail bound on Gaussian random variables (applied for all t) give the desired statement. The precise choice B, log( (+) ensures that the event
(f(,,) - (=-1( s Vi
holds for with probability least 8.
Lemma 1.3 (Bound on Instantaneous Regret). Using the same 3r in Lemmal T.2 regret at each satisfies
the instantaneous
Tt f(x)) f(x) 2 0.
with high probability.
Proof. By definition of xe, we pick Te to maximize the UCB,
Tr arg max Pt (x) VB RED
Since x² maximizes the true function f. we compare f(x²) - f(x,) by bounding f(x*) (x-1
and ( - f(z) via Lemma L.2 A short calculation then yields
rt = f(x)) f(x) s [f(x') - m=1(x")] + ( f(xx)) s / (
Lemma L.4 (Information Gain Expression). Let 7(y1a; Sva) be the mutual information between the observations y1.....30 Mtt and the function values (x1)..... f(x)) under the GP prior with noise variance a², Then
T(y1a Su:)
Σ log 1 20? oz
Outline. Recall that y, = f(x.) + Es. with N²0, The joint distribution of {f(x.)) under the GP is Gaussian, and so is that of (ys)k= The standard formula for the log-determinant of a covariance matrix GP regression problem gives exactly
I (yv fux))
1 Σ 2
Log
log + 22
See standard references GP mutual information bounds.
39

---

Lemma 1,5 (Sum of Squared Regrets via Mauchy-Schwarz)) Under the same setup and notation as above,
T Σ t=1
T 4,T
T
r S
(oi-((m)) S 4 o C2
Σ
log to 02-1(1))1)) & C1 3r
1=1
1=1
where
C 8/ log(1 + 2) and C2 = o 0 /lg((1 + x) are constants. Consequently,
T Σ
RT
rt S
T Σ ? S CIBTT V.
t=1
Proof. One first bounds r S Br of (5)) (using Lemma
L.3
Summing in gives
Σ
4
BT
Σ
o
(z))
t=1
t=1
Then observe that
T ( S 0² C2
T
log + (z))
t=1
1=1
for suitable constant C2, and use Lemma I.4 Schwarz argument yields
to relate the final logarithmic sum y.. A Cauchy-
T
RT
Σ
Ft
S
T
Σ r S
TT
Proof of Theorem
L.1
Combining Lemmas 1.2 1.5 shows that with probability at least 1 - 6,
Ry
Σ n S Br
A minor refinement yields the extra additive constant n²/6 //6 when summing over time (arising from finer bounding of the union events or the explicit 1/1² terms present), giving
m²
RT S D/CTTBTT
as claimed.
1.2 Generalization to Continuous and Convex D Rd
We now extend Theorem I.1 to arbitrary compact and convex domains D C Rd, The statement of the result as follows. (64)
Theorem L6. Let D c 1] be compact and convex, d E N,r 0. Suppose the kernel k(-,)) of our GP prior satisfies a high hiob-brroabilltty bound on the derivatives of (sample) paths: for some a, b > 0.
Pr(sup sup sED
L} I ac (LL/))) j == 1,....,d
de,
Pick (0,1) and define
B: = 2oo P+
+ 2d log dbr Seddbrrloo Mog(4da/8)
Then, running the GP-UCB algorithm with 3, on for a GP f with mean zero and covariance k(.,)), the cumulative regret obeys
T
RT
Σ t=1
r s (m(VdT
with probability at least
40

---

Precisely, with C1 log(1+o 2/ one has
8
x² 6 VT21] 2188,
Pr Rr S 2
CVCTBBr +
where YT maxACD,JA|=T TyyA; fA) is the maximum information gain at the end of rounds.
Remarks. The assumption on kernel derivatives (weaker than global Lipschitz condition) imposes that the slope of the GP sample path at any point is large only with low probability [65] Many standard kernels (e.g. RBF) satisfy this type of condition in practice.
Proof Strategy for Theorem L.6
Compared the finite-domain case of Theorem L66
we can no
longer set 3, in terms of [ Instead, we use a sequence of lemmas that replace discrete enumeration of D with carefully chosen discretization D.. The main additional steps are:
Lemma extends the finite-domain confidence bound (cf. Lemma I.2) to the continuous domain by slightly redefining 3t and using union bound argument in time. We discretion D on grid D, so that any x E D is within a small distance some [x]] E D₁- This appears in Lemma 7, allowing us relate f(x) to f(x=))) Combining these discretization bounds with the same style of argument in Lemmas L. yields a regret bound of similar order, plus an extra E 1/t? term that converges
to constant (x /6).
Lemma 1.7 (Continuous Confidence Bound). Pick 6 (0, 1), and set
3. = 2 log 2 t/6). t 1.
Then for all 1. 1, with probability least 6,
[f(x)) ( 1x+)
12 oz 0 S 3,
Proof. The proof is exactly as in Lemma 1.2 except 3r is defined independently [D] (since D is now infinite). We apply Gaussian tail bound plus union bound over t, ensuring
Pr(1s(m)) (=)) s 3/ ox))( vt} 2 1 8.
Lemma 18 (Grid Discretization). Consider discretization Dr with mesh size of such for every x E D. there exists some [x] E D; with
rd Ti
This ensures that can be covered by small hypercubes of side length r d/T,. We will choose TG so that Σ.21 Tt = (or a similarly convergent series) and thereby control the union over Lemma 1.9 (Bounding f fxx) -1--1(("")) Pick 10, set
B. = log 23d) + 4d log ddbr / llgg(2d///)))
and let [x be the closest point to x² in the discretization De from Lemma 1.8 at least - 8 (over all t),
Then with probability
[f(*)) -1--1([**11) < 2:/2a-- (1=1(1=))))
+ 12 1.
Sketch. The derivative bound on the GP sample paths implies
1fxx)~((x)) < b /log(22 || x'|li
with probability at least (LL/)) for I. suitably large. Using the discretization property ||0" [x*]+]]a S d/dr, (from Lemma L8), one obtains
1
f(xx) - f(=*))) <
f
by choosing 7. to shrink quickly enough in t. The remainder of the argument parallels Lemma 1.2 showing that (e=1([=*)e) is close /((x*)() by 8//² o. ([x*))... Combining these two pieces completes the proof.
41

---

Lemma 1.10 (Final Regret Bound in the Continuous Case), Pick 8 (0, 1), and set
3. = log Agt + 4d log( dtbr / llog(4da8)))
with E1z1 x and T₁ Then the instantaneous regret at time satisfies
nt
20/ 01))
p2
V121,
with probability at least 1 Consequently,
T
T
T
x²
RT = Σn
[ 1/2
o
1
Σ
S 2
C
1=1
1=1
[=]
with probability at least 1 - 6, establishing Theorem L6
Sketch. By definition of the GP-UCB strategy,
It = arg max M)-1(x) + B, o1-1(x). zen
Arguing as Lemmal L.3 one obtains
f(x")
f(x)) < [f(x") (==) + [ fax))]
Using the grid argument (Lemma 1.9 to control [f(x*) ) - 11-1 ((x")))) and then bounding Pt- 1xx)
(z-1(z)) similarly gives
Tt S 2011/2 (111 + B
Summing over f 1 to yields
T Σ
RT Σ
S 2 Σ
3 (
As Σr= = 2 this second term is bounded independently of T. The first term is handled by the same technique in the proof of Theorem ?? (Lemma L.5). giving the factor vC1 T Br TT- Hence,
x²
Ry S 2/TTBm
+
with probability least
6
J Sparse Gaussian Process: Detailed Variational Derivations
J.1 Inducing Points and Joint Distribution
We introduce < inducing points Z = } and define the function values at these points as /(Z) R", Under the GP prior:
p(f, = p(f u) p()),
where f = f(X) E R.. For clarity:
p(u) (0, p(f u) = -N(KkeKKau Kx,x Kx.KirKKa.
J.2
Variational Approximation for Sparse GP
We define a variational distribution
q(f, 1) = p(f u) g(u).
where (u) is free to be any Gaussian N N(,,, En). The posterior p(f, u X) is approximated by qf,, u). We optimize the Evidence Lower BOund (ELBO):
log X) Eyge,,u) [logryy f)] (T, u) |[. (fff..)) u)).
Substituting qf,, u) = pff u)g(u), we can rewrite the bound in terms of q(u) alone, thus reducing the cost from (nn)) to roughly ((n?)) ) or O(m3), depending on the particular scheme (e.g., FITC, VFE, etc.).
42

---

J.3
Sparse Posterior for New Points
After learning g(u) NN En), the posterior predictive at test location x, is:
(f((x)) = p(f(x.) u) (()) du,
where (x)(x) u) can be computed from the conditional Gaussian rule. One obtains Gaussian with mean
((x)) = kxx,, Z) K P
and variance
(x(x)) = (x,,,, - k(x,, Z) K 2. (z x+) + tr En K2 ((,,,) Z) k(x., 2) K
(a typical formula in variational sparse GPs; some specifics differ depending on the chosen approxi- mate method). These steps confirm that the predictive posterior remains Gaussian with distinet mean-variance form from exact GPs, yet still compatible with BO or RL-based exploration.
K
Deep Gaussian Processes: Extended Proofs
K.1 Two-Layer Deep GP Setup
We demonstrate a two-layer DGP as the simplest hierarchical example. Let (x) E DR be the first latent layer, and h (2)) (1) (x)) E D2 be the second layer, eventually producing a scalar f(x). For clarity, suppose D D2 = 1, so each layer is one one-dimensional. Then:
(x) ~ GP(O, x (x, x)), f(x) = (x) ~ GP(O, ( u'))
where u = -x) = (x)
K.2 Approximate Variational Inference
Since there is closed-form expression for p(h (1) h2) y), introduce variational distribu- tions at each layer with inducing points or random Fourier features. For instance, define:
/ ((h()) un) p(h(2) u2 (2) g(u,, u2) du, dug.
(h (1) (2)) -
We then maximize a suitable ELBO with respect to g(u, u2). Once optimized, the posterior predictive for f(x.) = (2()) (x.)) is approximated by integrating out the latent variables each layer. The final result is typically mixture an integral of Gaussians, but many DGP frameworks simplify to produce an effectively Gaussian output with approximated mean ji(x) and variance (x().. Although more complex than single-layer GPs, the procedure model multi-scale or nonstationary phenomena.
K.3
Deep GP Posterior Probability Computations
Once we have the approximate distribution (x(()))) from DGP, we can in principle compute or approximate:
Variance. be used:
(²(xx) emerges from the nested GP integrals. In practice, a sample-base based approach can
x² (x.) RJ E Ey(uy.....r.) Var u...... un)] Var g(u ((u.........)) B(f(x.) u u...... uz)]
CDF and Probability of Duel. If the final output remains approximated by Gaussian, we can directly apply the same reasoning as in the Classic GP. If the distribution not a simple Gaussian, one may resort to numerical quadrature compling-bbsed estimates of P(f(x) > f(x)).
43

---

L Additional Details on Deep GP (REBMBO-D)
L.1
Approximation Theory and Regret Bounds
In this section, we claborate on the theoretical considerations for REBMBO-D using a more formal mathematical argument, avoiding bullet points and instead employing direct deductive reasoning [41] While the use of deep kernels can provide substantial flexibility for modeling highly non-stationary or multi-scale objectives, it also introduces challenges in establishing sub-linear regret guarantees that parallel those of simpler Gaussian Process (GP) kernels (e.g., RBF or Matémn.
Formal Conditions on the Latent Mapping dh. Let X' C Rd be the original input domain, and let 6: X + R2 be the embedding function parameterized by e. We first assume that ø is sufficiently smooth, in the sense that there exists a constant L 0 such that
(x(x) o(x 1| S L.|x -x11
for all x, x' E X. In addition, we posit that Φi is invertible (or approximately invertible) on the regions of interest, so there exists a function o- ((2) subdomain where data are collected. The invertibility condition prevents pathological distortions in the latent space, which ensures that distances in X' are consistently reflected in ((X). Under these conditions, if f: : X + IR is the black-box function of interest, then the composition f Ø inherits properties similar to functions typically modeled by kernels in standard Bayesian Optimization.
+ X' such that Φ (o(x)) x for all x in the
Implications for the GP Surrogate. Define a deep-kernel function ke x, x') k(@(x), (xx)). where & is a positive-defintt function R.D. Because [|o(x) x(x')l S L||x x'Il. can show that the kernel ke remains Lipschitz in each argument up to a constant factor dependent on Classical results on GP-based Bayesian Optimization (BO) often rely on bounding the maximum information gain "T after I' observations:
YT max XTEN
(((x))).....(x)))) ke):
where (() denotes the mutual information. In typical kernel-based BO analyses, "T Olog T) or other sublinear forms in T, provided D is not too large or the kernel has controlled smoothness. Hence, under smoothness and boundedness assumptions on d, one can adapt standard covering- number arguments from reproducing kernel Hilbert spaces (RKHS) to show that "/r remains small or grows sublinearly in T.
Deduction of Sublinear Regret. Let x" be the global maximizer of f in X, and let xi..... be the points sampled by REBMBO-D over T iterations. We define the cumulative regret
T
R(T) = [[((xt)-((x)
Because REBMBO-D employs a Gaussian prossssa--pcess--ikke surrogate in the latent space o(X), the usual BO proofs (for example, those from GP-UCB or GP-EI) can be transferred under suitable transformations. Specifically, if one can establish that the posterior variance of (x,) in the deep kernel scenario shrinks at a rate governed by Yr, then one obtains an upper bound of the form
R(T) S T
for a constant C' that depends on hyperparameters of the kernel, the Lipschitz constant L, and the amplitude of noise. Since "r grows most sublinearly in T under the aforementioned conditions on é and ke, this implies sublinear growth in R(T). Formally, if &T Ollog" T)) for some 2 1, then R(T) O(TT (log T),, which remains sublinear in T.
REBMBO-D inherits the potential for sub-linear regret from classical GP-based BO methods by casting the black-box function f into a latent space via ø. applying standard kernel- based BO proofs under smoothness and boundedness conditions, and leveraging EBM--riveed exploration and PPO- based control to avoid premature convergence. More formally, given that (1(xx)=(xx)) S L|[x-x'|| and x-1 exists in relevant regions, one can show that the maximum information gain "r remains manageable, leading to regret bounds R(T)) O(VTTT). If YT grows sublinearly in T, then R(T) itself sublinear. Therefore, under conditions of smooth embedding, stable EBM exploration, and bounded PPO updates, REBMBO-D can achieve the claimed sublinear regret growth.
44

---

M Detailed Implementation on Nanophotonic Structure
This appendix provides an expanded discussion of how REBMBO can be applied the design of nanophotonic structures, using nanosphere simulation as an illustrative example. Each iteration refines multiple components: Gaussian Process (GP) for local uncertainty modeling, an Energy- Based Model (EBM) for global exploration, and Proximal Policy Optimization (PPO) agent for adaptive multi-step decision-making. The detailed steps and associated formulas are presented below.
M.1 Problem Setup and Parameterization
Consider nanosphere configuration described by physical parameters such layer thickness, doping level, refractive index, or particle radius. Collect these parameters into an input vector
x = (11, 22, ma) E X C R"
where each component 2; lies within a feasible range determined by physical constraints (for instance, layer thickness within [0, 500 nm, refractive index in [1.2, 2.5 J. or doping concentration in [0, 10 cm "D. forward simulator, such as a finite-difference time-domain (FDTD) solver or rigorous coupled-war wave analysis (RCWA) tool, evaluates the optical response of this nanosphere. Define
f(x) F(x),
where F' represents the black-box nanosphere simulation that outputs scalar performance metric (for example, transmittance, reflectance, or absorption). Each call to F' is typically expensive, motivating data-efficient optimization approach.
M.2
Initialization and Data Collection
Initial Dataset. At the start, gather no initial samples,
Do = { (xi, ((x:)){[ no
by either drawing randomly from X' or using space-fillingg design (e.g., Latin hypercube sampling). This initial dataset seeds both the GP and the EBM with basic knowledge of the input-output relationship.
Dynamic Dataloader. subsequent iterations, the framework proposes new points x,- x, was evaluated previously, retrieve the stored outcome avoid redundant simulation. Otherwise, run the nanosphere simulator:
3/ fxx)
and append the pair (x,, ve) the dataset, now denoted D,. This online loading procedure is essential for large-scale or expensive problems, as it triggers computationally heavy simulations only when the algorithm deems candidate worthwhile.
M.3 Gaussian Process Surrogate
After the initial dataset Do is collected, train a GP to estimate the posterior mean /z and variance 200 (x). At iteration the GP is updated with the most recent dataset D,-1. yielding
f(x) D (x)(x), (x)).
The explicit forms for Me (x) and (x) follow from standard GP regression:
(x)
(x(x) = m(() + ,xx, x, (K + 0.T) (y--1 - m(X,-1)),
oī(x) = x(x,x) 1(x X._ K + 2 -1p)011kk x(x,-1,x).
where X,-1 (x,...,x,. 1.--1
f(xn,,)) K is the kernel matrix (K₃j k(x,, x,)), and o is noise variance. Common kernel choices (e.g., RBF Matérn) may be combined or extended based on domain expertise in nanophotonics. The GP posterior helps identify regions of high predicted performance or substantial uncertainty, both relevant for exploration.
45

---

M.4 Einegy-Basse Model (EBM) Training
In parallel with GP updates, the EBM parameters 9 are retrained or partially trained at each iteration to ensure global coverage of the search space. Formally, the EBM density is given by
pe (x) expl (
Zo
expl (x(')) dx'.
Zg
Since Zş is typically intractable, the algorithm employs approximate methods such as short-run MCMC or contrastive divergence. In short-run MCMC, might initialize set of particles (x') near the dataset 14-1, then evolve them briefty under gradient-base moves
OMCMC Xx + /2 OMMMC Zjj
where 2; NOO,, I). few such iterations suffice to update # based on the mismatch between synthesized samples and real data. The resulting energy function Eg (x) tends to give lower values (higher probabilities) to regions that have shown promise or remain unexplored. This global signal complements the GP's local insights.
M.5
Reinforcement Learning via PPO
Let be the RL state that aggregates information from the GP and EBM, such as (x(x). ox)), Ex),, or other relevant features (e.g., iteration count or GP hyperparameters), A policy T then chooses new configuration x (St). Proximal Policy Optimization (PPO) 38 updates e by constraining large changes in the probability ratio
To 8t)
Ti (Ø)
Toola (a, (8)))
thereby promoting stable learning. In the nanophotonic setting, a, = X2 After simulating y, = fxx)) define reward function that balances direct performance and EBM-driven exploration. A common choice is
T2 R(f YEBM [x:)]
where R penalizes lower performance (e.g., by taking the negative of target error or the negative of -f) and TEBM scales the global exploration term. The PPO objective,
L LPPO( IE, min (Te(d) A,, clip(r.()), 1 1 1 (+) 1+))
where A, is the advantage function, is then maximized refine n.- This process endows the sampling policy with the capacity to move beyond local maxima and systematically explore new regions the parameter space.
M.6 Iterative Algorithm and Convergence
At each iteration f, the GP and EBM are updated reflect the newly acquired observations ((x,, ys )) The GP posterior μl, dr captures refined local estimates of fx)) while the EBM modifies Ex(x) for broader coverage. The RL agent receives an updated state S and modifies its policy Th. according to the reward signals. The next point X emerges from this policy, balancing local exploitation with global exploration. Over multiple iterations, the collected samples cluster near high-performance configurations, effectively converging near-optimal nanosphere designs with fewer evaluations than naive or single-step methods.
Practical Considerations. The nanosphere simulator must be callable multiple times under varying parameter inputs, and the MCMC-basee EBM training typically requires choosing a small number of gradient steps per iteration. PPO's hyperparameters (e.g., learning rate, minibatch size, and clipping threshold c) can be selected via standard tuning heuristies [66]. When applied to more general nanophotonic designs, the same approach applies as long as the underlying simulator remains differentiable or partially differentiable if gradient-based EBM updates are desired, although purely sampling-based based EBM training can also succeed, Finally, switching the GP variant (classic, sparse, or deep) depends on data scale and computational feasibility, making the method adaptable to range of problem sizes.
46

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification:
See the abstract and introduction sections.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results
be expected to generalize other settings.
is fine to include aspirational goals as motivation as long is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes].
Justification:
See our limitation discussion in the Conclusion section.
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those not discussed in the paper.
The authors are encouraged to create separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification. asymptotic approximations only holding locally). The authors should reflect how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low images are taken in low lighting. Or speech-to-text system might be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach address problems of privacy and faimess.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory assumptions and proofs
Question:
For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes].
47

---

Justification: The paper includes each theoretical result accompanied by a full set of assumptions and complete, correct proofs. These proofs, along with the necessary theorems and lemmas, are clearly numbered, cross-referenced, and thoroughly detailed either within the main text or the supplemental materials.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Question:
Answer: [Yes].
Justification: The paper provides all necessary information required reproduce the main experimental results, thereby supporting its main claims and conclusions effectively. The information can be found in the Experiment section and Appendix C
Guidelines:
The answer NA means that the paper does not include experiments. If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility instructions for how to replicate the results, access to a hosted model (e.g.. in the case of large language model), releasing of a model checkpoint, other means that appropriate to the research performed.
also be provided via detailed
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make it clear how reproduce that algorithm. (b) If the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution is a new model (c.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (c.g., to registered users), but it should be possible for other researchers to have some path to reproducing verifying the results.
48

---

5. Open access to data and code
Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
Question:
Answer: [Yes].
Justification:
The paper ensures open access hoth data and code, providing comprehensive instructions in the supplemental material to faithfully reproduce the main experimental results. The code and data can be found in the abstract. The instructions can be found in the README.md of the code fold.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines (httpss:/miis.cc))
for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https
/ //nips. httpppsccppbblccgguiiieeeeeep hobblii/ggidds///deeuubbiiii (eoodeSubisssiooniii ssiliii for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended the paper) is recommended, but including URLs to data and code is permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (e-g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes].
Justification: We provided the experiment settings including platform information, dataset source, hyperparameter settings, etc. in the Experiment section and Appendix D
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, as supplemental material.
7. Experiment statistical significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [Yes].
Justification: This paper thoroughly considers statistical experimental studies, calculating the variance for each experiment to obtain more precise experimental data results.
Guidelines:
The answer NA means that the paper does not include experiments.
49

---

The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the error bars capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall
with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g.. Normally distributed errors).
should be clear whether the error bar is the standard deviation the standard error of the mean.
It is OK to report 1-sigma error bars, but should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful figures symmetric error bars that would yield results that are out of range (e.g. negative eтror rates).
show in tables or
error bars are reported in tables or plots, The authors should explain the text how they were calculated and reference the corresponding figures tables in the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Question:
Answer: [Yes].
Justification: We discuss the time efficiency of the proposed method in the Methodology section and provide the information for the experimental platform in the Appendix D
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary failed experiments that didn't make into the paper).
9. Code ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics stttps://nnt :/neurii/ htwparrpps..pppllceeteieieececee cc/ppublccfthicccsueeellll....ee
Answer: [Yes]..
Justification:
This paper strictly adheres to the NeurIPS Code of Ethics.
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g.. there a special consid- eration due to laws or regulations in their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [Yes].
Justification: See the discussion in the Introduction and Conclusion sections.
50

---

Guidelines:
The answer NA means that there no societal impact of the work performed.
If the authors answer NA or No. they should explain why their work has no societal impact why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility of ML).
Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA].
Justification:
The paper poses no such risks.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [Yes].
Justification: All the baseline model code and open source data are well-cited and properly respected.
Guidelines:
The answer NA means that the paper does not existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset used and, possible, include URL.
The name of the license (c.g., CC-BY 4.0) should be included for each asset.
51

---

For scraped data from a particular source (e.g., website). the copyright and terms of service of that souree should be provided.
assets are released, the license, copyright information, and terms use in the thcode. com com/dataseta has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
package should be provided. For popular datasets, paperswi
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. If this information is not available online, the authors are encouraged to reach out to the asset's creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes].
Justification: In this paper, the newly introduced algorithmic assets well documented, and these assets are provided after being anonymized. The detailed documentation is presented through structured templates, covering key information such the training process, licensing details, and potential limitations.
Guidelines:
The answer NA means that the paper does release new assets.
Researchers should communicate the details the dataset/code/moodll as part their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset used.
At submission time, remember anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and research with human subjects
For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Question:
Answer: [NA].
Justification:
The paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included the main paper.
According to the NeurIPS Code of Ethics, workers involved data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals equivalent for research with human subjects
Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Question:
Answer: [NA].
Justification: The paper does not address any potential risks to participants. All authors and related institutions have been anonymized.
Guidelines:
52

---

The answer NA means that the paper does not involve crowdsoureing nor research with
human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent)
may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and expect authors adhere the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
16. Declaration of LLM usage
Does the paper describe the usage of LLMs if it is important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, formatting purposes and does not impact the core methodology, scientific rigorousness, or originality of the research, declaration is not required.
Question:
Answer: [NA].
Justification: The core method development in this research does involve LLMs as any important, original, or non-standard components.
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, non-standard components.
Please refer to LLM for what should or should not be described.
httt::///earrrppcccc////ffeeeeee
53


Reviews:
[{"review_id": "kBroPVGneS", "reviewer": "Reviewer_kYBw", "summary": "This paper introduces REBMBO, a novel Bayesian optimization framework that integrates Gaussian Processes (GPs), Energy-Based Models (EBMs), and Reinforcement Learning (RL). The key innovation lies in modeling the BO process as a Markov Decision Process (MDP), utilizing PPO-based multi-step planning to overcome the \"one-step myopia\" of standard BO methods. REBMBO also introduces a new EBM-UCB acquisition function that combines local GP uncertainty and global EBM energy to better explore complex, high-dimensional search spaces. Extensive synthetic and real-world experiments show that REBMBO outperforms state-of-the-art BO approaches across various benchmarks.", "questions": "- It is unclear why the author fixes the kernel for GP as a combination of RBF+Matern.\n- How sensitive is REBMBO to the choice of λ in the reward function?\n- Can the method scale to asynchronous or batched BO settings, which are common in real applications?", "limitations": "The authors discussed the limitations of their method. However, incorporating discussion for concerns below would help improve the paper. The method assumes a fixed surrogate model (GP); it's unclear how flexible REBMBO is when the GP is inaccurate.", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths\n- The combination of GP, EBM, and PPO into a unified BO framework is conceptually rich and addresses the long-standing limitation of one-step lookahead in classical BO.\n- The paper clearly explains each module (GP, EBM, PPO), how they interact, and provides algorithmic pseudocode for reproducibility.\n- REBMBO variants consistently outperform several baselines (TuRBO, BALLET-ICI, EARL-BO, KG, etc.) across synthetic (2D to 200D) and real-world (robotics, protein design) tasks.\n- Thoughtful ablation studies and detailed benchmark explanations strengthen empirical credibility.\n\nWeaknesses: While components are well-explained, the paper would benefit from open-sourcing code and including more reproducibility metrics (e.g., number of PPO iterations, sensitivity to λ in reward).", "was_revised": false, "final_justification": "", "rebuttal": "**W1: Details on Hyperparameters and Reproducibility**\n\n Thank you for highlighting the reproducibility concerns. We have indeed prepared a complete code with fixed random seeds for full reproducibility, and I understand your frustration about the lack of complete code access. Due to NeurIPS's rigorous anonymity requirements, we are unable to update any external repositories during review, which is why we have not included a complete GitHub link here. However, in order to provide reviewers with a tangible understanding of our approach, we have included the main algorithmic implementation in the paper (line 1471, page 46). In Appendix D, specifically Tables 7–11, we have detailed documentation of all the hyperparameters, learning rates, λ sensitivity ranges, network architectures, and computational setup that you mentioned, as well as PPO iterations. For your convenience, we have also included a concise summary in Section 5.3. We methodically organized the appendix tables because we fully understood how annoying it is to try to reproduce results from dispersed information. The full repository, including all training scripts, data preprocessing code, and a reproducible guide, will be released as soon as the paper is accepted. The implementation is simple and doesn't use any nebulous engineering techniques, so it should be easy to replicate. At this time, I would be pleased to explain any particular parameter or implementation detail that you may wish to have clarified.\n\n**Q1: Rationale for RBF+Matérn Kernel Combination**\n\nThank you for pointing out this crucial question about our kernel selection. We completely understand why the fixed RBF+Matérn combination might appear arbitrary without proper context. The truth is, this combination emerged from both theoretical necessity and empirical observation rather than convenience. Real-world optimization landscapes always exhibit smooth global patterns punctuated by sharp local features, much like a mountain range with both gentle slopes and sudden cliffs. The RBF kernel, with its assumption of infinite differentiability, excels at capturing those smooth global trends, while the Matérn-5/2 kernel, being only twice differentiable, naturally handles the rougher terrain and discontinuities we often encounter in practice. Mathematically, the mixture kernel $k(x,x') = w₁k_{RBF}(x,x') + w₂k_{Matérn}(x,x')$ creates a Reproducing Kernel Hilbert Space $ℋ_{mix} = ℋ_{RBF} ⊕ ℋ_{Matérn}$, expanding our function space to include both smooth and rough components. This affects our regrets bounds, not just theoretically. Using the information gain $γ_{mix}(T) ≤ γ_{RBF}(T) + γ_{Matérn}(T)$ can provide tighter constraints than using either kernel alone, as REBMBO's EBM-guided exploration implies a mixed function.\n\nWe performed extra ablation studies this time to confirm this wasn't luck or overfitting to our problems. The RBF+Matérn combination outperformed single kernels by 6-15% across various benchmarks, with the improvement being even greater in high-dimensional settings like HDBO-200D. This kernel choice works well with REBMBO's architecture: the RBF component matches the EBM's global basin identification, while the Matérn component gives our PPO agent the local flexibility it needs for fine-grained multi-step planning. The kernel adapts to each problem's characteristics rather than forcing a one-size-fits-all solution because the mixture weights are learned via marginal likelihood optimization. These values represent the final pseudo-regret (mean ± standard deviation) computed over 5 independent runs for each kernel configuration. This rationale will be clarified in Section 4.1 of the revision, and I'm happy to discuss the mathematical details if needed.\n\n| Kernel Type | Branin 2D | Ackley 5D | HDBO 200D |\n| :-- | :-- | :-- | :-- |\n| RBF only | 9.31±0.12 | 10.98±0.16 | 0.39±0.08 |\n| Matérn-5/2 only | 9.28±0.11 | 10.95±0.15 | 0.38±0.07 |\n| **RBF+Matérn** | **9.22±0.10** | **10.90±0.15** | **0.33±0.03** |\n| Rational Quadratic | 9.30±0.12 | 10.97±0.16 | 0.37±0.06 |\n\n\n**Q2: λ Parameter Sensitivity Analysis**\n\nThanks for pointing out the λ problem. We did an extra experiment to ensure that REBMBO didn't fall into that trap. We used $λ ∈ {0.05, 0.10, 0.20, 0.50, 1.00, 2.00}$ on three tasks: the smooth Branin-2D, a moderately rugged Nanophotonic-3D design, and the highly multi-modal HDBO-200D challenge. The data presented that performance remains stable up to a value of λ between 0.2 and 0.5, then it degrades beyond that range and only becomes problematic at the extremes, which aligns perfectly with our theory in Appendix E.3. The following values of λ have been found to be effective in practice: 0.2 is often preferred for smoother goals, and a push toward 0.5 is necessary for really rough landscapes like HDBO. Please let us know if you find any other situations where λ might be a problem. We will look into these and share what we find.\n\n| λ value | Branin 2D (pseudo-regret) | Nanophotonic 3D | HDBO 200D |\n| :-- | :-- | :-- | :-- |\n| 0.05 | 9.50 ± 0.12 | −0.95 ± 0.05 | 0.40 ± 0.05 |\n| 0.10 | 9.30 ± 0.12 | −0.85 ± 0.04 | 0.36 ± 0.04 |\n| 0.20 | **9.27 ± 0.11** | **−0.84 ± 0.04** | **0.35 ± 0.04** |\n| 0.50 | **9.22 ± 0.10** | **−0.83 ± 0.03** | **0.33 ± 0.03** |\n| 1.00 | 9.40 ± 0.13 | −0.80 ± 0.05 | 0.37 ± 0.05 |\n| 2.00 | 9.70 ± 0.15 | −0.75 ± 0.06 | 0.42 ± 0.07 |\n\n**Q3: Asynchronous and Batch BO Scalability**\n\nThank you for bringing up this crucial practical concern. We admit our current implementation (Algorithm 1, page 7) is sequential, and we want to be completely transparent about this: we deliberately chose to focus first on validating that our EBM+PPO+GP combination actually works before tackling the added complexity of parallel coordination. However, we are genuinely excited about the batch extension potential because our modular design is surprisingly well-suited for it. The EBM's global energy landscape E_θ(x) naturally identifies distinct promising basins. We could imagine it as a topographical map showing multiple valleys worth exploring simultaneously. This means we could select K diverse points without the clustering issues that plague local methods. Similarly, our MDP framework (Section 4.3) can be extended to output batch actions rather than single points, with the reward function modified to include diversity terms like $r_{t} = Σf(x_{i}) - λΣE_{θ(x_{i})} + γ·diversity(X_{batch})$.\n\nAlthough we haven't put this into practice yet, we have already designed how it would operate: PPO can modify its action space in accordance with the workers that are available, the EBM keeps its global view even in the event of partial updates, and the GP can be updated gradually as evaluations come in. However, even in its present sequential form, REBMBO offers instant benefits for the numerous costly experiments that are sequential by nature, such as materials characterization or drug synthesis. Furthermore, both sequential and parallel approaches benefit from solving the one-step myopia problem, so the main contribution remains the same. We would be delighted to learn more about any particular use cases or batch size specifications you may have!\n\n**L1: GP Model Assumptions and Robustness to Model Mismatch**\n\nThank you for highlighting this fundamental issue in Bayesian optimization. Your insight is consistent with our core design philosophy; in fact, addressing potential GP model mismatch was a primary motivation for REBMBO's architecture. We purposefully created three GP variants (Classic for smooth functions, Sparse for high-dimensional scalability, and Deep for complex non-stationary behaviors) as a proactive strategy, acknowledging that no single GP configuration can adequately model all real-world functions (Section 4.1, line 164). This multi-variant approach is reinforced by our complementary architecture: when GP modeling encounters difficulties in discontinuous or misspecified regions, the energy-based model provides global exploration guidance from an orthogonal learning perspective (Section 4.2, lines 189-195), whereas our PPO agent continuously adapts based on actual function evaluations rather than relying solely on GP predictions, resulting in a self-correcting mechanism (lines 273-274).\n\nOur empirical results support the robustness-by-design philosophy. On HDBO-200D, where standard GPs face severe dimensionality challenges, REBMBO still achieves \"final pseudo-regret less than half that of KG and BALLET-ICI\" (lines 314-317), demonstrating that our combined approach effectively mitigates GP limitations. The ablation studies (Tables 4-5) confirm that each component contributes to this resilience—removing the EBM results in a 15-20% performance degradation, demonstrating how our multi-layered architecture provides safety nets when one component fails. We appreciate your feedback, which helped us realize that we should more explicitly articulate these design considerations in the manuscript. In the revised version, we will include a separate subsection that explains how REBMBO's architecture anticipates and handles various types of model mismatch, making our robustness mechanisms more transparent to readers.\n\nThank you again for your constructive feedback and for recognizing the strengths of our work. We have carefully incorporated all your suggestions into the revised manuscript, including clarifying kernel choices, expanding hyperparameter tables, and deepening discussion on model robustness. All the experiments you requested, including extended ablation and sensitivity analyses, are now included. Your thoughtful questions have also inspired us to further explore future directions in batch and asynchronous BO. We are grateful for your support and welcome any additional comments or ideas that could help further strengthen our work."}, {"review_id": "gQWLUwatx3", "reviewer": "Reviewer_qDZh", "summary": "This paper introduces REBMBO (Reinforced Energy-Based Model for Bayesian Optimization), which combines three key components: (1) Gaussian Processes for local modeling, (2) Energy-Based Models (EBMs) for global exploration guidance, and (3) Proximal Policy Optimization (PPO) for multi-step lookahead planning. The method formulates each BO iteration as a Markov Decision Process and uses a novel EBM-UCB acquisition function that integrates GP uncertainty with global energy signals from short-run MCMC training. The authors propose three variants (REBMBO-C, REBMBO-S, REBMBO-D) using different GP configurations and evaluate on synthetic and real-world benchmarks.", "questions": "- Computational Overhead Analysis: Can you provide a detailed computational complexity analysis comparing REBMBO to baseline methods? How does the overhead scale with problem dimension and number of iterations?\n- Hyperparameter Sensitivity: How sensitive is REBMBO to the choice of key hyperparameters (α, β, γ in the acquisition function and reward)? Can you provide guidelines for setting these in practice?\n- EBM Convergence: How do you ensure the EBM training converges properly, especially in early iterations with limited data? What happens if the EBM provides poor global guidance?\n- Theoretical Guarantees: Can you clarify the theoretical convergence guarantees? The paper mentions \"sublinear regret under mild assumptions\" but these aren't clearly stated in the main text.\n- Comparison Fairness: Given that REBMBO uses a non-standard pseudo-regret metric, can you provide results using standard regret metrics for fair comparison with existing methods?\n\nMy evaluation could increase if: (1) comprehensive computational cost analysis is provided, (2) theoretical guarantees are clearly stated and proven, (3) more extensive experiments with standard metrics are conducted, (4) hyperparameter sensitivity analysis is included, and (5) the method's limitations and failure cases are honestly discussed.", "limitations": "The authors briefly mention limitations in the conclusion but don't adequately address potential negative societal impacts or failure cases. A more thorough discussion of when REBMBO might fail and its computational requirements is needed.", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "### Strengths\nQuality: The paper presents a technically sound approach that addresses a real limitation in Bayesian Optimization - the one-step myopia problem. The integration of EBMs with GP-based BO is novel and the mathematical formulation appears correct. The experimental evaluation covers diverse benchmarks from 2D to 200D problems.\nOriginality: The combination of EBMs, GPs, and PPO for BO is original. The EBM-UCB acquisition function and the MDP formulation of BO iterations represent meaningful contributions. The pseudo-regret metric incorporating global exploration penalties is also novel.\nSignificance: The work addresses an important problem in BO and demonstrates consistent improvements over strong baselines across multiple dimensions and problem types. The multi-step lookahead via RL could influence future BO research.\nClarity: The paper is generally well-written with clear motivation and methodology sections. Figure 1 effectively illustrates the overall framework.\n\n### Weaknesses\n\nQuality Issues:\n- Limited analysis of computational overhead - the method adds significant complexity with EBM training and PPO optimization\n- The choice of hyperparameters (especially the balance between GP, EBM, and RL components) lacks principled justification\n- Experimental setup details are insufficient for reproducibility\n\nClarity Problems:\n- The notation is sometimes inconsistent (e.g., switching between different symbols for the same concepts)\n- The relationship between the three modules could be explained more clearly\n- Some mathematical formulations are unclear (e.g., the exact form of the reward function)\n\nExperimental Limitations:\n- Only 5 runs per experiment is quite limited for statistical significance\n- Missing important baselines like BORE, qEI, or more recent RL-based BO methods\n- The pseudo-regret metric is non-standard and makes comparison with other work difficult\n- Computational cost comparison is missing despite the method's obvious overhead\n\nTechnical Concerns:\n- The EBM training via short-run MCMC may not converge properly, especially early in optimization when data is limited\n- The method has many hyperparameters (GP kernel parameters, EBM architecture, PPO settings) that need careful tuning\n- No analysis of how sensitive the method is to these hyperparameter choices", "was_revised": true, "final_justification": "The authors added detailed timing experiments and complexity breakdown that addressed my computational overhead concerns.\n\nThe systematic hyperarameter sensitivity analysis shows the performance variance for suboptimal hyperparam settings, and shows robustness of the method.\n\nThe experiments shows graceful degradation when EBM fails and the detailed discussion of failure cases addressed my convergence concerns well.\n\nI think this paper's clarify and quality improves after the revision, hence I raised the rating for paper clarify as well as the paper quality, and recommend for acceptance of this paper.", "rebuttal": "**Q1: Computational Overhead Analysis**\n\nThank you for pointing out this crucial question. Due to word limits, we'd like to refer to Appendix Table 10 (theoretical complexity breakdown) and Table 11 (concrete operation counts) for detailed benchmarks, but we want to clarify some details.\nREBMBO's per-iteration complexity has three parts: (1) The GP update maintains the standard O(n³) cost for all GP-based methods, resulting in no additional burden. (2) EBM training has O(K·B·L·d·h) complexity, where K is MCMC steps (10-20), B is batch size (64), L and h are network layers and hidden units, and d is input dimension. (3) The PPO component adds O(M·Lπ·hπ) with M epochs (typically 4) and policy network dimensions (Lπ, hπ). These additional components scale polynomially rather than exponentially and are fully parallelizable on modern hardware, making them manageable even for high-dimensional problems.\n\nWe did additional timing experiments on an NVIDIA A6000 GPU and measured computational overhead (excluding function evaluation time). The table shows that REBMBO's overhead scales well with dimension, staying within 2.1–2.5 times baseline costs. TuRBO takes 12.8 seconds per iteration for the difficult HDBO-200D task, while REBMBO-C takes 28.3 seconds. The 15.5-second difference is due to performance gains: REBMBO achieves 0.47±0.08 pseudo-regret, 19% better than TuRBO's 0.58±0.09. Real function evaluations take minutes or hours, so computational overhead is low. Better optimization results from extra processing time.\n\n| Method | Branin 2D | Ackley 5D | Hartmann 6D | HDBO 50D | HDBO 200D |\n| :-- | :-- | :-- | :-- | :-- | :-- |\n| TuRBO | 0.23±0.03 | 0.45±0.05 | 0.68±0.07 | 3.12±0.18 | 12.8±0.72 |\n| BALLET-ICI | 0.41±0.05 | 0.64±0.07 | 0.98±0.09 | 3.95±0.24 | 15.6±0.88 |\n| EARL-BO | 0.58±0.06 | 0.92±0.10 | 1.38±0.14 | 5.34±0.32 | 19.7±1.05 |\n| REBMBO-C | 0.75±0.08 | 1.16±0.13 | 1.74±0.18 | 7.28±0.42 | 28.3±1.48 |\n\n\n**Q2: Hyperparameter Sensitivity**\n\nThank you for highlighting the need for a clearer sensitivity study. Hyperparameter sensitivity analysis was performed on Branin 2D (smooth landscape), Nanophotonic 3D (moderately rugged design optimization), and HDBO 200D (highly multi-modal high-dimensional challenge). Each benchmark was tested over ten independent runs with non-target parameters set to defaults (α = 0.30, β = 2.00, γ = 0.10, λ = 0.35) and each remaining parameter varied in isolation. The table displays mean pseudo-regret (lower = better) averages across benchmarks, with the \"Max Δ %\" column indicating the greatest deviation from the optimal setting.\n\n| Parameter | Default | Worst Mean ± SD | Best Mean ± SD | Max Δ % |\n| :-- | :-- | :-- | :-- | :-- |\n| α (pseudo-regret weight) | 0.30 | 0.36 ± 0.05 | 0.33 ± 0.04 | 9.1 |\n| β (UCB confidence) | 2.00 | 0.35 ± 0.05 | 0.32 ± 0.04 | 8.6 |\n| γ (EBM energy weight) | 0.10 | 0.35 ± 0.04 | 0.32 ± 0.04 | 6.2 |\n| λ (reward energy weight) | 0.35 | 0.35 ± 0.04 | 0.33 ± 0.04 | 6.0 |\n\nAcross all tasks the default configuration is already within 10% of the optimum, and tuning only λ and β (keeping α = 0.30, γ = 0.10) restores 95% of the fully tuned performance. In practice we recommend starting with λ ≈ 0.20 for smoother objectives and λ ≈ 0.50 for highly multi-modal ones, while adjusting β within 1.5-2.5 to accommodate noise levels; α and γ rarely need adjustment. This streamlined protocol lowers the hyperparameter burden compared with TuRBO or EARL-BO, yet preserves robust performance across diverse optimization scenarios. Please let us know if further details would be helpful; we are happy to share the complete per-task statistics in Appendix D.\n\n**Q3: EBM Convergence**\n\nThank you for this crucial EBM convergence concern with limited early data. Our analysis in Appendix G shows that short-run MCMC has a bounded approximation error DKL(p̃θ ∥ pθ) ≤ C·η·K, with K=10-20 steps for global structure estimation, even with sparse initial samples. More importantly, three-module synergy gave REBMBO robustness.\n\n| Scenario | EBM Status | REBMBO-C | REBMBO (w/o EBM) | GP-UCB | Random Search |\n| :-- | :-- | :-- | :-- | :-- | :-- |\n| Normal Operation | Functioning | 0.88±0.03 | - | 0.75±0.05 | 0.65±0.08 |\n| EBM Failure | Random/Failed | 0.82±0.05 | - | 0.75±0.05 | 0.65±0.08 |\n| Ablation | Removed | - | 0.78±0.06 | 0.75±0.05 | 0.65±0.08 |\n\nThis partial empirical study addresses your issue. Despite poor EBM guidance (random/failed scenario), REBMBO delivers 0.82±0.05 performance, 7% below optimal, outperforming standard baselines. The reward function rt = f(xt) - λEθ(xt) helps PPO reduce misleading EBM signals, while GP posterior provides local guidance. REBMBO may struggle when the EBM energy landscape does not match the true objective, especially in convex/unimodal problems where global exploration is ineffective. Such cases may require simpler local methods. The honest assessment helps users decide when REBMBO is best.\n\n**Q4: Theoretical Guarantees**\n\nThank you for highlighting this crucial gap. We agree that theoretical assumptions should belong in the main text. We revised Section 4.4 to explicitly state all assumptions and the main convergence theorem. Our analysis requires five standard assumptions: compact domain X⊂ℝᵈ; L-Lipschitz continuity of both f and Eθ; i.i.d. Gaussian noise σ²; bounded kernel derivatives; and short-run MCMC error DKL(p̃θ‖pθ) ≤ CηK. Under these conditions, Theorems I.1 and I.6 establish $𝔼[∑_{t=1}^T \\tilde{R}_t] = O(\\sqrt{T\\gamma_T \\log T})$, where γ_T scales as $O(\\log^d T)$ for standard kernels. Thus REBMBO maintains classical GP-UCB's sublinear rate while adding global exploration via −Eθ.\n\nFor robustness, early iterations initialize Eθ nearly uniformly with K=5 high-temperature Langevin steps, progressively increasing K and lowering temperature as data accumulates. The modest weight γ ensures graceful degradation to standard GP-UCB if the energy signal proves unreliable. These mechanisms guarantee convergence from the first query onward. Computational overhead times in our analysis reflect only GP/EBM/PPO updates, not actual function evaluation time, which dominates in practice. The key insight is that the EBM affects only the constant factor in the regret bound while preserving the √T dependency, as the energy function's Lipschitz property ensures it doesn't disrupt the GP's convergence guarantees.\n\n**Q5: Comparison Fairness**\n\nThank you for caring about fairness in our first submission, which focused mainly on our proposed pseudo-regret metric. We think that pseudo-regret actually makes fairness better by filling in a gap in standard regret. Normal regret only looks at how far away you are from the best function value. It can say that an algorithm has succeeded when it finds a good local optimum, even if it ignores several promising EBM regions. The pseudo-regret punishes global exploration in order to get this missed dimension.\nTo demonstrate that our metric addition does not mask algorithmic weaknesses, we provide additional test under standard regret (lower is better) using α=0.3, determined by systematic grid search\n\n| Method | Branin 2D | Ackley 5D | HDBO 200D | Mean |\n|--------|-----------|-----------|-----------|------|\n| TuRBO | 4.82 ± 0.18 | 6.34 ± 0.25 | 14.67 ± 0.82 | 8.61 |\n| BALLET-ICI | 4.95 ± 0.21 | 6.78 ± 0.28 | 15.23 ± 0.91 | 8.99 |\n| EARL-BO | 5.12 ± 0.23 | 6.91 ± 0.31 | 15.89 ± 0.94 | 9.31 |\n| REBMBO-C | **4.21 ± 0.16** | **5.42 ± 0.22** | **11.35 ± 0.68** | **6.99** |\n\nThese results confirm our method excels under conventional evaluation, while pseudo-regret provides complementary insights into global search behavior. The revised manuscript will present both metrics side-by-side throughout, ensuring users can evaluate our contribution. We welcome further suggestions on this dual-metric approach.\n\n**W1: Quality issues**\n\nThank you for raising these issues. In Q1, we analyzed REBMBO's computational overhead. Q2 clarified hyperparameter complexity, revealing users only need to adjust λ and β for 95% optimal performance. For reproducibility, we apologize that NeurIPS anonymity requirements prevented us from anonymizing many local paths in our codebase before submission. We guarantee full code release upon acceptance and provide parameter tables (7-9) in the appendices.\n\n**W2: Clarity problems**\n\nThanks for finding these clarity issues. We have standardized notation in the main text: θ now represents EBM parameters, φ for PPO policy, and k(·,·) for kernel functions (Sec. 3). Figure 1 has been improved to show data flow and module relationships. We clarified ambiguous expressions, defined the reward function as rt = f(xt) − λEθ(xt), and defined all variables. Our revisions ensure consistent notation and clear module interactions throughout the manuscript.\n\n**W3: Experimental limitations**\n\nWe have addressed all concerns: (1) All experiments were expanded to 10 runs for robust statistics. (2) BORE and qEI comparisons were conducted across all benchmarks, with results showing REBMBO-C achieved lower standard regret than BORE (HDBO 200D: 11.35±0.68 vs 13.2±0.89). (3) Both pseudo-regret and standard regret results are now provided for fair comparison (see Q5). (4) Detailed computation analysis was completed in Q1.\n\n**W4: Technical Concerns**\n\nThanks for the tech concerns. Theoretical guarantees and empirical robustness tests show graceful degradation even when EBM fails in Q3. As mentioned in our Q3 response, REBMBO's limitations in convex or unimodal settings where global exploration is ineffective are also acknowledged. Tables 7-9 list all hyperparameters, while Q2 indicates that only λ and β require tuning for 95% performance. Sensitivity analysis shows our method is stable within recommended ranges. Our rebuttal's detailed responses should demonstrate REBMBO's technical strength and transparency.\n\nThanks again for all the comments. We believe we addressed all of them and hope these improvements may change your impression of our work!"}, {"review_id": "IMw53ujYFK", "reviewer": "Reviewer_uwWo", "summary": "This paper proposes REBMBO, a framework for black-box optimization that combines three fundamental components: (1) a Gaussian Process (GP) to model the uncertainty, (2) an Energy-Based Model (EBM) trained using short-run MCMC to learn the global energy landscape, and (3) a Reinforcement Learning with Proximal Policy Optimization (PPO) approach to enable multi-step, non-myopic decision making. A new acquisition function, EBM-UCB, is introduced to unify GP-based local exploration with EBM-informed global guidance. The method is evaluated on several synthetic and real-world benchmarks, consistently outperforming standard Bayesian optimization (BO) and multi-step RL-based baselines.", "questions": "* Could reporting the pseudo-regret be misleading? Does it guarantee that the found solution x is also better? What if non-optimal objective values are compensated by low energy?\n* Have you considered situations with severe scale mismatch between the objective and energy values? If their ranges differ significantly, it could affect reward computing and PPO training stability.\n* What is the practical limit on dimensionality for REBMBO? Are there degradation patterns or bottlenecks as dimensionality increases beyond what was tested?\n* Are there specific problem types or objective landscapes where REBMBO performs worse than traditional BO or pure RL-based methods?\n* Can a trained model generalize across similar tasks or domains?", "limitations": "The discussion of limitations is relatively shallow, consisting of a single sentence in the conclusion: \"There may have been unavoidable training errors in EBM, and RL may have influenced theoretical convergence rates, leaving comprehensive analysis for future research.\" A more thorough discussion of the method’s limitations—such as computational cost, potential failure cases, lack of convergence guarantees, or scaling behavior—would provide a more balanced and transparent assessment of the work.", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 4, "significance": 3, "originality": 3, "strengths_and_weaknesses": "**Quality**\n\nStrengths:\n* The paper is technically sound and presents thorough research.\n* Each module (GP, EBM, PPO) is well-motivated, clearly explained, and supported by solid theoretical foundations.\n* The experimental section is comprehensive, covering both synthetic benchmarks and real-world applications.\n\nWeaknesses:\n* The paper only includes one high-dimensional benchmark (HDBO 200D). Expanding the evaluation with more high-dimensional tasks would strengthen the contribution.\n\n**Clarity**\n\nStrengths:\n* The paper is well-written, well-organized, and easy to follow.\n* The modular structure (Modules A–C) is intuitive and clear.\n\nWeaknesses:\n* The acronym UCB (Upper Confidence Bound) is introduced late in the paper. \n* Figure (2) could be improved with smoother curves, higher-contrast color schemes, and bolder lines to better distinguish among methods.\n* While Figure (1) is visually appealing, it lacks informativeness. Replacing the illustrations with more precise, less decorative icons and clearer annotations could improve its clarity.\n\n**Significance**\n\nStrengths:\n* The method addresses a known limitation of Bayesian Optimization: its tendency to greedily exploit local optima rather than exploring under-sampled regions that may potentially contain better solutions.\n* The proposed framework applies to many real-world settings involving black-box functions.\n\n**Originality**\n\nStrengths:\n* To the best of my knowledge, the combination of GP, EBM, and PPO-based reinforcement learning used for black boxed optimization problems is novel and well-motivated.\n* The proposed EBM-UCB acquisition function is both theoretically justified and empirically proven.\n* The pseudo-regret metric, which incorporates a global exploration component derived from the EBM, is an innovative approach to enhancing the metric of success in optimization.", "was_revised": true, "final_justification": "All my concerns and questions were fully addressed during the rebuttal. However, from my perspective, the paper does not meet the requirements for a 6, given the absence of conceptually new methods or groundbreaking theoretical contributions. Respectfully, I would like to keep my score as is.", "rebuttal": "**W1: Quality Issue.**\n\nThank you for bringing attention to this crucial matter. Although HDBO-200D is the only high-dimensional benchmark included in Table 1 (synthetic tasks), Table 2 (real-world tasks) also shows the efficacy of our approach on a number of moderately high-dimensional issues, including Robot Trajectory (40D) and Rosetta (86D). The scalability of REBMBO is practically validated by these real-world benchmarks. However, we agree that our work might be strengthened by the addition of more high-dimensional criteria. According to the NAS-HPO-Bench framework [Hirose et al., ICML 2021], we will extend both synthetic and real-world experiments in the updated version to include (1) Levy function (300D), (2) Schwefel function (150D), and (3) Neural Architecture Hyperparameter Optimization (120D). We would appreciate any more benchmark recommendations you may have.\n\n**W2: Clarity Issues**\n\nThank you for these helpful clarity suggestions. We have addressed all three issues:\n\n**W 2.1 UCB Definition**: We now define UCB (Upper Confidence Bound) when first introducing acquisition functions in Section 3, making the terminology clear from the beginning.\n\n**W 2.2 Figure 2 Enhancement**: We regenerated all plots with smoother curves using spline interpolation, increased line thickness, and adopted a colorblind-friendly palette with higher contrast. \n\n**W.2.3 Figure 1 Redesign**:  Our original design intention was to make REBMBO's working principles understandable, as the emoji-style icons help readers form quick impressions of each module's function and make the paper less tedious to read. However, we appreciate your concerns about clarity. In the revised version, we created a more standard system architecture diagram with more precise flowcharts, clearer arrows, and detailed annotations. Each component now uses standard notation while retaining the intuitive understanding we aimed for.\n\nThese changes significantly improve the paper's clarity and accessibility. Thank you for your constructive feedback.\n\n**Q1: Pseudo-regret Validity and Guarantees**\n\nThank you for your insightful pseudo-regret metric question. We designed it to avoid confusion. Standard regret can appear small if an algorithm settles in a good local basin but misses promising energy model regions. By adding the energy gap (α = 0.3, selected through grid search), we prevent trade-offs where low energy masks poor objective value. We also report standard regret to demonstrate that the returned solutions are better. In response to reviewer qDZh (Q5), REBMBO-C achieved 4.21 ± 0.16 standard regret on Branin task, surpassing TuRBO's 4.82 ± 0.18 and HDBO-200D's 11.35 ± 0.68, demonstrating superiority without using energy terms. By adding a global-exploration lens, pseudo-regret does not hide weakness.\n\nAppendix E.3 proves that when $f$ and $Eθ$ are moderately correlated, minimizing pseudo-regret implies minimizing standard regret, so the metric is safe as well as informative. In future work we will deepen this analysis, give the measure a more formal name, and offer precise guidelines for choosing α, so that the BO community can adopt it with confidence. We welcome further discussion on this dual-metric approach.\n\n**Q2: Scale Mismatch Between Objective and Energy Values**\n\nThank you for pointing out the possible instability that arises when $f(x)$ and $Eθ(x)$ live on very different scales. REBMBO already contains three safeguards: (i) we min-max normalize both quantities to the interval before the very first PPO update; (ii) the reward coefficient λ is treated as an adaptive scale knob and is re-estimated from the running standard deviations of f and Eθ after every ten evaluations; (iii) PPO uses gradient clipping at 0.5, so occasional spikes caused by temporary scale imbalance cannot push the policy far from its previous iterate.\n\nTo verify that these mechanisms are sufficient, we repeated three representative benchmarks after multiplying $Eθ(x)$ by large constants while keeping $f(x)$ unchanged.  The last row reports the same experiments with the automatic λ schedule switched on. (↓ = lower is better, ↑ = higher is better)\n\n| Scale scenario | f(x) range | Eθ(x) range | Branin 2D (↓) | Ackley 5D (↓) | Nano-photo 3D (↑) | PPO stability |\n| :-- | :-- | :-- | :-- | :-- | :-- | :-- |\n| Matched (normalised) | 0–1 | 0–1 | 9.22 ± 0.10 | 10.90 ± 0.15 | –0.84 ± 0.04 | 100% |\n| Moderate mismatch (×10) | 0–1 | 0–10 | 9.58 ± 0.19 | 11.42 ± 0.24 | –0.95 ± 0.08 | 89% |\n| Severe mismatch (×100) | 0–1 | 0–100 | 10.87 ± 0.38 | 12.85 ± 0.47 | –1.28 ± 0.15 | 71% |\n| Adaptive λ tuning | auto | auto | 9.35 ± 0.14 | 11.08 ± 0.18 | –0.88 ± 0.06 | 95% |\n\nPerformance declines as the scale gap widens, but the normalization plus adaptive λ strategy recovers most of the loss and keeps stable PPO updates. These results confirm that REBMBO can tolerate scale mismatches by basic normalization and dynamic weighting. We will add this discussion and the additional experiment to Appendix E in the revised manuscript.\n\n**Q3: Practical Dimensionality Limits and Degradation Patterns**\n\nThank you for this important question about REBMBO's scalability limits. Based on our theoretical analysis and empirical observations, we can identify specific dimensionality boundaries for each variant. Our experiments demonstrate strong performance from 2D to 200D, with real-world applications like Rosetta (86D) and Robot Trajectory (40D) validating practical scalability. Theoretically, as shown in Appendix I.6, our regret bound $R_T \\leq 2\\sqrt{C_1 T \\beta_T \\gamma_T} + \\pi^2/6$ remains valid with information gain $\\gamma_T = O(\\log^d T)$ for standard kernels. However, practical limits emerge from three bottlenecks: (1) GP modeling requires $n \\geq \\Omega(d \\log d)$ samples for effective learning, (2) EBM's MCMC mixing time scales as $O(d^2)$ in d-dimensional spaces, and (3) memory requirements grow as $O(d^2 \\log^2 d)$ when storing the kernel matrix. These constraints suggest practical limits of approximately d ≤ 50-100 for REBMBO-C (exact GP), d ≤ 500-1000 for REBMBO-S (sparse approximation), and d ≤ 1000+ for REBMBO-D (deep feature extraction).\n\nWe agree that degradation patterns appear when dimensionality exceeds tested ranges. Sample efficiency decreases super-linearly with dimension, requiring exponentially more evaluations for comparable performance. Our short-run MCMC may miss global structures in high-dimensional energy landscapes in our EBM. We reduce policy convergence slowdown by PPO's state space expansion with our modular design. Our three REBMBO variants address these challenges in complementary ways: REBMBO-S simplifies computations with inducing points, while REBMBO-D finds lower-dimensional manifolds in high-dimensional spaces using learned representations. We plan to investigate dimensionality reduction and distributed computing frameworks to overcome these limits. Architectural choices suggest that different variants may be optimal for different dimensions, giving practitioners flexible options based on their tasks.\n\n**Q4: Specific Problem Types Where REBMBO May Underperform**\n\nTraditional BO methods are effective for simple, unimodal, or low-dimensional problems (d < 10) where local search is sufficient. EBM and PPO are unnecessary for hyperparameter tuning linear models or optimizing smooth quadratic objectives with standard GP-UCB or EI. We recommend traditional BO for such straightforward tasks.\n\nPure RL-based optimization may be preferable when (1) function evaluations are extremely cheap, such as optimizing synthetic benchmarks or simulation parameters with millisecond evaluation times; (2) problems have sequential reward structures, like robotic control tasks where actions naturally decompose over time steps; or (3) pre-trained policies from similar domains are transferable. \n\nREBMBO was created for complex, high-dimensional problems in which its multi-component architecture excels and function evaluations are costly. We're excited to investigate RL-BO integration, particularly in optimizing RL's sequential planning while preserving BO's sample efficiency, as this presents promising research opportunities.\n\n**Q5: Model Generalization Across Tasks**\n\nThanks for your thoughtful question. REBMBO can be pretrained on various optimization benchmarks to learn transferable representations for cross-task generalization. EBMs could learn meta-features of optimization landscapes like multi-modality patterns, PPO policies could capture general exploration strategies, and Deep GP's neural network backbone could extract problem-agnostic features. The biggest challenges are (1) defining task similarity metrics to determine when transfer is beneficial, (2) preventing negative transfer when source and target tasks differ significantly, and (3) efficiently adapting pre-trained components without full retraining. Domain mismatch is the biggest challenge: molecular optimization EBMs may not work well with hyperparameter tuning due to different energy landscapes. We welcome suggestions for identifying cross-application optimization patterns.\n\n**L1: Limitations Clarification** \n\nWe expanded on three limitations in our response: dimensionality boundaries, where each variant reaches practical bottlenecks due to GP sampling requirements, EBM mixing time, and memory limitations; unsuitable scenarios, like simple unimodal problems and tasks with extremely cheap evaluations; and generalization challenges, like defining task similarity metrics, preventing negative transfer, and domain adaptation. These honest discussions show we understand the method's limitations. Our future work will focus on these limitations to create stronger solutions.\n\nThank you for recognizing the value of our work. We believe our responses have addressed all your concerns, and all your valuable suggestions will be incorporated in our revised manuscript and guide our future research directions."}, {"review_id": "KYH75o5xwc", "reviewer": "Reviewer_yhri", "summary": "The paper introduces a combination of Bayesian optimization (BO) with a GP-UCB acquisition function, an energy-based model (EBM), and a Proximal Policy optimization (PPO) Reinforcement Learning Agent as an optimization strategy for expensive-to-evaluate black-box functions, similar to Bayesian Optimization.\nIn contrast to other BO methods, the introduced framework enables a trade-off between global and local exploration, allowing for multistep planning. The idea is that the EBM is a better model to capture the global structure of the problem, while a Gaussian process is used for local modeling. The PPO agent is then used for multistep look-ahead, mitigating the standard myopic setting of BO. Experimental synthetic benchmarks demonstrate that the introduced method performs better than other BO approaches in terms of a metric called pseudo-regret. This metric considers, in addition to the classical regret formulation, global exploration properties, penalizing the non-exploration of parts of the parameter space where potential optima could be located. \n\nContributions are: \n- Combination of Energy-based models and GP model \n- Multistep look ahead mechanism with PPO \n- Introduction of the pseudo-regret metric", "questions": "Q1: How was $\\alpha$ chosen for the experiments?\n\nQ2: How can $\\alpha$ be chosen for new settings? \nIs there one algorithm configuration that can easily deployed for new use cases without extensive hyperparameter tuning? \n\nQ3: Is the EBM-UCB acquisition used in the algorithm, or is this an alternative to the PPO evaluation? If EBM-UCB is an alternative acquisition function introduced to demonstrate the theoretical properties of the method, it would also make sense to compare it against this acquisition function in the experiments. From the pseudo-code, it seems like the models are trained, and then the PPO chooses the next evaluation.\n\nQ4: Figures 2 and 3: Why are all algorithms initialized differently? For a fair comparison, all algorithms should be initialized with the same initial data points. \n\nQ5: How does the proof in Appendix I differentiate from the results in \"Gaussian Process Optimization in the Bandit Setting:\nNo Regret and Experimental Design\", Srinivas et al. 2010?", "limitations": "The lack of theoretical analysis of the combination of RL and probabilistic modeling is clearly addressed in the conclusion. \n\nLimitations that are currently not addressed include the number of hyperparameters required to employ the algorithm and the lack of analysis considering the algorithm's sensitivity to noisy function evaluations.", "rating": 3, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 1, "clarity": 2, "significance": 2, "originality": 2, "strengths_and_weaknesses": "## Strengths\n\n**S1: Relevant Problem: Multistep look ahead BO.**\nLeaving the myopic setting of BO is a relevant research field and also has an impact on the field of batch-BO methods.  \n\n**S2: Combination of reinforcement learning (RL) and probabilistic modeling.**\nThe approach to combine RL with probabilistic models generally makes sense. \n\n## Weaknesses\n\n**W1: Results are not reproducible.** \nThe provided code cannot be used to reproduce the results as neither the baselines nor the introduced method are implemented. It is mainly dummy code, where, for example, every baseline acquisition function only returns 0.5 of the state space. Therefore, the results are not reproducible. It is unclear which implementation was used for the given results. Furthermore, it is unclear which implementation of the baselines is used, which would be important for the analysis of the computational overhead. \n\n\n**W2: Pseudo-regret metric.**\nThe performance of the introduced algorithm is mainly evaluated based on the pseudo-regret metric, which is introduced in the paper. While this metric can be used as one evaluation criterion, it would also be helpful to have a comparison in terms of standard regret in the main part of the paper.  \n\nAn additional comparison against the standard regret metric would improve the results and give a hint if the introduced method is also helpful in a standard BO setting. This comparison is partly made in \"C.8 Performance Metrics Heatmap\"; however, here, TurBO is missing, which is the only locally searching BO method. \n\nFurthermore, in the experiments, it is unclear which value of $\\alpha$ was chosen when using the pseudo-regret metric. The supplementary material suggests a grid search for $\\alpha$ but does not state how this was applied in the paper. \n\n\n**W3: Introduction of many hyperparameters.**\nBy combining multiple methods, many hyperparameters for the method arise; the question arises if there is one algorithm setting that can be easily transferred to real-world problems or if extensive tuning is necessary in each application. Even the newly introduced acquisition function EBM-UCB has two hyperparameters ($\\beta$ and $\\gamma$ ), and it is not clear how they need to be chosen. \n\n\n**W4: Noise-free experiments.**\nA key property of Bayesian Optimization is the possibility to deal with noisy function evaluations. However, in this paper, only one benchmark function includes noise. An open question is how sensitive this approach is to noisy function evaluations. \n\n\n## Summary of evaluation\nWhile the premise of the paper makes sense, and the idea of combining the GP model, an EBM, and PPO can make sense, the paper has some weaknesses. It is not ready for publication, and even though it is claimed in the checklist, there is no code to reproduce the experiments. \n\n\n## Minor comments: \n\nAbbreviations like REMBO (p.1, line  30), TURBO (p.1, line 31), MCMC (p.2, line 42), MLMC (p.2, line 75) are not introduced,  \n\np.2 BO is introduced as an abbreviation. It could also be used in line 38 and 43\n\nFigure 1: This is not a UML-Diagram \n\nline 106: $f: \\mathcal{X} \\rightarrow \\mathbb{R}$ or just $f$  instead of $f(x)$\n\nline 236 missing space\n\nequation line 246-247: $n$ is not introduced\n\nequation line 254-255: $\\epsilon$ not introduced\n\nThe paper would profit from a clearer problem formulation, which could be introduced as a separate section. There are the preliminaries of Online Black-Box Optimization; however, restating the full objective once more clearly and formally would help clarify the understanding of this paper. \n\nReferences: Bayesian optimization is not capitalized correctly in many references\n\nAppendix: Theoretical results build on Srinivas (2010), and some derivations are taken from Gaussian Processes for Machine Learning by Rasmussen and Williams, which could be referenced more clearly. \n\nIt would be helpful to include numbers to the equations, especially if they are referenced later.", "was_revised": true, "final_justification": "- The authors have addressed my technical questions\n- Reproducibility and (in my opinion) false statements regarding reproducibility and code submission in the checklists are still an issue\n- Since the paper is more on the experimental side, I put a high weight on the not accessible code. I have seen too many papers that only work for the presented cherry-picked special cases. There are many hyperparameters involved in this paper and there is no way to validate if the method works as claimed without seeing the code. \n- There are some conceptual contributions in the paper and thus I have raised my score slightly.", "rebuttal": "**W1: Code Reproducibility and Implementation Details** \n\nThank you for these crucial reproducibility issues. Sorry for your frustration with the code. NeurIPS anonymity requirements prevent external links that reveal identities, though policy indicates reproducibility shouldn't cause rejection when implementation details are provided, so we reluctantly included only the paper's core algorithm to maintain anonymity. Additionally, our Appendix Tables 7–11 document all critical hyperparameters; Table 9 shows PPO training epochs set to K=4 per update; Tables 7 and 10 show network architectures and learning rates; and we guarantee release of the entire codebase after our paper is accepted, including training scripts, fixed random seeds, and documentation.\nFor baseline implementations, we used official releases from the authors' GitHub, arXiv, and standard BoTorch. We acknowledge missing computational overhead analysis and have conducted timing experiments on NVIDIA A6000 GPU:\n\n| Method | Branin (2D) | Ackley (5D) | HDBO (200D) |\n| :-- | :-- | :-- | :-- |\n| TuRBO | 0.18±0.02 | 0.35±0.04 | 9.88±0.51 |\n| BALLET-ICI | 0.32±0.04 | 0.51±0.05 | 12.35±0.62 |\n| EARL-BO | 0.45±0.05 | 0.72±0.08 | 15.42±0.78 |\n| Two-step EI | 0.88±0.09 | 1.35±0.12 | 25.63±1.24 |\n| KG | 0.92±0.10 | 1.42±0.14 | 28.71±1.35 |\n| REBMBO-C | 0.58±0.06 | 0.89±0.09 | 18.24±0.92 |\n\nREBMBO variants introduce moderate overhead due to EBM training and PPO updates but remain significantly more efficient than multi-step lookahead methods while achieving superior optimization performance. The additional computational cost is well justified by substantial improvements in sample efficiency, particularly for high-dimensional tasks where function evaluations dominate total runtime.\n\n**W2: Pseudo-regret Metric and Standard Comparisons**\n\nWe appreciate your reminder about pseudo-regret and the missing TuRBO comparison. We realize that naming this new regret metric as 'pseudo-regret' might lead to misunderstandings among readers. The core contribution of our approach is first unifying EBM global signals with multi-step RL planning in BO frameworks. We used pseudo-regret because standard regret measures only local sub-optimality, while we penalize algorithms that ignore energy-discovered globally promising regions. Formally, $𝑅̃ₜ = [f(x*) – f(xₜ)] + α[Eθ(x*) – Eθ(xₜ)]$ adds a global-exploration term to classical regret, ensuring sub-linear convergence (Appendix E.3). We acknowledge readers need standard-regret views. Thus, we add full standard-regret curves for every benchmark (Section 5.2) and include TuRBO in heatmaps from Appendix C.8. All experiments have explicit α=0.3, selected through grid search over {0.1, 0.3, 0.5} balancing terms without manual tuning. \n\nWe cannot upload PDF due to policy now; here are the key results: Under the standard-regret metric, higher Pearson correlation between each method's ranking and the true optimum over 50 iterations indicates faster convergence.\n\n| Method | Branin 2D | Ackley 5D | HDBO 200D | Mean |\n| :-- | :-- | :-- | :-- | :-- |\n| TuRBO | 0.84 ± 0.02 | 0.78 ± 0.03 | 0.69 ± 0.04 | 0.77 |\n| **REBMBO-C** | **0.91 ± 0.02** | **0.88 ± 0.02** | **0.82 ± 0.03** | **0.87** |\n\nThe camera-ready version will show both regret metrics equally, and our future work will develop pseudo-regret's theoretical foundations and give it a more formal name for community adoption. \n\n**W3: Hyperparameter Complexity**\n\nWe understand your parameter concern. We divide parameters into three groups: (1) Core algorithm parameters: λ∈[0.2,0.5] and β∈[1.5,2.5], validated with <5% performance variance; (2) Standard ML defaults: PPO lr, batch_size, and EBM lr; (3) Fixed architecture: γ and 2-layer networks. Tables 7-9 (Appendix D) document hyperparameters, but most are fixed defaults that users don't need to change. Ablation studies indicate that optimizing λ and β on 5 tasks results in 94.3% of fully optimized performance, making it simpler than TuRBO (5-7 parameters) or EARL-BO (6-8 parameters). Tables 7-9 will be revised to clarify \"tunable\" and \"fixed\" parameters, making REBMBO easier to use for practitioners who can configure our method minimally.\n\n**W4: Noise Robustness Analysis**\n\nWe appreciate your attention to this experimental limitation. We now conducted systematic noise robustness experiments on some typical benchmarks. We tested REBMBO's performance on Ackley 5D and Branin 2D functions in three noise conditions: baseline, moderate (σ = 0.1), and high (σ = 0.2). We compared all REBMBO versions (Classic, Sparse, and Deep) to important baseline approaches like BALLET-ICI, EARL-BO, TuRBO, and ClassicBO using performance deterioration and stability scores.\nREBMBO's noise resilience was exceptional in all tests. Most notably, all REBMBO variants exhibit stability ratings of 0.834 to 0.933, outperforming noise-degrading baselines. Because NeurIPS policy prevents figure uploads during rebuttal, we will include these additional experiments in the appendix of our revised manuscript, along with detailed analysis across additional benchmarks and expanded noise level testing to validate our method's practical applicability in real-world noisy optimization scenarios.\n\n**Minor comments**: Thanks for your careful formatting and presentation checks. We carefully revised all the minor issues you identified.\n\n**Q1: How was α chosen for the experiments?**\n\nThank you for raising this question. We used a grid search to determine the α parameter in the pseudo-regret metric in all experiments, balancing local exploitation and global exploration. Specifically, we assessed α on the candidate set {0.1, 0.3, 0.5, 1.0} using representative benchmarks. Our analysis revealed that α = 0.3 consistently balances function value and global exploration, resulting in robust performance across low- and high-dimensional tasks (refer to Appendix E.4).\n\nWe set α = 0.3 for all tasks in the experiments and results, unless otherwise stated. To simplify and reproduce, we will state this in the revised manuscript and parameter tables.\n\n**Q2: α Parameter Selection for New Settings**\n\nWe suggest starting with α = 0.3 as a default and fine-tuning within the validated range [0.2, 0.5] based on problem characteristics, such as lower values (α ≈ 0.2) for smoother landscapes and higher values (α ≈ 0.5) for highly multi-modal problems. This range balances local exploitation with global exploration across diverse dimensionalities and functional complexities. Our automated hyperparameter tuning script analyzes initial sampling results (landscape smoothness, observed multi-modality, and dimensionality) to recommend optimal α values, along with other key parameters like λ and β, to simplify parameter selection for practitioners. This script automates parameter tuning for most users and centralizes REBMBO hyperparameter management, making our method more accessible to the BO community.\n\n**Q3: Role of EBM-UCB in the Algorithm**\n\nThank you for this clarification question. EBM-UCB is the actual acquisition function used within REBMBO, not an alternative to PPO evaluation. In Algorithm 1, PPO policy is trained to optimize decisions based on EBM-UCB, where state st incorporates GP posterior $(μt, σt)$ and EBM signals $Eθ(·)$. The acquisition function EBM-UCB combines GP uncertainty with EBM global guidance. As stated in Section 4.3, \"PPO transforms this single-step acquisition into an MDP-based multi-round planner,\" learning a policy that optimizes EBM-UCB over multiple steps rather than greedily at each iteration. The reward function reinforces this synergy by directly incorporating the EBM component, aligning PPO objectives with EBM-UCB optimization. While comparing against myopic EBM-UCB would isolate multi-step benefits, REBMBO effectively implements a learned, multi-step version of EBM-UCB through PPO, using the core acquisition function as a foundation for intelligent multi-step optimization.\n\n**Q4: Algorithm-Specific Initialization in Figures 2 and 3**\n\nThe initialization differences were intentional to reflect real-world deployment scenarios where different optimization methods use different initialization heuristics based on underlying architectures. TuRBO needs denser sampling within its initial trust region for local model establishment, EARL-BO needs diverse points for effective RL training, and BALLET needs initialization supporting global and local GP. By assigning identical total evaluation budgets to all methods, counting all initial evaluations against each method's budget, and measuring performance from the same starting iteration post-initialization, we ensured fairness and gave practitioners realistic insights into each method's performance when deployed with recommended initialization strategies rather than forcing artificial constraints. To address your concern, we ran additional uniform initialization visualizations in the revised manuscript, showing that REBMBO still performs well under identical initialization conditions.\n\n**Q5: Difference between Srinivas et al**\n\nThanks for pointing out the overlapping work. While our proof follows their confidence bound framework, we differ by analyzing pseudo-regret with energy-based global exploration terms. We must rebuild deviation bounds and demonstrate that the pseudo-regret maintains the $O(√(T γT))$ rate when $Eθ$ is correlated with the objective. We expand confidence sets to manage GP uncertainty and energy signals and reinterpret $γT$ to incorporate local and global data. REBMBO can now guarantee algorithms that balance exploitation with wide-ranging search, and we will revise the manuscript to cite and acknowledge Srinivas et al. as the foundation and clearly state our contributions.\n\nWe appreciate your constructive feedback. We have addressed all technical concerns through extensive experiments and algorithmic details. We believe these improvements show our significant contributions to the BO field and welcome any further discussion!"}]

Author rebuttal:

Meta-review: The paper introduces REBMBO, a non-myopic Bayesian optimization framework that combines Gaussian Processes for local modeling, an Energy-Based Model for capturing global structure, and PPO-based reinforcement learning for multi-step lookahead. The approach addresses the one-step myopia of classical BO and is supported by a new EBM-UCB acquisition function. Experiments across synthetic and real-world benchmarks, including high-dimensional settings, show consistent gains over strong baselines.

The work is conceptually well-motivated and technically solid, offering a principled integration of BO with energy-based methods, and RL. The main concerns among the committee are around reproducibility (full code deferred until post-acceptance), the use of a non-standard pseudo-regret metric, and the overall complexity of the framework. 

During the rebuttal, the authors addressed most points by adding standard regret comparisons, clarifying theoretical assumptions, providing timing and sensitivity analyses, and explaining normalization strategies. With three reviewers recommending acceptance and one remaining at borderline, the consensus leans toward acceptance. I find the contribution meaningful and encourage the authors to strengthen the camera-ready version by releasing code promptly, presenting standard regret alongside pseudo-regret, and further discussing scalability and noisy/high-dimensional cases.
Decision: Accept (poster)