Transition Constrained Bayesian Optimization via Markov Decision Processes

Authors: jose pablo folch
calvin  tsay
robert m lee
andreas  krause eth
zurich  zurich
mark  van der wilk
ruth  misener
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 16:05:49

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.

Main paper: 
 introduction: 
 Many areas in the natural sciences and engineering deal with optimizing expensive black-box functions. Bayesian optimization (BayesOpt) [1] [2] [3] , a method to optimize these problems using a probabilistic surrogate, has been successfully applied to a myriad of examples, e.g. hyper-parameter selection [4] , robotics [5] , battery design [6] , laboratory equipment tuning [7] , and drug discovery [8] . However, state-of-the-art algorithms are often ill-suited when physical sciences interact with potentially dynamic systems [9] . 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-decision 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 [10] [11] [12] , environmental monitoring [13] [14] [15] [16] [17] , lake surveillance with drones [18] [19] [20] , energy systems [21] , vapor compression systems [22] , electron-laser tuning [23] and seabed identification [24] . For example, Figure 1 depicts an application in environmental monitoring where autonomous sensing vehicles must avoid obstacles (similar to Hitz et al. [18] ). Our main focus application are transient flow reactors [25] [26] [27] . Such reactors allow efficient data collection by obtaining semi-continuous time-series data rather than a 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 quasi-steady-state 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 to identify the optimal configuration of a physical system governed by a black box function f , namely, x ⋆ = arg max x∈X f (x). 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) + ϵ(x), where ϵ has a known Gaussian likelihood, which is possibly heteroscedastic. We assume that f can be modeled probabilistically using a Gaussian process prior that we introduce later. Importantly, the order of the evaluations is dictated by known, potentially stochastic, dynamics modeled by a Markov chain that limits our choices of x ∈ 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:EQUATIONdepending on all the past data at iteration t, X t . Eq. ( 1) 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 a utility as,EQUATIONwhere U is the greedy utility of the policy π and X t 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 it to admit policies depending on the full optimization history.Contributions. We present a BayesOpt framework that tractably plans over the complete experimentation horizon and respects Markov transition constraints, building on active exploration in Markov chains [17] . Our key contributions include:• We identify a novel utility function for maximum identification as a function of policies, and greedily optimize it. The optimization is tractable, and does not scale exponentially in the policy horizon. In many cases, the problem is convex in the natural representation.• We provide exact solutions to 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 nonconvex objective. Interestingly, in both cases, the resulting policies are history-dependent (non-Markovian).• We analyze the scheme theoretically and empirically demonstrate its practicality on problems in physical systems, such as electron laser calibration and chemical reactor optimization.The previous discussion also allows us to 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-Mirrison-Woodbury lemma.We now work backwards from (26) , and first write the objective in terms of features of arbitrarily large size. Using the shorthand, σ2 = σ 2 /T H, let us define a diagonal matrix that describes the state-action distribution D = diag({d x : x ∈ X }) of the size |X |×|X |, and Φ(X ) which corresponds to the unique (possibly infinite-dimensional) embeddings of each element in X , ordered in the same way as D.σ2 x∈X d(x)Φ(x)Φ(x) + Iσ 2 -1 = σ2 Φ(X ) T DΦ(X ) + Iσ 2 -1 = I -Φ(X ) T (σ 2 D -1 + Φ(X )Φ(X ) ⊤ ) -1 Φ(X )If we then pre-multiply by Φ(z) ⊤ and Φ(z ′ ) we obtain:EQUATION) which allows us to calculate the objective for general kernel methods at the cost of an |X | × |X | inversion. Upon identifying the z, z ′ that maximize the above, we can use them in an optimization procedure. This holds irrespective of whether the state space is discrete or continuous. In continuous settings however, we again require a parametrization of the infinite dimensional probability distribution 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. 
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 a visual overview of how we selected the individual components for tractability 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 a GP, e.g., for inputs (x 1 , .., x p ), the values {f (x j )} p j=1 , are normally distributed. We adopt a Bayesian approach and assume f is a sample from a GP prior with a known covariance kernel, k, and zero mean function, f ∼ GP(0, k). Under these assumptions, the posterior of f , given a Gaussian likelihood of data, is a GP that is analytically tractable. Classical BayesOpt is naturally myopic in its definition as a greedy one-step update (see (1) ), 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-exploitation trade-off. In this paper we follow similar ideas; however, we do not focus on regret but instead on gathering information to maximize our chances to identify x ⋆ , the maximizer of f . Maximum Identification via Hypothesis testing. Maximum identification can naturally be expressed as a multiple hypothesis testing problem, where we need to determine which of the elements in X is the maximizer. To do so, we require good estimates of the differences (or at least their signs) between individual queries f (x i ) -f (x j ); x i , x j ∈ X . For example, if f (x i ) -f (x j ) ≤ 0, then x i cannot be a maximizer. Given the current evidence, the set of arms which we cannot rule out are all potential maximizers, Z ⊂ X . At termination we report our best guess for the maximizer as:x T = arg max x∈Z µ T (x), where µ T is the predictive mean at termination time T . Suppose we are in step t out of T , then let X t be the set of previous queries, we seek to identify new X new that when evaluated minimize the probability of returning a sub-optimal arm at the end. For a given function draw f , the probability of returning a wrong maximizer z ̸ = x ⋆ f is P (µ T (z) -µ T (x ⋆ f ) ≥ 0|f ). We can then consider the worst-case probability across potential maximizers, and taking expectation over f we obtain a utility through an asymptotic upper-bound on the log-probability, indeed for large T we obtain:EQUATIONThe expectation is on the current prior (posterior up to X t ), the kernel k is the posterior kernel given observations X t ∪ X new . 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 if the values of f (z) and f (x ⋆ f ) are far apart. The non-trivial distribution of f (x ⋆ ) [29] renders the utility intractable; therefore we employ a simple and tractable upper bound on the objective (3) which can be optimized by minimizing the uncertainty among all pairs in Z:EQUATIONSuch objectives can be solved greedily in a similar way as acquisition functions in Eq. ( 1) by minimizing U over X new . Note that Fiez et al. [30] derive this objective for the same problem with linear bandits, albeit they consider the frequentist setting and (surprisingly) a 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 Eq.( 4) can be found in Appendix C.1-C.2.Utility with kernel embeddings. For illustrative purposes, consider a special case where the kernel k has a low rank due to existence of embeddings Φ(x) ∈ R m , i.e., k(x, y) = Φ(x) ⊤ Φ(y). 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,EQUATIONwhere θ are weights with a Gaussian prior. We can then rewrite the objective Eq. ( 4) as:EQUATIONThis 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 × m matrix instead of |X | × |X | (see Sec. 4). Appendix D.1 provides a utility without the finite rank-assumptions that is more involved symbolically and computationally.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 a policy for H time-steps (horizon) and obtain a trajectory, τ = (x 0 , a 0 , x 1 , a 1 , ..., x H-1 , a H-1 ). From the trajectory, we obtain a sequence of noisy observations y(τ) := {y(x 0 , a 0 ), ..., y(x H-1 , a H-1 )} s.t. y(x h ) = f (x h , a h ) + ϵ(x h , a h ),where ϵ(x h , a h ) is zero-mean Gaussian with known variance which is potentially state and action dependent. The trajectory is generated using a known transition operator P (x h+1 |x h , a h ). A Markov policy π(a h |x h ) is a mapping that dictates the probability of action a h in state x h . Hence, the state-to-state transitions are P (x h+1 , x h ) = a∈A π h (a|x h )P (x h+1 |x h , a). In fact, an equivalent description of any Markov policy π is the corresponding distribution giving us the probability of visiting a state-action pair under the policy, which we denote d π ∈ D, whereD := ∀h ∈ [H] d h | d h (x, a) ≥ 0, a,x d h (x, a) = 1, a d h (x ′ , a) = x,a d h-1 (x, a)p(x ′ |x, a)We will use this polytope to reformulate our optimization problem over trajectories. Any d ∈ 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 T times (episodes) to obtain the final dataset of the form X T = {τ i } T i=1 .Notice that the utility U in Eq. 6 depends on the states visited and hence states of the trajectory. In our notation, X t 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| H ). In fact, for stochastic transitions, we cannot pick the trajectories directly, so instead we work in the space of distributions. For a given policy, through sampling, we are able to create an empirical distribution of all the state-action pairs visited during policy executions, dπ (x, 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,EQUATIONThis 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 a maximum, as 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: min dπ∈D U(d π ), and then construct policy π via marginalization. We refer to this optimization as the planning problem. The constraint D encodes the dynamics of the MDP.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] , Zahavy 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 [11, 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. [44] tackles a 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 is 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 et al. [30] 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] [49] [50] [51] [52] [53] [54] seeks to 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, 56] ). 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-myopic analysis is non-Markovian, in Sec. 4, we show that adaptive resampling iteratively approximates this non-myoptic optimal policies in a numerically tractable way via receeding horizon planning. In our experiments we comfortably plan for over a hundred steps. 
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 + 1 given we previously queried state x h . This mean that the transition operator is P (x h+1 |x h , a) = 0 for any transition x h → x h+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 a mix of the previous two, where we obtain observations with unspecific a delay.In section 2.1 we introduced the utility for maximum identification. Using the same simplifying assumption (finite rank approximation of GPs in Sec. 2.1, Eq. ( 4)), we can show that the expected utility U can be rewritten in terms of the state-action distribution induced by X new :EQUATIONwhere V(d π ) =x,a∈X ×A dπ(x,a)Φ(x,a)Φ(x,a) ⊤ σ 2 (x,a)+ 1 T H I . The variable d π (x, a) is a state-action 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.1 in Appendix D, a condition required for the expression as a function of state-action visitations [17] . Additionally, by rewriting the objective in this form, the dependence and convexity with respect to the state-action density d π becomes clear as it is only composition of a linear function with an inverse operator. Also, notice that the constraint set is a convex polytope. Therefore we are able to use convex optimization to solve the planning problem (see Sec. 4).Input: Procedure for estimating sets of maximizers, initial point x 0 , initial set of maximizer candidates Z 0 Initialize the empirical state-action distribution d0 = 0 Set of potential maximizers Z. The definition of the objective requires the use of a set of maximizers. In the ideal case, we can say a particular input x, is not the optimum if there exists x ′ such that f (x ′ ) > f (x) with high confidence. We formalize this using the GP credible sets (Bayesian confidence sets) and define:for t = 0 to T -1 do for h = 0 to H -1 do U t,h (d π ) ← U(d π ⊕ dt,h |Z t,EQUATIONwhere 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 t .Until this point, our formulation focused on discrete S and A for ease of exposition. However, the framework is compatible with continuous state-action spaces. The probabilistic reformulation of the objective in Eq. ( 7) is possible irrespective of whether X (or A) is a discrete or continuous subset of R d . In fact, the convexity of the objective in the space of distributions is still maintained. The difference is that the visitations d are no longer probability mass functions but have to be expressed as probability density functions d c (x, a). To recover probabilities in the definition of V, we need to replace sums with integrals i.e. x∈X ,a∈A d(x) Φ(x,a)Φ(x,a) ⊤ σ(x,a) 2→ x∈X ,a∈A d c (x, a) Φ(x,a)Φ(x,a) ⊤ σ(x,a) 2 .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 a non-trivial constrained optimization problem when X is continuous. As an alternative, we propose approximating the set Z 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.1). In Appendix E.5, we numerically benchmark reasonable choices of K, and show that the performance is not significantly affected by them.The general algorithm combines the ideas introduced so far. We present it in Algorithm 1. Notice that apart from constructing the current utility, keeping track of the visited states and updating our GP model, an essential step is planning, where we need to find a policy that maximizes the utility. As this forms the core challenge of the algorithm, we devote Sec. 4 to it. In short, it solves 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.1 in Appendix C.4.Question: For each theoretical result, does the paper provide the full set of assumptions and a 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 crossreferenced.• 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 if they appear in the supplemental material, the authors are encouraged to provide a 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. 
solving the planning problem: 
 The planning problem, defined as min dπ∈D U(d π ), can be 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 5 in Appendix A for a breakdown of the different components of our solution. Following developments in Hazan et al. [34] and Mutný et al. [17] , we use the classical Frank-Wolfe algorithm [59] . It proceeds by decomposing the problem into a series of linear optimization sub-problems. Each linearization results in a policy, and we build a mixture policy consisting of optimal policies for each linearization π mix,n = {(α i , π i )} n i=1 , and α i step-sizes of Frank-Wolfe. Conveniently, after the linearization of U the subproblem on the polytope D corresponds to an RL problem with reward ∇U for which many efficient solvers exist. Namely, for a single mixture component we have,EQUATION)Due to convexity, the state-action distribution follows the convex combination, d πmix,n = n i=1 α i d πi . The optimization produces a Markovian policy due to the subproblem in Eq. ( 10) being optimized by one. We now detail how to construct a non-Markovian policies by adaptive resampling.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 h 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 t and time-point h we follow a Markov policy for a single step, the overall policy is a history-dependent non-Markov policy.We define the empirical state-action visitation distribution,dt,h = 1 tH + h ( t j=1 x,a∈τj δ x,avisited states in past trajectoriesEQUATIONwhere δ x,a denotes a delta mass at state-action (x, a). Instead of solving the objective U(d) as in Eq. ( 10), we seek to find a correction to the empirical distribution by minimizing,EQUATIONWe use the same Frank-Wolfe machinery to optimize this objective:d π t,h = arg min dπ∈ D U t,h (d π ).The distribution d π t,h 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-problem can be solved extremely efficiently to exactness using dynamic programming. As can be seen in Appendix B.4, our solving times are just a few seconds, even if planning for very long horizons. The resulting policy π can be found by marginalization π h (a|x) = d π,h (x, a)/ a d π,h (x, a), a basic property of MDPs [60] .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 {a h } T h=1 with the linear state-space model, x h+1 = Ax h + Ba h . While this formalism is not as general as it could be, it gives us a tractable sub-problem formulation common to control science scenario [61] that is practical for our experiments and captures a vast array of problems. The optimal set of actions is solved with the following problem, where we state it for the full horizon H:EQUATIONsuch that ||a h || ≤ a max , x h ∈ X , and x h+1 = Ax h + Ba h , where the known dynamics serves as constraints. Notice that instead of optimizing over the policy d π , we directly optimize over the parameterizations of the policy {a h } H h=1 . 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 [62] , our objective is non-convex and tends to focus on gathering information rather than stability. Due to the non-convexity in this parameterization, we need to solve it heuristically. We identify a number of useful heuristics to solve this problem in Appendix G. The line plots denote the best prediction regret, while the bar charts denote the percentage of 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 transition constraints.Question: Does the paper fully disclose all the information needed to reproduce the main experimental results of the paper to the extent that it 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, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.• If the contribution is a dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.• Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to 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 submissions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example (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, it 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 to have some path to reproducing or verifying the results. 
experiments: 
 Sections 5.1 -5.3 showcase real-world applications under physical transitions constraints, using the discrete version of the algorithm. Section 5.4 benchmarks against other algorithms in the continuous setting, where we consider the additive transition model of Section 4.2 with A = B = I. We include additional results in Appendix B. For each benchmark, we selected reasonable GP hyper-parameters and fixed them during the optimization. These are summarized in Appendix E.2. As we are interested maximizer identification, in discrete problems, we report the proportion of reruns that succeed at identifying the true maximum. For continuous benchmarks, we report inference regret at each iteration: Regret t = f (x * ) -f (x µ,t ), where x µ,t = arg max x∈X µ t (x). All statistics reported are over 25 different runs.Baselines. We include a naive baseline that greedily optimizes the immediate reward to showcase a method with no planning (Greedy-UCB). Likewise, we create a baseline that replaces the gradient in Eq. ( 10) with Expected Improvement [63] (MDP-EI), a weak version of planning. In the continuous settings, we 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, one using Thompson Sampling (MDP-BO-TS) and one using Upper Confidence Bound (MDP-BO-UCB). 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 see that the paths taken by MDP-BO-TS are more consistent, and the final performance is comparable to TrSnAKe. While in the asynchronous setting, we found creating the maximization set using Thompson Sampling gave a stronger performance, in the synchronous setting, UCB is preferred. LSR gives a very strong performance, comparable to MDP-BO-UCB in almost all benchmarks.Our chemical reactor benchmark synthetizes Knorr pyrzole in a transient flow reactor. In this experiment, we can control the flow-rate (residence time) τ and ratio of reactants B in the reactor. We observe product concentration at discrete time intervals and we can also change inputs at these intervals. Our goal is to find the best parameters of the reaction subject to natural movement constraints on B, and τ . In addition, we assume decreasing the flow rate of a reactor can be easily 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.The kernel. Schrecker et al. [27] indicate the reaction can be approximately represented by simple kinetics via a differential equation model. We use this information along with techniques for representing linear ODE as constraints in GP fitting [65, 66] to create an approximate ODE kernel k ode through the featurization:Φ ode (τ, B) = (1 -S(B))y (1) (τ, B) + S(B)y (2) (τ, B)where y (i) (τ, B) are equal to:γ i (B) λ (i) 2 λ (i) 1 -λ (i) 2 e λ (i) 1 τ - λ (i) 1 λ (i) 1 -λ (i) 2 e λ (i) 2 τ + 1 for i = 1, 2, where λ (i) 1 and λ (i)2 are eigenvalues of the linearized ODE at different stationary points, γ 1 (B) = B, γ 2 (B) = 1 -B, and S(x) := (1 + e -αsig(x-0.5) ) -1 is a sigmoid function. Appendix H holds the details and derivations which may be of independent interest. As the above kernel is only an approximation of the true ODE kernel, which itself is imperfect, we must account for the model mismatch. Therefore, we add a squared exponential term to the kernel to ensure a non-parametric correction, i.e.:k(τ, B) = α ode k ode (τ, B) + α rbf (τ, 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-filling and avoid sub-optimal areas because of our choice of 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 (τ R , B) = (0, 0). Figure 2 reports quantitative results and shows that the best-performing algorithm is MDP-BO.Samaniego et al. [20] investigated automatic monitoring of Lake Ypacarai, and Folch et al. [11] and Yang et al. [40] benchmarked different BayesOpt algorithms for the task of finding the largest contamination 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 3a . 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.Apart from hard constraints, we can apply our framework to state-dependent BayesOpt problems involving transitions. For example, the magnitude of noise ϵ may depend on the transition. This occurs in systems observing equilibration constraints such as a free-electron 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, x ′ ) = s(1 + w||x -x ′ || 2 ). By choosing A = X , we rewrite the problem as σ(s, a) = s(1 + w||x -a|| 2 ). 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 3b reports the comparison between worst-case stateless BO and our algorithm. Our approach substantially improves performance.V(d π ) = x,a∈X d π (x, a) 1 σ 2 (x,a) Φ(x)Φ(x) ⊤ + 1 T H I. FigureWe benchmark on a variety of classical BayesOpt problems while imposing local movement 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 [68] , which we treat as asynchronous as per Folch et al. [11] . 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 in 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 it 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.Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental material?Answer: [Yes]Justification: The code will be made public upon acceptance, and an anonymized version is included for the review process.Guidelines:• The answer NA means that paper does not include experiments requiring code.• Please see the NeurIPS code and data submission guidelines (https://nips.cc/ public/guides/CodeSubmissionPolicy) for more details.• While we encourage the release of code and data, we understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for a 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: //nips.cc/public/guides/CodeSubmissionPolicy) for more details.• The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.• The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only a subset of experiments are reproducible, they should state which ones are 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 to the paper) is recommended, but including URLs to data and code is permitted. 
conclusion: 
 We considered transition-constrained BayesOpt problems arising in physical sciences, such as chemical 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.Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?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.• The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them.• The full details can be provided either with the code, in appendix, or as supplemental material. 
experiment statistical significance: 
 Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?Answer: [Yes] Justification: We provide regret quantiles in our experiments, reproduce on 25 random seeds and use a 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, confidence intervals, or statistical significance tests, at least for the experiments that support the 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).• It should be clear whether the error bar is the standard deviation or the standard error of the mean.• Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.• The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point 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 or 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 to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
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)? 
new assets: 
 Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?Answer: [Yes] Justification: The only new asset will be the code implementation which will be published upon acceptance.Guidelines:• The answer NA means that the paper does not release new assets.• Researchers should communicate the details of the dataset/code/model as part of 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 is used.• At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file. 
crowdsourcing and research with human subjects: 
 Question: 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 as details about compensation (if any)?Answer: [NA] Justification: No crowd-sourcing or human subjected were used.Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.• Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should be included in the main paper.• According to the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector. 
institutional review board (irb) approvals or equivalent for research with human subjects: 
 Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?Answer: [NA]  
 
 Appendix: In Figure 5 we summarize how our algorithm creates non-Markovian policies for maximizer identification and the corresponding connections to other works in the literature. In black we 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, we take inspiration from hypothesis testing and focus on (b) the variance reduction in a set of maximizers, which leads to our (c) acquisition function. The objective is the same as Fiez et al. [30] introduced in the linear bandits literature from a frequentist perspective. To optimize it, we follow developments in Mutný et 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. 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. We also run the Ypacarai experiment with immediate feedback. To increase the difficulty, we used large observation noise, σ 2 = 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 it struggles to identify which of the two optima is the global one.We also include results for the Knorr pyrazole synthesis with immediate feedback. In this case we observe very strong early performance from MDP-BO, but by the end MDP-EI is comparable. The greedy method performs very poorly.Finally, we also include additional results on more synthetic benchmarks for both synchronous and asynchronous feedback. The results are shown in Figures 8 and 9 . The results back the conclusions in the main body. All benchmarks do well in 2-dimensions while highlighting further that MDP-BO-UCB and LSR can be much stronger in the synchronous setting than Thompson Sampling planning-based approaches (with the one exception of the Levy function). 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, We call this the selection the recommendation rule. We focus on this recommendation rule as this rule is interpretable to the facilitator of the analysis and experimenters. In this derivation we use that f = θ ⊤ Φ(x). More commonly, the notation ⟨θ, Φ(x)⟩ is used, where the inner product is potentially infinite dimensional. We use use ⊤ notation for simplicity for both cases. Same is true for any other functional estimates, e.g., for the posterior mean estimate, we use µ t (x) = Φ(x) ⊤ µ t . The inner product is in the reproducing kernel Hilbert space associated with the kernel k.x T = arg max x∈Z µ T (x).Now, suppose there is a given f (we will take expectation over it later), then there is an x ∈ X achieving optimum value, denoted x ⋆ f (suppose unique for this development here). Hence, we would like to model the risk associated with predicting a fixed z ̸ = x ⋆ f , which is still in Z at time T . Suppose we are at time t, we develop the utility to gather additional data X new on top of the already acquired data X t . These should improve the discrepancy of the true answer, and the reported value the most.Suppose there are two elements in Z simple = {z, x ⋆ f }. We will generalize to a composite hypothesis later. In two-element case, the probability of the error in incurred due to selecting z is:P (µ T (z) -µ T (x ⋆ f ) ≥ 0|f ). The randomness here is due to the observations y = f (X new ) + ϵ that are used to fit the estimator µ T (x). Namely due to ϵ ∼ N (0, σ 2 ). Given f (equivalently θ), the distribution of our estimator (namely the posterior mean) is Gaussian. Hence, given f :µ T ∼ N ((V T + I H ) -1 V T θ, σ 2 (V T + I H ) -1 V T (V T + I H ) -1 ),where V = T i=1 1 σ 2 Φ(x i )Φ(x i )⊤ is an operator on the reproducing kernel Hilbert space due to k as H → H, and I H the identity operator on the same space. This is the posterior over the posterior mean as a function.A posterior over the specific evaluation isµ T (z) -µ T (x ⋆ f ) ∼ N (θ ⊤ (V T + I H ) -1 V T (Φ(z) -Φ(x ⋆ f )) a , σ 2 (Φ(z) -Φ(x ⋆ f ) ⊤ (V T + I H ) -1 V T (V T + I H ) -1 (Φ(z) -Φ(x ⋆ f )) b 2).We can now bound the probability of making an error using a Gaussian tail bound inequality:P(µ T (z) -µ T (x * ) ≥ 0) = P(µ T (z) -µ T (x * ) ≥ a z + (-a z )) ≤ e -a 2 2b 2with the caveat that the inequality only holds when the a z is negative. However note that a z → f (z) -f (x * ) < 0 as T → ∞ therefore it will hold once T is large enough. From this we can take logarithms and then the expectation across the randomness in the GP:E f ∼GP [log P(µ T (z) -µ T (x * )|f )] ≤ - 1 2 E f ∼GP a 2 z b 2 zwhich is called the log Bayes' factor and is expected log failure rate for the set of potential maximizers Z simple . The expectation is over the posterior including the evaluations X t (or prior at the very beginning of the procedure). In fact, we can think of the posterior as being the new prior for the future 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:EQUATIONFor moderate to large T ≫ 0, we can upper bound this objective via elegant argument to yield a very transparent objective:EQUATIONwhere we have used an lower and upper bound on the a z and b z , respectively as follows:a 2 z = (θ ⊤ (V T + I H ) -1 V T (Φ(z) -Φ(x ⋆ f ))) 2 = (θ ⊤ (V T + I H ) -1 (V T + I H -I H )(Φ(z) -Φ(x ⋆ f ))) 2 = (θ ⊤ (Φ(z) -Φ(x ⋆ f )) -θ ⊤ (V T + I H ) -1 (Φ(z) -Φ(x ⋆ f ))) 2 T ≫0 ≈ (θ ⊤ (Φ(z) -Φ(x ⋆ f ))) 2 = (f (z) -f (x ⋆ f )) 2 b 2 z = σ 2 (Φ(z) -Φ(x ⋆ f ) ⊤ (V T + I H ) -1 V T (V T + I H ) -1 (Φ(z) -Φ(x ⋆ f )) ≤ σ 2 (Φ(z) -Φ(x ⋆ f )) ⊤ (V T + I H ) -1 (Φ(z) -Φ(x ⋆ f )) = k X (z, x ⋆ f). 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 E f for large T . The objective Eq. ( 15) is intractable due to the expectation of the prior and which involves expectation over the maximum f (x ⋆ f ), which is known to be very difficult to estimate. Interestingly, the denominator is independent of f if we adopt the worst-case perspective over the x ⋆ f , 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 to the smallest gap gap(f ) between two arms in X . Hence, we can upper bound the objective as:-E f sup z∈Z\{x ⋆ f } (f (z) -f (x ⋆ f )) 2 k X∪Xnew (z, x ⋆ f ) ≤ -sup z∈Z\{x ⋆ f } E f [gap(f )] k X∪Xnew (z, x ⋆ f ) ≤ -E f [gap(f )] sup z∈Z\{x ⋆ f } 1 Var f (z) -f (x ⋆ f )|X t ∪ X new .As the constant in front of the objective does not influence the optimization problem, we do not need to consider it when defining the utility. Furthermore, in order to minimise the probability of an error we can just minimise the variance in the denominator instead (since arg min x -g(x) is equivalent to arg min x 1 g(x) 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 Z:EQUATIONSurprisingly, 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 as they try to minimize T for a 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 z ) can be neglected.Let us now briefly summarize the Nyström approximation [31, 69] . Given a kernel k(•, •), and a data-set X, we can choose a sub-sample of the data x1 , ..., xm . Using this sample, we can create a low r-rank approximation of the full kernel matrixKr = K b K † K b where K b = [k(x i , xj )] N ×m , K = [k(x i , xj )] m×m and K † denotes the pseudo-inverse operation. We can then define the Nyström features as:EQUATION)where Dr is the diagonal matrix of non-zero eigenvalues of Kr and Vr the corresponding matrix of eigenvectors. It follows that we obtain a finite-dimensional estimate of the GP:EQUATIONwhere Φ(x) = (ϕ 1 (x), . . . ϕ m (x)) T , and θ are weights with a Gaussian prior.The fact that our objective is derived using Bayesian decision theory makes it well-rooted in theory.In addition to the derivation of Section C.1, 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 to Z as Z t . Let us contemplate for a second what could the optimal policy. As the set of Z t is changing, we follow the line of work of started by Russo [70] and introduce an optimal algorithm that knows the true x ⋆ f for each possible realization of the prior f . In other words, its an algorithm that any time t, would follow:d ⋆ t = min d∈D E f max z∈Zt\{x ⋆ f } k dt⊕d (z, x ⋆ f ) ,where in the above dt ⊕ d represents the weighted sum as in the main Algorithm 1 that scales the distributions properly according to t 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 z ′ ∈ Z, but fixes it to the value x ⋆ f that the hypothetical algorithm has privileged access to. To eliminate the cumbersome notation, we will refer to the objectives as U(d|Z t , Z t ) as the objective used by our algorithm (real execution) and U(d|Z t , {x ⋆ f }), as the objective that the privileged algorithm is optimizing which serves as theoretical baseline.The visitation of d ⋆ t represents the best possible investment of the resources (of the size T -t) to execute at time t had we known the x ⋆ f instead of only Z t . This is interpreted as if the modeler knows x ⋆ f , and sets up an optimal curriculum that is being shown to an observer in order to convince him/her of that x ⋆ f 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 t . Hence, our goal is to show that we are closely tracking the performance of these optimal policies in time t, and eventually there is little difference between our sequence of executed policies (visitations) dt and the algorithm optimal d ⋆ t . In order to prove the theorem formally, we need to assume that Z t is decreasing. The rate at which this set is decreasing determines the performance of the algorithm to a large extent. Namely, we assume that given two points in time, having the same empirical information dt . Given, f , supposesup d∈D |d ⊤ (∇U( dt |Z t , Z t ) -∇U( dt |Z t , {x ⊤ f }))| ≤ C t . (⋆)As we gather information in our procedure the, {x ⋆ f } ⊂ Z t ⊆ Z t-1 , but the exact decrease depends on how Z t is constructed. We leave the particular choice for C t to make the above hold for future work. We conjecture that this is decreasing as C t ≈ γt Proposition C.1. Assuming episodic feedback, and suppose that for any Z,EQUATION)for α ∈ (0, 1) and η, h ∈ ∆ p , L := max η,α L η,α 4. condition in (⋆) holds with Bayesian posterior inference, we can show that the Algorithm 1 satisfied for the sequences of iterates { dt } T t=1 :1 T T -1 t=1 U(d t |Z t , {x ⋆ f }) -U(d ⋆ t |Z t , {x ⋆ f }) ≤ O 1 T T -1 t=1 C t + L log T T + B √ T log 1 δ ,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 of the experimental campaign T is increasing, we are increasingly converging to the optimal allocation of the experimental resources on 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 to the best possible allocation of the resources. Now, we prove the Proposition. The proof is an extension of the Theorem 3 in [17] . Whether the objective satisfied the above conditions depends on the set X . Should the objective not satisfy smoothness, it can be easily extended by using the Nesterov smoothing technique as explained in 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. 3. in Mutný et al. [17] .Let us start by notation. We will use the notation that U t is the privileged objective U(d|Z t , {x f t }), while the original objective will be specified as U(d|Z t , Z t ).First, what we follow in the algorithm:EQUATIONThe executed visitation is simply generated via sampling a trajectory from q t . Let us denote the empirical visiation of the trajectory as δ t ,EQUATIONFor the analysis, we also need the best greedy step for the unknown (privileged) objective U asEQUATIONLet us start by considering the one step update:U t ( dt+1 ) = U t dt + 1 t + 1 (δ t -dt ) L-smooth ≤ U t ( dt ) + 1 t + 1 ∇U t ( dt ) ⊤ (δ t -dt ) + L 2(1 + t) 2 δ t -dt 2 U( dt+1 ) bounded ≤ U t ( dt ) + 1 t + 1 ∇U t ( dt ) ⊤ (δ t -dt ) + L (1 + t) 2 = U t ( dt ) + 1 t + 1 ∇U t ( dt ) ⊤ (q t -dt ) + 1 t + 1 ∇U t ( dt ) ⊤ (-q t + δ t ) ϵt + L (1 + t) 2We will now carefully insert and subtract two set of terms depending on the real objective so that we can bound them using (⋆):= U t ( dt ) + 1 t + 1 ∇(U t ( dt ) ⊤ -∇U t ( dt |Z t , Z t ) ⊤ )(q t -z t ) + 1 1 + t U t ( dt ) ⊤ (z t -dt ) + 1 1 + t ϵ t + L (1 + t) 2 Using ⋆ ≤ U t ( dt ) + 2 1 1 + t C t + 1 t + 1 U t ( dt ) ⊤ (z t -dt ) + 1 1 + t ϵ t + L (1 + t) 2Carrying on,U t ( dt+1 ) (22) ≤ U t ( dt ) + 1 t + 1 ∇U t ( dt ) ⊤ (d ⋆ t -dt ) + 1 1 + t ϵ t + L (1 + t) 2 convexity ≤ U t ( dt ) - 1 t + 1 (U t ( dt ) -U t (η ⋆ t )) + 1 1 + t ϵ t + L (1 + t) 2 + 1 1 + t C t U t ( dt+1 ) -U t (d ⋆ t ) ≤ U t ( dt ) -U t (d ⋆ t ) - 1 t + 1 (U t ( dt ) -U t (d ⋆ t )) + 1 1 + t ϵ t + L (1 + t) 2 + 1 1 + t C t ≤ t 1 + t U t ( dt ) -U t (d ⋆ t ) + 1 1 + t ϵ t + L (1 + t) 2 + 1 1 + t C t = t 1 + t U t ( dt ) -U t (d ⋆ t ) + 1 1 + t ϵ t + L (1 + t) 2 + 1 1 + t C tNow multiplying by t + 1 both sides, and summing on 1T -1 T -1 t=1 . Using the shorthand ρ t ( dt ) = U t ( dt ) -U t (d ⋆ t )we get:1 T T -1 t=1 (t + 1)ρ t+1 ( dt+1 ) ≤ 1 T T -1 t=1 tρ t ( dt ) + 1 T T -1 t=1 (ϵ t + C t + L/(1 + t))First notice that T -1 t=1 1 1+t ≤ L log TT . The rest is eliminated by the reccurence of the terms, and using that U(d|Z t , {x ⋆ f } ≤ U(d|Z t-1 , {x ⋆ f }) for any d. This is due to set Z t decreasing over time. We report the result in asymptotic notation as function of T and log(1/δ).For the main objective we try to optimize over a subset of T trajectories X = {τ i ∈ X H } T i=1 . Let X H be the set of sequences of inputs τ = (x 1 , ..., x H ) where they consist of states in the search space X . Furthermore, assume there exists, in the deterministic environment, a constraint such that x h+1 ∈ C(x h ) 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:EQUATIONWe define the objective as:EQUATION)Our goal is to show that optimization over sequences can be simplified to state-action visitations as in Mutný et al. [17] . For this, we require that the objective depends additively involving terms x, a separately. We formalize this in the next result. In order to prove the result, we utilize the theory of reproducing kernel Hilbert spaces [72] .Lemma D.1 (Additivity of Best-arm Objective). Let X be a collection of t trajectories of length H. Assuming that f ∼ GP(0, k). Assuming that k has Mercer decomposition as k(x, y)= k λ k ϕ k (x)ϕ k (y). f (x) = k ϕ k (x)θ k θ k ∼ N (0, λ k ).Let d X be the visitation of the states-action in the trajectories in X, asd X = 1 T H T t=1x,a∈τt δ x,a , where the δ x,a represent delta function supported on x, a. Then optimization of the objective Eq. ( 23) can be rewritten as:U (d X ) = 1 T H max z,z ′ ∈Z ||Φ(z) -Φ(z ′ )|| 2 V(d X ) -1 ,whereV(d) = i x,a∈τi d(x, a)Φ(x)Φ(x) ⊤ + Iσ 2 /(T H) is a operator V(d) : H k → H k , the norm is RKHS norm, and Φ(z) k = ϕ k (z).Proof. Notice that the posterior GP of any two points z, z′ is (f (z), f (z ′ )) = N ((µ(z), µ(z ′ )), K z,z ′ ), where K z,z ′ is posterior kernel (consult Rasmussen and Williams [28] for details) defined via a posterior kernel kX (z, z ′ ) = k(z, z ′ ) -k(z, X)(K(X, X) + σ 2 I) -1 k(X, z ′ ).Utilizing k(z, z ′ ) = Φ(z) ⊤ Φ(z) (RKHS inner product) with the Mercer decomposition we know that k t (z) = Φ(X)ϕ(z). Applying the matrix inversion lemma, the above can be written as usingV = T t=1 x∈τt Φ(x)Φ(x) ⊤ + σ 2 I H . k X (z, z ′ ) = k(z, z ′ ) -k t (z) ⊤ (K X,X + σ 2 I) -1 k t (z ′ ) Mercer = Φ(z) ⊤ Φ(z ′ ) -Φ(z) ⊤ Φ(X) ⊤ (Φ(X)Φ(X) ⊤ + σ 2 I) -1 Φ(X)Φ(z ′ ) Lemma D.3 = Φ(z) ⊤ Φ(z ′ ) -Φ(z) ⊤ V -1 (V -Iσ 2 )Φ(z ′ ) = Φ(z) ⊤ V -1 VΦ(z ′ ) -Φ(z) ⊤ V -1 (V -Iσ 2 )Φ(z ′ ) = Φ(z) ⊤ V -1 (V -V + Iσ 2 )Φ(z ′ )Leading finally to:EQUATION) Let us calculate Var[f (z) -f (z ′ )|X].The variance does not depend on the mean. Hence,Var [f (z) -f (z ′ )|X] = Var(f (z)) -Var(f (z ′ )) -2Cov(f (z), f (z ′ )) = k X (z, z) + k X (z ′ , z ′ ) -2k X (z, z ′ ) (25) = (Φ(z) -Φ(z ′ )) T t=1 x∈τt Φ(x)Φ(x) ⊤ + σ 2 I H k -1 (Φ(z) -Φ(z ′ )) = (Φ(z) -Φ(z ′ )) T H T H T t=1 x∈X #(x ∈ τ t )Φ(x)Φ(x) ⊤ + σ 2 I H k -1 (Φ(z) -Φ(z ′ ))to arrive at:Var [f (z) -f (z ′ )|X] = (Φ(z) -Φ(z ′ )) T H x∈X d(X)Φ(x)Φ(x) ⊤ + σ 2 T H I H k -1 (Φ(z) -Φ(z ′ )) (26)The symbol # counts the number of occurrences. Notice that we have been able to show that the objective decomposes over state-action visitations as d X decomposes over their visitations Note that the objective equivalence does not imply that optimization problem in Eq. ( 23) is equivalent to finding,EQUATIONIn other words, optimization over trajectories and optimization over d π ∈ D 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 a 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 as 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, we state the auxiliary lemma. We make use of the Sherman-Morrison-Woodbury (SMW) formula, [76] :EQUATION)Here we do the opposite, and invert an n × n matrix instead of a q × q one. to show the following.EQUATION) Note that instead of inverting n × n matrix, we can invert a q × q matrix. Proof.A ⊤ (AA ⊤ + ρ 2 I) -1 SMW = A ⊤ (ρ -2 I -ρ -2 A(ρ 2 I + A ⊤ A) -1 A ⊤ ) = (ρ -2 I -ρ -2 A ⊤ A(ρ 2 I + A ⊤ A) -1 )A ⊤ = (ρ -2 (ρ 2 I + A ⊤ A) -ρ -2 A ⊤ A)(ρ 2 I + A ⊤ A) -1 A ⊤ = (A ⊤ A + ρ 2 I) -1 A ⊤ D.To apply our method, we find ourselves having to frequently solve RL sub-problems where we try to maximize x,a d(x, a)∇F (x, a). To approximately solve this problem in higher dimensions, it becomes very important to understand what the linearized functional looks like.Remark D.4. Assume the same black-box model as in Lemma D.1, and further assume that we have a mixture of policies π mix with density d πmix , such that there exists a set X mix satisfying d πmix = 1 N x∈Xmix δ x for some integer N . Then:∇F (d πmix )(x, a) ∝ -(Cov[f (z * ), f (x)|X mix ] -Cov[f (z ′ * ), f (x)|X mix ]) 2 where z * , z ′ * = arg max z,z ′ ∈Z Var[f (z) -f (z ′ )].Proof. To show this, we begin by defining:Σ θ,d = x∈X Φ(x)Φ(x) T d(x) + σ 2 I -1 z * , z ′ * = arg max z,z ′ ∈Z ||Φ(z) -Φ(z ′ )|| 2 Σ θ,d z * = Φ(z * ) -Φ(z ′ * )In the definition above we 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 Theorem that:∇U (d)(x) = ∇z T * Σ θ,d z * = ∇Tr z * zT * Σ θ,d = Tr z * zT * ∇Σ θ,d = -Tr z * zT * Σ θ,d Φ(x)Φ(x) T Σ θ,d (as ∂K -1 = -K -1 (∂K)K -1 )) = -Tr zT * Σ θ,d Φ(x)Φ(x) T Σ θ,d z * = -zT * Σ θ,d Φ(x) Φ(x) T Σ θ,d z * ∝ -(Cov[f (z * ), f (x)] -Cov[f (z ′ * ), f (x)]) (Cov[f (x), f (z * )] -Cov[f (x), f (z ′ * )]) = -(Cov[f (z * ), f (x)] -Cov[f (z ′ * ), f (x)]) 2In this section we provide implementation details, and show some studies into the effects of specific hyper-parameters. We note that the implementation code will be made public after public review.We give details of the two methods used for approximating the set of potential maximizers. In particular, we first focus on Thompson Sampling [58] :Z (T S) cont = arg max x∈Xc f i (x) : f i ∼ GP(µ t , σ t ) K i=1where K is a 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 set using BayesOpt through the UCB acquisition function [71] :Z (U CB) cont = arg max x∈Xc µ t (x) + β i σ t (x) : β i ∈ B K i=1where B = linspace(0, 2.5, K) which serves as scaling for the size of set Z(U CB) cont. Both cases reduce optimization over Z to enumeration as with discrete cases.For all benchmarks, aside from the knorr pyrazole synthesis example, we use a standard squared exponential kernel for the surrogate Gaussian Process:k rbf (x, x ′ ) = σ 2 rbf exp - ||x -x ′ || 2 2 2ℓ rbfwhere σ 2 rbf is the prior variance of the kernel, and ℓ rbf the kernel. We fix the values of all the hyperparameters a priori and use the same for all algorithms. The hyper-parameters for each benchmarks are included in Table 2 .For the knorr pyrazole synthesis example, we further set α ode = 0.6, α rbf = 0.001, k 1 = 10, k 2 = 874, k 3 = 19200, α sig = 5. Recall we are using a finite dimensional estimate of a GP such that:EQUATIONin this case we set a prior to the ODE weight such that ω ode ∼ N (0.6, 0.0225). This is incorporating two key pieces of prior knowledge that (a) the product concentration should be positive, and (b) we expect a maximum product concentration between 0.15 and 0.45.The number of features for each experiment, M , is set to be M = |X | in the discrete cases and M = min 2 5+d , 512 where d is the problem dimensionality.In the case of Local Search Region BayesOpt (LSR) [22] we set the exploration hyper-parameter to be γ = 0.01 in all benchmarks.Table 2 : Benchmark and hyper-parameter information. ∆ max represents the size of the box constraints in the traditional benchmarks. For the synchronous benchmarks and for SnAr we used a noise level of σ 2 = 0.001. For the asynchronous benchmarks, and the knorr pyrazole example we used σ 2 = 0.0001. For the Ypacarai example we used σ 2 = 0.001 and σ 2 = 0.01 for the episodic and immediate feedback respectively. Samaniego et al. [20] investigated the use of Bayesian Optimization for monitoring the lake quality in Lake Ypacarai in Paraguay. We extend the benchmark to include additional transition constraints, as well as initial and end-point constraints. These are all shown in Figure 12 .We use the simulator from Mutný et al. [67] that optimizes quadrupole magnet orientations for our experiment with varying noise levels. We use a 2-dimensional variant of the simulator. We discretize the system on 10 × 10 grid and assume that the planning horizon H = 100. The simulator itself is a GP fit with γ = 0.4, hence we use this value. Then we make a choice that the noise variance is proportional to the change made as σ 2 (x, a) = s(1 + w||x -a|| 2 ), 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 desired. We showcase the difference to classical BayesOpt, which uses the worst-case variance σ = sup x,a s(1 + w||x -a|| 2 ) for modeling as it does not take into account the state in which the system is. We see that the absence of state modeling leads to a dramatic decrease in performance as indicated by much higher inference regret in Figure 3b . 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 the Ypacari example (see Figure 14 ) we see very little difference in the results, while in the Knorr pyrazole synthesis (see Figure 13 ) we observe a much bigger difference. A single component gives a much stronger performance than multiple ones -we conjecture this is because the optimum is on the edge of the search space, and adding more components makes the policy stochastic and less likely to reach the boarder (given episodes are of length ten and ten right-steps are required to reach the boarder).Overall, it seems the performance of a single component is better or at worst comparable as using multiple components. This is most likely due to 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. . Overall we found the performance of the algorithm to be fairly robust to the size of the set in all benchmarks, with a higher K generally leading to a little less spread in the performance. Our objective is motivated by hypothesis testing between different arms (options) z and z ′ . In particular,EQUATIONOne could maximize the information of the location of the optimum, as it has a 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 via information-theoretic acquisition functions [29, [77] [78] [79] . 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 [80] . This is chiefly due to problem of dealing with the distribution of f (x ⋆ ). Defining a posterior value for f (z) is easy. Using, this and the worst-case perspective, an alternative way to approximate the best-arm objective, could be:EQUATIONWhat are we losing by not considering the differences? The original objective corresponds to the X Y-allocation in the bandits literature. The modified objective will, in turn, correspond to the Gallocation, which has been argued can perform arbitrarily worse as it 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 to be as expected: performance was very similar in general, however in some cases not considering the differences led to much poorer performance. As an example, see Figure 17 for results on the synchronous Branin2D benchmark.From remark D.4 it becomes clear that for decreasing covariance functions, such as the squared exponential, ∇F will consist of two modes around z * and z ′ * The sub-problem seems to find a 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-problem objective in (13) for the shortest paths from x t-1 → z * and x t-1 → z ′ * , 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 a 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 a good and quick approximation.The kernel is based on the following ODE model, which is well known in the chemistry literature and given in [27] .EQUATION)and then:dy 1 dt = R 2 dy 2 dt = -R 1 dy 3 dt = -R 1 dy 4 dt = R 1 -R 2 dy 5 dt = R 1 + R 2Our main goal is to optimize the product concentration of the reaction, which is given by y 1 . We do this by sequentially querying the reaction, where we select the residence time, and the initial conditions of the ODE, in the form y 0 = [0, A, B, 0, 0], where A = 1 -B.Due to the non-linearity in Eq. (34) we are unable to fit a 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 = {y 1 = a 1 , y 2 = b 1 , y 3 = 0, y 4 = 0, y 5 = c 1 |a 1 , b 1 , c 1 ∈ R} S 2 = {y 1 = a 2 , y 2 = 0, y 3 = b 2 , y 4 = 0, y 5 = c 2 |a 1 , b 1 , c 1 ∈ R}And the Jacobian of the system is:J =      0 0 0 k 3 0 0 -k 1 y 3 -k 1 y 2 k 2 y 5 k 2 y 4 0 -k 1 y 3 -k 1 y 2 k 2 y 5 k 2 y 4 0 k 1 y 3 k 1 y 2 -k 2 y 5 -k 3 -k 2 y 4 0 k 1 y 3 k 1 y 2 -k 2 y 5 + k 3 -k 2 y 4     Giving:J 1 = J| S1 =      0 0 0 k 3 0 0 0 -k 1 b 1 k 2 c 1 0 0 0 -k 1 b 1 k 2 c 1 0 0 0 k 1 b 1 -k 2 c 1 -k 3 0 0 0 k 1 b 1 -k 2 c 1 + k 3 0      J 2 = J| S2 =      0 0 0 k 3 0 0 -k 1 b 2 0 k 2 c 2 0 0 -k 1 b 2 0 k 2 c 2 0 0 k 1 b 2 0 -k 2 c 2 -k 3 0 0 k 1 b 2 0 -k 2 c 2 + k 3 0     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:d⃗ y dt = J 1 ⃗ y d⃗ y dt = J 2 ⃗ yWe focus on the first system for now. The matrix has the following eigenvalues:λ 1,2 = - 1 2 b 1 k 1 + c 1 k 2 + k 3 ± b 2 1 k 2 1 + c 2 1 k 2 2 + k 2 3 + 2b 1 c 1 k 1 k 2 -2k 3 (b 1 k 1 -c 1 k 2 )λ 3,4,5 = 0Note that the three eigenvalues give us the corresponding solution based on their (linearly separable) eigenvectors:v 3 = [1 0 0 0 0] , v 4 = [0 1 0 0 0] , v 5 = [0 0 0 0 1] ⃗ y(t) = p 3 v 3 + p 4 v 4 + p 5 v 5where p i 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:b 2 1 k 2 1 + c 2 1 k 2 2 + k 2 3 + 2b 1 c 1 k 1 k 2 -2k 3 (b 1 k 1 -c 1 k 2 ) ≥ b 2 1 k 2 1 + k 2 3 -2b 1 k 1 k 3 = (b 1 k 1 -k 3 )2 ≥ 0 and therefore all eigenvalues will always be real. Therefore we can write down the solution as:⃗ y(t) = p 1 v 1 e λ1t + p 2 v 2 e λ2t + p 3 v 3 + p 4 v 4 + p 5 v 5where we ignore the case of repeated eigenvalues for simplicity (this is the case where b2 1 k 2 1 + 2b 1 c 1 k 1 k 2 + c 2 1 k 2 2 -2k 3 (b 1 k 1 -c 1 k 2 ) + k 23 is exactly equal to zero). We further note that the eigenvalues will be non-negative as:b 1 k 1 + c 1 k 2 + k 3 = (bk 1 + c 1 k 2 + k 3 ) 2 = b 2 1 k 2 1 + c 2 k 2 2 + k 2 3 + 2b 1 c 1 k 1 k 2 + 2c 1 k 2 k 2 + 2b 1 k 1 k 3 ≥ b 2 1 k 2 1 + c 2 1 k 2 2 + k 2 3 + 2b 1 c 1 k 1 k 2 + 2c 1 k 2 k 2 -2b 1 k 1 k 3 therefore: λ 1 ≤ λ 2 ≤ 0which means the solutions will always be a linear combination of exponentially decaying functions of time plus constants.The eigenvectors have the closed form:v 1 =          1, 1 2 b 1 k 1 + c 1 k 2 -k 3 + b 2 1 k 2 1 + 2b 1 c 1 k 1 k 2 + c 2 1 k 2 2 -2(b 1 k 1 -c 1 k 2 )k 3 + k 2 3 /k 3 , 1 2 b 1 k 1 + c 1 k 2 -k 3 + b 2 1 k 2 1 + 2b 1 c 1 k 1 k 2 + c 2 1 k 2 2 -2(b 1 k 1 -c 1 k 2 )k 3 + k 2 3 /k 3 , -1 2 b 1 k 1 + c 1 k 2 + k 3 + b 2 1 k 2 1 + 2b 1 c 1 k 1 k 2 + c 2 1 k 2 2 -2(b 1 k 1 -c 1 k 2 )k 3 + k 2 3 /k 3 , -1 2 b 1 k 1 + c 1 k 2 -3k 3 + b 2 1 k 2 1 + 2b 1 c 1 k 1 k 2 + c 2 1 k 2 2 -2(b 1 k 1 -c 1 k 2 )k 3 + k 2 3 /k 3          =      1, -λ 1 /k 3 -1, -λ 1 /k 3 -1, λ 1 /k 3 , λ 1 /k 3 + 2      v 2 =          1, 1 2 b 1 k 1 + c 1 k 2 -k 3 -b 2 1 k 2 1 + 2b 1 c 1 k 1 k 2 + c 2 1 k 2 2 -2(b 1 k 1 -c 1 k 2 )k 3 + k 2 3 /k 3 , 1 2 b 1 k 1 + c 1 k 2 -k 3 -b 2 1 k 2 1 + 2b 1 c 1 k 1 k 2 + c 2 1 k 2 2 -2(b 1 k 1 -c 1 k 2 )k 3 + k 2 3 /k 3 , -1 2 b 1 k 1 + c 1 k 2 + k 3 -b 2 1 k 2 1 + 2b 1 c 1 k 1 k 2 + c 2 1 k 2 2 -2(b 1 k 1 -c 1 k 2 )k 3 + k 2 3 /k 3 , -1 2 b 1 k 1 + c 1 k 2 -3k 3 -b 2 1 k 2 1 + 2b 1 c 1 k 1 k 2 + c 2 1 k 2 2 -2(b 1 k 1 -c 1 k 2 )k 3 + k 2 3 /k 3          =      1, -λ 2 /k 3 -1, -λ 2 /k 3 -1, λ 2 /k 3 , λ 2 /k 3 + 2     We are optimizing over initial set of conditions y 0 = [0, A, B, 0, 0], so solving for the specific values of the constants gives:p 1 = λ 2 λ 1 -λ 2 B p 2 = - λ 1 λ 1 -λ 2 B p 3 = B p 4 = A -B p 5 = 2BFinally, since we are setting A = 1 -B and we are only optimizing the first component of ⃗ y we can obtain it in closed form: The second ODE is very similar to the first, recall it depends has the following matrix:EQUATIONJ 2 = J| S2 =      0 0 0 k 3 0 0 -k 1 b 2 0 k 2 c 2 0 0 -k 1 b 2 0 k 2 c 2 0 0 k 1 b 2 0 -k 2 c 2 -k 3 0 0 k 1 b 2 0 -k 2 c 2 + k 3 0     The resulting ODE is symmetric to the alternate linearization giving the same solution:y(t) = p 1 v 1 e λ1t + p 2 v 2 e λ2t + p 3 v 3 + p 4 v 4 + p 5 v 5with the only difference being the eigenvectors now are:v 3 = [1 0 0 0 0] , v 4 = [0 0 1 0 0] , v 5 = [0 0 0 0 1]which in turn leads to solutions of the form:y 1 (t, B) = λ 2 λ 1 -λ 2 Ae λ1t - λ 1 λ 1 -λ 2 Ae λ2t + A = A λ 2 λ 1 -λ 2 e λ1t - λ 1 λ 1 -λ 2 e λ2t + 1where A = 1-B. Note that we now have four different eigenvalues, which depend on the linearization points:λ (1) 1,2 = - 1 2 b 1 k 1 + c 1 k 2 + k 3 ± b 2 1 k 2 1 + c 2 1 k 2 2 + k 2 3 + 2b 1 c 1 k 1 k 2 -2k 3 (b 1 k 1 -c 1 k 2 ) λ (2) 1,2 = - 1 2 b 2 k 1 + c 2 k 2 + k 3 ± b 2 2 k 2 1 + c 2 2 k 2 2 + k 2 3 + 2b 2 c 2 k 1 k 2 -2k 3 (b 2 k 1 -c 2 k 2 )Giving solutions:y (1) 1 (t, B) = B λ (1) 2 λ (1) 1 -λ (1) 2 e λ (1) 1 t - λ (1) 1 λ (1) 1 -λ (1) 2 e λ (1) 2 t + 1 y (2) 1 (t, B) = A λ (2) 2 λ (2) 1 -λ (2) 2 e λ (2) 1 t - λ (2) 1 λ (2) 1 -λ (2) 2 e λ (2) 2 t + 1Due to the length of the derivation, we confirm that our analysis is correct by comparing the numerical solution of the ODE to the exact solution we found in Figure 18 . Finally, we look at interpolating between the two solutions; so given the solutions y (1) (t, B) and y (2) (t, B) corresponding to the linearization with stationary point in S 1 and S 2 respectively, we consider a solution of the form:EQUATIONwhere S(x) := (1 + e -αsig(x-0.5) ) -1 is a sigmoid function centered at B = 0.5 and where we have introduced a new hyper-parameter α sig . Finally, given Eq. ( 37) we can obtain the kernel. In particular, we want (37) to be a feature we are predicting on; therefore the kernel is simply the (dot) product of the features therefore:k ode ((t, B), (t ′ , B ′ )) = y(t, B|k 1 , k 2 , k 3 , α) × y(t ′ , B ′ |k 1 , k 2 , k 3 , α)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:k joint ((t, B), (t ′ , B ′ )) = α ode k ode ((t, B), (t ′ , B ′ )) + α rbf k rbf ((t, B), (t ′ , B ′ ))where α ode and α rbf are parameters we can learn, e.g. using the marginal likelihood.Justification: We discuss limitations in the conclusion: we are restricted to a 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 a 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 on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a 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, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.• The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.• If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.• While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.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 a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.• Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.• We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets Question: 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?Answer: [Yes] Justification: All exisiting assests are cited and the corresponding licenses respected.Guidelines:• The answer NA means that the paper does not use 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 a URL.• The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.• If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset.• For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.• If this information is not available online, the authors are encouraged to reach out to the asset's creators.√t , where γ t is the information gain due to Srinivas et al.[71]. We are now ready to state the formal theorem along with its assumptions. 


Decision: 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.

Reviews: review: paper_summary: paper_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.
strength_weakness: strength_weakness: 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.

2. 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.

3. 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.

4. 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.

***References***

1. 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.
2. 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.
3. Astudillo, Raul, and Peter Frazier. "Bayesian optimization of function networks." Advances in Neural Information Processing Systems 34 (2021): 14463-14475.
questions: questions: 1. Could the author comment on the practice of limiting the extension of MDP-BO to either TS or UCB?
2. One minor question: in figure 3(a), is the legend "MDP-BO" written as "MDP-B0"?
3. 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.

4. 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. 

***References***

1. 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.
2. 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.
3. 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.
4. 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: limitations: Discussed above.
paper_summary: paper_summary: &nbsp;

The 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.

&nbsp;
strength_weakness: strength_weakness: &nbsp;

__MAJOR POINTS__

&nbsp;

1. 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.

2. The code would benefit from a README describing how to reproduce the results of the paper.

&nbsp;

__MINOR POINTS__

&nbsp;

1. It would be great if the references appeared in numbered order.

2. There is some missing capitalization in the references e.g. "bayesian" in place of "Bayesian".

3. Line 90, typo, "we seek to identify".

4. Line 566, typo, "due to selecting".

5. Line 566/567, missing full stop at the end of the equation.

6. On line 567 epsilon is homoscedastic but in the main paper the authors state that epsilon can be heteroscedastic?

7. Line 584, typo, "hypotheses".

8. Line 680, typo, extraneous "be".

9. 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.

10. 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?

11. 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.

12. In the proof of statement 25, it might be worth giving the acronym for SMW, namely Sherman-Morrison-Woodbury, somewhere in the text.

13. On line 136, have the authors defined S? This is the state space presumably.

14. In the related work on look-ahead BayesOpt, it may also be worth mentioning [1].

15. Line 179, typo, "a mix of the previous two".

16. In Section 2.2, the variable d could be explained in more detail.

17. When referencing Thompson sampling, it would be worth citing the original paper [2].

18. 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.

19. 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?

20. 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.

21. It would be great if the authors could provide the number of random trials over which the error bars were computed in Figure 4.

22. Line 622, "scales them properly". What is "them"?

23. Line 670, typo, "using".

&nbsp;

__REFERENCES__

&nbsp;

[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.

[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.

[3] VR Saltenis (1971). One Method of Multiextremum Optimization. Avtomatika i Vychislitel’naya Tekhnika (Automatic Control and Computer Sciences) 5(3):33–38.

[4] Garnett, R., 2023. Bayesian optimization. Cambridge University Press.

&nbsp;
questions: questions: &nbsp;

1. 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.

2. 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?

3. 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?

4. In Section 5, what was the motivation for the decision to fix the GP hyperparameters?

5. For Figures 2 and 3 why are there no error bars for the average regret curves?

&nbsp;
limitations: limitations: &nbsp;

1. 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.


&nbsp;
paper_summary: paper_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.
strength_weakness: strength_weakness: The paper may be improved by clarifying several questions below.
questions: 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?

2. Line 575: Please explain why the proprotionality holds. Monotonicity does not imply proportionality (same ratio).

3. Line 587: Please explain the two inequalities.

4. 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).

5. 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).

6. 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?
   
7. Would the authors please clarify equation (12) using the correct subscript of d (d_{t,h})?
   
8. Some minor remarks: In equation (11), should j be from 1 to t-1? Line 122: y(xh) should be y(xh,ah).
limitations: limitations: There are no potential negative societal impacts.
paper_summary: paper_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.
strength_weakness: strength_weakness: General Presentation: 

Generally, 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.
What 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.

Background Section: 

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.

Related Work (Section 2.4): 

I believe one of the most important related concepts, Bayesian Reinforcement Learning or Dual Control, is not discussed at all.

Acquisition Function Analysis: 

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?
questions: questions: - Is the utility function presented in the paper a contribution? Could the authors elaborate a bit on this?
- I am a bit confused about the notation in line 132. Could the authors please define $\tau_i$?
 - In line 136 and 199 is $S$ a switch in notation?
-	Could the authors elaborate a bit on line 209. What exactly is sampled using Thompson Sampling?
-	In line 259, why is this specific contstraint of 0.5 chosen?
-	I 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: limitations: -

Comments: title: Response to reviewer
Response to reviewer 2 / 2
Response to reviewer 1 / 2
Many Thanks to the Authors for their Rebuttal
comment: We disagree with the reviewer’s assessment, we do show that $\mathcal{U}(d_\pi)$ (as defined in the Appendix, and after fixing the typo in Eq. (8)) is proportional to Eq. (4). We stand by correctness of our work and the theoretical and empirical results we presented.

We would like to thank the reviewer for their time.
Thank you for the clarification. As I expected, $\mathcal{U}(d_\pi)$ is not proportional to the expression derived from Eq. (4) due to the noise term in $V$. This issue affects the main objective of the paper, and without knowing whether $1 - \alpha$ term is the same as $1 / TH$, I am unable to verify its correctness. Considering other errors in the proof and notations, my assessment of the paper remains unchanged.
We would like to thank the reviewer for their prompt response.

> The corrected proof resolves my concern

We are happy to hear that. Please consider raising your score, especially if you find the following clarifications useful.

> However, the explanation regarding...

The reviewer's intuition is correct, let us explain in detail what this "constant" is and how we use it (in essence the constant is the total length of the trajectories):

Firstly, we would like to clarify that empirical state visitations are also normalized. This can be seen, for example, in Eq.(11) where we introduced them.

Let us also note that the posterior variance of a Gaussian process does not depend depend on observed values, but only the points $x$, hence we can forward plan reduction in uncertainty. This is why Eq. (4) is equal to Eq. (6). Now, we show that the relaxation of the objective to state distributions holds in the simplest setting, for notation we will use $X$ to refer to the trajectory and $S$ to refer to the state space (in the paper we use $\mathcal{X}$):

For simplicity assume we have a single trajectory of length $H$ and no observation noise:
\begin{equation}
    V(X) = \sum_{x \in {X}} \Phi(X) \Phi(X)^T = H \sum_{x \in S} \hat{d_X} (x) \Phi(X) \Phi(X)^T:= H V(\hat{d_X})
\end{equation}
Now we note that:
\begin{equation}
   ||\Phi(z) - \Phi(z')|| \_{V(X)^{-1}} = (\Phi(z) - \Phi(z'))^T V(X)^{-1} (\Phi(z) - \Phi(z'))
\end{equation}
\begin{equation}
= \Phi(z)^T V(X)^{-1} \Phi(z) + \Phi(z')^T V(X)^{-1} \Phi(z') - 2 \Phi(z)^T V(X)^{-1} \Phi(z')
\end{equation}
\begin{equation}
= \text{Var}[f(z) | X] + \text{Var}[f(z') | X] - 2 \text{Cov}[f(z), f(z') | X] \\
\end{equation}
\begin{equation}
= \text{Var}[f(z) - f(z') | X]
\end{equation}
Note that since $V(X) \propto V(\hat{d}\_X)$, then:
\begin{equation*}
    \text{Var}[f(z) - f(z') | X] \propto ||\Phi(z) - \Phi(z')||\_{V(\hat{d}\_X)^{-1}}
\end{equation*}
so the optimization of the objectives is equivalent, note that the proportionality constant is exactly the normalization constant. In Lemma D.1 we formalize this more carefully in the setting with many trajectories and observation noise. We have to change the regularization of the objective appropriately (e.g. we must scale the contribution of the prior by $1 / H$, when planning for one episode, to counteract the normalization of the state-distributions), i.e., for the objectives to be equivalent we must use:
\begin{equation*}
    V(d_\pi) = \left(\sum_{x \in S} d_\pi(x) \frac{\Phi(x)^T \Phi(x)}{\sigma^2} + \frac{1}{H} I \right) = \frac{1}{\sigma^2}\left(\sum_{x \in S} d_\pi(x) \Phi(x)^T \Phi(x)+ \frac{\sigma ^2}{H} I \right).
\end{equation*}
For future reference, there are still two minor typos in the paper that we found, Eq. (8) it should read:
\begin{equation*}
    V(d_\pi) = \left(\sum_{x \in S} d_\pi(x,a) \frac{\Phi(x)^T \Phi(x)}{\sigma^2(x,a)} + \frac{1}{TH} I \right).
\end{equation*}
where we plan for $T$ episodes (trajectories) to make it consistent with the correct result in Lemma D.1. And in the proof of Lemma D.1, in the last equation of L705 there should not be a $TH$ before the sum. We wanted to manipulate $TH$ out but forgot it inside the inverse -- an honest typo. We will fix the typos and add a clarification that we need the correct regularization for the objectives to be equivalent.

To be fully clear, our implementation of the objective in the code does include the correct regularization procedure (see file mdpexplore/functionals/bandit\_functionals.py L242 where we multiply the identity by $1 - \alpha$).
Thank you for the clarification. The corrected proof resolves my concern. However, the explanation regarding $\mathcal{U}(d_{\pi})$, where $d_{\pi}$ is from the visitation polytope (normalized), and the original definition of $U(X)$ in Eq. (4)  and Lemma D.1 ($X$ is considered as "count") is still unclear to me.

As I mentioned in my previous comment, $U(X)$ represents the posterior variance of $f(z) - f(z') | X$. How do we then compute the posterior variance of $f(z) - f(z') | d_{\pi}$ (given a state visitation)? If I understand correctly, the expressions seem not equivalent. Specifically, $V(d_{\pi})$ involves a normalized $d_{\pi}$, while $V(X)$ involves an unnormalized $\hat{d}_X$. For them to be equivalent (or proportional for optimization purposes), we would need to multiply $d\_{\\pi}$ by some constant and check if $\\mathcal{U}(d\_{\\pi})$ is proportional to $U(d\_{\\pi} \\times \\text{constant})$ (such that $d\_{\\pi} \\times \\text{constant}$ is un-normalized), where $U(d\_{\\pi} \\times \\text{constant})$ is evaluated using Eq. (4).

Could the author provide further explanation on this or derive $\\mathcal{U}(d\_\pi)$ from Eq.(4)?
> From the author’s response, I now understand that what the author mean may be the empirical state distribution as described in Mutti et al. (2023). I believe the paper's notation should be revised to reflect this.

We happily agree, this will improve clarity.

> However, this also raises a concern...

There seems to be confusion between empirical visitation and a visitation as a variable. For simplicity, we will only talk about state visitations, but the more complete treatment with state-action visitations is analogous.

There is a $d_\pi$ a variable from the (state) visitation polytope, and then there is empirical visitation of the states after a policy has been executed $\hat{d}$. The notation we follow is that $\hat{d}$ are empirical ones. All objectives are functions of $d_\pi$ and we optimize over it, and it leads to convex programs in this variable. Sometimes the objective features $\hat{d}$ serving as an offset for the objective, not a variable. 

For experiment design we propose the optimization of Eq.(8) given by:
\begin{eqnarray}
      \mathcal{U}(d_\pi) = \max_{z, z' \in \mathcal{Z}} || \Phi(z) - \Phi(z')||^2_{V(d_\pi)^{-1}}  
\end{eqnarray}
where $V(d_\pi) = \left( \sum_{x \in \mathcal{X}} \frac{d_\pi(x)\Phi(x) \Phi(x)^T}{\sigma^2(x)} + I \right)$. This is an objective with a state distribution $d_\pi$ as a variable. Now assume we sample a trajectory $X$. In the objective introduced in Eq.(4), we showed in Lemma D.1 that we can equivalently write it as:
\begin{equation}
    U(X) = \max_{z, z' \in \mathcal{Z}} \text{Var}[f(z) - f(z') | X] = \max_{z, z' \in \mathcal{Z}} || \Phi(z) - \Phi(z')||^2_{V(X)^{-1}}  
\end{equation}

where $V(X) = \left( \sum_{x \in \mathcal{X}} \frac{\hat{d}_X(x)\Phi(x) \Phi(x)^T}{\sigma^2} + I \right)$. In particular, note that we can write $V(X)$ (and hence $U(X)$) solely as a function of the state-distribution of the trajectory $\hat{d}_X$. Most importantly, the value is exactly the same as the objective in Eq.(8) where we evaluated $\mathcal{U}(\hat{d}_X)$, i.e. setting d_pi to the empirical distribution of trajectory $X$. As the reviewer points out, in this case, not every state distribution translates to a trajectory, but this is not an issue for our objective.

Secondly, when we optimize over state-action polytope $d$, we are in essence doing a continuous relaxation of a discrete optimization over empirical visitations with a fixed $T$. Once we have an optimal state density, we are able to build a policy from it. And once we have a policy we can deploy it in the environment to obtain a new trajectory. If it helps, it is similar to solving an integer optimization problem by considering a continuous relaxation of a function. The "rounding" in the relaxation proceeds via sampling from the policy for our case. It is an analog of convex relaxations in the experiment design literatures. For more examples, see the book on convex optimization (Boyd and Vandenberghe , 2003) or for something more relevant to our work, Mutny et. al. (2023), where they study different methods of rounding.
Thank you the response and for carefully analyzing our derivation. We want to again mention that the objective used in our paper has been used in works of Fiesz et. al. (2019) published at NeurIPS. We hope there is no doubt that this is a meaningful objective. We only provide a Bayesian motivation for it. 

Upon a more meticulous look of our derivation, we found we had an unnecessarily complicated path that led to the sign blunder and wrong fraction order that the reviewer points out (which ended up cancelling each other out). Allow us instead to provide a more elegant and simple derivation (which will be added to the revised paper). It follows from the Gaussian tail bound that states for any $X \sim N(\mu, \sigma^2)$ and $t > 0$ then:
$$ \mathbb{P}(X \geq \mu + t) \leq e^{-\frac{t^2}{2 \sigma^2}} $$
We can now bound the probability of making an error, indeed $\mu_T(z) - \mu_T(z^*) \sim N(a_z, b_z^2)$:
$$ \mathbb{P}(\mu_T(z) - \mu_T(z^*) \geq 0) = \mathbb{P}(\mu_T(z) - \mu_T(z^*) \geq a_z + (-a_z) ) \leq e^{-\frac{a^2}{2 b^2}}$$
with the caveat that the bound only holds asymptotically as $T \rightarrow \infty$, as in this regime $a_z \approx f(z) - f(x_f^*)$, due to asymptotic argument we made previously. Notice that by definition of $x_f^*$, $f(z) - f(x_f^*)\leq 0$, hence substituting for $t = -a_z$,
\begin{equation}
    E_{f\sim GP}[\log P(\mu_T(z) - \mu_T(x_f^\star) \geq 0 | f)] \leq - \frac{1}{2} E_{f\sim GP}\left[\frac{a_z^2}{b_z^2}\right]
\end{equation}
by taking first the log, and second the expectation on both sides. From here the arguments follow in a similar fashion to the original paper, with the caveat that originally we had an error sign and the numerator and denominator the wrong way around and so the errors cancelled out. For full clarity here is the argument:

Our aim is now to minimize making an error, so minimizing $-E_{f \sim GP}\left[\frac{a_z^2}{b_z^2}\right]$. It then follows that in the large $T$ limit (from the previously discussed bounds on $b_z$ and $a_z \approx f(z) - f(x_f^\star)$): 
\begin{equation}
    - E_{f\sim GP}\left[\frac{a_z^2}{b_z^2}\right] \leq - E_{f\sim GP} \left[ \frac{(f(z) - f(x_f^\star))^2}{k_{\mathbf{X}}(z, x_f^\star)} \right]
\end{equation}
combining the two equations we obtain:
\begin{equation}
    E_{f\sim GP}[\log P(\mu_T(z) - \mu_T(x_f^\star) \geq 0 | f)] \leq - \frac{1}{2} E_{f\sim GP} \left[ \frac{(f(z) - f(x_f^\star))^2}{k_{\mathbf{X}}(z, x_f^\star)} \right]
\end{equation}
note the correct interpretation of the bound: the probability of an error is small if (i) the uncertainty is small, or (ii) the difference in values is large. Then we simply follow the arguments in Section C.3 to arrive to the objective:
\begin{equation}
    - E_{f\sim GP} \left[ \frac{(f(z) - f(x_f^\star))^2}{k_{\mathbf{X}}(z, x_f^\star)} \right] \leq - \frac{E_{f\sim GP}[\text{gap(f)}]}{\text{Var}[f(z) - f(z') | \mathbf{X}]}
\end{equation}
And finally, we note that due to the negative sign, and the denominator being constant, minimizing the RHS is equivalent to minimizing the objective introduced in Eq.(4), i.e. $\arg\min_{x}-g(x)$  is equivalent to $\arg\min_{x} \frac{1}{g(x)}$ when $g(x) > 0$.
Thank you for your clarification. I still have a few questions, and I would appreciate your assistance in understanding the following points.

**On the derivation of the objective**
Let us go through your derivation step by step.

$P(\mu_T(z) - \mu_T(x^*_f) \ge 0|f) \propto \exp (- \frac{a_z^2}{b_z^2}) \frac{1}{b_z} \Rightarrow P(\mu_T(z) - \mu_T(x^*_f) \ge 0|f) = C_1 \exp (- \frac{a_z^2}{b_z^2}) \frac{1}{b_z}$

for some constant $C_1 \neq 0$. By taking the logarithm of both sides,

$\log P(\mu_T(z) - \mu_T(x^*_f) \ge 0|f) = \log C_1 - \frac{a_z^2}{b_z^2} + \log \frac{1}{b_z} = \log C_1 - \frac{a_z^2}{b_z^2} - \log b_z$

Then, how does it imply

$\log P(\mu_T(z) - \mu_T(x^*_f) \ge 0|f) \propto \frac{b_z^2}{a_z^2} + \log b_z$

and

$E_{f \sim GP} [...] \le \frac{1}{C} E_{f \sim GP} [\frac{b_z^2}{a_z^2}]$ (note that we still have a constant $\log C_1$)


**On expected trajectory**
I have just skimmed through the work of Mutti et al. (2023) mentioned in the author response, so please correct me if I'm wrong. From what I read, they first review the original RL formulation that maximizes the expected total reward (not the reward of the expected trajectories) (Equation (1) in Mutti et al. 2023), which can be rewritten in terms of the state distribution induced by the policy. They then discuss a generalization of the RL formulation (Equation (2) in Mutti et al. 2023) by directly introducing a utility function based on the state distribution induced by a policy. It is important to note that this utility function is defined for a state distribution (or its empirical distribution), rather than for a specific trajectory (the expected trajectory).

From the author’s response, I now understand that what the author mean may be the empirical state distribution as described in Mutti et al. (2023). I believe the paper's notation should be revised to reflect this.

However, this also raises a concern: As the author responded,

$U(d_\pi) = U(\sum_{\tau} \eta_{\pi}(\tau) \delta_\tau)$

In Equation (4), $U$ can only take trajectories (or a set of inputs) as it represents the posterior variance of $f(z) - f(z')$ given the observations. When considering $U(\sum_{\tau} \eta_{\pi}(\tau) \delta_\tau)$, it is unclear how to evaluate $U$ for $\sum_{\tau} \eta_{\pi}(\tau) \delta_\tau$​, which is a weighted average (empirical distribution) of the inputs. This raises the question of how to evaluate the posterior variance of $f$ given an empirical distribution of the inputs. Furthermore, Lemma D.1 also pertains to the posterior variance given a set of inputs $X$, rather than an empirical distribution of the inputs.
Reviewer mJAD, could you please respond do the additional comments made by the authors? Currently there is a rather large discrepancy between your score and the other reviewers, and your review didn't point to any major weaknesses other than clarity (given your questions). I would like to understand whether the authors' additional responses helped clarify your confusion. Thanks.
&nbsp;

Many thanks to the authors for their clarifications. For figures 2 and 3 it would have been nice to see the errorbars in the 1-page rebuttal document that could have been attached to the rebuttal comment.

The major point of clarification stems from the problem settings under which the authors' method can be applied relative to LSR and SnAKe.  Given that this point has been fully addressed, I am raising my score.

&nbsp;
**On eq. (7)**

Our final objective will be a function of all executed trajectories. At each episode $t$, we deploy a policy $\pi_t$. $\pi$ refers to the aggregation of all policies $\pi: \\{1, ..., T\\} \times X \times A \rightarrow [0, 1]$.

**Q: Why not take the expectation of the utility over the trajectories ... how should we conceptualize the "expected trajectory"?**

Taking the expectation outside as $\mathbb{E}_{\tau\sim \pi}[F(\tau)]$ is intractable in nearly all cases, see Mutti et. al. (2023) for discussion in the context of reinforcement learning. Hence, all works studying remotely similar problem in formal fashion (providing tractable certificates) focus on the case where expectation is inside. 

By considering utilities of the expectation, we make the problem tractable (not NP-hard; even not NP-hard to approximate as in Prajapat et. al. (2024)). This follows the classical approach in the field of experiment design (Fedorov, 1997), convex reinforcement learning (Zahavy et. al., 2019) and prior work of Mutny et al. (2023), where the expectation is inside. 

The expected trajectory is not well conceptualized in our manuscript. As we state in L139, the concept borrows on the "expected trajectory" from Mutny et. al. (2023) on which this work builds up. We do not utilize the formalization of this anywhere later in the paper, hence we felt no need to reproduce it. 

Let us summarize here what we mean by $U$ (non-caligraphic), and $U(\mathbb{E}[X])$ in this context. The conceptualization begins by introducing the space of probability distributions over trajectories $\mathcal{P}$, here, $\eta_\pi(\tau) \in \mathcal{P}$ is a distribution describing the probability of obtaining trajectory $\tau$ upon deployment of policy $\pi$. Each policy $\pi$ induces $\eta_\pi(\tau)$. Similarly any finite execution of policy $\pi$ (say $n$ times) leading to observed trajectories $\mathbf{X}$ as in  Eq. (4) induces an empirical distribution $\hat\eta_\pi(\tau)$; this is the distribution where all our *observed* trajectories are given equal mass. The function $U$ is a function over probability distributions over trajectories $U:\mathcal{P}\rightarrow \mathbb{R}$, and for **all** empirical measures they coincide with $U(\mathbf{X}) =U(\hat\eta_\pi)$ from Eq. (4). Expectation of the empirical measure leads to $\mathbf{E}[\hat\eta_\pi] = \eta_\pi$, and this is the way to think about it. In particular, for a single episode, we write:
$$
\mathcal{U}(d_\pi) \stackrel{!} = U(\mathbb{E}[\mathbf{X}]) := U\left(E_{\tau \sim \pi} \left[ \delta_{\tau}\right]\right) = U \left(\sum_{\tau} \eta_{\pi}(\tau) \delta_{\tau} \right)
$$

Note that the sum is over all possible trajectories, making it intractable. For our objective, however, we can show it can be equivalently written in the variables of space of state-action distributions (!), we only have to sum over all states and actions which makes solving the problem solvable ($(|\mathcal{X}||\mathcal{A}|)^H$ to $|\mathcal{X}||\mathcal{A}|H$ variables). For completeness, for more episodes:
$$
\mathcal{U}(d_\pi) = U\left(E_{\tau_1 \sim \pi_1, ..., \tau_t \sim \pi_t}\left[\frac{1}{t}\sum_{i=1}^t \delta_{\tau_i}\right]\right) = U \left(\frac{1}{t} \sum_{i=1}^t\sum_{\tau} \eta_{\pi_i}(\tau) \delta_{\tau} \right)
$$
We will refer the readers to Appendix for the clarification of these concepts, but the most natural way to eliminate the confusion is most likely to introduce expected empirical visitations $\mathbb{E}[\hat{d}(x,a)]$ instead of $\mathbb{E}[\mathbf{X}]$. We have used $\mathbb{E}[\mathbf{X}]$ to make the connection to prior works, but its most likely causing confusion.

References:

- Mutti e.t al. (2023), Convex Reinforcement Learning in Finite Trials, JMLR

- Zahavy et. al. (2021), Reward is enough for convex MDPs, NeurIPS 2021

- Fedorov (1997), Model-Oriented Design of Experiments, Springer

- Prajapat et. al. (2024), Submodular Reinforcement Learning, ICLR 2024
Thanks to the reviewer for reading our responses carefully, and bringing up the following points of discussion. We will split our answer in two parts, first carefully providing details of the derivation, and secondly clarifications on the notation and concepts of eq. (7).

**On the derivation of the objective:**

Based on the discussion, L575 now reads:
$$ P(\mu_T(z) - \mu_T(x_f^\star) \geq 0 | f) \propto \exp\left(-\frac{a_{z}^2}{b_{z}^2}\right)\frac{1}{b_{z}}\implies \log  P(\mu_T(z) - \mu_T(x_f^\star) \geq 0 | f) \propto \frac{b_z^2}{a_z^2} + \log(b_{z}) $$
and L578:
$$ E_{f\sim GP}[\log P(\mu_T(z) - \mu_T(x_f^\star) \geq 0 | f)] \leq \frac{1}{C} E_{f\sim GP}\left[\frac{b_z^2}{a_z^2}\right] $$
where $C$ is a constant of proportionality, and the upper bound holds asymptotically as $T \rightarrow \infty$ or in fact for a sufficiently large $T$, (where $\log b_z \leq 0$ as long as $b_z \leq 1$, which has to eventually happen at finite $T$).

Then the following lines (587) will read:
\begin{eqnarray*}
a_z^2 & = & (\theta^\top(V_T+I_{H})^{-1} V_T (\Phi(z) - \Phi(x^\star_f)))^2 \\ & \stackrel{T\gg 0}\approx & (\theta^\top (\Phi(z) - \Phi(x^\star_f)))^2 = (f(z) -  f(x^\star_f))^2\\
\end{eqnarray*}
which holds in the asymptotic regime as $(V_T + I_H)^{-1}V_T \approx I$. In fact, we can give order dependence for the specific instance where $V_T$ is made up from unit vectors. In this case, $a_z^2 = (f(z) -  f(x^\star_f))^2(1 - \mathcal{O}\left(\frac{1}{T}\right))$, so the gap rapidly goes down.

The rest of the derivation remains the same. As mentioned in our general rebuttal, the objective we derived has been studied before. We chose the objective because of its asymptotic optimality properties and good empirical performance. We developed a new, Bayesian, motivation for this existing objective rather than provide a formalized theorem. The derivation's reliance on asymptotic arguments ($T\rightarrow\infty$) is in line with the frequentist studies (Fiez et al. 2019).

We agree that the approximation quality remains unclear for small $T$. Claiming a "tight upper bound" is perhaps overselling, by this we meant that there is an instance where the above holds with equality. This is in fact true for any instance in the large $T$ limit.  We will add this note to the camera-ready, and remove the words "tight upper bound". However note, we do not use this upper bound for the algorithm. We clearly state we make yet another upper bound which is perfectly formal, and then use this in the algorithm. This second gap we formally analyze in Appendix C.5.
Reviewer J8j4, could you please review the authors' response and see whether it addresses your open questions? Please acknowledge having done so in a comment together with a discussion of whether / how it has affected your assessment of the submission. Thanks.
I appreciate the authors' response. I maintain my original evaluation that the setting and algorithm are novel and valid. Though my fellow reviewers have raised clarity issues, I don't foresee significant soundness problems. Hence, I'll keep my score.
Thank you for your response. However, I remain unclear about some of the notations and derivations presented in the paper.

After reviewing the suggested lines 574-578, it appears that there is some confusion between proportionality and monotonicity within the proof. It occurs in several places, making it challenging to verify the correctness of the arguments, for instance, the expression following line 575, the subsequent omission of the log term, and the expectation taken over $f$ after line 578. Even if we replace "proportional to" with "approximate to," the quality of this approximation remains ambiguous. Moreover, the proof of Equation (14) is intended to establish a "tight upper bound" (as noted in line 94), rather than an approximation.

The author posits that: A trajectory is a random set, and therefore we must define a utility with respect to the expected trajectory when using a policy $\pi$." However, this raises the question: why not take the expectation of the utility over the trajectories, similar to reinforcement learning approaches that maximize expected total reward over trajectories rather than the reward of the expected trajectory? Additionally, in the context of a discrete state space, how should we conceptualize the "expected trajectory"?

Regarding Equation (7), $\pi$ is still unclear to me, does it refer to $\pi_1$​ or $\pi_2$ or something else​?
I have read the rebutal and thank the authors for their answers. I will keep my score and still vote for acceptance of the paper.

I want to mention one minor clarification that non-convex MPC is far from being exotic. If you design a controller for a non-linear dynamical system, lets say for a toy example like a cart-pole, the MPC problem is solved using a direct method. When the dynamics are not linearized, this leads to a non-linear dynamics constraints, which is non-convex.