Transition Constrained Bayesian Optimization via Markov Decision Processes
Authors: Jose Pablo Folch
Calvin Tsay
Robert Matthew Lee
Behrang Shafei
Weronika Ormaniec
Andreas Krause
Mark van der Wilk
Ruth Misener
Mojmir Mutny
Venue: NeurIPS 2024 poster (2024)
Keywords: Bayesian Optimization
Transition Constrained
Markov Decision Process
Linear Bandits
Convex Reinforcement Learning
Primary area: optimization
TL;DR: We do Bayesian Optimization under transition constraints by creating and solving tractable long-term planning problems in Markov Decision Processes.

Abstract: Bayesian optimization is a methodology to optimize black-box functions. Traditionally, it focuses on the setting where you can arbitrarily query the search space. However, many real-life problems do not offer this flexibility; in particular, the search space of the next query may depend on previous ones. Example challenges arise in the physical sciences in the form of local movement constraints, required monotonicity in certain variables, and transitions influencing the accuracy of measurements. Altogether, such *transition constraints* necessitate a form of planning. This work extends classical Bayesian optimization via the framework of Markov Decision Processes. We iteratively solve a tractable linearization of our utility function using reinforcement learning to obtain a policy that plans ahead for the entire horizon. This is a parallel to the optimization of an *acquisition function in policy space*. The resulting policy is potentially history-dependent and non-Markovian. We showcase applications in chemical reactor optimization, informative path planning, machine calibration, and other synthetic examples.

Transition Constrained Bayesian Optimization via Markov Decision Processes
Jose Pablo Folch Imperial College London London, UK
Calvin Tsay Imperial College London London, UK
Robert M Lee BASF SE Ludwigshen, Germany
Behrang Shafei
Weronika Ormaniec
Andreas Krause
BASF SE Ludwigshen, Germany
ETH Zurich Zurich, Switzerland
ETH Zurich Zurich, Switzerland
Mark van der Wilk Imperial College London London, UK
Ruth Misener
Mojmír Mutný
Imperial College London London, UK
ETH Zurich Zurich, Switzerland
Abstract
Bayesian optimization is methodology optimize black-box functions. Tra- ditionally, focuses on the setting where you can arbitrarily query the search space. However, many real-life problems do offer this flexibility; in particular, the search space of the next query may depend on previous ones. Example chal- lenges arise in the physical sciences in the form of local movement constraints, required monotonicity in certain variables, and transitions influencing the accuracy of measurements. Altogether, such transition constraints necessitate a form of planning. This work extends classical Bayesian optimization via the framework of Markov Decision Processes. We iteratively solve a tractable linearization of our utility function using reinforcement learning to obtain policy that plans ahead for the entire horizon. This parallel to the optimization of an acquisition function in policy space. The resulting policy is potentially history-dependent and non-Markovian. We showcase applications chemical reactor optimization, informative path planning, machine calibration, and other synthetic examples.
Introduction
Many areas in the natural sciences and engineering deal with optimizing expensive black-box functions. Bayesian optimization (BayesOpt) B. method to optimize these problems using probabilistic surrogate, has been successfully applied to myriad of examples, e.g. prper-parameter selection
4 robotics [5], battery design 66. laboratory equipment tuning [7], and drug discovery []. However, state-of-the-an algorithms are often ill-suited when physical sciences interact with potentially dynamic systems ( In such circumstances. real-life constraints limit our future decisions while depending on the prior state of our interaction with the system. This work focuses on transition constraints influencing future choices depending on the current state of the experiment. In other words, reaching certain parts of the decision space (search space) requires long-term planning in our optimization campaign. This effectively means we address a general sequential-decisoo problem akin to those studied in reinforcement learning (RL) or optimal control for the task of optimization. We assume the transition constraints are known a priori to the optimizer.
Applications with transition constraints include chemical reaction optimization [10H12), environ- mental monitoring [3]17] lake surveillance with drones [8), energy systems 21). vapor
compression systems 22 electron-laser tuning 24 and seabed identification (24] For example,
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

Figure depicts an application in environmental monitoring where autonomous sensing vehicles must avoid obstacles (similar to Hitz al. [8]. Our main focus application are transient flow reactors 27 Such reactors allow efficient data collection by obtaining semi-continuous time-series data rather than single measurement after reaching the steady state of the reactor. As we can only change the inputs of the reactor continuously and slowly to maintain guasi-steady-statee operation, allowing arbitrary changes, as in conventional BayesOpt, would result in measurement sequences which are not possible due to physical limitations.
Problem Statement, More formally, we design an algorithm identify the optimal configuration of a physical system governed by a black box function f, namely, x²' = arg max fz). The set X' summarizes all possible system configurations, the so called search space. We assume that we can sequentially evaluate the unknown function at specific points x in the search space and obtain noisy observations, y = f()) + <(x), where has a known Gaussian likelihood, which possibly heteroscedastic. We assume that If can be modeled probabilistically using Gaussian process prior that we introduce later. Importantly, the order of the evaluations is dictated by known, potentially stochastic, dynamics modeled by Markov chain that limits our choices of E X.
BayesOpt with a Markov Decision Processes. The problem of maximizing an unknown function could be addressed by BayesOpt, which typically chooses to query f(x) by sequentially maximizing an acquisition function, u.
x++1 = arg max u(x[X,), zEX'
(I)
depending on all the past data iteration 4, X, Eq. TT arises as a greedy one-step approximation whose overall goal is to minimize e.g. cumulative regret, and assumes that any choice of point in the search space X' is available. However, given transition constraints, we must traverse the search space according to the system dynamics. This work extends the BayesOpt framework and provides a method that constructs a potentially non-Markovian policy by myopically optimizing utility as,
1 arg max (1==((/()) (T (=)..),
where U is the greedy utility of the policy т and X₁ encodes past trajectories through the search space. In the following sections, we will show how to tractably formulate the overall utility, how to greedily maximize it, and how to adapt
admit policies depending on the full optimization history.
Contributions. We present a BayesOpt framework that tractably plans over the complete experimen- tation horizon and respects Markov transition constraints, building on active exploration in Markov chains [7]. Our key contributions include:
We identify novel utility function for maximum identification greedily optimize it. The optimization tractable, and does not scale exponentially in the policy horizon. In many cases, the problem is convex in the natural representation.
function of policies, and
We provide exact solutions the optimization problems using convex optimization for discrete Markov chains. For continuous Markov chains, we propose a reparameterization by viewing our problem as an instance of model predictive control (MPC) with a non- convex objective. Interestingly. in both cases, the resulting policies are history-dependent Markovarrni..).
We analyze the scheme theoretically and empirically demonstrate its practicality on problems in physical systems, such electron laser calibration and chemical reactor optimization.
2 Background
Because our contributions address experimental design of real-life systems by intersecting design of experiments, BayesOpt and RL, we review each of these of components. Refer to Figure 5 in Appendix A for visual overview of how we selected the individual components for tractabilityy of the entire problem.
Gaussian Processes To model the unknown function f. we use Gaussian processes (GPs) (28). GPs are probabilistic models that capture nonlinear relationships and offer well-calibrated uncertainty estimates. Any finite marginal of GP, e.g., for inputs (#1..., xp)) the values (fxx)) J are normally distributed. We adopt Bayesian approach and assume f is sample from GP prior with a known covariance kernel, k and zero mean function, f ~ GP k). Under these assumptions, the posterior of f. given a Gaussian likelihood of data, is a GP that is analytically tractable.

---

Figure 1: Representative task of finding pollution in river while following the current. (a) Problem formulation: The star represents the maximizer and the arrows the Markov dynamics. (b) Objective formulation: Orange balls represent potential maximizers, with size corresponding to model uncer- lainty. (c) Optimization: Deploy potentially stochastic policy that minimizes our objective.
2.1 Maximum Identification: Experiment Design Goal
Classical BayesOpt is naturally myopic in its definition greedy one-step update (see (T). but has the overall goal to minimize, e-g., the cumulative regret. Therefore u needs to chosen such that overall non-myopic goals can be achieved, usually defined as balancing an exploration-explotanee trade-off. In this paper we follow similar ideas; however, we do not focus on regret but instead on gathering information to maximize chances identify x*, the maximizer of f.
Maximum Identification via Hypothesis testing. Maximum identification can naturally be expressed as multiple hypothesis testing problem, where we need to determine which of the elements in N' is the maximizer. To do so, we require good estimates of the differences (or least their signs) between individual queries (x;;) - fzzz); Jy E X. For example, if fz;) f(z)) < 0, then z; cannot be maximizer. Given the current evidence, the of arms which we cannot rule out are potential maximizers, Z X. At termination we report our best guess for the maximizer as:
IT = arg max maxxmaxaaaa (T(x),
where pT is the predictive mean at termination time T.
ZEZ
Suppose in step out of T, then let X, be the of previous queries, seek to identify new Xnew that when evaluated minimize the probability of returning sub-optimal arm at the end. For given function draw f. the probability of returning wrong maximizer 2 # ry is f) We can then consider the worst-case probability across potential maximizers, and taking expectation over we obtain utility through an asymptotic upper- upper-bound bound on the log-pro bog-probability, indeed for large T' we obtain:
P(ur(z) pr(x))
(f(z) (12)-f(zt)) f(z))) :E2((05) Xx,,Xx.ev (z, x;,
E, sup log (( - uT(z)) -r =EZ\(#}}
Olf)
$ -
sup
The expectation on the current prior (posterior up to X,), the kernel is the posterior kernel given observations X, U Since we consider the probability of an error, it is more appropriate to talk about minimizing instead of 'maximizing the utility" but the treatment is analogous. Further, note the intuitive interpretation of the bound: the probability of an error will be minimized if the uncertainty is small or the values (=) and [(x;)) are apart. The non-trivial distribution of f(**) 29 renders the utility intractable; therefore we employ a simple and tractable upper bound on the objective 33) which can be optimized by minimizing the uncertainty among all pairs in Z:
U(Xnew)
max 2',
Var[f(z)
{x Xaew]-
Such objectives can be solved greedily in similar way acquisition functions in Eq. ! by minimizing over X Note that Fiez et al. [30] derive this objective for the same problem with linear bandits, albeit they consider the frequentist setting and (surprisingly) different optimality criterion: minimizing T for a given failure rate. For their setting, the authors prove that it is an asymptotically optimal objective to follow. They do not consider any Markov chain structure. Derivation of the Bayesian utility and its upper bound in (4) can found in Appendix C.I
Utility with kernel embeddings. For illustrative purposes, consider special case where the kernel has a low rank due existence of embeddings ((x) E R i.e., k(x, y) +(x) + Such
embeddings can be, e.g., Nyström features 31 or Fourier features 32 33
While not necessary,

---

these formulations make the objectives considered in this work more tractable and easier to expose to the reader. With the finite rank assumption, the random function f becomes,
f(x) = #(x)² 0 and 0 (0, Imxm)
(5)
where 9 are weights with Gaussian prior. We can then rewrite the objective Eq. (4 as:
nemm)) = max 2,
(+(2 - (
(6)
X
This reveals an essential observation that the utility depends only on the visited states; not their order. This suggests a vast simplification, where we do not to model whole trajectories, and Markov decision processes sufficiently describe our problem. Additionally, numerically, the objective involves the inversion of an m x m matrix instead of |X] x | (see Sec. 4 Appendix provvides a utility without the finite rank-assumptions that is more involved symbolically and computationally.
2.2 Markov Decision Processes
To model the transition constraints, we use the versatile model of Markov Decision processes (MDPs). We assume an environment with state space X' and action space A, where we interact with an unknown function f X A - R by rolling out policy for H time-steps (horizon) and obtain a trajectory, 7 (Zg. 0g. 21, 01 ..... IH-1, -1) From the trajectory, we obtain a sequence of noisy observations y(T) {y(2o, an),
-2 s.t. y(za) = f(z, as) + c(z), an). where e(T,, an) is zero-mean Gaussian with known variance which is potentially state and action dependent. The trajectory is generated using a known transition operator P( |2n, an). A Markov policy x (a, (ar) is a mapping that dictates the probability of action an in state The Hence, the state-to-state transitions are P(æR+1, =T.) EaeA Th (al (a)(ra)) P( ( a). In fact, an equivalent description of any Markov policy и is the corresponding distribution giving us the probability of visiting state-action pair under the policy, which we denote d₂ E D, where
D = vh E [H] da (x, a) 2 O, Σ (z,za))=1 1, Edh(s',a) 0,K
x,a
We will use this polytope to reformulate our optimization problem over trajectories. Any d E D can be realized by a Markov policy π and vice-versa. We work with non-stationary policies, meaning the policies depend on horizon count h. The execution of deterministic trajectories is only possible for deterministic transitions. Otherwise, the resulting trajectories are random. In our setup, we repeat interactions 7 times (episodes) obtain the final dataset of the form X,- (T))
2.3 Experiment Design in Markov Chains
Notice that the utility U in Eq. deeeends on the states visited and hence states of the trajectory. In our notation, X, will now form a set of executed trajectories. With deterministic dynamics, we could optimize over trajectories, but this would lead to an exponential blowup (i.e. (X|#). In fact, for stochastic transitions, we cannot pick the trajectories directly, so instead we work in the space of distributions. For given policy, through sampling, of all the state-action pairs visited during policy executions, dz (2, a), which assigns equal mass to each state-action visited during our trajectories. This allows us to focus on the expected utility over the randomness of the policy and the environment, namely,
are able to create an empirical distribution
u(d+)) ==
(d
(7)
This formulation stems from Mutný et al. 17 who try to tractably solve such objectives that arise in experiment design by performing planning in MDPs. They focus on learning linear operators of an unknown function, unlike identifying maximum, we do here. The key observation they make is that any policy æ induces a distribution over the state-action visitations, dę. Therefore we can reformulate the problem of finding the optimal policy, into finding the optimal distribution over state-action visitations as: mind, ED Udz)), ). and then construct policy и via marginalization. We refer to this optimization the planning problem. The constraint D encodes the dynamics of the MDP.
2.4 Additional Related Works
The most relevant prior work to ours is exploration in reinforcement learning through the use of Markov decision processes as in Mutný et al. [17 and convex reinforcement learning of Hazan et al. 34 Zabavy et al. [35] which we will use to optimize the objective. Other related works are:

---

Pure exploration bandits objectives. Similar objectives to ours have been explored for BayesOpt. Li and Scarlett [36] use the G-allocation variant of our objective for batch BayesOpt. achieving good theoretical bounds. Zhang et al. [37 and recently Han et al. 38 take advantage of possible maximizer sets to train localized models, while Salgia et al. [39 show that considering adaptive maximization sets yields good regret bounds under random sampling. Contrary to them, motivation and derivation in terms of a Bayesian decision rule do not appear elsewhere according to our best knowledge. We also recognize that we can relax the objective and optimize it in the space of policies.
Optimizing over sequences. Previous work has focused on planning experimental sequences for minimizing switching costs [1I 21 40 41 however they are only able adhere to strict constraints under truncation heuristics (20) 22 42 Recently, Qing et al. [43] also tackle Bayesian optimization within dynamical systems, with the focus of optimizing initial conditions. Concurrent work of Che et al. [4 tackles constrained variant of a similar problem using model predictive control with a different goal.
Regret vs Best-arm identification. Most algorithms in BayesOpt focus on regret minimization. This work focuses on maximizer identification directly, i.e., to identify the maximum after a certain number of iterations with the highest confidence. This branch of BayesOpt mostly addressed in the bandit literature [45]. Our work builds upon prior works of Soare et al. [46], Yu et al. [47], and specifically upon the seminal approach of Fiez al. [0] to design an optimal objective via hypothesis
testing. Novel to our setting is the added difficulty of transition constraints necessitating planning.
Non-myopic Bayesian Optimization. Look-ahead BayesOpt (48-44 seeks improve the greedy aspect of BayesOpt. Such works also use an MDP problem formulation, however, they define the state space to include all past observations (e.g. (55] (6)). This comes at the cost of simulating expensive integrals, and the complexity grows exponentially with the number of look-ahead steps (usually less than three steps). Our work follows a different path, we maintain the greedy approach to control computational efficiency (i.e. by optimizing over the space of Markovian policies), and maintain provable and state-of-art performance. Even though the optimal policy through non-mi monmic analysis is non-Markovian in Sec. 4 we show that adaptive rexampling iteratively approximates this non-myoptic optimal policies a numerically tractable way via receeding horizon planning. In our experiments we confortably plan for over a hundred steps.
3 Transition Constrained BayesOpt
This section introduces BayesOpt with transition constraints. We use MDPs to encode constraints. Namely, the Markov dynamics dictates which inputs we are allowed to query at time-step h + given we previously queried state FA. This mean that the transition operator is P(r 1 a) 0 for any transition Tú + T4+1 not allowed by the physical constraints.
Motivated by our practical experiments with chemical reactors, we distinguish two different types of feedback. With episodic feedback we can be split the optimization into episodes. At the end of each episode of length H, we obtain the whole set of noisy observations. On the other hand, instant feedback is the setting where we obtain a noisy observation immediately after querying the function. Asynchronous feedback describes mix of the previous two, where we obtain observations with unspecific a delay.
3.1 Expected Utility for Maximizer Identification
In section 2.T we introduced the utility for maximum identification. Using the same simplifying assumption (finite rank approximation of GPs in Sec. Sec.2.1 Eq. 4 we can show that the expected utility 24 can be rewritten in terms of the vate-action distribution induced by Xnew:
U(d=) max [1((2) +b> (b==))[((()>>-
(8)
2,2'EZ
where V(d+) - -E de (a. + TH 1). The variable di (2., a) is a state-action (=,0))
visitation, ((x) are e.g. Nyström features of the GP. We prove that the function is additive in terms of state-action pairs in Lemma D.I in Appendix [ a condition required for the expression as function of state-action visitations [7]. Additionally, by rewriting the objective in this form, the dependence and convexity with respect to the state-action density dz becomes clear as it is only composition of a linear function with an inverse operator. Also, notice that the constraint set a convex polytope.
Therefore we are able to use convex optimization to solve the planning problem (see Sec.
5

---

Algorithm Transition Constrained BayesOpt via MDPs
Input: Procedure for estimating sets of maximizers, initial point 20, initial set of maximizer
candidates Z0
Initialize the empirical state-action distribution do for t. = 0 to 7
do
for h 0 to H 1 do
Ue. ) - U(d @ de, Ze., Lt,. Tt.h)
II define the objective, see eq. (
Wr, arg mid,,,,....... U
Il solve MDP planning problem deploy policy
2t,++1 = Te. (21,6)
if feedback is immediate then
|| obtain observation Il update model and maximizer candidate set 1 (2) + de + ((21R+112)) updaate empirical state-action distribution, see eq. (II
Mt.h+1 = (I Et
GPi,, Ze, f- Update(X1,h, Y,,)
if feedback is episodic then
obtain observations // update model and maximizer candidate set Return: Estimate of the maximum using the GP posterior's s mean S. arg max (T(
Y.H f(X,,)) C.
GP++1.: 1++1: f Update(X1,r) Y..
The definition of the objective requires the use of set of maximizers, In the ideal case, we can say a particular input x, is not the optimum there exists x' such that [a") > f(x) with high confidence. We formalize this using the GP credible sets (Bayesian confidence sets) and define:
Set of potential maximizers 2.
2: = {s X (CB(
2 sup LCB(J(x') s'EX
where UCB and LCB correspond to the upper and lower confidence bounds of the GP surrogate with a user specified confidence level defined via the posterior GP with data up to X
3.2 Discrete vs Continuous MDPs.
Until this point, our formulation focused on discrete S and A for of exposition. However, the framework is compatible with continuous state-action spaces. The probabilistic reformulation of the objective Eq. is possible irrespective of whether X' (or .A) discrete or continuous subset of R² In fact, the convexity of the objective in the space of distributions is still maintained. The difference is that the visitations are longer probability mass functions but have to be expressed as probability density functions dr( a). To recover probabilities in the definition of V, we need to ((x.a)((z.a) replace sums with integrals i.e. ExEx,aEA d(x)) ((.)) - SSEX de(x,, 0) d(2,a)"
In the Eq. 8) we need to approximate a maximum over all input pairs in Z. While this can be enumerated in the discrete case without issues, it poses non-trivial constrained optimization problem when X' is continuous. As an alternative, we propose approximating the set using a finite approximation of size K which can be built using Thompson Sampling [57 58 or through maximization of different UCBs for higher exploitation (see Appendix[E.] [..)p. In Appendix E..] we numerically benchmark reasonable choices of K, and show that the performance is not significantly affected by them.
3.3 General algorithm and Theory
The general algorithm combines the ideas introduced so far. We present Algorithm I Notice that apart from constructing the current utility, keeping track of the visited states and updating our GP model, an essential step is plamning, where we need to find a policy that maximizes the utility. As this forms the core challenge of the algorithm, a sequence of dynamic programming problems defined by the steps of the Frank- Wolfe algorithm. From a theoretical point of view, under the assumption of episodic feedback, the algorithm provably minimizes the utility as we show in Proposition C.Iin in Appendix C.4
devote Sec. to it. In short, it solves
4 Solving the planning problem
The planning problem, defined as mina, ED (d=,,, can thought of as analogous to optimizing an acquisition function in traditional BayesOpt, with the added difficulty of doing it in the space of

---

policies. See the bottom half of Figure in Appendix A for a breakdown of the different components of our solution. Following developments Hazan al. [34 and Mutný al. [17], we use the classical Frank-Wolfe algorithm [59]. It proceeds by decomposing the problem into series of linear optimization sub-problems. Each linearization results in a policy, and we build mixture policy consisting of optimal policies for each linearization m {(a,, mT) j== and a, step-sizes of Frank-Wolfe. Conveniently, after the linearization of U the subproblem on the polytope D corresponds to an RL problem with reward Vu for which many efficient solvers exist. Namely, for a single mixture component we have,
d
arg min
Σ FU(d=mii...)a a)da(z a).
(10)
dED
v,a_h
Due convexity, the state-action distribution follows the convex combination, dmmix- =E== a.d The optimization produces a Markovian policy due to the subproblem in Eq. (TO) being optimized by one. We now detail how construct a non-Markovian policies by adaptive resampling.
4.1 Adaptive Resampling: Non-Markovian policies.
A core contribution of our paper is receding horizon re-planning. This means that we keep track of the past states visited in the current and past trajectories and adjust the policy at every step of the horizon H in each trajectory indexed by t. At h, we construct a Markov policy for a reward that depends on all past visited states. This makes the resulting policy history dependent, While in episode and time-poin h we follow a Markov policy for single step, the overall policy is history-dependent non-Marken policy.
We define the empirical state-action visitation distribution,
deh tH +H
1
Σ ΣΣ 6m
Σ δ2,a)
(11)
J=1 IEEETT
2.0ETs/1
visited states in past trajectories
states at ep. i up
where δz,a denotes a delta mass at state-action (x,, a). Instead of solving the objective U(d) as in Eq. I0). we seek to find correction to the empirical distribution by minimizing,
1 H h (c) = 4 и 1+t d+
(H=
tH + h
de.1
(12)
1
We use the same Frank-Wolfe machinery to optimize this objective: dr,,, arg min Ue (d+). The distribution d=,,, represents the density of the policy to be deployed at trajectory t and horizon counter h. We now need to solve multiple (n due to FW) RL problems at each horizon counter h. Despite this, for discrete MDPs, the sub- sub-problem can be solved extremely efficiently to exactness using dynamic programming. As can be seen in Appendix B.4, our solving times are just few seconds, even if planning for very long horizons. The resulting policy can be found by marginalization TA (a|z)) = d (x, a)/ Σ. dz,h(x, a), a basic property of MDPs (60
4.2 Continuous MDPs: Model Predictive Control
With continuous search space, the sub-problem can be solved using continuous RL solvers. However, this can be difficult. The intractable part of the problem is that the distribution d, needs to be represented in a computable fashion. We represent the distribution by the sequence of actions taken \T= with the linear state-space model, zh-1 Ax + Ban- While this formalism is not as general as could be, it gives us a tractable sub-problem formulation common control science scenario [61] that is practical for our experiments and captures vast array of problems. The optimal set of actions is solved with the following problem, where we state it for the full horizon H:
H
arg miin o(d. (In, (zh)an), an),
(13)
A=
such that ||q, < Omax Th X' and Øh+1 Ax, + Ban, where the known dynamics serves as constraints. Notice that instead of optimizing over the policy dn, we directly optimize over the parameterizations of the policy {ax 12 In fact, this formulation is reminiscent of the model predictive control (MPC) optimization problem. Conceptually, these are the same. The only caveat in our case is that unlike in MPC (2), our objective is non-convex and tends to focus gathering information rather than stability. Due to the non-convyxity in this parameterization, we need to solve it heuristically. We identify number of useful heuristics to solve this problem in Appendix G
7

---

0.05
GUBB
ratio of initial reactants (B)
1.0
Tue maximizer
Best Arm Identified (%)
Potential maximizers
Average Regret
MDP-8
0.75
Discarded states
0.5
0.5
0.25
0.0
0.00
0.0
20
30
Episode
residence time (T)
Figure 2: The Knorr pyrazole synthesis experiment. On the left, we show the quantitative results. The line plots denote the best prediction regret, while the bar charts denote the percentage runs that correctly identify the best arm at the end of each episode. On the right, we show ten paths in different colours chosen by the algorithm. The underlying black-box function is shown as the contours, and we can see the discretization as dots. We can see four remaining potential maximizers (in orange), which includes the true one (star). Notice all paths are non-decreasing in residence time, following the fransition constraints.
5 Experiments
Sections 5.1 5.3 showcase real-world applications under physical transitions constraints, using the discrete version of the algorithm. Section benchmarks against other algorithms in the continuous setting, where we consider the additive transition model of Section 4.2 with additional results in Appendix and fixed them during the optimization. These are summarized in Appendix [.2 As we are interested maximizer identification, in discrete problems, we report the proportion of reruns that succeed identifying the true maximum. For continuous benchmarks, we report inference regret at each iteration: Regret, f(x.) (xr,,), where Ftt arg maxzex (c All statisties reported over 25 different runs.
= I. We include B For each benchmark, we selected reasonable GP hyper-parameters
Baselines. We include a naive baseline that greedily optimizes the immediate reward to showcase a method with no planning (Greedy-UCB). Likewise, we create baseline that replaces the gradient in Eq. (TO) with Expected Improvement [63 (MDP-EI), weak version of planning. the continuous settings, compare against truncated SnAKe (TrSnaKe) [42]. which minimizes movement distance, and against local search region-constrained BayesOpt or LSR [22] for the same task. We compare two variants for approximating the set of maximizers, using Thompson Sampling (MDP-BO-TS) and one using Upper Confidence Bound (MDP-BO-UCB).
5.1 Knorr pyrazole synthesis
0.037
Greedy-UUch
1.0
0.739
NEP-80
1.0
State-less
Best Arm Identifeed (36)
Average Regret
MDP:2
Average Regret
2.5
Identifieed (%6
0.5
40
60
100
Episode
Iteration
(a) Monitoring Lake Ypacarai.
(b) Free-electron laser tuning.
Figure 3: Results for Ypacarai and free electron-laser tuning experiments. On the left, the line plots denote the best prediction regret, while the bar charts denote the percentage of runs that correctly identify the best arm the end of each episode. On the right, We plot the regret and compare against standard BO without accounting for movement-dependeen noise.

---

Reque
MOP
Moddonn Regret
Median Rognut
0.0
0.5-
6.0
10
40
50
(a) Michalewicz 3D. (synch)
(b) Hartmann 3D, (synch)
Hartmann 6D, (synch)
1.5
Modlun degpet
Nostan Nogred
3.0
64
0.
20
30
40
50
70
80
90
100
20
30
60
30
00
70
10
90
100
20
30
40
50
70
80
90
100
Ite
Poration
Baration
(d) Michalewicz 2D. (asynch)
(e) Hartmann 3D. (asynch)
(l) SnAr 4D. (asynch)
Figure 4: Results of experiments on the asynchronous and synchronous benchmarks. We plot the median predictive regret and the 10% and 90% quantiles. For the asynchronous experiments, we can that the paths taken by MDP-BO-TS to TrSnAKe. While in the asynchronous setting, we found creating the maximization set using Thompson Sampling gave stronger performance, in the synchronous setting. UCB is preferred. LSR gives very strong performance, comparable to MDP-BO-UCB in almost all benchmarks.
more consistent, and the final performance is comparable
Our chemical reactor benchmark synthetizes Knorr pyrzole in transient flow reactor. In this experiment, we can control the flow-rate (residence time) t and ratio of reactants B in the reactor. We observe product concentration discrete time intervals and intervals. Our goal is to find the best parameters of the reaction subject to natural movement constraints on B, and In addition, achieved. However, increasing the flow rate can lead to inaccurate readings [64]. A lower flow rate leads to higher residence time, so we impose that must be non-decreasing.
also change inputs these
assume decreasing the flow rate of reactor be easily
The kernel. Schrecker et al. 27 indicate the reaction can be approximately represented by simple kinetics via differential equation model. We this information along with techniques for representing linear ODE as constraints in GP fitting (65] [66] to create an approximate ODE kernel kode through the featurization:
odde (r, B) = (1
where (T, B) are equal to:
1
A(B)
1
2
for 1,2, where ( and X/2 (i)
eigenvalues of the linearized ODE at different stationary points, (B) (2) 1- B, and S(x) (1 e=0n)((--0.3))- is sigmoid function. Appendix Holld the details and derivations which may be of independent interest. As the above kernel is only an approximation the true ODE kernel, which itself imperfect, must account for the model mismatch. Therefore, we add squared exponential term to the kernel to ensure a non-parametric correction, i.e.: k(T, B) aode Kode (T, B) + arty (T, B).
We report the examples of the trajectories in the search space in Figure 2 Notice that all satisfy the transition constraints. The paths are not space- falceffilling and avoid sub-optimal arcas because of our choice non-isotropic kernel based on the ODE considerations. We run the experiment with episodic feedback, for 10 episodes of length 10 each, starting each episode with (Tn, B) (0, 0). Figure reports quantitative results and shows that the best-performing algorithm MDP-BO.
5.2 Monitoring Lake Ypacarai
Samaniego al. (20 investigated automatic monitoring of Lake Ypacarai, and Folch et al. and Yang et al. 40 benchmarked different BayesOpt algorithms for the task of finding the largest con- tamination source in the lake. We introduce local transition constraints to this benchmark by creating the lake containing obstacles that limit movement (see Figure 12 in the Appendix). Such obstacles

---

in environmental monitoring may include islands or protected areas for animals. We add an initial and final state constraint with the goal of modeling that the boat has to finish at a maintenance port.
We focus on episodic feedback, where each episode consists of 50 iterations. Results can be seen in Figure [a] MDP-EI struggles to identify the maximum contamination for the first few episodes. On the other hand, our method correctly identifies the maximum in approximately 50% of the runs by episode two and achieves better regret.
5.3 Free-electron laser: Transition-driver corruption
Apart from hard constraints, can apply our framework to sate-dependent BayesOpt problems involving transitions. For example, the magnitude of noise e may depend on the transition. This occurs in systems observing equilibration constraints such free-clectron laser [23]. Using the simplified simulator of this laser [67, we use our framework to model heteroscedastic noise depending on the difference between the current and next state, (2, x²) s(1 + w||æ - x' 12112). By choosing A = X, we rewrite the problem as (8, a) = s(1 + w|| all2). The larger the move, the more noisy the observation. This creates a problem, where the BayesOpt needs to balance between informative actions and movement, which can be directly implemented in the objective (8) via the matrix (dz)) Ex.ex (x, a) a) [azaa
(xx) Th Figure 3b reports the comparison between worst-case stateless BO and our algorithm. Our approach substantially improves performance.
5.4 Synthetic Benchmarks
We benchmark constraints and considering both immediate and asynchronous feedback (by introducing an observation delay of 25 iterations). We also include the chemistry SnAr benchmark, from Summit 168). which we treat as asynchronous as per Folch et al. [II]. Results are in Figure 4 In the synchronous setting, we found using the UCB maximizer criteria for MDP-BO yields the best results (c.f. Appendix for details of this variant). We also found that LSR performs very competitively on many benchmarks, frequently matching the performance of MDP-BO. In the asynchronous settings we achieved better results using MDP-BO with Thompson sampling. TrSnAKe baseline appears to be competitive all synthetic benchmarks as well. However, MDP-BO is more robust having less variance in the chosen paths as seen in the quantiles. It is important to highlight that SnAKe and LSR are specialist heuristic algorithms for local box-constraints, and therefore is not surprising they perform strongly. Our method can be applied to more general settings and therefore it is very encouraging that MDP-BO is able to match these SOTA algorithms in their specialist domain.
variety of classical BayesOpt problems while imposing local movement
6 Conclusion
We considered conistioon--ntttrineee BayesOpt problems arising in physical sciences, such as chemi- cal reactor optimization, that require careful planning to reach any system configuration. Focusing on maximizer identification, we formulated the problem with transition constraints using the framework of Markov decision processes and constructed a tractable algorithm for provably and efficiently solving these problems using dynamic programming or model predictive control sub-routines. We showcased strong empirical performance in a large variety of problems with physical transitions, and achieve state-of-the-art results in classical BayesOpt benchmarks under local movement constraints. This work takes an important step towards the larger application of Bayesian Optimization to real-world problems. Further work could address the continuous variant of the framework to deal with more general transition dynamics, or explore the performance of new objective functions.
Acknowledgments
JPF is funded by EPSRC through the Modern Statisties and Statistical Machine Learning (StatML) CDT (grant no. EP/S023151/I) and by BASF SE, Ludwigshafen am Rhein. RM acknowledges support from the BASF/ Royal Academy of Engineering Research Chair Data-Driven Optimisation. This publication was created as part of NCCR Catalysis (grant number 180544), National Centre of Competence in Research funded by the Swiss National Science Foundation, and was partially supported by the European Research Council (ERC) under the European Union's Horizon 2020 research and Innovation Program Grant agreement no. 815943. We would also like to thank Lindee Schrecker, Ruby Sedgwick, and Daniel Lengyel for providing valuable feedback the project.
10

---

References
[1] Peter Frazier. A tutorial on Bayesian optimization. arXiv preprint arXiv: 180..0281, 2018.
[2] Donald R Jones, Matthias Schonlau, and William Welch. Efficient global optimization of expensive black-box f'unctions. Journal of Global Optimization, 13(4):455-492, 1998.
[3] Bobak Shahriari, Kevin Swersky, Ziyu Wang, Ryan P. Adams, and Nando de Freitas, Taking the human out of the loop: A review of Bayesian optimization. Proceedings of the IEEE, 104 (1):148-175, 2016.
[4] Frank Hutter, Lars Kotthoff, and Joaquin Vanschoren, editors. Automated Machine Learning Methods, Systems, Challenges. Springer, 2019.
[5] Alonso Marco, Felix Berkenkamp, Philipp Hennig, Angela Schoellig, Andreas Krause, Stefan Schaal, and Sebastian Trimpe. Virtual real: Trading off simulations and physical experiments reinforcement learning with Bayesian optimization. In 2017 IEFE International Conference on Robotics and Automation (ICRA), pages 1557-1563. IEEE, 2017.
[6] Jose Pablo Folch, Robert M Lee, Behrang Shafei, David Walz, Calvin Tsay, Mark van der Wilk, and Ruth Misener. Combining multi-fidelity modelling and asynchronous batch Bayesian optimization. Computers & Chemical Engineering, 172:108194, 2023.
[7] Thomas M Dixon, Jeanine Williams, Maximilian Besenhard, Roger M Howard, James MacGre- gor, Philiip Peach, Adam D Clayton, Nicholas Warren, and Richard Bourne. Operator-free HPLC automated method development guided by Bayesian optimization. Digital Discovery, 3 (8):1591-1601, 2024.
[8] Jennifer Brennan, Lalit Jain, Sofia Garman, Ann E Donnelly, Erik Scott Wright, and Kevin Jamieson. Sample-efficient identification of high-dimensional antibiotic synergy with a normal- ized diagonal sampling design. PLOS Computational Biology, 1807)=1010311, 2022.
[9] Alexander Thebelt, Johannes Wiebe, Jan Kronqvist, Calvin Tsay, and Ruth Misener. Maximizing information from chemical engineering data sets: Applications to machine learning. Chemical Engineering Science, 252: 17469, 2022.
[10] Pratibha Vellanki, Santu Rana, Sunil Gupta, David Rubin, Alessandra Sutti, Thomas Dorin, Murray Height, Paul Sanders, and Svetha Venkatesh. Process-constraing batch Bayesian optimisation. In Advances in Neural Information Processing Systems. Curran Associates, Inc., 2017.
[11] Jose Pablo Folch, Shiqiang Zhang, Robert Lee, Behrang Shafei, David Walz, Calvin Tsay, Mark van der Wilk, and Ruth Misener. SnAKe: Bayesian optimization with pathwise exploration. In Advances in Neural Information Processing Systems, volume 35, pages 252299,,,,,, 2022.
[12) Sarah L. Boyall, Holly Clarke, Thomas Dixon, Robert WM Davidson, Kevin Leslie, Graeme Clemens, Frans Muller, Adam Clayton, Richard Bourne, and Thomas W Chamberlain. Automated optimization of multistep. multiphase continuous flow process for pharmaceutical synthesis. ACS Sustainable Chemistry & Engineering, 2024.
[13] 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.
[14] Sigrid Passano Hellan, Christopher G Lucas, and Nigel H Goddard. Bayesian optimisation for active monitoring of air pollution. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 36, pages 1198--1916. 2022.
[15] Sigrid Passano Hellan, Christopher G Lucas, and Nigel H Goddard. Bayesian optimisation against climate change: Applications and benchmarks. Data-centric Machine Learning Research (DMLR) Workshop at the 40th International Conference on Machine Learning, 2023,
11

---

[16] Clara Stoddart, Lauren Shrack, Richard Sserunjogi, Usman Abdul-Ganiy Engineer Baino- mugisha, Deo Okure, Ruth Misener, Jose Pablo Folch, and Ruby Sedgwick. Gaussian processes for monitoring air-quality in Kampala. NeurlPS 2023 Workshop: Tackling Climate Change with
Machine Learning. 2023.
[17] Mojmir Mutný, Tadeusz Janik, and Andreas Krause. Active Exploration via Experiment Design in Markov Chains. In Proceedings of The 26th International Conference Artificial Intelligence and Statistics, volume 206, pages 7349-7374, 2023.
[18] Gregory Hitz, Alkis Gotovos, Marie-Eve Gameau, Cédric Pradalier, Andreas Krause, Roland Y Siegwart, et al. Fully autonomous focused exploration for robotic environmental monitoring. In 2014 IEEE International Conference on Robotics and Automation (ICRA), pages 2658-2664. IEEE, 2014.
[19] Alkis Gotovos, Nathalie Casati, Gregory Hitz, and Andreas Krause. Active learning for level set estimation. In IJCAI 2013, 2013.
[20] Federico Peralta Samaniego, Daniel Gutiérrez Reina, Sergio L Toral Marín, Mario Arzamendia, and Derlis Gregor. A Bayesian optimization approach for water resources monitoring through an autonomous surface vehicle: The Ypacraii lake case study. IEEE Access, 09916391799 2021.
[21] Shyam Sundhar Ramesh, Pier Giuseppe Sessa, Andreas Krause, and Ilija Bogunovic. Movement penalized Bayesian optimization with application to wind energy systems. In Advances in Neural Information Processing Systems, volume 35, pages 27036-27048 2022.
[22] Joel A Paulson, Farshud Sorouifar, Christopher Laughman, and Ankush Chakrabarty. LSR- BO: Local search region constrained Bayesian optimization for performance optimization of vapor compression systems. 2023 American Control Conference (ACC), pages 576-582. IEEE, 2023.
[23] Johannes Kirschner, Mojmír Mutný, Andreas Krause, Jaime Coello de Portugal, Nicole Hiller, and Jochem Snuverink. Tuning particle accelerators with safety constraints using Bayesian optimization. Phys. Rev. Accel. Beams, 25:062802, 2022.
[24] Matthew Sullivan, John Gebbie, and John Lipor. Adaptive sampling for seabed identification from ambient acoustic noise. IEEE CAMSAP, 2023.
[25] Sergey Mozharov, Alison Nordon, David Littlejohn, Charlotte Wiles, Paul Watts, Paul Dallin, and John M Girkin. Improved method for kinetic studies in microreactors using flow manipula- tion and noninvasive Raman spectrometry. Journal of the American Chemical Society, 133(10); 3601-3608, 2011.
[26] Linden Schrecker, Joachim Dickhaut, Christian Holtze, Philipp Staehle, Andy Wieja, Klaus Hellgardt, and King Kuok Mimi Hii. An efficient multiparameter method for the collection of chemical reaction data via 'one-pot" transient flow. Reaction Chemistry & Engineering, 8(12): 3196-3202, 2023.
[27] Linden Schrecker, Joachim Dickhaut, Christian Holtze, Philipp Staehle, Marcel Vranceanu, Klaus Hellgardt, and King Kuok Mimi Hi. Discovery of unexpectedly complex reaction pathways for the Knorr pyrazole synthesis via transient flow. Reaction Chemistry Engineering, 8(1):41-46, 2023.
[28] Carl Edward Rasmussen and Christopher K. I. Williams. Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning). The MIT Press, 2005.
[29] Philipp Hennig and Christian J Schuler. Entropy search for information-efficient global opti- mization. Journal of Machine Learning Research, 13(6), 2012.
[30] Tanner Fiez, Lalit Jain, Kevin G Jamieson, and Lillian Ratliff. Sequential experimental design for transductive linear bandits. Advances in neural information processing systems, 32, 2019.
[31] Christopher Williams and Matthias Seeger. Using the Nyström method to speed up kernel machines. Advances in neural information processing systems, 13, 2000.
12

---

[32] Mojmír Mutný and Andreas Krause. Efficient high dimensional Bayesian optimization with additivity and quadrature fourier features. Advances in Neural Information Processing Systems, 31, 2018.
[33] Ali Rahimi and Benjamin Recht. Random features for large- scale kernel machines. Advances in neural information processing systems, 20, 2007.
[34] Elad Hazan, Sham Kakade, Karan Singh, and Abby Van Soest. Provably efficient maximum entropy exploration. In International Conference on Machine Learning, pages 2681-2691. PMLR, 2019.
[35] Tom Zahavy, Brendan O'Donoghue, Guillaume Desjardins, and Satinder Singh. Reward is enough for convex MDPs. Advances in Neural Information Processing Systems, 34:25746- 25759, 2021.
[36] Zihan Li and Jonathan Scarlett. Gaussian process bandit optimization with few batches. In International Conference on Artificial Intelligence and Statistics, pages 92-107, 2022.
[37] Fengxue Zhang, Jialin Song. James 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 on Machine Learning, pages 41579 41595. PMLR, 2023.
[38] Minbiao Han, Fengxue Zhang, and Yuxin Chen. No-regret learning of Nash equilibrium for black-box games via Gaussian processes. arXiv preprint arXiv:2405.0833 083318 2024.
[39] Sudeep Salgia, Sattar Vakili, and Qing Zhao. Random exploration in Bayesian optimization: Order-optimal regret and computational efficiency. In Forty-first International Conference Machine Learning. 2024.
[40] Adam X Yang, Laurence Aitchison, and Henry Moss. MONGOOSE: Path-wise smooth
Bayesian optimisation via meta-learning. arXiv preprint arXiv:2302.11533, 2023.
[41] Qiyuan Chen and Raed Al Kontar. The traveling bandit: A framework for Bayesian optimization with movement costs. arXiv preprint arXiv:2410. 145333,..., 2024.
[42] Jose Pablo Folch, James Odgers, Shiqiang Zhang, Robert M Lee, Behrang Shafei, David Walz, Calvin Tsay, Mark van der Wilk, and Ruth Misener. Practical path-based NeurlPS 2023 Workshop on Adaptive Experimental Design and Active Learning in the Real World, 2023.
based Bayesian optimization.
[43] Jixiang Qing, Becky D Langdon, Robert M Lee, Behrang Shafei, Mark van der Wilk, Calvin Tsay, and Ruth Misener. System-awar neural ODE processes for few-shot Bayesian optimiza- tion. arXiv preprint 0rXiv:2406.02352, 2024.
[44] Ethan Che, Jimmy Wang, and Hongseok Namkoong. Planning contextual adaptive experiments with model predictive control. NeurIPS 2023 Workshop on Adaptive Experimental Design and Active Learning in the Real World, 2023.
[45] Jvan-Yves Audibert and Sébastien Bubeck. Best arm identification in multi-armed bandits. In Conference on Learning Theory, pages 13-p, 2010.
[46] Marta Soare, Alessandro Lazaric, and Rémi Munos. Best-arm identification in linear bandits. Advances in Neural Information Processing Systems, 27, 2014.
[47] Kai Yu, Jinbo Bi, and Volker Tresp. Active learning via transductive experimental design. In Proceedings of the 23rd international conference Machine learning, pages 1081-1088. 2006.
[48] David Ginsbourger and Rodolphe Le Riche. Towards Gaussian process-based optimization with finite time horizon. In Advances in Model-Oriented Design and Analysis: Proceedings of the 9th International Workshop in Model-Oriented Design and Analysis, pages 89-96. Springer, 2010.
[49] Roman Marchant, Fabio Ramos, Scott Sanner, et Sequential Bayesian optimisation for spatial-temporal monitoring. In UAI, pages 553-562, 2014.
13

---

[50] Remi Lam, Karen Willcox, and David H Wolpert. Bayesian optimization with a finite budget: An approximate dynamic programming approach. Advances in Neural Information Processing Systems, 29, 2016.
[51] Shali Jiang, Henry Chai, Javier Gonzalez, and Roman Garnett. Binoculars for efficient, nonmy- opic sequential experimental design. In International Conference on Machine Learning, pages 4794 4803. PMLR, 2020.
[52] Eric Hans Lee, David Eriksson, Valerio Perrone, and Matthias Seeger. A nonmyopic approach to cost-constrained Bayesian optimization. In Uncertainty Artificial Intelligence, pages 368-577. PMLR, 2021.
[53] Joel A Paulson, Farshud Sorouifar, and Ankush Chakrabarty. Efficient multi-step lookahead Bayesian optimization with local search constraints. In 2022 IEEE 6 Conference on Decision and Control (CDC), pages 123-129. IEEE, 2022.
[54] Mujin Cheon, Haeun Byun, and Jay H Lee. Reinforcement learning based multi-step look-ahead bayesian optimization. PWAC-PapessOnl.ing 55(7):100-105, 2022.
[55] Shali Jiang, Daniel Jiang, Maximilian Balandat, Brian Karrer, Jacob Gardner, and Roman Garnett. Efficient nonmyopic bayesian optimization via one-shot multi-step trees. Advances in Neural Information Processing Systems, 33:18039-18049, 2020.
[56] Raul Astudillo, Daniel Jiang, Maximilian Balandat, Eytan Bakshy, and Peter Frazier. Multi-step budgeted bayesian optimization with unknown evaluation costs. Advances in Neural Information Processing Systems, 34:20197-20209, 2021.
[57] William R Thompson. On the likelihood that one unknown probability exceeds another in view of the evidence of two samples. Biometrika, 25(3-4):285-294, 1933.
[58] Kirthevasan Kandasamy, Akshay Krishnamurthy. Jeff Schneider, and Barnabás Póczos. Paral- lelised Bayesian optimisation via Thompson sampling. In International Conference on Artificial Intelligence and Statistics, pages 133-142. PMLR, 2018.
[59] Martin Jaggi. Revisiting Frank- Wolfe: Projection-Free Free Sparse Convex Optimization. In Sanjoy Dasgupta and David McAllester, editors, Proceedings of the 30th International Conference on Machine Learning, volume 28, pages 427-435. PMLR, 17-19 Jun 2013.
[60] Martin L Puterman. Markov decision processes: discrete stochastic dynamic programming. John Wiley Sons, 2014.
[61] James Blake Rawlings, David Q Mayne, and Moritz Diehl. Model predictive control: theory, computation, and design, volume 2. Nob Hill Publishing, 2017.
[62] Carlos E. García, David M. Prett, and Manfred Morari. Model predictive control: Theory and practice-A survey, Automatica, 25(3)333--34,, 1989.
[63] Vydu-nas Šaltenis. One method of multiextremum optimization. Avtomatika i Vychislitel naya Tekhnika (Automatic Control and Computer Sciences), 5(3):33 38, 1971.
[64] Linden Schrecker, Joachim Dickhaut, Christian Holtze, Philipp Staehle, Marcel Vranceanu, Andy Wieja, Klaus Hellgardt, and King Kuok Hii. A comparative study of transient flow rate steps and ramps for the efficient collection of kinetic data. Reaction Chemistry & Engineering, 2024.
[65] Mojmár Mutný and Andreas Krause. Experimental design for linear functionals in reproducing kernel Hilbert spaces. Advances in Neural Information Processing Systems, 35:20175 7501011--1111100 88, 2022.
[66] Andreas Besginow and Markus Lange-Hegrrmann. Constraining Gaussian processes to systems of linear ordinary differential equations. Advances Neural Information Processing Systems, 35:29386-2998 29399,, 2022.
14

---

[67] Mojmir Mutný, Johannes Kirschner, and Andreas Krause. Experimental design for optimization of orthogonal projection pursuit models. In Proceedings of the 34th AAAI Conference on Artificial Intelligence (AAAT), 2020.
[68] Kobi Felton, Jan Rittig, and Alexei Lapkin. Summit: Benchmarking Machine Learning Methods for Reaction Optimisation. Chemistry Methods, February 2021.
[69] Tianbao Yang, Yu-Feng Li, Mehrdad Mahdavi, Rong Jin, and Zhi-Hua Zhou. Nyström method random Fourier features: A theoretical and empirical comparison. Advances in neural information processing systems, 25, 2012.
[70] Daniel Russo. Simple bayesian algorithms for best arm identification. In Conference on Learning Theory, pages 1417-1418. PMLR, 2016.
[71] Niranjan Srinivas, Andreas Krause, Sham M Kakade, and Matthias Seeger. Gaussian pro- cess optimization in the bandit setting: No regret and experimental design. arXiv preprint arXiv:0912.3995, 2009.
[72] Felipe Cucker and Steve Smale. On the mathematical foundations of learning. Bulletin of the American mathematical society, 39(1): 2002.
[73] Kathryn Chaloner and Isabella Verdinelli. Bayesian Experimental Design: A Review. Statist. Sci., 10(3):273-304, 08 1995.
[74] Friedrich Pukelsheim. Optimal Design of Experiments (Classics in Applied Mathematics) (Clas- sics in Applied Mathematics, 50). Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 2006. ISBN 0898716047.
[75] Stephen Boyd and Lieven Vandenberghe. Convex optimization. Cambridge university press, 2004.
[76] Max A Woodbury. Inverting modified matrices. Memorandun Rept. 42, Statistical Research Group, page 4. Princeton Univ., 1950.
[77] José Miguel Hernándezz-Lobtto Matthew W Hoffman, and Zoubin Ghahramani. Predictive entropy search for efficient global optimization of black-box functions. Advances in neural information processing systems, 27, 2014.
[78] Ben Tu, Axel Gandy, Nikolas Kantas, and Behrang Shafei. Joint entropy search for multi- objective Bayesian optimization. Advances in Neural Information Processing Systems, 35: 9922 9938, 2022.
[79] Carl Hvarfner, Frank Hutter, and Luigi Nardi. Joint entropy search for maximally-informed Bayesian optimization. Advances in Neural Information Processing Systems, 35:1 1494-11506, 2022.
[80] Zi Wang and Stefanie Jegelka. Max-value entropy search for efficient Bayesian optimization. In International Conference Machine Learning, pages 3627-3635. PMLR, 2017.
15

---

A Visual abstract of the algorithm
In Figure | we summarize how our algorithm creates -an-Mvaoviian policies for maximizer identifi- cation and the corresponding connections to other works in the literature.
Transition Constrained BO Problem Formulation
(a) Need plan ahead
Hypothesis Testing
Objective Formulation
(b) Variance reduction of trajectory
Unconstrained Linear Bandits (Soare et al. [46]) (frequentist)
Look-ahead Utilities and averaging out uncertainty (intractable)
Adaptive Objective with optimality guarantees for Linear Bandits (Fiez al. 30) (frequentist)
(c) Acquisition function
Acquisition function
Direct optimization the space of trajectories (intractable)
(d) Relaxation to space of state-action distributions
optimization
I (e) Solvable by Frank-Wolfe Algorithm
Mutaý et al. Hazan et al. 14
1
Efficient RL solvers, e.g. Dynamic Pro- gramming
Convergence to optimal Markovian policy (Mutný et al. (17)
(f) Iterative Reinforcement
Learning sub-problem
1
Adaptive resampling at every iteration
(g) Non-Markovian Policies
Figure 5: Visual abstract of the work. In black show the method presented in this paper, with literature connections shown in blue. In red we show solutions which we did not pursue due to intractability. The problem creates the (a) need to plan ahead. To do this, from hypothesis testing and focus on (b) the variance reduction in a set of maximizers, which leads our (c) acquisition function. The objective is the same as Fiez al. (30) introduced in the linear bandits literature from a frequentist perspective. To optimize it, we follow developments in Mutný al. [17), Hazan et al. 34 by (d) relaxing the acquisition function to the space of state-action distributions and (e) solving the planning problem using the Frank-Wolfe algorithm.
take inspiration
This consists of iteratively solving tractable (f) reinforcement learning sub-problems which give us optimal Markov policies. We then apply adaptive resampling to obtain (g) non-Markovian policies.
16

---

1.0
Greedy-UCB
Best Arm Identified (%))
MDRE
MDP-30
Log Regret
0
0.2
-12
0.0
Iteration
Figure High noise constrained Ypacari experiment with immediate feedback.
1.0
UCB
MDP MDP-80
Average Log Regret
0.8
Best Amm Identified (%)
-0.
12
0.0
100
Iteration
Figure 7: Knorr pyrazole synthesis with immediate feedback
B
Additional Empirical Results
B.1 Constrained Ypacarai
We also run the Ypacarai experiment with immediate feedback. To increase the difficulty, we used large observation noise, o² 0.01. The results can be seen in Figure 6 The early performance of MDP-EI is much stronger, however, it gets overtaken by our algorithm from episode three onwards, and gives the worst result at the end, as struggles to identify which the two optima is the global one.
B.2 Knorr pyrazole synthesis
We also include results for the Knorr pyrazole synthesis with immediate feedback. this case we observe very strong early performance from MDP-BO, but by the end MDP-EI comparable. The greedy method performs very poorly.
B.3 Additional synthetic benchmarks
Finally, we also include additional results on more synthetic benchmarks for both synchronous and asynchronous feedback. The results are shown in Figures (nnur The results back the conclusions the main body. All benchmarks BO-UCB and LSR can be much stronger in the synchronous setting than Thompson Sampling planning-baa based approaches (with the one exception of the Levy function).
well in 2-dimensions while highlighting further that MDP-
Table 1: Average acquisition function solving times for each practical benchmark. We give the solving times the nearest second, and provide the size of the state-space, (S), the maximum number of actions one can take from specific state, [A()),, and the planning horizon. In all benchmarks we are able to solve the problem few seconds.
Benchmark Solve Time
S
Maximum A(S)
Planning horizon
Knorr pyrazole
Is
100
6
10
Ypacarai Electron laser
3s
100
8
50
15s
100
100
100
17

---

Nedien Beger
MDP BO-TS
HDP -0-05
MOM SC-UCB
NDP-80-UCB
MDP 6O LCB
T/S
TSSnAK
Mosiun
Modian
05
Iteration
Raration
koration
(a) Branin 2D. (asynch)
(b) Michalewicz 3D. (asynch)
(c) Hartmann 6D. (asynch)
Figure 8: Additional asynchronous results.
6.5
0.20
Reque
C4.
MDN-80-UCB
Wedan Begee
HDP BC-UCB
MIP--OO-TS
0.00
MENNIO-IC
NDN-8O-1T
PSSAMic
Lbit
0.16
6.04
0.05
0.02
0.00
10
40
20
10
30
60
20
80
Barsstio
Bhe Berratoo
(a) Branin 2D.
(h) Michalewiez2D.
(c) Levy4D.
Figure 9: Additional synchronous results.
B.4 Computational study
We include the average acquisition function solving time for each of the discrete problems. For the continuous case the rumning time was comparable to Truncated SnAKe [42) since most of the computational load was create the set of maximizers using Thompson Sampling. The times were obtained in a simple 2015 MacBook Pro 2.5 GHz Quad-Cor Intel Core i7. The bulk of the experiments was ran in parallel on High Performance Computing cluser, equipped with AMD EPYC 7742 processors and 16GB of RAM.
B.5 Median plots for Ypacarai and reactor experiments
In Figures [0 11
[[ we give the median and quantile plots for the Knorr pyrazole synthesis and the Ypacarai experiment, which were included in the main paper avoid cluttering the graphics.
Greedy-UCB
Greedy-UCB
0.05
0.06
MDP-BO
MDP-B0
0.05
Regret
Regret
0.04
0.03
0.03
0.02
0.01
0.01
0.00
0.00
2
5
6
10
20
30
40
50
60
70
80
90
100
Episode
Iteration
(a) Episodic feedback results
(b) Immediate feedback results
Figure 10: Median and 10th/90th quantile plots for Knorr pyrazole synthesis experiment.
C Utility function: Additional Info
We describe the utility function in complete detail using the kernelized variant that allows to extend the utility beyond the low-rank assumption in the main text.
C.1 Derivation of the Bayesian utility
Suppose that our decision rule is to report the best guess of the maximizer after the T steps as,
FT = arg max (x:)
18

---

0.12
Greedy-UCB
Greedy-UCB
MDP-EI
MDP-EI
MDP-80
0.12
MDP-80
0.08
Regret
Regret
0.10
0.06
0.06
0.04
0.06
0.04
0.02
100
200
500
Episode
Iteration
Episodic feedback results
(b) Immediate feedback results
Figure 11: Median and 10th/90th quantile plots for Ypacarai experiment.
We call this the selection the recommendation rule. We focus on this recommendation rule this rule is interpretable the facilitator of the analysis and experimenters. In this derivation we use that f f= ((x). More commonly, the notation (0, ((x) used, where the inner product is potentially infinite dimensional. We use use T notation for simplicity for both cases. Same is true for any other functional estimates, e.g., for the posterior mean estimate, we use (x(x)) - ((x) Pe. The inner product is in the reproducing kemell Hilbert space associated with the kemel k.
Now, suppose there given (we will take expectation over later), then there is an æ X achieving optimum value, denoted zz (suppose unique for this development here). Hence, we would like to model the risk associated with predicting a fixed 2 * z²,, which still in 2 at time T. Suppose we are at time I, we develop the utility to gather additional data Xnew on top of the already acquired data X,, These should improve the discrepancy of the true answer, and the reported value the most.
Suppose there are two elements in simple = {z, 1. We will generalize to a composite hypothesis later. In two-element case, the probability of the error incurred due to selecting 2 is:
P(ur(=) pt (w+) 0lʃ)
The randomness here is due the observations y = f (Xmmew) that are used to the estimator ur(x), Namely due to (N, o²). Given / (equivalently 0), the distribution of our estimator (namely the posterior mean) Gaussian. Hence, given f:
PT N(VT V00, + Ix) -a))"VV((Vr Ix)
where V= ET (Z.)) is an operator on the reproducing kernel Hilbert space due to
as > H. and In the identity operator on the same space.
This is the posterior over the posterior mean as a function. the specific evaluation is (T(2)
A
μT NO (x) Vr ( ( ((()))) ((x))) (Vr + IH) V In) (1(2)
posterior
(Q( -xxx)))
We can now bound the probability of making an error using a Gaussian tail bound inequality:
P(uT(2) (r(x*) 0) P(uT(2) (T(x*)) 2 a (-a=))
20
with the caveat that the inequality only holds when the 0₂ is negative. However note that f(=)-((** < as - xo therefore will hold once is large enough. From this we can take logarithms and then the expectation across the randomness in the GP:
1 E Esar [log P(ur(=) pr (z)))) E
which is called the log Bayes' factor and is expected log failure rate for the set of potential maximizers Zimple- The expectation is over the posterior including the evaluations X, (or prior at the very beginning of the procedure). In fact, can think of the posterior being the new prior for the future
19

---

at any time point, Now assuming that Z has more than one additional element, we want to ensure the failure rate is small for all other failure modes, all other hypothesis. Technically this means, we have an alternate hypothesis, which is composite composed of multiple point hypotheses. We take the worst-case perspective as its common with composite hypotheses. In expectation over the prior, we want to minimize:
log ((ur(z) JuT (x 2 (ule
min E, X
sup (x))(x))
(14)
For moderate to large T >> O, we can upper bound this objective via elegant argument yield very transparent objective:
(J(=) - (())))
min Ey Xura
sup =E2((m7)
log ((uT(z)--TtT) pT(z7) 2 (l
< min Es,
sup sE2\(x3)
2 X..
1xx, (z, 137)
(15)
where we have used an lower and upper bound on the u2 and b. respectively as follows:
a²
((2))))
((z])) -n)
(0T - (((=)))(rr))) (V. + I2) 1x)=1(((()))))) - ((z)))))
T
(0'()())-((4)))
- -(())))()))))))))) (f(z) - ((x7))²
62
(0 -
+ In] ( (0(2) b(2)))
o² (9(2) {(2))) (V. (+()) (<)+))
lx (x, x)).
In the last line we have used the same identity as in Eq. 25 We will explain how to eliminate the expectation in Section C.2
C.2 Upper-bounding the objective: Eliminating B, for large T.
The objective Eq. (15) is intractable due to the expectation of the prior and which involves expec- tation over the maximum (x which is known be very difficult estimate. Interestingly, the denominator is independent of f if we adopt the worst-case perspective over the æý, and hence the only dependence is through the set Z as well as the denominator. Given all current prior information, we can determine Z, and hence split the expectation. Let us now express
At any time point, we can upper- bound the denominator by the minimum as done by Fiez et al. [30] Even if Z decreases, as we get more information, the worst-case bound is always proportional the smallest gap gap(f) between two arms in X. Hence, we can upper bound the objective as:
(f(=)- [(r3)) =EZ\(az) (2,
E, [gap(f)]
-E, sup
sup =EZ\(a7)
kxux (2,
Ey [gap(f)] zEZ)aa)) Var [(=) - f(x))
sup
Xt [x
New
As the constant in front of the objective does not influence the optimization problem, we do not need consider it when defining the utility. Furthermore, in order minimise the probability of an error we can just minimise the variance in the denominator instead (since arg min,
-g(x) is equivalent to
arg min, ş(z) when g(x) > 0). However, the non-trivial distribution of f(x*) (29 renders the utility
intractable; therefore we employ a simple and tractable upper bound on the objective by minimising the uncertainty among all pairs in 2:
U U(Xmew)
max 2E,
Varlf()) f(2') U X.ew]
(16)
Surprisingly. this objective coincides with the objective from Fiez et al. [30] which has been derived as lower bound to the best-arm identification problem (maximum identification) with linear bandits. Their perspective is however slightly different they try to minimize T for fixed å failure rate. Perhaps it should not be surprising that the dual variant, consider here, for fixed T and trying to minimize the failure rate leads to the same decision for large T when log(b.) can be neglected.
20

---

C.3
Approximation of Gaussian Processes
Let us now briefly summarize the Nyström approximation 31 data-set X, we can choose sub-sample of the data 21, ................................ d'm Using this sample, we can create low r-rank approximation of the full kernel matrix
69
Given a kernel :), and a
K, K K₁
where K₁- [x; N K [k +5)] mm and K"t denotes the pseudo-inverse operation. We can then define the Nyström features as:
0n(2) - D, V (k(r, ((x, xmm))²
(17)
where D, is the diagonal matrix of non-zero eigenvalues of K, and V. the corresponding matrix of eigenvectors. It follows that we obtain tnite-dimensional estimate of the GP:
f(x) x D(z)"0
(18)
where ((x) = (I)
and are weights with a Gaussian prior.
C.4 Theory: convergence to the optimal policy
The fact that our objective is derived using Bayesian decision theory makes well-rooted in theory. In addition to the derivation of Section we can prove that our scheme is able to converge in terms of the utility.
Notice that the set of potential maximizers is changing over time, and hence we add a time subscript as Z₁. Let us contemplate for second what could the optimal policy. As the set of Z, is changing. we follow the line of work of started by Russo [70] and introduce an optimal algorithm that knows the true æ for each possible realization of the prior f. other words, an algorithm that any time t, would follow:
di
min EG dED
max (E2()()))
where in the above de + d represents the weighted sum as in the main Algorithm that scales the distributions properly according to and T, so to make the sum of them a valid distribution. Notice that in contrast to our objective, it does not take the maximum over 2' Z. but fixes it to the value z that the hypothetical algorithm has privileged access to. To eliminate the cumbersome notation, we will refer to the objectives (dd] 2,) as the objective used by our algorithm (real execution) and {2y }) the objective that the privileged algorithm is optimizing which serves as theoretical baseline.
The visitation of d; represents the best possible investment of the resources (of the size T - () to execute at time f had we known the zÿj instead of only Z, This is interpreted as if the modeler knows r.. and sets up an optimal curriculum that is heing shown to an observer in order to convince him/her of that æt is the optimal value. He or she is using statistical testing to elucidate it from execution of the policy. Like the algorithm, the optimal policy changes along the optimization procedure due to changes in Z, Hence, our goal is to show that we are closely tracking the performance of these optimal policies in time I, and eventually there is little difference between our sequence of executed policies (visitations) d, and the algorithm optimal d;.
In order to prove the theorem formally, we need to assume that Z, is decreasing. The rate at which this set is decreasing determines the performance of algorithm to a large extent. Namely, we assume that given two points in time, having the same empirical information de. Given, f. suppose
sup Id (VU(d, (2 2,) - (ud Ze, (x (z7)))1 S Ce.
dED
As we gather information in our procedure the, {a}} Ze S Z,-1. but the exact decrease depends on how Z, is constructed. We leave the particular choice for C. to make the above hold for future work. We conjecture that this is decreasing as C2 & In where is the information gain due to Srinivas et al. (71). We are now ready to state the formal theorem along with its assumptions.
21

---

Proposition C.1. Assuming episodic feedback, feedbac,, and suppose that for any Z,
1. U is convex D
2. B-locally Lipschitz continuous under || 11.0 norm
3. locally smooth with constant L, i.e,
Im.a
U(n + aa) S s(()) + vu(n) +
11bl/2
(19)
2
for a E (0,1) 1) and n, h E Ap, L == maxn,a Lyga
4. condition in ( holds with Bayesian posterior inference,
we can show that the Algorithm 1 satisfied for the sequences of iterates {di)i=1²
L Llog T
B VT
u(dolZe, (x7}) (xd))
log
with probability 1 on the sampling from the Markov chain. The randomness on the confidence set is captured by Assumption in Eq.
The previous proposition shows that as the budget
experimental campaign T is increasing, we are increasingly converging to the optimal allocation of the experimental resources average also on the objective that is unknown to us. In other words, our algorithm is becoming approximately optimal also under the privileged information setting representing the best possible algorithm. Despite having a limited understanding of potential maximizers at the beginning by following our procedure, we show that we are competitive the best possible allocation of the resources. Now, we prove the Proposition. The proof is an extension of the Theorem 3 in [7]. Whether the objective satisfied the above conditions depends the set X. Should the objective not satisfy smoothness, it can be easily extended by using the Nesterov smoothing technique as explained the same priorly cited work.
Proof of Proposition C.1 The proof is based on the proof of Frank-Wolfe convergence that appears Appendix B.4 in Thm. in Mutný al. [17.
Let us start by notation. We will the notation that U4 the privileged objective Udd {af }). while the original objective will be specified as (d] Z,, Z,).
First, what we follow in the algorithm:
q: arg min (Ud,,,2,,,)) dED
(20)
The executed visitation is simply generated via sampling a trajectory from qe- Let us denote the empirical visiation of the trajectory as 82
8 gc.
(21)
For the analysis, also need the best greedy step for the unknown (privileged) objective 2/ as
2t arg min vu, (ds) d.
(22)
deD
Let us start by considering the one step update:
Ue (d++z)
Ui de +
t (5: de))
1
=
L-smooth
ui(d:)
(de) [6, - -de) +
L
6, de 2
S
+
1. 1.+1
2(1 + 1)2
u(d(41)
boundedd S
Ui(de) + vu, (d.) T(6. - -d,) +
t+1
(1 t)z
u(dd)) + (d)) T (g. d.) 1
1
1
Vu T( ( q: + 6z)
L
++1
(1 + t)²
22

---

We will now carefully insert and subtract two set of terms depending on the real objective so that we can bound them using +
(((de)
- 22) + 1
vu(ddzz,
1 (idd) (2e de)
+ t t+1
+ 1+ 61 +
Using Un(ds) + 1 +
1
C2
+
1
Ui Un(d)) (z--d) ds) + 1 ct +
1
L
1
++
(1 t)²
Carrying on,
Uc
ug(ds) 1
(d) - dd)
1
L
t+1
1+t
ce
(1+
convexity
(,dd)) i+1 (d.ss) - (((())) + 1+t ct
1
1
I.
+
C2
(1+ t)3
1+t
1c(de+1) - ((d(d))
ue(d)) Ui({)) [+1 (Ur(dr)) - Ucd()) +
1
1
Et +
L
I+1 Ct
1
l
(1 (1+ +1z
t
(ui(a)) uu((t))
1
L
1
s 1+1
+ 1+
ct +
+
1+1 Ct 1
(1
+12
1 (ui(da)- - -u(d))
t
1
L
+ 1+t
+
C.
(1+
1+t
Now multiplying by ++ both sides, and summing on イー. 1 Uc(de) (cdt) get:
Using the shorthand pe(ds)
T-1 Σ tt (d+z22)))))
T-1 Σ + L/(++)) Σ (de + T
T-1
S
T
T
First notice that - Ei= & martingale difference sequence. The other term is the on --1 Σ Ct which appears in the main result. The sum on E(T 1 and using that (dd Z,, {z7} s u(d {z+}) for any d. This is due to set Z, decreasing over time. We report the result in asymptotic notation as function of T and log(1/5).
B log(1/5) by Lemma in Mutný et al. due to e being vT
1-1 < L To The rest is eliminated by the reccurence of the terms,
Li
D Objective reformulation and linearization
For the main objective we try to optimize over a subset of T trajectories X = {T; E x Let x be the set of sequences of inputs T = (I....., where they consist of states in the search space X. Furthermore, assume there exists, the deterministic environment, a constraint such that ah+1 E C(xA) for all h = 1, ..... H 1. Then we seek to find the set X., consisting of T trajectories (possibly repeated), such that we solve the constrained optimization problem:
X, = arg min max Varlf(=) f(=)) XEXTH 2,2'2
s.t.
Tht1 E C(xn)
VI =-=1 1,...,h 1-1
(23)
We define the objective as:
U(X) max Var 2,2'E2
- [ [(=))\x
(24)
Our goal is to show that optimization over sequences can simplified to state-action visitations as in Mutný et al. 17.. For this, we require that the objective depends additively involving terms ær, separately. We formalize this in the next result. In order to prove the result, we utilize the theory of reproducing kernel Hilbert spaces (72).
23

---

Lemma D.1 (Additivity of Best-arm Objective), Let X be collection of trajectories of length H. Assuming that f (). Assuming that has Mercer decomposition as k(z,y) y)
Σx dR (v).
f(x) Σ ou(1)) 02 Ak).
Let dx be the visitation of the states-uction in the trajectories in X. as dx - TH Σt where the Sz represent delta function supported on I. a. Then optimization of the objective Eq. 23 can be rewritten as:
TH
Ex.oen 8g.0.
1
U(dx)= = TH
max (10(2) -(=2')l 1
2.e'E2
where V(d)= =
- d(z, (z)x))(x)) T + Io (TH) is a operator V(d) : H4 - Hr. the norm is RKHS norm, and ((2)x = $1(2).
Proof. Notice that the posterior GP of any two points 2, is f())) (((((=), p(=')) :,,; where =.,,' is posterior kernel (consult Rasmussen and Williams (8] for details) defined via posterior kernel xx (z, =²) k(z,z') k(z, X) + - z').
Utilizing k(2, 2²) +(2) (()) (RKHS inner product) with the Mercer decomposition we know that (2) ((X)(()) Applying the matrix inversion lemma, the above can be written as using
Ei= Exer. (xx) I
kx(2,
-z,z2))k((z)) +xxx
Mercer +12) (12)) 4(2) ((X) (x(X) -
Lemana
D3
- o())) VV Io=)0(==) -bbz)=
((:) V--V====) 0(z)) V-1(V -V -I²))((=')
((=) -V-
Leading finally to:
T
(x(2,2))
Σ
Σ (())
+ IH
$(2").
(25)
WET
Let us calculate Var[f(z) f(2') ()). The variance does not depend on the mean. Hence,
Var [2) f(x')|X]
[ar(f(z)) - (fr(f(z')) - 2Coov((().
(x(222) 2) + 1x (x','' 2²) 2xx(2,2)
(Φ(=) -((')) Σ (x) + +001111 (0(2) +b(1))
1= mEn
TH ΣΣ #(= (())))) + toelH TH
(4(2)- - 4('')
(0(2) (((2))
[=1 tE
to arrive at:
Var [f(] (2(2) ((²'))
O² + TH Ins
Σ (x))) ZEX
( - +(=')) (26)
TH
The symbol # counts the number of occurrences. Notice that we have been able to show that the objective decomposes over state-action visitations as dx decomposes over their visitations D
24

---

Note that the objective equivalence does not imply that optimization problem in Eq. finding,
23 is equivalent
d_ = arg min max d,ED 2,2'Ez
[1((2) - + <<<))\lv(aa>>-
(27)
In other words, optimization over trajectories and optimization over de E DD is not equivalent. The latter is merely a continuous relaxation of discrete optimization problems to the space of Markov policies. It is in line with the classical relaxation approach addressed in experiment design literature with rich history. e.g., Chaloner and Verdinelli [73]. For introductory texts on the topic, consider Pukelsheim [74] for the statistical perspective and Boyd and Vandenberghe [75) for the optimization perspective. However, as Mutný et al. [17] points out, reducing the relaxed objective does reduce the objective Eq. (23) as well. In other words, by optimizing the relaxation with a larger budget of trajectories or horizons, we are able to decrease Eq. (23) as well.
For completeness, state the auxiliary lemma. We make use of the Sherman- Morrison-Woodbury Woodbury (SMW) formula, [76]:
Lemma D.2 (Sherman-Morrison-Wee
(SMW)). Let A E Ruxy and D E R2x9 then: (A DA + p²11)= p --- I- p p= - p² AAT)-'A.
(28)
Here we do the opposite, and invert an n X n. matrix instead of a x q one. to show the following.
Lemma D.3 (Matrix Inversion Lemma). Let A E Raxq then A A (AA + +p21)= = (AJ + p²1)-1 AA ATT
(29)
Note that instead of inverting n x matrix, we can invert a X Q matrix.
Proof.
A (AAT +p-1)
SMW
AT -I - - A AAA)
I - "A'A(P++ATA ATA))AA A
( -ATA))))AAA) "AAA)+++A'A) (ATA+
D.1 Objective formulation for general kernel methods
The previous discussion also allows us write the objective in terms of the general kernel matrix instead of relying on finite dimensional embeddings. The modification is very similar and relies again on Sherman-Mirrisee Mirrison- -ooodbrry lemma.
We now work backwards from (26). and first write the objective in terms of features of arbitrarily large size. Using the shorthand, /² a² R/HH, let define diagonal matrix that describes the state-action distribution D = diag( (dg : x E X}) of the size | x | and ((2) which corresponds the unique (possibly infinie-dimmensionaa) embeddings of each element in X, ordered in the same way as D.
Σ (()))) 18²) )o²) IE
= =
-I-
If we then pre-multiply by (22) and 4(2') we obtain: kx(2,z')= 4(2) T (xz)) (22) (22))
(ō'D - + +(x) - (())((2))
Finally giving:
kx
(30)
which allows us to calculate the objective for general kernel methods at the cost of an |X| x [X] inversion. Upon identifying the x' that maximize the above, we can use them in an optimization procedure. This holds irrespective of whether the state space discrete or continuous. continuous settings however, we again require parametrization of the infinite dimensional probability distribu- tion by some finite means such as claiming that D, contains some finite dimensional simplicity. This is what we do with the linear system example in Section 4.2
25

---

D.2 Linearizing the objective
To apply our method, we find ourselves having to frequently solve RL sub-problems where we try to maximize Ex.a d(z, FV, a). To approximately solve this problem in higher dimensions, becomes very important understand what the linearized functional looks like.
Remark D.4. Assume the same black-box model as in Lemma D.I and further assume that have a mixture of policies Tiix with density dxmix' such that there exists a set Xmix satisfying Taii ExeXmix Óz for some integer N. Then:
(F(dtmna where Zo, 2444444
[() Cov(S(E)), Xmix))
arg 2a Var[f(=) [(')).
Proof. To show this, we begin by defining:
Eg,d Σ ((z) + +
Ze; 2 = arg max
{(() + (1)
2,4'EZ
2, (2.)) (
In the definition above dropped the constant pre-factors since they do not influence the maximizer of the gradient as they are related by a constant multiplicative factor.
It then follows, by applying Danskin's 's Theorem that:
(z) V Le,d
vTr &
Tr { ( V
-Trr (2. E (z)) Lo,a)
(as -KK K-11)
=-Tr (1))))))))))))) Lp,
(= (0(x)
(Co] (ov
E Implementation details and Ablation study
In this section provide implementation details, and show some studies into the effeets of specific hyper-parameters. We note that the implementation code will be made public after public review.
E.1 Approximating the set of maximizers using Batch BayesOpt
We give details of the two methods used for approximating the set of potential maximizers. In particular, we first focus on Thompson Sampling [58]:
(TS)
arg max [1(x) GP(H,, 5-1 ZEX.
cont
where K is new hyper-parameter influencing the accuracy of the approximation of Z. We found that the algorithm could be too exploratory in certain scenarios. Therefore, we also propose an alternative that encourages exploitation by guiding the maximization using BayesOpt through the UCB acquisition function [71]
arg max (4(x) (,0(()) ((0(()) E K
2 (UCB) cont
zex.
i=1
where B linspace (0, 2.5, which serves as scaling for the size of ZUCBB) Both cases reduce optimization over Z to enumeration with discrete cases.
26

---

E.2
Benchmark Details
For all benchmarks, aside from the knorr pyrazole synthesis example, we use standard squared exponential kernel for the surrogate Gaussian Process:
kapf x - exp (-11e 101-2113 22 )
where dres is the prior variance of the kernel, and bnb the kernel. We fix the values of all the hyper- parameters a priori and use the same for algorithms. The hyper-parameters each benchmarks are included in T
For the knorr pyrazole synthesis example, we further set ande 0.6, ardy 0.001, 10, k2 = 874, 03 = 19200, Olsig = 5. Recall we are using finite dimensional estimate of GP such that:
M
(xx) s Wode T ode (x)
+ Σ Wrbf, (rbf(2)
(31)
in this case we set prior to the ODE weight such that Wode ~ (0.6, 0.0225). This is incorporating two key pieces of prior knowledge that (a) the product concentration should positive, and (b) we expect maximum product concentration between 0.15 and 0.45.
The number of features for each experiment, M, M min (25+d 512) where d is the problem dimensionality.
be M |X in the discrete cases and
In the case of Local Search Region BayesOpt (L.SR) we set the exploration hyper-parameter to be 0.01 in all benchmarks.
Table 2: Benchmark and hyper-parameter information. Amez represents the size of the box constraints in the traditional benchmarks. For the synchronous benchmarks and for SnAr we used a noise level of x² 0.001. For the asynchronous benchmarks, and the knorr pyrazole example x² = 0.0001. For the Ypacarai example we used o² = 0.001 and o² = 0.01 for the episodic and immediate feedback respectively.
used
Benchmark Name
A Amaz
Variance Jrbf
Lengthscale brb
Knorr pyrazole
0.001
0.1
Constrained Ypacarai
1
0.2
Branin2D
0.05
0.6
0.15
Hartmann3D
0.1
2.0
0. 0.13849
Hartmann6D
0,2
1.7
0.22
Michalewiez2D
0.05
0.35
0. 0.179485
Michalewiez3D
0.1
0.85
0.179485
Levy4D
0.1
0.6
0.14175
SnAr
0.1
0.8
0.2
E.3 Ypacarai Lake
Samaniego et al. (20) investigated the use Bayesian Optimization for monitoring the lake quality in Lake Ypacarai in Paraguay. We extend the benchmark to include additional transition constraints, as well initial and end-point constraints. These are all shown in Figure 12
E.4 Free-electron Laser
We use the simulator from Mutný experiment with varying noise levels. We use 2-dimensional variant of the simulator. We discretize the system on 10 10 grid and assume that the planning horizon Hl 100. The simulator itself is GP fit with y 0.4, hence we use this value. Then we make choice that the noise variance is proportional to the change made as (² a) s(1 + w||x - all2), where s 0.01 and w 20. Note that x |-0.5, 0.5]2 in this modeling setup. This means that local steps are indeed very difference to classical BayesOpt, which uses the worst-case variance sup s(1 + w||l = for modeling system is. We see that the absence of state modeling leads to a dramatic decrease performance as indicated by much higher inference regret in Figure 3b
[67] that optimizes quadrupole magnet orientations for our
desired. We showcase
it does not take into account the state in which the
27

---

global optima
local optima chosen trajectory
Figure 12: Lake Ypacari with the added movement constraints. We show one local optimum and one global one. The constraints of problem requiring beginning and ending the optimization dark square.
E.5 Ablation Study
E.5.1 Number of mixture components
We investigate the effect number of components used in the mixture policy when optimizing the Frank- Wolfe algorithm. We tested on the four real-world problem using N = 1, 10 and 25. In very little difference in the results, while in the Knorr pyrazole synthesis (see Figure [3) we observe much bigger difference. single component gives a much stronger performance than multiple ones we conjecture this is because the optimum on the edge of the search space, and adding more components makes the policy stochastic and less likely reach the boarder (given episodes of length ten and ten right-steps required reach the boarder).
the Ypacari example (see Figure 14)
Overall, it seems the performance of a single component is better or at worst comparable as using multiple components. This is most likely due the fact that we only follow the Markovian policies for a single time-step before recalculating, making the overall impact of mixture policies smaller. Based on this, we only present the single-component variant in the main paper.
1.0
1.0
0.037
MDP BD (1)
(10)
MDP 80 (25)
SeentFiie A fast
Average Log Regret
MDP ADP-B0 [25)
Best Arm Identifiie (%)
Average Regret
0,5
0
090 Pap
0.000
10
0.0
-12
10
20
30
40
50
80
70
80
50
100
Episode
Iteration
(a) Episodic feedback results
(b) Immediate feedback results
Figure 13: Ablation study on the number of mixture components the Knorr pyrazole synthesis benchmark
E.5.2 Size batch for approximating the set of maximizers
We explore the effect of the number of maximizers, K, in the maximization sets (T'S) and cont Overall we found the performance the algorithm be fairly robust to the size of the set in all benchmarks, with higher K generally leading little less spread in the performance.
(UCB)
28

---

100
1.0
0.033
MDP (10)
MUP-B (1)
Best Arm Identified
MDP-80
(10)
Best Arm Identified (%6)
Average Regret
(25)
MDP 80 (23)
Log Regret
0.0
-12
3
100
200
300
400
500
Episode
Iteration
Episodic feedback results
(b) Immediate feedback results
Figure 14: Ablation study on the number of mixture components on the Ypacarai benchmark
3.5
MDP-BO-TS (25)
MDP-BO-TS (50)
MDP-BO-TS (100)
2.5
Median Regret
1.5
1.0-
0.5
0.0
20
30
40
50
60
70
80
90
100
Iteration
Figure 15: Ablation study into the size of the Thompson Sampling maximization set in the asyn- chronous Hartmann3D function. We can that the performance of the algorithm is very similar for all values of = 25, 50, 100.
MDP SD UCE (:0)
3.5
HDP
0.10
HDP 80 UCE 1101
3.0-
MDP 80 UCB (100)
MDR-88 80 UCB (25)
(100)
0.08
HDP 90 uCB 1100)
HDP 80 uCB 1251
2.5
0.06
2,0
Medon
0.04
2.5
3.0
0.5
0.5
20
30
40
50
60
70
100
90
100
0
20
60
100
80
100
(a) Hartmann6D (asynch).
(b) Hartmann3D.
Levy4D.
Figure 16: Ablation study into the size of the UCB maximization set in variety of benchmarks. We can see that the performance the algorithm is very similar for all values of K 10, 25, 100.
F X3 alloacacn vs G-allocation
Our objective is motivated by hypothesis testing between different arms (options) and z'. In particular,
(d) = max Varlf(z) (()()))))))(dd))))) [xx)) 2.:22
332)
One could maximize the information the location of the optimum, has Bayesian interpretation.
This is at odds in frequentist setting, where such interpretation does not exists. Optimization of information about the maximum has been explored before, in particular
ieformationnthhorree
acquisition functions 29 77479. However, good results (in terms of regret) have been achieved by
focusing only on yet another surrogate to this, namely, the value of the maximum [8). This is chiefly
due to problem of dealing with the distribution of f(**). Defining posterior value for f(2) is easy.
29

---

0.5
C-allocation
0.4
0.3
Modian
0.2
0.1
0.0
Iteration
Figure 17: Comparison of using X'V-allocation against G-allocation as the basis for the objective. In both cases the maximization sets were created using Thompson Sampling. Overall the performances were often similar, however in few examples, such as Branin2D which we showcase here, G- allocation performed very poorly. This is consistent with what we can expect from the bandits literature.
Using, this and the worst-case perspective, an alternative way to approximate the best-arm objective, could be:
(d) max (=)l =EZ
(33)
What are we losing by not considering the differences? The original objective corresponds the X'' 2'laallocation in the bandits literature. The modified objective will, in turn, correspond to the g- allocation, which has been argued can perform arbitrarily worse as does not consider the differences, e.g. see Appendix A in Soare et al. [46]. We nonetheless implemented the algorithm with objective (33, and found the results be as expected: performance was very similar general, however in some cases not considering the differences led to much poorer performance. As an example, see Figure T7 results on the synchronous Branin2D benchmark.
G Practical Planning for Continuous MDPs
From remark D.4 it becomes clear that for decreasing covariance functions, such as the squared exponential, VF will consist of two modes around z, and 2 The sub-problem seems find sequence that maximizes the sum of gradients, therefore the optimal solution will try to reach one of the two modes as quickly as possible. For shorter time horizons, the path will reach whichever mode is closest, and for large enough horizons, the sum will be maximized by reaching the larger of the two modes.
Therefore we can approximately solve the problem by checking the value of the sub sub-problem objective in (3) for the shortest paths from #z-1 - 2. and 0t-1 =!, which are trivial to find under the constraints in (13). Note that the paths might not necessarily be optimal, as they may be improved by small perturbations, e.g., there might be small deviation that allows us to visit the smaller mode on the way to the larger mode increasing the overall value of the sum of gradients, however, they give us good and quick approximation.
H
Kernel for ODE Knorr pyrazole synthesis
The kernel is based on the following ODE model, which is well known in the chemistry literature and given in (27).
R
k1 k234495
(34)
R2 = k3y4
(35)
30

---

and then:
dy1 R
dt
dys R
dt
dy'
R1
d/
dy4 R1
R2
dt
dyss
R1
R2
Our main goal is to optimize the product concentration of the reaction, which is given by y1. We do this by sequentially querying the reaction, where we select the residence time, and the initial conditions of the ODE, in the form yo [0, A, B, 0J, where A
B.
Due to the non- -ineaiiy in Eq. (34) we are unable to GP to the process directly. Instead, we first linearize the ODE around two equilibrium points. The set of points of equilibrium are given by:
S = {1/1 = 01 1/2 = b1 y3 = 0, 11
0, y/s = C1 a1 b1, c1 E R
S {u 0gg 1/2 O, y/s
b2. Ua
0, y/s
C2 a1 b1, c1 R.
And the Jacobian of the system is:
K3
0
k1y3
1/1222
k2
k234
J
K19s
k:2l/4
O
Kiys
k1y2
k2ys
lss
k2y4
k y/
k192
k28s + ls
k284
Giving:
T
0
07
O
0
kig d2
lzzc1
J1 J|ss
O
0
ki b1
k2c1
O
0
k b1
K2C1
ka
0
ki bi
K201 + K3
1o
O
- k b2
k202
J2 = J|sz
Ki b2
K222
0
k k1b2
0
K202
K3
k1b2
K202
+ k3
Unfortunately, since the matrices are singular, we do not get theoretical results on the quality of the linearization. However, linearization is still possible, with the linear systems given by:
dy
dý
J
J2j
dt
di
We focus on the first system for now. The matrix has the following eigenvalues:
1 (1ki + cike + hs +
\1.2
Gi t +
2kg(biki 21kc2
c1k2)
A5,4,5 = 0
Note that the three eigenvalues give us the corresponding solution based on their (linearly separable) cigenvectors:
03
84 [0
0
7(t) PaUs + P4v4 + P5U6
31

---

where p; are constants. The behaviour of the ODE when this is not the case will depend on whether the remaining eigenvalues will be real or not. However, note:
ke + c + k3 + 26 2k c1 k2) 2 b2 + kg 2b,k1 = - - k3
0
and therefore all eigenvalues will always be real. Therefore can write down the solution as:
jit) Pivie + Pauge P3 P1 25
where we ignore the case repeated eigenvalues for simplicity (this is the case where bā 2 + c 2kg(61 c2k2) + kg is exactly equal to zero). We further note that the eigenvalues will be non-negative as:
h1 + c1ke + kz
v (bk) + c1 + ka)²
v DiKk + c hig + 2b k2k2 + 2b kn kig
U929 + cRg + kg 2b1 201 20122 20, k3
therefore:
S A2 0
which means the solutions will always be linear combination of exponentially decaying functions of time plus constants.
The eigenvectors have the closed form:
1 b + cake ks + by + 2bicikike + cke - - - + kg
/k3,
1.
-h1/k 1,
2
1 biki + cıki2 - k v Bi + 2b + t 2(b - k2k2)kss + hig
/hi3
=
-1/ka 1,
(ikii + t / / 2b c 2(bik) - ei k2) + kzz
/kig.
A A//k,,
A1 + 2
/ (b11) + cckk2 3k3 + k 201ck o - -1 k2) + kiz
//3
1.
1 11 + c2k2 - K3 02kk + 2b - -
/k3
1. -2/ks 1,
2
02
61k1 + c1k2 - k3 / Big + 2b + chg 2(b - k2 + hig
/k3-
d2/ li's
1,
(ykk + +ck22 + kss V 6&7 2b ci 2(b2 e₁ k2) + kzz
/lis.
Ag X2/k3,
A2 +
2 (bika +ak 3k3 Dig 2bici cgk 2(b,k) k2)k3 + hi /ba
We are optimizing over initial set of conditions yo = [0, A, B, 0, 0). solving for the specific values of the constants gives:
X2 B
P1 = A1 X2
A1
P2 X1 Av2
B
P3 = B
P4 A B
Ps = 2B
Finally, since setting A = 1
B and only optimizing the first component of ğ can
obtain it in closed form:
A
A1
(t,
B)
Betst +
A1
A2
\1
A2
X
B
Axt
chat +
(36)
\1
X2
A1
x2
32

---

0.40
numerta
0.35
0.30
0.20
0.15
0.10
0.00
30
Time
Figure 18: Comparing the numerical solution against the solutions found in equation (36).
The second ODE is very similar to the first, recall it depends has the following matrix:
0
ky b2
k222
0
J2 Js,
kibb b2
6262
k1b2
0
-
K202
k3
0
61b2
k202
+ kis
0
The resulting ODE is symmetric to the alternate linearization giving the same solution:
y(t) Proje + P2u2ee2t + P332 + P4t44 + Ps
with the only difference being the eigenvectors now are:
U/ which in turn leads to solutions of the form:
24 0
0
0 83 = 00001
(t, = A2
A1
Ac
Aet + A
\ A2
A1
A2
A2
A
Aat
e
+1
\
A2
A2
where A = 1-B. Note that we now have four different eigenvalues, which depend on the linearization points:
(1)
1
b1 + c2k2 + his + / b9k7 + c713 + 26,c - (b, C1 c2k2)
1.2
2
1
X
byki + c2k2 + h - b2At + cg13 + A3
2b - 2ka(b2kk c2k22
1,2
2
Giving solutions:
(1)
1
(1) ( B) BB
+1
A
(1)
1 (1)
2
(2)
/ (t, =
A
+1
(2)
Due to the length of the derivation, confirm that analysis is correct by comparing the numerical solution of the ODE to the exact solution we found in Figure [8] Finally, we look at interpolating between the two solutions; so given the solutions y (t)) B) and (t) B) corresponding to the consider solution of the form:
linearization with stationary point in S, and S2 respectively,
Blk,, K2 K3 a) (1 S(B)) ( (t, S(B))
(37)
33

---

where S(x) == (1 + e (= is a sigmoid function centered at B = 0.5 and where we have introduced a new syper--arameter axig. Finally, given Eq. (37) we can obtain the kernel. In particular, we want (37) to be feature we are predicting on; therefore the kernel is simply the (dot) product the features therefore:
Kode ((1, B), (V, B\)) v(t, Blk,, A2, a3 x u(t', B'] k2, K3, or)
And because we know we are simply approximating the data we can simply correct the model by adding an Gaussian Process correction; giving us the final kernel:
kjoint (tt, B), (t', B)))
Code (tt, B), (t', B3)) + arby Eebykroy((t, B), (t', BB))
where odee and Orbb are parameters we can learn, e.g. using the marginal likelihood.
34

---

NeurIPS Paper Checklist
The checklist is designed to encourage best practices for responsible machine learning research, addressing issues of reproducibility. transparency, research ethics, and societal impact. Do not remove the checklist: The papers not including the checklist will be desk rejected. The checklist should follow the references and follow the (optional) supplemental material. The checklist does NOT count towards the page limit.
Please read the checklist guidelines carefully for information on how answer these questions. For each question in the checklist:
You should answer [Yes], [No] [NA].
[NA] means either that the question Not Applicable for that particular paper or the relevant information is Not Available.
Please provide short (1-2 sentence) justification right after your answer (even for NA).
The checklist answers are an integral part of your paper submission. They are visible to the reviewers, area chairs, senior area chairs, and ethics reviewers. You will be asked to also include (after eventual revisions) with the final version of your paper, and its final version will be published with the paper.
The reviewers of your paper will be asked to use the checklist one of the factors in their evaluation.
While "[Yes] " generally preferable to "[No] ", is perfectly acceptable to answer "[No] " "rovidedd a proper justification is given (e.g., "error bars are not reported because it would be too computationally expensive" or "we were unable to find the license for the dataset we used"). "[No] " or "[NA] is not grounds for rejection. While the questions are phrased in binary way, acknowledge that the true answer often more nuanced, so please just use your best judgment and write justification elaborate. All supporting evidence can appear either in the main paper the supplemental material, provided in appendix. If you answer [Yes] to a question, in the justification please point to the section(s) where related material for the question can be found.
general, answering
IMPORTANT, please:
Delete this instruction block, but keep the section heading "NeurIPS paper checklist",
Keep the checklist subsection headings, questions/answers and guidelines below.
Do not modify the questions and only use the provided macros for your answers.
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification:
All claims in the abstract are backed up by and developed in the paper.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
It is fine to include aspirational goals as motivation long it is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations the work performed by the authors?
Answer: [Yes]
35

---

Justification: We discuss limitations in the conclusion: we are restricted to specific parameterization in the continuous state-space case.
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are 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 on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect the scope of the claims made, e.g., if the approach was only tested on a few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or speech-to-text system might not be used reliably provide closed captions for online lectures because it fails handle technical jargon.
The authors should discuss the computational efficiency the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to 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, 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 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]
Justification; All proofs are in the Appendix, where we provide full assumptions and proofs,
or cite the relevant results. Specifically see Appendix C D
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 in 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 to 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
36

---

Question: 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)?
Answer: [Yes]
Justification:
We include all the method and computational details in G
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, No answer 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 make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, provide access the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how replicate the results, access a hosted model (e.g.. in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
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) the contribution is primarily a new algorithm, the paper should make clear how reproduce that algorithm.
(b) If the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) If the contribution is new model (e.g., a large language model), then there should either way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset instructions for how to construct the dataset).
(d)
We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to describe the particular way they provide for reproducibility. In the case of closed-source models, may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers have some path reproducing or verifying the results.
5. Open access to data and code
Question: 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?
Answer: [Yes]
Justification: The code will be made public upon acceptance, and anonymized version is included for the review process.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines (https://nips.ce)
/tool////iiiie///eddsssssssooool for more details.
37

---

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 NeurIPS code and data submission guidelines (https for more details.
//nips
The authors should provide instructions 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
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?
Question:
Answer: [Yes]
Justification:
We include all the method and computational details in
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core paper to a level detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, 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: We provide regret quantiles in our experiments, reproduce on 25 random seeds and use variety of different benchmarks.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, statistical significance tests, at least for the experiments that support main claims of the paper.
The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run 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 or the standard error of the mean.
38

---

It is OK to report 1-sigma error bars, but one 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 emrors is not verified.
For asymmetric distributions, the authors should be careful to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates).
If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables 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 to reproduce the experiments?
Question:
Answer: [Yes]
Justification: See section B.4
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or 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 as estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported the paper (e.g., preliminary or failed experiments that didn't 't make into the paper).
9. Code Of Ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics
Answer: [Yes]
Justification: Upon reading the guidelines, the paper seems to conform to every point in the NeurlPS 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 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).
Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification: There are no broader impacts to discuss, outside of those already established by research into the area of design of experiments.
Guidelines:
The answer NA means that there is 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.
39

---

Examples of negative societal impacts include potential malicious unintended uses (e.g., disinformation, generating fake profiles, surveillance), faireess considerations (e.g., deploymeet 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 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 out 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 unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Question: 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)?
Answer: [NA]
Justification:
There are no risks for misuse of our work.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have 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 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 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 exisiting assests are cited and the corresponding licenses 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 or dataset.
The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (c.g., CC-BY 4.0) should be included for each asset.
40

---

For scraped data from a particular source (e.g., website). the copyright and terms of service of that souree should be provided.
If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswi thcode. com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re- eepackaged, both the original license and the license of the derived asset (if 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
Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Question:
Answer: [Yes]
Justification: The only new asset will be the code implementation which will published upon acceptance.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details the dataset/codeemmodll 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 to 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:
No crowd-sourcing human subjected were used.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing 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 in the main paper.
According to NeurIPS Code Ethics, workers involved data collection, curation, other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional Review Board (IRB) Approvals or 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 equivalent approval/review based on the requirements of your country or institution) were obtained?
Question:
Answer: [NA]
41

---

Justification:
No human subjects were used in this research.
Guidelines:
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 to adhere to 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 the institution conducting the review.
42


Reviews:
[{"review_id": "F2HmQ9LDZl", "reviewer": "Reviewer_MNRf", "summary": "The paper studies the best arm identification task in the transition-constrained setting where the MDP captures the transition. To solve the resulting planning problem defined as maximizing the acquisition propogated through the MDP, the paper studies using the RL solvers and heuritics in the discrete and continuous cases.", "questions": "1. Could the author comment on the practice of limiting the extension of MDP-BO to either TS or UCB?\n2. One minor question: in figure 3(a), is the legend \"MDP-BO\" written as \"MDP-B0\"?\n3. The objective defined in Section 2.1 and the general BO framework presented in Algorithm 1, beyond the MDP-based planning solver, are commonly seen in previous works. Li and Scarlett (2022) studied the batched version of this formulation, Zhang et al. (2023) studied the high-dimensional treatment, and Han et al. (2024) studied the formulation in game theory, among others. There seems to be potential to generalize the proposed MDP-BO framework into broader applications. It would be a bonus point to additionally discuss the corresponding potentials of MDP-BO.\n\n4. Moreover, recent advancements by Salgia et al. (2024) suggest that random exploration could be sufficient on $Z_t$. I encourage the authors to explore this if it is of interest. \n\n***References***\n\n1. Li, Zihan, and Jonathan Scarlett. \"Gaussian process bandit optimization with few batches.\" In International Conference on Artificial Intelligence and Statistics, pp. 92-107. PMLR, 2022.\n2. Zhang, Fengxue, 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 on Machine Learning, pp. 41579-41595. PMLR, 2023.\n3. Han, Minbiao, Fengxue Zhang, and Yuxin Chen. \"No-Regret Learning of Nash Equilibrium for Black-Box Games via Gaussian Processes.\" In The 40th Conference on Uncertainty in Artificial Intelligence.\n4. Salgia, Sudeep, Sattar Vakili, and Qing Zhao. \"Random Exploration in Bayesian Optimization: Order-Optimal Regret and Computational Efficiency.\" In Forty-first International Conference on Machine Learning.", "limitations": "Discussed above.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "1. Valid scenario, well-motivated and well-illustrated algorithm.\n2. Comprehensive empirical evidence substantiating the effectiveness of the proposed algorithm, along with verification of parameter choices.", "weaknesses": "1. Although employing an RL solver, as admitted by the authors, the objective is non-convex and differs from conventional optimization in RL. This demands heuristics, which are not sufficiently discussed.\n\n2. In the experiments, the GP hyperparameters are fixed, which is questionable due to potential misspecification and manually introduced bias. A typical practice might include either a full Bayesian treatment of the GP hyperparameters or kernel learning. This is particularly relevant for the ODE kernel due to its higher complexity.\n\n3. The results in Figure 4 appear to lack statistical significance. Additionally, there is no clear preference among the variants of MDP-BO, resonating with the challenge in conventional BO.\n\n4. There is a series of works on applying BO to graphical structures. Although the MDP formulation is different, adding them to the discussion would enhance completeness as some design choices are conceptually related.\n\n***References***\n\n1. Kusakawa, Shunya, Shion Takeno, Yu Inatsu, Kentaro Kutsukake, Shogo Iwazaki, Takashi Nakano, Toru Ujihara, Masayuki Karasuyama, and Ichiro Takeuchi. \"Bayesian optimization for cascade-type multistage processes.\" Neural Computation 34, no. 12 (2022): 2408-243.\n2. Aglietti, Virginia, Xiaoyu Lu, Andrei Paleyes, and Javier González. \"Causal Bayesian optimization.\" In International Conference on Artificial Intelligence and Statistics, pp. 3155-3164. PMLR, 2020.\n3. Astudillo, Raul, and Peter Frazier. \"Bayesian optimization of function networks.\" Advances in Neural Information Processing Systems 34 (2021): 14463-14475.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "**W: The objective is non-convex and differs from conventional optimization in RL. This demands heuristics, which are not sufficiently discussed.**\n\nThis is true for the continuous case and common to most setups in control theory. However, special to our setting, for all discrete environments in the experiments we are able to solve the sub-problems provably with convex solvers, see Frank-Wolfe updates.\n\n**W: In the experiments, the GP hyperparameters are fixed...**\n\nPlease see the general rebuttal.\n\n**W: The results in Figure 4 appear to lack statistical significance. Additionally, there is no clear preference among the variants of MDP-BO, resonating with the challenge in conventional BO.**\n\nPlease see the general rebuttal.\n\n**W: There is a series of works on applying BO to graphical structures...**\n\nCausal BO uses graphs that relate different variables in the optimization process according to a Bayesian causal model. Our graphs relate transitions between states. This is different in spirit, but we will add to the discussion.\n\n**Q: Could the author comment on the practice of limiting the extension of MDP-BO to either TS or UCB?**\n\nWe obtained good results with both, so we did not try to use any other batch algorithms for estimating the maximizer candidate set. However, there is flexibility in the algorithm to use different estimation methods and could be the focus of future work.\n\n**Q: One minor question: in figure 3(a), is the legend \"MDP-BO\" written as \"MDP-B0\"?**\n\nCorrect, this is an error, thanks for catching it!\n\n**Q: The objective defined in Section 2.1 and the general BO framework presented in Algorithm 1, beyond the MDP-based planning solver, are commonly seen in previous works. Li and Scarlett (2022) studied the batched version of this formulation, Zhang et al. (2023) studied the high-dimensional treatment, and Han et al. (2024) studied the formulation in game theory, among others. There seems to be potential to generalize the proposed MDP-BO framework into broader applications. It would be a bonus point to additionally discuss the corresponding potentials of MDP-BO.**\n\nThanks for all the references, they are indeed relevant. We will include a discussion of them in the paper.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "gQv15kolfO", "reviewer": "Reviewer_J8j4", "summary": "&nbsp;\n\nThe authors address the problem of transition-constrained Bayesian optimization, modeling the transition constraints as a Markov decision process. They take a novel approach in deriving their utility function for this problem setting based on maximum identification and bounding the probability of choosing an incorrect maximizer. While the method is well-described and has the potential to yield theoretical insight into the structure of transition-constrained Bayesian optimization, the empirical performance does not seem hugely compelling compared to existing methods such as SnAKe and LSR. As such, my recommendation is borderline with potential to increase my score if the authors can address the main criticism of the advantages of the method compared to competitor approaches.\n\n&nbsp;", "questions": "&nbsp;\n\n1. In Equation 3, do the authors have any intuition as to why k(z, x_f*) / (f(z) - f(x_f*))^2 is related to the probability of returning an wrong maximizer? Beyond the derivation presented in Section C.1 of the appendix.\n\n2. In Equation 4, the authors state their objective as a maximization problem. However, in the text, the authors refer to minimizing the uncertainty among all pairs in \\mathcal{Z}. Could the authors explain this?\n\n3. In Section D.1 of the appendix, the authors give the formulation of the objective for general kernel matrices. The authors state that the objective may be computed by inverting a |\\Chi| x |\\Chi| matrix. How restrictive is this for the scale of problems that can be addressed with this formulation of the objective?\n\n4. In Section 5, what was the motivation for the decision to fix the GP hyperparameters?\n\n5. For Figures 2 and 3 why are there no error bars for the average regret curves?\n\n&nbsp;", "limitations": "&nbsp;\n\n1. The main limitation as described above would appear to be the lack of compelling evidence for problem settings under which the authors' algorithm outperforms other transition-constrained Bayesian optimization methods such as LSR and SnAKe.\n\n\n&nbsp;", "rating": 7, "confidence": 3, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "&nbsp;\n\nThe authors derive a novel transition-constrained Bayesian optimization algorithm, leveraging the literature on maximum identification via hypothesis testing to derive their utility function.\n\n&nbsp;", "weaknesses": "&nbsp;\n\n__MAJOR POINTS__\n\n&nbsp;\n\n1. The main weakness of the work would appear to be the lack of compelling evidence for why the authors' method should be chosen over competitor methods such as LSR or SnAKe on practical problems. I think the paper would benefit substantially if the authors could highlight problem settings in which they expect their algorithm to perform better in relation to at least one practical performance metric.\n\n2. The code would benefit from a README describing how to reproduce the results of the paper.\n\n&nbsp;\n\n__MINOR POINTS__\n\n&nbsp;\n\n1. It would be great if the references appeared in numbered order.\n\n2. There is some missing capitalization in the references e.g. \"bayesian\" in place of \"Bayesian\".\n\n3. Line 90, typo, \"we seek to identify\".\n\n4. Line 566, typo, \"due to selecting\".\n\n5. Line 566/567, missing full stop at the end of the equation.\n\n6. On line 567 epsilon is homoscedastic but in the main paper the authors state that epsilon can be heteroscedastic?\n\n7. Line 584, typo, \"hypotheses\".\n\n8. Line 680, typo, extraneous \"be\".\n\n9. Line 680, in the notation should there be some kind of indication that X corresponds to T trajectories e.g. an indication of the dimensions of the space X corresponds to? The same on line 690.\n\n10. In Equation 4, the conditioning bar is difficult to read, could the authors introduce some more space to the left and right of the bar?\n\n11. In the Subsection on \"Utility with Embeddings\", the authors may wish to specify \"Kernel Embeddings\" so as to disambiguate from embeddings over the input space X.\n\n12. In the proof of statement 25, it might be worth giving the acronym for SMW, namely Sherman-Morrison-Woodbury, somewhere in the text.\n\n13. On line 136, have the authors defined S? This is the state space presumably.\n\n14. In the related work on look-ahead BayesOpt, it may also be worth mentioning [1].\n\n15. Line 179, typo, \"a mix of the previous two\".\n\n16. In Section 2.2, the variable d could be explained in more detail.\n\n17. When referencing Thompson sampling, it would be worth citing the original paper [2].\n\n18. In Section E.2, the expression for the RBF kernel implicitly assumes x is 1-dimensional. The authors could replace the squared distance with the squared 2-norm to generalize the expression.\n\n19. In Algorithm 1, would it be more descriptive to refer to Z_0 as the initial set of maximizer candidates rather than the initial set of maximizers?\n\n20. It would be worth citing [3] (as per the discussion the historical section by Garnett [4]) for the Expected Improvement acquisition function given that it is used.\n\n21. It would be great if the authors could provide the number of random trials over which the error bars were computed in Figure 4.\n\n22. Line 622, \"scales them properly\". What is \"them\"?\n\n23. Line 670, typo, \"using\".\n\n&nbsp;\n\n__REFERENCES__\n\n&nbsp;\n\n[1] Roman Marchant, Fabio Ramos, and Scott Sanner. Sequential Bayesian optimisation for spatial-temporal monitoring. In Proceedings of the International Conference on Uncertainty in Artificial Intelligence, 2014.\n\n[2] W. R. Thompson. On the likelihood that one unknown probability exceeds another in view of the evidence of two samples. Biometrika, 25(3-4):285–294, 1933.\n\n[3] VR Saltenis (1971). One Method of Multiextremum Optimization. Avtomatika i Vychislitel’naya Tekhnika (Automatic Control and Computer Sciences) 5(3):33–38.\n\n[4] Garnett, R., 2023. Bayesian optimization. Cambridge University Press.\n\n&nbsp;", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "**W: The main weakness of the work would appear to be the lack of compelling evidence for why the authors' method should be chosen over competitor methods such as LSR or SnAKe on practical problems...**\n\nPlease see the general rebuttal.\n\n**W: The code would benefit from a README describing how to reproduce the results of the paper.**\n\nWe agree and we will add this after peer review. We did have a README but hastily deleted it to preserve anonymity!\n\n**W: On line 567 epsilon is homoscedastic but in the main paper the authors state that epsilon can be heteroscedastic?**\n\nFor simplicity, the derivation is done for the homoskedastic setting, however the same derivation holds in the heteroskedastic setting by simply changing $\\mathbf{V} = \\sum_{i=1}^T \\frac{1}{\\sigma_i^2}\\Phi(x_i)\\Phi(x_i)^\\top$ instead of for heteroskedastic least squares $\\mathbf{V} = \\sum_{i=1}^T \\frac{1}{\\sigma^2}\\Phi(x_i)\\Phi(x_i)^\\top =: \\mathbf{X} \\mathbf{\\Sigma}\\mathbf{X}$. We can add this comment to the discussion.\n\n**W: Other minor points**\n\nThanks for pointing out minor issues which we will address.\n\n**Q: In Equation 3, do the authors have any intuition...**\n\nIt allows an interpretation as noise to signal ratio. In the numerator, we have uncertainty (covariance) between $f(z)$ and $f(x_f^*)$, while in the denominator we have the level of signal that is the gap between $f(z)$ and $f(x_f^*)$. The smaller the gap, the harder it is to know which of the two is maximal, increasing the probability of error. On the other hand, if the two states $x_f^*$ and $z$ are similar to each other, we do not need to put additional effort as they are very correlated, hence balancing these two is the right metric. We will add this intuition to the paper to improve clarity.\n\n**Q: In Equation 4, the authors state their objective as a maximization problem. However, in the text, the authors refer to minimizing the uncertainty among all pairs in $\\mathcal{Z}$. Could the authors explain this?**\n\nThe ultimate goal is to maximize an unknown function $f$. For exposition, as is common in BayesOpt, we talked about utility or acquisition function maximization, but then deriving the utility we find we need to minimize the uncertainty among all pairs. This lends itself to explain as minimization better. We could equally talk about maximization of utility $=$ 1/loss.\n\n**Q: In Section D.1 of the appendix, the authors give the formulation of the objective for general kernel matrices. The authors state that the objective may be computed by inverting a $|\\mathcal{X}| \\times |\\mathcal{X}|$ matrix. How restrictive is this for the scale of problems that can be addressed with this formulation of the objective?**\n\nIndeed, this is a limitation seen also with some classical BO approaches such as Thompson sampling which need inversion of $|\\mathcal{X}| \\times |\\mathcal{X}|$ matrix -- same complexity as our method. \n\nCommon remedies are some finite dimensional approximations of the kernels, leading to $m \\times m$ inversions or other reformulations improving the efficiency of matrix inversion that we cite in our work. Akin to Thompson sampling, even combinations of different techniques such as feature approximation and different sampling, namely Matheron rule for sampling, could be a very practical remedy (Wilson et al. 2020). This issue is a fundamental difficulty associated with reasoning (planning) in non-parametric function spaces, not necessarily something specific to our method. Also notice that in continuous domains the reparametrization of $d$ eliminates this inversion but leads to potentially non-convex program in the reparametrization.\n\n**Q: In Section 5, what was the motivation for the decision to fix the GP hyperparameters?**\n\nPlease see the general rebuttal.\n\n**Q: For Figures 2 and 3 why are there no error bars for the average regret curves?**\n\nWe did not include them to avoid cluttering the graphics. The percentage of runs we identify the maximizer allows us to see the robustness of each method. Nonetheless, we are happy to include complete regret plots in the appendix.\n\n**L: The main limitation as described above would appear to be the lack of compelling evidence for problem settings under which the authors' algorithm outperforms other transition-constrained Bayesian optimization methods such as LSR and SnAKe.**\n\nPlease see the general rebuttal.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "W5AFoxQYxK", "reviewer": "Reviewer_mJAD", "summary": "The paper introduces a new Bayesian optimization problem that finds the optimal policy to optimize a black-box function subject to transition constraints on the query. The method works for both discrete and continuous Markov chains. The paper empirically demonstrates several practical applications in physical systems, electron laser calibration, and chemical reactor optimization.", "questions": "1. Section 2 is the background but there are no citations in section 2.1. Is this a new solution proposed by the paper, or is it based on existing work?\n\n2. Line 575: Please explain why the proprotionality holds. Monotonicity does not imply proportionality (same ratio).\n\n3. Line 587: Please explain the two inequalities.\n\n4. Would the authors please clarify equation (7)? What does it mean by the expectation of a state $X_t$? What is $\\pi$ in $d_{\\pi}$? Z_t and GP posterior may change after every time step, so I have a hard time understand why U remains the same and the expectation is only taken over X_t in equation (7). Furthermore, as the utility depends on the time step (the GP posterior and Z_t), two visitations to the same state at different time steps have two different utilities. Then, how do we combine all visitations into $d_\\pi$ (lines 143-144).\n\n5. The pseudocode in Algorithm 1 is unclear. What is d_pi (there is no initialization of it)? Why does the update of $\\hat{d}\\_{t,h+1}(x)$ match equation (11)? Why are GP_{t+1} and Z_{t+1} only updated with X_{t,H} and Y_{t,H} (the observations from episode t, not all episodes from 1 to t)? What are X_{t,H} and Y_{t,H} (the time step is only until H-1).\n\n6. The policy depends on the time step, so pi_{t,0} is only deployed at the first time step. How is the state visitation frequency of such a policy estimated well if it is only deployed at the 1st time step of each episode?\n   \n7. Would the authors please clarify equation (12) using the correct subscript of d (d_{t,h})?\n   \n8. Some minor remarks: In equation (11), should j be from 1 to t-1? Line 122: y(xh) should be y(xh,ah).", "limitations": "There are no potential negative societal impacts.", "rating": 3, "confidence": 2, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "The problem formulation is new, and the paper effectively motivates it using several practical applications. The solution is based on principled methods, with several experiments supporting the practicality of the problem.", "weaknesses": "The paper may be improved by clarifying several questions below.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "**Q: Section 2 is the background but there are no citations in section 2.1. Is this a new solution proposed by the paper, or is it based on existing work?**\n\nPlease see the general rebuttal.\n\n**Q: Line 575: Please explain why the proportionality holds. Monotonicity does not imply proportionality (same ratio).**\n\nThis is a typo. With the first $A \\propto B$ we refer to there existing a finite constant such that $A(x) = C B(x)$. Our goal is to find $\\arg\\min_x A(x)$, so we can equivalently minimize $\\min \\log B(x)$ to identify the minimizer of $A(x)$ as well, and this is what the second $\\propto$ refers to, but it should be replaced with $\\stackrel{!}{=}$ and explained in the text for more clarity. We apologize for the confusion. In addition, later in line 578, there is a typo; instead of $\\propto$, it should read $\\approx$, see the preceding line for an explanation.\n\n**Q: Line 587: Please explain the two inequalities.**\n\nThe first inequality is missing a square (a typo, see 574 for the correct definition). The step is incorrectly explained and in fact holds only for large $T \\gg 0$, or when $\\mathbf{V}_T \\gg \\mathbf{I}$, where then $(\\mathbf{V}_T + \\mathbf{I})^{-1}\\mathbf{V} \\approx \\mathbf{I}$.\n\nThe second inequality follows by $\\mathbf{V}_T = \\mathbf{X}^\\top \\mathbf{X}$, push-forward identity, $\\mathbf{V}_T \\preceq \\mathbf{V}_T + \\mathbf{I}$ and monotonicity of the matrix inverse under Löwener ordering. The two are then substituted back to the original objective. Notice that the objective is as such only when $T$ is large (asymptotic regime).\n\n**Q: Would the authors please clarify equation (7)? What does it mean by the expectation of a state $X_t$? What is $\\pi$ in $d_{\\pi}$? $Z_t$ and GP posterior may change after every time step, so I have a hard time understanding why $U$ remains the same and the expectation is only taken over $X_t$ in equation (7). Furthermore, as the utility depends on the time step (the GP posterior and $Z_t$), two visitations to the same state at different time steps have two different utilities. Then, how do we combine all visitations into $d_\\pi$ (lines 143-144).**\n\nIn the discussion around Eq. (7), $X_t$ refers to the set of trajectories (see L132). We analyze and derive the algorithm in the *episodic* setting (trajectory feedback), where we execute a trajectory in full before updating the posterior (GP and $Z_t$). A trajectory is a random set, and therefore we must define a utility with respect to the expected trajectory when using a policy $\\pi$. Each policy $\\pi$ induces a state-action distribution, which tells us how likely we are to visit a state in a realization of a trajectory, which we denote $d_\\pi$. To our benefit, as we show in the derivations, the objective depends only on $d_\\pi$.\n\nFor trajectory feedback, $Z_t$ does not change at every time-step, only after the whole trajectory is executed. When planning, we keep $Z_t$ fixed. Later we introduce *instant feedback*, where GP and policy are recalculated at each time step. However, we do not provide a formal derivation of optimality for this setup. The utility is still defined properly where we optimize the remaining part of the trajectory and by extension $\\pi$ or $d_\\pi$ or the time-steps from $h$ to $H$. This draws a parallel to the MPC literature of using receding horizon re-planning.\n\n**Q: The pseudocode in Algorithm 1 is unclear. What is $d_\\pi$ (there is no initialization of it)? Why does the update of $\\hat{d}\\_{t,h+1}(x)$ match equation (11)? Why are $GP_{t+1}$ and $Z_{t+1}$ only updated with $X_{t,H}$ and $Y_{t,H}$ (the observations from episode $t$, not all episodes from 1 to $t$)? What are $X_{t,H}$ and $Y_{t,H}$ (the time step is only until $H-1$).**\n\nNote $d_\\pi$ does not require initialization as it is the variable of the utility; it denotes the state-action distribution of a policy. In the second line of the algorithm you can see we take an argmax over this variable.\n\nThroughout the algorithm, we want to keep track of the visited states. We do this via the empirical state-action distribution, defined by the states we have visited so far. After visiting new state(s), we want to add them to the empirical state-action distribution which we can do by following the normalization procedure defined in eq. (11). For clarity, the empirical state-action distribution $\\hat{d}\\_{t,h}$ denotes _past_ states, while the variable $d_\\pi$ is representing _future_ states we must choose.\n\nThe GP and $Z$ are updated using observations from all previous trajectories. In the paper, we use lower-case variables to denote individual observations, and upper-case variables to denote all trajectories up to the time-step specified by the index.\n\n**Q: The policy depends on the time step, so $\\pi_{t,0}$ is only deployed at the first time step. How is the state visitation frequency of such a policy estimated well if it is only deployed at the 1st time step of each episode?**\n\nThe state-action frequency of a given policy can be calculated exactly; there is no need to deploy the policy itself. We assume the transition dynamics are known.\n\n**Q: Would the authors please clarify equation (12) using the correct subscript of $d (d_{t,h})$?**\n\nNotice that $\\hat{d}_{t,h}$ refers to the empirical visitation. In other words, the frequency of visited states in trajectory $t$ at time horizon $h$. This visitation changes as we execute the policies for a longer time and more trajectories—this is updated depending on the trajectory we have followed up to time-step $t, h$. The empirical visitation $\\hat{d}$ and also the variable $d$ in optimization Eq. (12) is $d:\\mathcal{S} \\times \\mathcal{A} \\times H \\rightarrow [0,1]$.\n\n**Q: Some minor remarks: In equation (11), should $j$ be from 1 to $t-1$? Line 122: $y(xh)$ should be $y(xh,ah)$.**\n\nYes, thanks for pointing out the mistakes.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "OS7RCI79JN", "reviewer": "Reviewer_2XmZ", "summary": "This paper explores extending Bayesian optimization to incorporate transition constraints through state dynamics. The approach treats the problem as an optimal control problem. For planning purposes, known state dynamics model the transition constraints, while the unknown objective is represented using a Gaussian process. To address this, the paper introduces a tractable acquisition function—an upper bound on the maximum probability of selecting the wrong optimizer. This acquisition function serves as a cost metric in a planning scenario, which is solved using the Frank-Wolfe algorithm. The paper’s effectiveness is evaluated across various benchmark tasks.", "questions": "- Is the utility function presented in the paper a contribution? Could the authors elaborate a bit on this?\n- I am a bit confused about the notation in line 132. Could the authors please define $\\tau_i$?\n - In line 136 and 199 is $S$ a switch in notation?\n-\tCould the authors elaborate a bit on line 209. What exactly is sampled using Thompson Sampling?\n-\tIn line 259, why is this specific contstraint of 0.5 chosen?\n-\tI am a bit confused about line 263. I do not think MPC is nescessarly convex. There is a lot of literature on non-linear MPC. Could the authors elaborate on this?", "limitations": "-", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper presents a framework for optimizing challenging-to-evaluate black-box functions with transition constraints, which could have practical applications. Its overall structure is well-organized.\n\nThe proposed solution strategy for solving the planning problem using the Frank-Wolfe algorithm is novel and represents a valuable contribution to the research community. Additionally, considering both discrete and continuous state problems is commendable.\n\nWhile the experiments are fair, including a theoretical synthetic example would have further strengthened the paper.", "weaknesses": "General Presentation: \n\nGenerally, I believe the paper’s presentation could be enhanced. Specifically, I am not very satisfied with how the generative model is derived and presented. In my opinion, the problem could be expressed as a stochastic optimal control with an unknown cost function. However, the paper takes a different approach, starting from Bayesian optimization, which typically does not incorporate control or dynamics. This choice leads to confusion in notation and extends the Bayesian optimization framework. For instance, the notation between sections 2.1 and 2.2 changes due to the introduction of the action variable a. Additionally, in Section 2.1, the authors discuss a set of arms, drawing from the bandit literature.\nWhat also confuses me is that the actions in the model are never formally introduced. Especially in Section 2.3, a formal definition of actions would have been helpful.\n\nBackground Section: \n\nIn the background section, it is not very obvious what the paper’s contribution is versus what constitutes background information. The citations are not very transparent.\n\nRelated Work (Section 2.4): \n\nI believe one of the most important related concepts, Bayesian Reinforcement Learning or Dual Control, is not discussed at all.\n\nAcquisition Function Analysis: \n\nAdditionally, I would have appreciated more analysis on why the specific acquisition function was chosen over others. How would the algorithm perform with a different acquisition function?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "**W: Generally, I believe the paper’s presentation could be enhanced...**\n\nIndeed, the reviewer is right we could start by presenting the problem as a control problem with an unknown objective, and then derive a surrogate for the unknown objective that is refined over episodes. Our main audience is practitioners of Bayesian optimization; hence, we chose the route in the paper. We believe this comment reflects preferences in exposition, not on the correctness or soundness of the work.\n\nThe actions are *formally* defined in Sec. 2.2. We could have introduced them already in Sec. 2.1, but it would clutter the notation to no benefit in exposition. We will make sure to improve clarity in Sec. 2.1.\n\n**W: In the background section, it is not very obvious what the paper’s contribution is versus what constitutes background information. The citations are not very transparent.**\n\nPlease see the general rebuttal.\n\n**W: I believe one of the most important related concepts, Bayesian Reinforcement Learning or Dual Control, is not discussed at all.**\n\nIndeed, these are related but also crucially different. For example, Bayesian RL uses a reward function that is linear in the state-action visitation, while our utility is non-linear but convex. Dual control goes more in the direction of an unknown transition system, which we assume to be known here but indeed constitutes a very reasonable further work. We can add a small discussion of both concepts.\n\n**W: Additionally, I would have appreciated more analysis on why the specific acquisition function was chosen over others. How would the algorithm perform with a different acquisition function?**\n\nThe objective was motivated by a goal to minimize the probability of a mistake. It was then refined to be computationally tractable with only one inequality in the derivation that makes it not tight (see general rebuttal). Our framework with the acquisition function as we derive it has the advantage that the induced state-action distribution of a trajectory allows for adaptive resampling and is theoretically grounded.\n\nIf we had a different goal to maximizer identification, a different acquisition function would perhaps be more appropriate and could be used within our framework provided that it depends only on the state-action visitations akin to described in Mutny et. al. (2023) in Lemma 1. We leave this to future work.\n\n**Q: Is the utility function presented in the paper a contribution? Could the authors elaborate a bit on this?**\n\nPlease see the general rebuttal.\n\n**Q: I am a bit confused about the notation in line 132. Could the authors please define $\\tau_i$?**\n\n$\\tau_i$ is the deployed trajectory at the $i$th episode. We will clarify this in the paper.\n\n**Q: In line 136 and 199 is $S$ a switch in notation?**\n\nThis was a typo, thanks for pointing it out.\n\n**Q: Could the authors elaborate a bit on line 209. What exactly is sampled using Thompson Sampling?**\n\nWe use Thompson Sampling to sample the elements of the candidate set of maximizers. In other words, it is used as a finite approximation of the uncertainty quantification provided by the continuous set $\\mathcal{Z}_t$. That is, we create $K$ samples of the GP and use the optimum of each sample to create a set of $K$ potential maximizers.\n\n**Q: In line 259, why is this specific constraint of 0.5 chosen?**\n\nThis is an arbitrary choice related to the code implementation, as we work in the search space $[-0.5, 0.5]$ (and normalize functions and data to those bounds). In reality, any bounds could be chosen, in which case the constraint becomes $x \\in [a, b]$. Thanks for pointing this out; we will clarify this in the paper.\n\n**Q: I am a bit confused about line 263. I do not think MPC is necessarily convex. There is a lot of literature on non-linear MPC. Could the authors elaborate on this?**\n\nOne thing is non-linear, another is non-convex. In fact, there is a lot of literature on both. Most MPC is studied with quadratic costs (non-linear), and with linear dynamics, then it is convex in the action variables. This is what is taught in graduate MPC courses. Of course, there is also non-convex non-linear MPC (common with mixed integer, for example), but this is more exotic and needs to be heuristically solved. We fall into this category as our cost is the *utility* which is non-convex in the action variables, but it can be heuristically solved.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:
We would like to thank all reviewers for the time spent reading and evaluating our paper. We are happy to have the strengths of our work recognized, and we are also pleased to have received so much constructive feedback. Here we address points common to multiple reviewers.

**Comparison with SnAKe and LSR:** Both SnAKe and LSR are algorithms specializing in a very specific type of transition constraint, when we have a restriction on how far we can move between experiments. The goal of our work is to deal with transition constraints in a more general and formalized setting. In fact, all the real-world experiments we consider are examples that SnAKe and LSR cannot handle. When considering local search constraints, it is not surprising that SnAKe and LSR perform strongly since they are specialist heuristics in this setting; however, we find it very encouraging that MDP-BO is able to match these SOTA algorithms.

**Section 2.1:** We included this section as background because similar objectives have been derived and used in the literature. Indeed, Fiez et al. [1] used the same objective for transductive linear bandits and there are other examples as reviewer MNRf pointed out. However, the motivation and derivation in terms of a Bayesian decision rule do not appear elsewhere according to our best knowledge. Furthermore, recognizing that we can relax the objective and optimize it in the space of policies is also a contribution. In the final version of the paper, we will carefully restate our contributions to make this clearer.

We had an omission of an assumption in the derivation of the objective. The first inequality in L587 only holds for large $T \rightarrow \infty$. This points to a similarity to Fiez et al. [1] being optimal in the asymptotic regime only as we explained in the Appendix. A better utility would require calculating the expectation of $ \mathrm{E}_{f \sim \text{GP}}[1/a_z^2]$ (see L574) which is unfortunately intractable. Were it possible to evaluate it would most certainly lead to a better utility. However, the gains are only to be seen for a very small regime of $T$ and quickly diminish. Based on this, the Eq. (3) in the paper should read:

$$
\min_{X_{\text{new}}}E_{f}\left[\sup_{z \in Z \setminus \{x_f^*\} }\log P(\mu_T(z) - \mu_T(x_f^*) \geq 0 | f)\right] \stackrel{\approx}{\leq}  \min_{X_{\text{new}}} E_f\left[ \sup_{z \in Z \setminus \{x_f^*\}} \frac{k_{X_{t} \cup X_{new}}(z,x_f^*)}{(f(z) - f(x_f^*))^2}\right]
$$

where we added the $\stackrel{\approx}{\leq}$ to denote that this is an approximate bound that will hold for large $T$, i.e., asymptotically. Note that the bound is on the log probability (this was a typo previously). The corresponding utility, the rest of the paper, and all empirical results are unchanged.

**Fixing the hyper-parameters for the experiments:** Good performance of Bayesian optimization is strongly dependent on such hyper-parameters, therefore there is a need to use previous data or expert knowledge to either fix the hyper-parameters (we use empirical Bayes, by giving all algorithms the same hyper-parameters learnt via marginal-likelihood on a small random hold-out set) or put a good prior while estimating them as we go (full Bayes). We fix them to compare without GP hyper-parameter choice being a co-founding effect since we are interested in quantifying the effect of the algorithm improvement. Note that all competing algorithms benefit or suffer from this. Investigating the best way of learning the GP's parameters would be interesting future work.

**References**

[1] Fiez, Tanner, et al. "Sequential experimental design for transductive linear bandits." Advances in Neural Information Processing Systems 32 (2019).

Meta-review: This paper presents an approach to incorporating transition constraints through Markovian state dynamics into a Bayesian Optimization framework. Reviewers generally agree that the problem formulation is novel, well-motivated, and well-illustrated by the empirical results. Reviewers also found the solution approach to be principled and novel. However, as it (necessarily, due to the problem complexity) relies on a number of approximations and relaxations, I found keeping track of the various components somewhat challenging and therefore ask the authors to consider ways of helping the reader with this. In addition, reviewers raised some valid points about discussing additional related work that the authors should take into account. Finally, while there were some questions about the technical details of the algorithm that were not fully resolved in the rebuttal discussion, I believe that those do not invalidate the results and thus the value of the paper.
Decision: Accept (poster)