BayeSQP: Bayesian Optimization through Sequential Quadratic Programming
Authors: Paul Brunzema
Sebastian Trimpe
Venue: NeurIPS 2025 spotlight ()
Keywords: Bayesian optimization
sequential quadratic programming
constrained optimization
Primary area: probabilistic_methods
TL;DR: We introduce BayeSQP, a novel black-box optimization algorithm that combines sequential quadratic programming with Bayesian optimization for high-dimensional constrained problems..

Abstract: We introduce BayeSQP, a novel algorithm for general black-box optimization that merges the structure of sequential quadratic programming with concepts from Bayesian optimization. BayeSQP employs second-order Gaussian process surrogates for both the objective and constraints to jointly model the function values, gradients, and Hessian from only zero-order information. At each iteration, a local subproblem is constructed using the GP posterior estimates and solved to obtain a search direction. Crucially, the formulation of the subproblem explicitly incorporates uncertainty in both the function and derivative estimates, resulting in a tractable second-order cone program for high probability improvements under model uncertainty. A subsequent one-dimensional line search via constrained Thompson sampling selects the next evaluation point. Empirical results show that BayeSQP outperforms state-of-the-art methods in specific high-dimensional settings. Our algorithm offers a principled and flexible framework that bridges classical optimization techniques with modern approaches to black-box optimization.

BayeSQP: Bayesian Optimization through Sequential Quadratic Programming
Paul Brunzema Sebastian Trimpe
Institute for Data Science in Mechanical Engineering
RWTH Aachen University
Aachen, Germany
(brunzema, tr trimpe (ddmppeeeedde...e wwth wwth-aachen..nnmmmmm
Abstract
We introduce BayeSQP, a novel algorithm for general black-box optimization that merges the structure of sequential quadratic programming with concepts from Bayesian optimization. BayeSQP employs second-order Gaussian process surrogates for both the objective and constraints to jointly model the function values, gradients, and Hessian from only zero-order information. At each iteration, local subproblem is constructed using the GP posterior estimates and solved obtain a search direction. Crucially, the formulation of the subproblem explicitly incorporates uncertainty in both the function and derivative estimates, resulting a tractable second-order cone program for high probability improvements under model uncertainty. A subsequent one-dimensional line search via constrained Thompson sampling selects the next evaluation point. Empirical results show that BayeSQP outperforms tate-offthee-a methods in specific high- bingh-iimnsiioall settings. Our algorithm offers a principled and flexible framework that bridges classical optimization techniques with modern approaches to black-box optimization.
Introduction
In recent years, Bayesian optimization (BO) has emerged as a powerful framework for black-box optimization ranging from applications in robotics [9, 6, 44] to hyperparameter tuning [55, 12] and drug discovery [21, 40, 10]. To address high-dimensional problems emerging in these fields, variety
of digh-dimensional BO (HDBO) approaches have been proposed, including the use of local BO (LBO) methods [15, 43] or methods that exploit specific structure in the objective [13]. Recently, growing debate has emerged over whether such HDBO methods are truly necessary, given that appropriate scaling of the prior can already yield strong performance on certain high-dimensional benchmarks [30, 67]. However, as shown by Papenmeier et al. [48], these approaches solve numerical issues in the hyperparameter optimization of the Gaussian process (GP) surrogate but their success can still be attributed to emerging local search behavior. We argue that it is not matter of choosing either HDBO approaches or standard approaches, but rather of leveraging recent advances how to achieve numerical stability also for HDBO methods.
Building on this perspective, we aim to integrate the strengths of established classical optimization techniques within the HDBO framework. Specifically, we extend the widely-adopted local method for HDBO GIBO [43, 45, 65, 16, 23ich can be interpreted as combining BO with first-order optimization methods to LBO with second-order methods. We introduce BayeSQP, novel algorithm for black-box optimization that merges the structure of sequential quadratic programming (SQP) with concepts from BO. BayeSQP employs GP surrogates for both the objective and constraints that jointly model the function values, gradients, and Hessians from only zero-order information (Figure I). At each iteration, local subproblem is constructed using the GP posterior estimates and solved to
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

BayeSQP
build sub-problem
solve sub-problem
Sequential Quadratic
Fie Vile Joil
Bayesian
Optimization
learn surrogale
models
Thompson sampling
as line scarch
Figure 1: Overview. BayeSQP combines ideas from sequential quadratic programming and Bayesian optimization for efficient high-dimensional black-box optimization.
obtain a search direction. Through constrained Thompson sampling, we select the point for the next iteration. In summary, the key contributions of this paper are:
C1 A novel algorithm BayeSQP leveraging GP surrogates to utilize the structure of classic SQP within BO for efficient high-dimensional black-box optimization with constraints.
C2 An uncertainty-aware subproblem for BayeSQP that accounts for the variance and covariance in function and gradient estimates, resulting a tractable second-orde cone program.
Empirical experiments demonstrating that BayeSQP outperforms state-of-the-art BO methods in specific high- ligh-dimensional constrained settings.
2 Problem formulation
consider the problem of finding an optimizer to the general non-convex optimization problem
arg min f(x ØEX
subject to cxxx) > O, Vi E Im == (1....,m. m)
(I)
where x R and constraints
X - R for all i E Im are black-box functions defined over the compact set X C Rd. At each iteration E Ir where T the total budget for the optimization, an algorithm selects query point E X and receives noisy zeroth-order feedback following the standard observation model in BO as f: = fax+) + e for the objective, and C₁ = es (xt)) Ec, for all E I for the constraints where f and Ec; are independent realizations from zero- mean Gaussian distribution with possibly different noise variances. From these observations, we construct independent datasets for the objective function Dj
and for each constraint Dr {( ( for all Im. which any zero-order method can leverage solve (1).
3 Preliminaries
3.1 Sequential quadratic programming
SQP represents powerful framework for solving nonlinear constrained optimization problems by iteratively solving quadratic subproblems. This method has become one of the most effective techniques for handling wide range of optimization problems. The foundation of constrained opti- mization rests on the Lagrangian function. defined as ((x,{) = f(x) Σn (xc((z). It combines the objective with the constraints, where each constraint is weighted by its Lagrange multiplier & Solving the optimization problem involves satisfying the Karush- Kuhn- Tucker conditions for this Lagrangian. To achieve this, at each iteration t, SQP constructs a quadratic approximation of the Lagrangian using the Hessian H. V2 œm L(x,, {) (or an appropriate approximation thereof) and linearizes the constraints around the current point 2;. This generates the following subproblem:
p: = arg min peRd
5P H.P Vf(x)) T + f(x))
(2)
subject
Va(x;) p -c;(x,), E 1m

---

The solution Pt provides search direction, and the next iterate is typically computed as a++1 2; + a,P: where a, a step size determined by an appropriate line search procedure that ensures adequate progress toward the optimum. For an overview of classical SQP methods, see [46].
Under standard assumptions, SQP exhibits local superlinear convergence when using exact Hessian information, and various quasi-Newton approximation schemes (such as BFGS or SRI updates, cf. [46]) can maintain good convergence properties while reducing computational overhead. This fast local convergence also makes it interesting for HDBO. However, the key challenge here is that usually only zero-order information on the objective and constraints is available.
3.2 Gaussian processes
GPs are a powerful and flexible framework for modeling functions in a non non-parametric way. A GP defined collection of random variables, any finite number of which have joint Gaussian distribution. Formally, a GP defined by its mean function m(x) = [=)] and kernel k( x') = Covl ( f(x')] [53]. In BayeSQP, we GPs to model the objective function f and the constraints as standard in BO [18]. Contrary to standard BO, we aim to leverage the following property of GPs: They closed under linear operations, i.e., the derivative of GP is again GP given that the kernel is sufficiently smooth [53]. This enables us to derive a distribution for the gradient and Hessian. We formulate the following joint prior distribution:
y
m(X)
Th (X + I21
(())
k(x, X) Vk(x, X) ()
k(2, x) Vk(a, x) (2((xx))
N
(3)
Vm
V2k(o, V*((o, V*k(z, x)
V m
where E IR are the n function observations, X = [i....., xnl gdxa is the matrix of all training inputs with æ; E Rd and æ E Ie is the test point. Here and the following, we use to denote symmetric entries for improved readability. The joint conditional distribution is then:"
()(x)) ( [a)))(x)
03(()) Ev)((x) (v)),(())
))
x. N,
Evy(a) Evs,,yy())
V2
Exe((x)
Following standard conditioning of multivariate normal distributions, we can directly compute the mean and covariance functions of the marginals of the posterior
(x(x) m(x) + k X) (y m(X)) ER,
Marginal GP of f:
(5)
07(x)) k(x, k( K --1K----x
ER
( = Vm(()) Vk (v. X - (1 m(X))
E Rd.
Marginal GP of Vf:
(6)
(x)(x) = V2k(z, Vh V VK(X,
E Rdxd
Here, we defined the Gram matrix as K = k( X) o with entries [K(X, x)l; k( xj) for i.j E In, and use the notation that k(X, x) E Rnx1 is the vector of kernel evaluations between each training point and the test point, with entries [k:(X, .));; (+(x,=))
k. Similarly, can obtain the covariance term E (x) as well as the mean estimate of the Hessian. As noted in Müller al. [43], must perform the inversion of the Gram matrix K only once. So, while calculating the gradient distribution and the mean of the Hessian is not for free, the additional computational overhead is limited with increasing data size.
3.3 Related work
Scalable Bayesian optimization For long, BO has been considered challenging for high- dimensional input spaces leading to the development of tailored algorithms for this setting. Such
We write the joint distribution over f. V and V in block matrix form convey intuition, though this abuse of notation. Formally, all components are vectorized and stacked into single multivariate normal vector. Specifically, we have [g. V veo (v f) T]
g=+1+d+d" Similar for the mean covariance. FFrom hereon, we will only consider the mean of the Hessian as storing the variance over all terms well as all covariances is very computationally intensive: Let 2 f(x) dddd then Cov f(x)] E dxdxd
3

---

approaches include LBO methods, which we will discuss in more detail in the following, as well as methods that aim to leverage a potential underlying structure or lower dimensional effective dimensionality [33, 62, 13, 47]. Recent results show that some of the core challenges in this high- dimensional setting are due to a numerical issues when optimizing the hyperparameters, which can be in part addressed by enforcing larger lengthscales [30, 67, 48). These developments do not make scalable approaches obsolete. Rather, we see them as a tool to further improve the modeling also for scalable BO approaches. To address scalability, in the sense of scaling with data, alternative surrogates for BO such as neural networks [56, 35, 8] or sparse GPs [42, 41] have been discussed; addressing these scalability issues, however, is not the focus of this work.
Local Bayesian optimization LBO methods aim to improve the efficiency of the optimization process by focusing on local regions of the search space. Approaches such as TuRBO [15] and SCBO [14] can be classified as pseudo-local methods: their trust-region approach still allows for the exploration of multiple local areas and only over time collapses to one local region. On the contrary, Müller et al. [43] introduced with GIBO a new paradigm of LBO combining gradient-based approaches with BO. Since then, the algorithm has been modified with different acquisition functions to actively learn the gradient [45, 58, 23, 16], theoretically investigated [65], and extended with crash constraints [61]. This class of algorithms operates fully locally. Our algorithm BayeSQP can also be classified as such a local method. In this sense, BayeSQP extends GIBO to second-order optimization by using a Hessian approximation from GP. Similar ideas have been leveraged in a quasi-Newton methods [11]. However, by incorporating ideas from SQP, BayeSQP is directly applicable to both unconstrained and constrained optimization problems-something which is not possible with GIBO.
Bayesian optimization and Gaussian processes in classical optimization There have been various papers integrating BO with first-order optimization, e.g., for line search [38, 57]. GPs have been successfully applied and leveraged in optimization-bott for local optimization [25, 24] and global optimization (essentially BO) [32, 18]. All of these can be classified numerics [27, 28J. Similar to our approach, Gramacy et al. [20] merged classical methods with BO by lifting the constraints into the objective using an augmented Lagrangian approach which later got extended to slacked [49] and recently a relaxed version [4]. These approaches are based on expected improvement (EI) and, crucially, Eriksson and Poloczek [14] showed that these approaches do not scale well to high-dimensional problems, BQyeSQP differs in the type of acquisition function for the line search as well as the framework as it builds on SQP. To our knowledge, we are the first to leverage a joint GP model of the function, its gradient and Hessian in a classical framework.
subfield of probabilistic
4 BayeSQP: Merging classic SQP and Bayesian optimization
This paper proposes the LBO approach BayeSQP. As described above, the main objects of this approach are GP models of the objective and possible constraints that jointly model the function value, the gradient as a well as the Hessian in single model. BayeSQp then leverages this model at each iteration construct quadratic uncertainty- improvement with high probability. In the following. we will first discuss our modeling approach. Based on this, we will construct the subproblem, followed by discussion line search. In the end, we touch on further practical extensions and give intuition on the optimization behavior.
subproblem for a search direction that yields
4.1 Second-order Gaussian processes as surrogate models for BayeSQP
In BayeSQP, we aim to leverage ideas from both SQP and BO to solve constrained black-box optimization problems as in (1). For this, we will model the objective and all constraints using second-order GP models introduced in Section 3.2 here stated for the objective:
xc, x x N (I
(x(x)) / (x)
oz(x) Ev)((x) x
[vec(Hp>)(=)) (
Exy(æ)
X
(7)
vec V
We use surrogate models of the same form for each constraint (x()). We do not compute the covariance of the Hessian (x) due the scaling issues with dimensions discussed in Section 3.2. Figure 2 demonstrates the effectiveness of such a joint GP model. We can estimate the gradient, identify local optima, and estimate curvature all from only zeroth-orde information.

---

Zero-order prediction
Gradient field
Hessian determinant
Training data
Critical point boundary
Test points
P3
P3
P3
P1
P1
20
TD
Mean curvalure and estinates
P2
P3
f(z)
1.003
f(x)
0.073
(())
-0965
0.85
f 33.15
-0.10
f
f
0.20
-2.61
0.26
A1
22.079
À 0.316
À1 17.689
17.391
3.567
21.351
Figure 2: The power of Gaussian processes. Although we only have zeroth-order information about the function, the differentiability of the GP allows us to estimate both the gradient and curvature. All estimates provided in expectation; the associated uncertainties are not shown.
Crucially, it not required to always evaluate the full posterior distribution for each test point. In SQP framework, approximate the Hessian of the Lagrangian once at our current iterate
H. (xx)) EE <(t
Mr (xx),
where e((-1) are the Lagrange multipliers from the solution the last subproblem, but for the subsequent line search, we directly work with the cheap marginal GP (N(r((z),
4.2 Deriving the subproblem for BayeSQP
Standard SQP approaches typically require exact knowledge of the objective function, constraints, and their respective gradients. our case, we only have access zero-order feedback and the question arises how to formulate a suitable subproblem given our choice of surrogate model.
Expected value SQP subproblem A straightforward approach is simply formulate a subproblem using expectations, leading to the following expected value subproblem:
Pt arg min PERe
E
H.p V (ax.) T f(x)) P
subject to E (c(x)) Vc,(x,) pl 2 0, Vi Im
While intuitive, this formulation fails to account for the inherent uncertainty in the estimates. discussed by Nguyen et al. [45] and He et al. [23]. taking into account the uncertainty of, e.g., the gradient, can crucial for improving with high probability for LBO approaches.
Uncertainty-aware SQP subproblem To address this limitation, we reformulate the standard QP subproblem into a robust version that explicitly accounts for uncertainty in the estimates:
1 1 H.P (x(x,) p f(x))
E arg min pERa
VaR -,
Objective value-at-risk with confidence level 1-6,
(10)
subject to
(c,(2))) + (.
2 1 8e, Vi E T
Constraint satisfaction with confidence 6_
This formulation accounts for uncertainty through two mechanisms: employing value-at-risk (VaR) for the objective function and enforcing probabilistic fcasibility for the constraints. The resulting
5

---

search direction minimizes the worst-case objective value while ensuring the constraints are satisfied with high probability.
Tractability through joint Gaussian process The robust formulation in (10) remains intractable without distributional assumptions. By modeling the objective and constraints jointly Gaussian with their gradients, we can transform (10) into a deterministic second-order cone program. Next, we derive this tractable reformulation for the constraints; the objective follows analogously.
For the constraints, we aim to ensure that P (c:(x))) + Vei(2;, T 2 0) 2 1- 8.. Since we have v ~ M( u, v v) for a multivariate Gaussian random variable z ~ M(u,, Σ₂) and v is deterministic vector, we know that (( + (2, T is also normal distributed with moments
=T
E (x)) + Vc((x) P = Me, (x:) + M (It) Var [<x (x,) p = - (a,))
(11)
( 2p'
(12)
where the last term accounts for the covariance between the function and its gradient. In the following, we drop the explicit evaluation æ; for all moments for notational convenience, i.e., Me; Me, (2t)
For Gaussian random variable to remain non-negative with probability at least - 8, we require its mean to exceed its standard deviation multiplied by the corresponding quantile. This yields:
T Me. + eve: Φ1 - 8 oZ. + p + 2p Vo,Ver
(13)
where 91-5 = Φ Φ- 8) denotes the (1 8) quantiie of the standard normal distribution. Rearrang- ing the terms and introducing an auxiliary variable fe; upper-bound the square root term allows the constraint to be reformulated as set of two inequalities:
uve. p + 91-sbe; S Me: and
V Voa. + +p Eve, p Ee, S be:
(14)
To express the square root term more compactly, we consider the full covariance matrix associated with the joint Gaussian distribution of c, and its gradient Vc, Specifically, we can state
Le +p Eve: p +p2p" 2p Ee, I
(15)
[
Ee,
Eve
P
By Cholesky decomposition of the covariance matrix, can express the square term second-order cone constraint:
V o + + p 2p Dc.
L
c2
p
=
L c2
P
S
be:
(16)
Using the same reasoning, we can reformulate the objective function by introducing the auxiliary variable by. In the end, we obtain the following formulation that we refer to as B-SUB.
The uncertainty-awan subproblem of BayeSQP
(B-SUB)
Pt E
arg min p.by, Jeetm
5P H.p + M fP + Ms + t1 bz
(17)
subject to
LJ
S by.
L
$ be.
Vi E Im:
P
P
2
Mvv.P P 91-8. b.b.s be, S Meg
Vi E Lm-
where b and L., are Cholesky factorizations as
of E EVi. Evi L.. LT
o2 Evars
Vi E Im.
(18)
E
We omitted the explicit dependency on x: for clarity but all moments are evaluated at Tt.
This formulation also naturally incorporates the subproblem formulation in (9).
Corollary (Recovering the expected value formulation). The solution for the search direction of B-SUB is equivalent to solution of (9) for on 0.5 and & 0.5.
(Proof in Appendix C) Remark practice, the numerical solver will have an influence on the obtained results. So while the cones no longer restrict the search direction, cone solver might still return different solution.

---

Query points
Feasible
Best
Direction ft
Path
Pay f = 0.5
QQP
Qa
C-logEI
f 0.2
6y 0.05
0.05
c(x)
c(x) <0
c(2) 0
c(w) < 0
300
To
Figure 3: Infuition on optimization behavior of BayeSQP. Disregarding uncertainty (6,, 6 0.5, left)
results in directions tangential to the circular constraint, while a for a very conservative configuration (6y, 6, = 0.05, center right), the constraint acts repellent to ensure feasibility. Values in between (center left) yield a desirable convergence path to the optimum. On the right, we see the space-filling behavior of constrained logEl which is fundamentally different compared to the local BayeSQP.
4.3 Line search through constrained posterior sampling
With the search direction given as the solution of the B-SUB subproblem, the next step is to decide on a step size a which which we can update the current iterate as 2t+11 = 2; + a.Pr. To implicitly decide on the step size, perform constrained posterior sampling [14] on the one-dimensional line segment spanned by pt- Specifically, we aim to solve
arg min
f(x subject to cx(x) 2 0, Vi E Im-
(19)
tapt oE(0,,])
This is similar to LineBO [34] but for an objective under potentially multiple constraints. However, in contrast to LineBO, our approach does not attempt global convergence along the line. Instead, we aim to select a sufficiently promising a, that yields progress given a limited evaluation budget M for the line search which set to in all experiments. Similar to [14]. either choose the next point to be the index of the best feasible point, or, none of the points are feasible, as the point with the least amount of constraint violations
arg min ( EF f(x("))
if
0t++ +
arg min ISISM Even max (0, -c (z, (z) otherwise.
(20)
{ (x(()) C (z(()) 0, Vi E lm] denotes the of feasible points among the M samples.
where F
4.4 Practical considerations and intuition optimization behavior
Local sub-sampling Unlike GIBO-style methods [43, 45, 23], we decide against adaptive sub- sampling which would require optimizing over the uncertainty of the Hessian which is computationally very expensive. Instead, to approximate local curvature after each line search, we sample K points from d-dimensional ball of radius centered at æ Rd. For this,
first draw Sobol sequence from the hypercube [0, 1]d+1 Each Sobol point (æ, u) E [0, 1 x [o, 1] is then transformed such that approximates a standard normal vector to yield a unit direction R, and determines the individual radius uul//d The final sample is then Zz + r
Slack variable fallback strategy The subproblem B-SUB may become infeasible due to constraint linearization high uncertainty in gradient estimates. To address this, we implement slack variable version of B-SUB as a fallback, which guarantees feasibility by design (cf. Appendix E). This approach aligns with established practices in classical SQP methods [46]. While the resulting search direction may not provide optimal robustness against uncertainty, the constrained posterior sampling along this direction will still seek to improve upon the current iterate.
Intuition on optimization behavior To gain intuition about the parameters of and 8, and their influence on the optimization process, study BayeSQP small toy example. We generate
7

---

TuRBO
logE1
SAASBO
MPD
Bay
Best Value (t)
-2
8
66
16
50
100
50
100
100
200
Bost Value (t)
6
8
6
10
d d= 32
10
64
12
d
0
200
0
250
500
500
1000
Iteration
Iteration
Iteration
Figure 4: Unconstrained within-model comparison. As the dimensions grow, the benefit of local search increases, with BayeSQP significantly outperforming the other baselines. Note that for SAASBO, no runs completed within the 24-hour time cap when the dimensionality exceeded 32.
two-dimensional within-model objective function (cf. Appendix A) with quadratic constraint, resulting in only small feasible region in the center. Figure 3 illustrates the optimization paths for different parameterizations. The initial step from the bottom left appears identical for all parameter settings. Subsequently, however, their behaviors differ significantly. In the expected value formulation (5y.6, 6 = 0.5), the linearization the quadratic constraint results tangential directions Pk. leading to limited or no improvement. We observe that incorporating uncertainty into the subproblem pushes the search direction toward the feasible set. Additionally, selecting very low value for &, effectively robustifies the constraints, as shown by the resulting directions PE.
5 Empirical evaluations
We next quantitatively evaluate our proposed method BayeSQP. Our evaluation first considers uncon- strained and then constrained optimization problems using BoTorch [5]. We benchmark against four baselines: logarithmic EI (logEI) [1, 32], TuRBO [15]. SAASBO [13], and MPD [45]. These baselines are widely used [40, 29, 51. 66] and represent complementary approaches--ogts employs classic global optimization strategy, TuRBO implements pseudo-local approach, SAASBO aims to automati- cally identify and exploit low-dimensional structure within high-dimensional search spaces through hierarchical sparsity prior, and MPD is fully local BO approach. Additionally, logEI and TuRBO can be readily adapted for constrained optimization through their respective variants: C-logEI [1, 17, 19] and SCBO [14] to which we compare on the constrained optimization problems.
In all subsequent plots, we present the median alongside the 5th 95th percentile range (90% inner quantiles) computed across 32 independent random seeds. For BayeSQP, we set the hyperparameters ō 8. 0.2 (unless stated otherwise) and d + 1, following Wu et al. [65, Corollary 1].
Unconstrained optimization We first consider unconstrained within-model problems [26] for which we adapt B-SUB accordingly. We generate the functions using random Fourier features fol- lowing [50, 64] (cf. Appendix A for all details). Optimizing such functions has gained relevance with recent advances latent space BO [60, 22, 40], where GP priors are enforced the latent space [52]. Figure 4 summarizes the results. BayeSQP outperforms the other baselines from dimension 16 onward. Furthermore, we can observe the siep-like behavior of Saye resulting from the subsampling followed by solving B-SUB and the subsequent line search which yields the improvement.
3 All simulations were performed on the same HPC cluster with Intel Xeon 8468 Sapphire at 2.1 GHz. Ressutt computed across that successfully found feasible solutions.
8

---

SCBO
C-logEl Pa
SCBO
C-logEI Bay
Best Value (+)
Optim. time
109
Parameerr 8y
Best Value (c-)
-5
10a
-10
10
4816
32
64
96
4816
32
64
96
0.01
0.1
0.5
Dimension d
Dimension
Parameter
(a) Performance.
(b) Optimization time."
(c) Hyperparrmeeerr sensitivity.
Figure 5: Constrained within-model comparison. BayeSQP demonstrates superior performance at high dimensions, fast optimization times, as well as low sensitivity to parameter choice.
Table 1: Results on popular BO benchmarks with multiple optima [14, 37]. BayeSQP's local search sometimes results in worse performance but crucially it always finds feasible solutions.
Method
Ackley5D
Hartmann
AoDley20D
Ac (constr.)
Hartmann (constr.)
A0Dley (constr.)
(C-)logEI
2.47
3.32
2.832
2.51 40) (feas. 32/ 32)
3.26 2:58
(feas. 32 32)
3. 3.41 (Teas. 15/32)
TURBO/
0.77.
-3.32
2.20 30
0.51 38 (feas. 32/
33.32 3.00 (feas. 32/ 32)
2. 2.03 32 (feas. 32/32)
SMS80
1.8
-3.2
76
YPD
12.57
0.01
13.36 AM
BayeSQ?
8.95
-3.3023.32
4M
10.66 48
6,257 (leas. 32/ 32)
52
-3.32 (leas. 32/ 32) 3.90 63 (feas. 32/ 32)
63
Constrained optimization Similarly, we can perform within- comparisons for the con- strained case. Here, also the constraint function c(x) sample from an GP. Again, all details provided in Appendix A. Figure 5 summarizes the constrained within-model results. As in the unconstrained case, BayeSQP outperforms the baselines at high dimensions (Figure 5a), while remaining orders of magnitude faster than SCBO and C-logEI despite computing full Hessians per B-SUB (Figure 5b). However, as we keep increasing dimensions, computing the Hessians of size d x d will results in computational overhead. Here, low-rank approximations might be useful for balancing the trade-off between computational efficiency and required accuracy of the subproblem-it is likely that especially in the context of BO, the accuracy of the Hessian is not of utmost importance. For a detailed runtime breakdown and discussion we refer to Appendix F. Lastly, in Figure 5c we can observe the influence of the parameters 6y and 8 of B-SUB the performance for d = 64. We can observe as visualized Figure 3, not considering uncertainty especially in the constraints (6, = 0.5) will result in suboptimal performance for such highly non-convex constraints. Including uncertainty results in a small buffer to the boundary, allowing the algorithm to escape local optima with a small region of attraction. The figure further highlights that beyond the decisive factor of taking uncertainty into account the overall sensitivity on how much uncertainty should be incorporated is rather low. The optimal values of these parameters may vary depending on the specific application.
Performance on standard benchmarks Lastly, we also evaluate BayeSQP on standard BO bench- marks. Here, we follow recent best practices and initialize lengthscales with v'd for all baselines 67.48), 48). The results are summarized in Table 1. We can clearly observe that BayeSQP sensitive to initialization highlighted by the large 90% quantile especially for Ackley. This is to be expected as the algorithm local and Ackley is very multi-modal. Still, importantly, BayeSQP is able to find feasible solutions for all seeds in all benchmarks contrary to C-logEI.
[30, 67,
To demonstrate the real-world applicability of BayeSQP, we compare constrained opti- mization baselines on the 7-dimensional Speed Reducer benchmark [36], which minimizes the weight of a speed reducer subject to 11 mechanical design -ineir constraints (more details Appendix A.4). The results are summarized in Table 2. All baselines are able to find feasible solutions for all seeds. C-logEI and BayeSQP show the best performance. In line with previous experiments,
Table 2: Performance on Speed Reducer [36].
Method
Performance
Avg. runtime (s)
SCBO
3006.85 33 13.25
(feas. 32 32)
286.46
c-logEI
3002.81
2006 67
(feas. 32)
3464.59
BayeSQP
3001.10 2996.97 97 (feas. 32 32)
3000.30
91.83

---

BayeSQP demonstrates a clear runtime advantage even in the presence of 11 constraints each requiring separate Hessian evaluations- and substantially larger B-SUB.
6 Discussion on limitations
While BayeSQP provides a novel framework combining classic optimization methods with BO, there several limitations and addressing them will be interesting future research.
Initialization matters As with any local approach, the initialization of BayeSQP will directly influence its per- formance (cf. Table 1). This further becomes clear when looking at the flow field of BayeSQP generated from 1000 different initial conditions on the Gramacy benchmark [20] in Figure 6 (details in Appendix A.5). Depending the initialization, the algorithm converges to a different local optimum of the constrained problem.
(a(3) < D
Optimization flow
1.0
0.4
Although global approaches can also exhibit sensitivity to initialization, this sensitivity is amplified in LBO ap- proaches, particularly in constrained optimization. How- ever, this sensitivity provides practitioners with the option to incorporate some expert knowledge into the optimiza- tion by choosing the initial guess; especially in engineering fields such as robotics, a feasible yet non-optimal solution is often known a-priori. An algorithm like BayeSQP will then become an automatic tool for fine-tuning.
0.2
0.0
Figure 6: Optimization behavior Gra- macy [20]. Depending on the initializa- tion, the BayeSQP will converge to a dif- ferent local optimum. Lighter colors and thicker lines indicate larger ||P:|/2.
Computational considerations We show that for up to 96 dimensions even with the additional cost of computing the Hessian, BayeSQP demonstrates as very low total runtime. Still, at very large dimensions or high number of constraints, computing as well as storing the Hessian of all constraints will become problematic. In principle, one could also incorporate Hessian uncertainty into B-SUB, for example following efficient schemes such as [2. 11]; whether this would lead to empirical performance improvements remains open question. Future work could focus evaluating the joint GP over only the most informative Hessian entries, adaptively selected during optimization, or on constructing the Lagrangian Hessian directly from gradient histories using BFGS-type update scheme.
Dependency on the kernel and model assumptions The performance of BayeSQP strongly depends on the choice of kernel and, more generally, on the modeling assumptions underlying the GP surrogate. Since the construction of the B-SUB directly relies on the accuracy of both gradient and Hessian estimates, poorly chosen kernel can lead to unreliable curvature information and ultimately to suboptimal search directions. While standard kernels such as the suwared-xxponential kernel perform well for smooth problems, they may struggle in settings with sharp nonlinearities or discontinuous constraints unless handled with additional care. Furthermore, kernel hyperparameters influence the scale and conditioning of the estimated Hessian, which can significantly affect the resulting search direction. Advances in GP modeling and training practices for BO (e.g., [30, 67]) are expected to directly improve the robustness and effectiveness of BayeSQP.
In Appendix B, we list possible extensions BayeSQp which in part address the limitations mentioned above as well as further interesting directions for future work.
7 Conclusion
In this paper, we presented BayeSQP as a bridge between classic optimization methods and BO. BayeSQP uses GP surrogates that jointly model the function, its gradient and its Hessian, which are then used to construct subproblems in an SQP-like fashion. Our results show that BayeSQP can outperform state-of-the-art methods in high-dimensional constrained optimization problems. We believe that BayeSQP provides a promising framework for integrating well-established classical optimization principles with modern black-box optimization techniques.
10

---

Acknowledgements
The authors thank David Stenger, Alexander von Rohr, Johanna Menn, Tamme Emunds, and Henrik Hose for various discussions on BO and optimization in general. Paul Brunzema is partially funded by the Deutsche Forschungsgemeinschatt (DFG, German Research Foundation)-Ree 2236/2 (Un- RAVeL). Simulations were performed in part with computing resources granted by RWTH Aachen University under projects rwth1579, p0022034, and p0021919.
References
[1] Sebastian Ament, Samuel Daulton, David Eriksson, Maximilian Balandat, and Eytan Bakshy. Unexpected improvements to expected improvement for Bayesian optimization. In Advances in Neural Information Processing Systems (NeurIPS), 2023.
[2] Sebastian E Ament and Carla P Gomes. Scalable first-order Bayesian optimization via structured automatic differentiation. In International Conference on Machine Learning (ICML), 2022.
[3] M. Andersen, J. Dahl, and L. Vandenberghe. CVXOPT: a Python package for convex optimiza- 2013.
tion (version 1.3.2). https:://cvooo/
[4] Setareh Ariafar, Jaume Coll-Font, Dana Brooks, and Jennifer Dy. ADMMBO: Bayesian optimization with unknown constraints using ADMM. Journal of Machine Learning Research (JMLR), 2019.
[5] Maximilian Balandat, Brian Karrer, Daniel R. Jiang, Samuel Daulton, Benjamin Letham, Andrew Gordon Wilson, and Eytan Bakshy. BoToreh: framework for efficient Monte-Carlo Bayesian optimization. In Advances in Neural Information Processing Systems (NeurIPS), 2020.
[6] Felix Berkenkamp, Angela Schoellig, and Andreas Krause. Safe controller optimization for quadrotors with Gaussian processes. In IEEE International Conference Robotics and Automation (ICRA), 2016.
[7] Eli Bingham, Jonathan P Chen, Martin Jankowiak, Fritz Obermeyer, Neeraj Pradhan, Theofanis Karaletsos, Rohit Singh, Paul Szerlip, Paul Horsfall, and Noah D Goodman. Pyro: Deep universal probabilistic programming. Journal of Machine Learning Research, 2019.
[8] Paul Brunzema, Mikkel Jordahn, John Willes, Sebastian Trimpe, Jasper Snoek, and James Harrison. Bayesian optimization via continual variational last layer training. In International Conference on Learning Representations (ICLR), 2025.
[9] Roberto Calandra, André Seyfarth, Jan Peters, and Mare Peter Deisenroth. Bayesian optimiza- tion for learning gaits under uncertainty: An experimental comparison on a dynamic bipedal walker. Annals of Mathematics and Artificial Intelligence, 2016.
[10] Lionel Colliandre and Christophe Muller. Bayesian optimization in drug discovery. High Performance Computing for Drug Discovery and Biomedicine, 2023.
[11] Filip De Roos, Alexandra Gessner, and Philipp Hennig. High-dimensional Gaussian process inference with derivatives. In International Conference on Machine Learning (ICML), 2021.
[12] Thomas Elsken, Jan Hendrik Metzen, and Frank Hutter. Neural architecture search: A survey. Journal of Machine Learning Research (JMLR), 2019.
[13] David Eriksson and Martin Jankowiak. High-dimensional Bayesian optimization with sparse axis-aligned subspaces. In Uncertainty in Artificial Intelligence (UAI), 2021.
[14] David Eriksson and Matthias Poloczek. Scalable constrained Bayesian optimization. In International Conference on Artificial Intelligence and Statistics (AISTATS), 2021.
[15] David Eriksson, Michael Pearce, Jacob Gardner, Ryan D Turner, and Matthias Poloczek. Scalable global optimization via local Bayesian optimization, In Advances in Neural Information Processing Systems (NeurIPS), 2019.

---

[16] Zheyi Fan, Wenyu Wang, Szu H Ng, and Qingpei Hu. Minimizing UCB: a better local search strategy in local Bayesian optimization. In Advances in Neural Information Processing Systems (NeurIPS), 2024.
[17] Jacob R Gardner, Matt Kusner, Zhixiang Eddie Xu, Kilian Weinberger, and John P Cun- ningham. Bayesian optimization with inequality constraints. In International Conference on Machine Learning (ICML), volume 2014, 2014.
[18] Roman Garnett. Bayesian optimization. Cambridge University Press, 2023.
[19] Michael A Gelbart, Jasper Snoek, and Ryan P Adams. Bayesian optimization with unknown constraints. Conference on Uncertainty in Artificial Intelligence (UAT), 2014.
[20] Robert B Gramacy, Genetha A Gray, Sébastien Le Digabel, Herbert KH Lee, Pritam Ranjan, Garth Wells, and Stefan M Wild. Modeling an augmented lagrangian for blackbox constrained optimization. Technometrics, 2016.
[21] Ryan-Rhys Griffiths and José Miguel Hernández-Lobato. Constrained Bayesian optimization for automatic chemical design using variational autoencoders. Chemical Science, 2020.
[22] Antoine Grosnit, Rasul Tutunov, Alexandre Max Maraval, Ryan-Rhys Griffiths, Alexander I Cowen-Rivers, Lin Yang. Lin Zhu, Wenlong Lyu, Zhitang Chen, Jun Wang, et al. High- dimensional Bayesian optimisation with variational autoencoders and deep metric learning. arXiv preprint arXiv:2 100.06696, 2021.
[23] Shiming He, Alexander von Rohr, Dominik Baumann, Ji Xiang, and Sebastian Trimpe. Simulation-aided aiddd policy tuning for black-box robot learning. IEEE Transactions on Robotics (TRO), 2025.
[24] Philipp Hennig. Fast probabilistic optimization from noisy gradients. In International Confer- ence on Machine Learning (ICML), 2013.
[25] Philipp Hennig and Martin Kiefel. Quasi-Newton methods: A new direction. Journal of Machine Learning Research (JMLR), 2013.
[26] Philipp Hennig and Christian J Schuler. Entropy search for information-efffcciet global opti- mization. Journal of Machine Learning Research (JMLR), 2012.
[27] Philipp Hennig, Michael A Osborne, and Mark Girolami. Probabilistic numerics and uncertainty in computations. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 2015.
(28) Philipp Hennig, Michael A Oshome,, and Hans P Kersting. Probabilistic Numerics: Computation as Machine Learning. Cambridge University Press, 2022.
[29] Henrik Hose, Paul Brunzema, Alexander von Rohr, Alexander Gräfe, Angela P Schoellig, and Sebastian Trimpe. Fine-tuning of neural network approximate MPC without retraining via Bayesian optimization. In CoRL Workshop on Safe and Robust Robot Learning for Operation in the Real World, 2024.
[30] Carl Hvarfner, Erik Orm Hellsten, and Luigi Nardi. Vanilla Bayesian optimization performs great in high dimensions. International Conference on Machine Learning (ICML), 2024.
[31] Hideaki Ishibashi, Masayuki Karasuyama, Ichiro Takeuchi, and Hideitsu Hino. A stopping crite- rion for Bayesian optimization by the gap of expected minimum simple regrets. In International Conference on Artificial Intelligence and Statistics (AISTATS), 2023.
[32] Donald R Jones, Matthias Schonlau, and William Welch. Efficient global optimization of expensive black-box functions. Journal of Global optimization, 1998.
[33] Kirthevasan Kandasamy, Jeff Schneider, and Barnbáss P6ezos. High dimensional Bayesian optimisation and bandits via additive models. International Conference on Machine Learning (ICML), 2015.
12

---

[34] Johannes Kirschner, Mojmir Mutny, Nicole Hiller, Rasmus Ischebeck, and Andreas Krause. Adaptive and safe Bayesian optimization in high dimensions via one-dimensional subspaces. In International Conference on Machine Learning (ICMI., 2019.
[35] Agustinus Kristiadi, Alexander Immer, Runa Eschenhagen, and Vincent Fortuin. Promises and pitfalls of the linearized laplace in bayesian optimization. Symposium on Advances in Approximate Bayesian Inference (AABI), 2023.
[36] Afonso CC Lemonge, Helio JC Barbosa, Carlos CH Borges, and Francilene BS Silva. Con- strained optimization problems in mechanical engineering design using a real-coded steady-state genetic algorithm. Mecánica Computacional, 29(95):9287-9) -903, 2010.
[37] Benjamin Letham, Brian Karrer, Guilherme Ottoni, and Eytan Bakshy. Constrained Bayesian optimization with noisy experiments. Bayesian Analysis, 2019.
[38] Maren Mahsereci and Philipp Hennig. Probabilistic line searches for stochastic optimization. Journal of Machine Learning Research (JMLR), 2017.
[39] Anastasia Makarova, Huibin Shen, Valerio Perrone, Aaron Klein, Jean Baptiste Faddoul, Andreas Krause, Matthias Seeger, and Cedric Archambeau. Automatic termination for hy- perparameter optimization. In International Conference on Automated Machine Learning, 2022.
[40] Natalic Maus, Haydn Jones, Juston Moore, Matt J Kusner, John Bradshaw, and Jacob Gard- ner. Local latent space Bayesian optimization over structured inputs. In Advances in Neural Information Processing Systems (NeurIPS), 2022.
[41] Natalic Maus, Kyurae Kim, Geoff Pleiss, David Eriksson, John P Cunningham, and Jacob R Gardner. Approximation-awwan aware Bayesian optimization. Advances in Neural Information Processing Systems (NeurIPS), 2024.
[42] Mitchell McIntire, Daniel Ratner, and Stefano Ermon. Sparse Gaussian processes for Bayesian optimization. In Conference on Uncertainty in Artificial Intelligence (UAI), 2016.
[43] Sarah Müller, Alexander von Rohr, and Sebastian Trimpe. Local policy search with Bayesian optimization. In Advances in Neural Information Processing Systems (NeurIPS), 2021.
[44] Matthias Neumann-Brosig, Alonso Marco, Dieter Schwarzmann, and Sebastian Trimpe. Data- efficient autotuning with Bayesian optimization: An industrial control study, IEEE Transactions on Control Systems Technology, 2019.
[45] Quan Nguyen, Kaiwen Wu, Jacob Gardner, and Roman Garnett. Local Bayesian optimization via maximizing probability of descent. In Advances in Neural Information Processing Systems (NeurIPS), volume 35, 2022.
[46] Jorge Nocedal and Stephen J Wright, Numerical optimization. Springer, 2006.
[47] Leonard Papenmeier, Luigi Nardi, and Matthias Poloczek. Increasing the scope as you learn: Adaptive Bayesian optimization in nested subspaces. Advances in Neural Information Process- ing Systems (NeurIPS), 2022.
[48] Leonard Papenmeier, Matthias Poloczek, and Luigi Nardi. Understanding high-dimensional 2025.
Bayesian optimization. arXiv preprint
[49] Victor Picheny, Robert B Gramacy, Stefan Wild, and Sebastien Le Digabel. Bayesian opti- mization under mixed constraints with a slack-variable augmented Lagrangian. In Advances in Neural Information Processing Systems (NeurIPS), 2016.
[50] Ali Rahimi and Benjamin Recht. Random features for large-scale kernel machines. In Advances in Neural Information Processing Systems (NeurIPS), 2007.
[51] Pantaleo Raimondi, Chamseddine Benabderrahmane, Paul Berkvens, Jean Claude Biasci, Pawel Borowiec, Dean-Francois Bouteille, Thierry Brochard, Nicholas B Brookes, Nicola Carmignani, Lee Carver, et al. The extremely brilliant source storage ring of the european synchrotron radiation facility. Communications Physics, 2023.
13

---

[52] Siddharth Ramchandran, Manuel Haussmann, and Harri Lähdesmäki. High-dimensiona bayesian optimisation with gaussian process prior variational autoencoders. In International Conference Learning Representations (ICLR), 2025.
[53] Carl Edward Rasmussen and Christopher K. I. Williams. Gaussian processes for machine learning. MIT Press, 2006.
[54] Rommel G Regis and Christine A Shoemaker. Combining radial basis function surrogates and dynamic coordinate search in -igh-dimensional expensive black-box optimization. Engineering Optimization, 2013.
[55] Jasper Snock, Hugo Larochelle, and Ryan P Adams. Practical Bayesian optimization of machine learning algorithms. In Advances in Neural Information Processing Systems (NeurIPS), 2012.
[56) Jasper Snoek, Oren Rippel, Kevin Swersky, Ryan Kiros, Nadathur Satish, Narayanan Sundaram, Mostofa Patwary, Mr Prabhat, and Ryan Adams. Scalable Bayesian optimization using deep neural networks. International Conference on Machine Learning. 2015.
[57] Shiro Tamiya and Hayata Yamasaki. Stochastic gradient line Bayesian optimization for efficient noise-robust optimization of parameterized quantum circuits. npj Quantum Information, 2022.
[58] Wei-Ting Tang and Joel A Paulson. CAGES: Cost-aware gradient entropy search for efficient local multi-fidelity Bayesian optimization. arXiv preprini SxXiv:2405.07766 2024.
[59] Wei-Ting Tang, Akshay Kudva, and Joel A Paulson. NeST-BO: Fast local Bayesian op- timization via Newton-step targeting of gradient and Hessian information. arXiv preprint arXiv:2510.05516. 2025.
[60] Austin Tripp, Erik Daxberger, and José Miguel Hernández-Lobato. Sample-efficient optimiza- tion in the latent space of deep generative models via weighted retraining. In Advances in Neural Information Processing Systems (NeurIPS), 2020.
[61] Alexander von Rohr, David Stenger, Dominik Scheurenberg, and Sebastian Trimpe. Local Bayesian optimization for controller tuning with crash constraints. al-Autonatisiernngsttiiiikk 2024.
[62] Ziyu Wang, Frank Hutter, Masrour Zoghi, David Matheson, and Nando De Feitas. Bayesian optimization in billion dimensions via random embeddings. Journal of Artificial Intelligence Research, 2016.
[63] James Wilson. Stopping Bayesian optimization with probabilistic regret bounds. In Advances in Neural Information Processing Systems (NeurIPS), 2024.
[64] James T Wilson, Viacheslav Borovitskiy, Alexander Terenin, Peter Mostowsky, and Marc Peter Deisenroth. Pathwise conditioning of Gaussian processes. Journal of Machine Learning Research (JMLR), 22(105), 2021.
[65] Kaiwen Wu, Kyurae Kim, Roman Gamett,, and Jacob Gardner. The behavior and convergence of local Bayesian optimization. In Advances in Neural Information Processing Systems (NeurIPS), 2023.
[66] Xiankui Wu, Xinyu Gu, and KW See. ADNNet: Attention-based deep neural network for air quality index prediction. Expert Systems with Applications, 2024.
[67] Zhitong Xu, Haitao Wang, Jeff M Phillips, and Shandian Zhe. Standard Gaussian process is all you need for dineeiiins Bayesian optimization. In International Conference on Learning Representations (ICLR), 2024.
14

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: The introduction and abstract were written highlight the main contributions of the paper. Thus, they are reflective of the paper's contributions and scope.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
It is fine to include aspirational goals as motivation long is clear that these goals are not attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: We explicitly discuss limitations of out method in separate section (Section 6). Furthermore, we provide empirical evidence of the limitations of our method (e.g.. Table 1).
Guidelines:
The answer NA means that the paper has 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 (c.g.. independence assumptions, noiseless settings, model sell-speciiiatiinn 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 few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low or images are taken in low lighting. Or speech-to-tex system might not be used reliably to provide closed captions for online lectures because it fails 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 faimesss.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community, Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory Assumptions and Proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
15

---

Answer: [Yes]
Justification: The paper does not include any major theoretical results. All derivations are clearly stated.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross referenced.
All assumptions should be clearly stated or referenced the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but they appear in the supplemental material, the authors encouraged provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental Result Reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Question:
Answer: [Yes]
Justification: In Appendix A, we provide all details on the experimental setup, including hyperparameters for all basclines.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility instructions for how to replicate the results, access to a hosted model (e.g.. in the case of large language model), releasing of a model checkpoint, other means that appropriate to the research performed.
also be provided via detailed
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make it clear how to reproduce that algorithm.
(b) the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) If the contribution is a new model (e.g., a large language model), then there should either way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reprodueibility may be tricky in some cases, in which 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 have some path reproducing or verifying the results.
5. Open access to data and code
16

---

reproducible, they Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Question:
Answer: [Yes]
Justification: All benchmarks are publicly available and we provide instructions to reproduce all experiments. The implementation provided under https: //t:::///tttbbb com/brunzema/ bayesqp and casily accessible via PyPI.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines (https://nips.co jublic/guides/Cooee
ssionPolicyy for more details.
While we encourage the release of code and data, we understand that this might not be possible, "No" an acceptable answer. Papers cannot be rejected simply for including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https:
/ //nips. cc/pccccpccccccccccc http:pulllggguiiii wtt.c//ge////cieecooocc Seoieeee///o cesuuubmisioonni asionPoolcyy for more details.
The authors should provide instructions on data access and preparation, including how access the raw data, preprocessed data, intermediate data, and generated data,
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only subset of experiments 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) recommended, but including URLs to data and code permirted.
6. Experimental Setting/Details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: In Appendix A, provide all details on the experimental setup, including hyperparameters for all basclines.
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 appreciate the results and make sense of them.
The full details can provided either with the code, appendix, or as supplemental material.
7. 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: In all our quantitative results, report median and 90% quantiles over 20 random seeds to best communicate average as well as tail performance.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by emoor bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
17

---

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 the standard error of the mean.
It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables figures symmetric error bars that would yield results that are out range (e.g. negative error rates).
If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments Compute Resources
Question:
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Answer: [Yes]
Justification: the main body of the paper, we provide the type of compute hardware used produce the results and in Figure 5b provide wall clock times.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g.. preliminary or failed experiments that didn't make into the paper).
9. Code Of Ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics https: http://www.p.p....a/.tttttttttta
Answer: [Yes]
Justification: The paper does not violate any of the guidelines in the NeurlPS Code of Ethics.
Guidelines:
The answer NA means that the authors have not reviewed the NeurlPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., there a special consid- eration due to laws regulations in their jurisdiction).
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [No]
We not discuss societal impacts. The proposed algorithm can help to better solve time-varying optimization problems. As this is very generic, we did not include discussion on this topic. As listed the guidelines below, such discussion not necessary for generic optimization algorithms.
Justification:
18

---

Guidelines:
The answer NA means that there no societal impact of the work performed.
If the authors answer NA or No. they should explain why their work has no societal impact why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious 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 considcrations.
The conference expects that many papers will be foundational research and 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 conrectly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper deseribe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification: The paper proposed general optimization algorithm. As such, the paper poses no such risks.
Guidelines:
The answer NA means that the paper poses such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model. for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [Yes]
Justification: The paper clearly cedits the main code base on which the empirical evaluation is based. The main package BoTorch under the MIT license.
Guidelines:
The answer NA means that the paper does not existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset used and, possible, include URL.
19

---

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.
assets are released, the license, copyright information, and terms use in the package should be provided. For popular datasets, paper thcode. com/ com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that 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 's creators.
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
Justification: We will provide the code upon acceptance of the paper. The code will be released under the MIT license.
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 used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and Research with Human Subjects
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 details about compensation (if any)?
Answer: [NA]
Justification: Our paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does involve crowdsourcing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According to the NeurlPS 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.
15. 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 equivalent approval/review based the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: Our paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsoureing nor research with human subjects.
20

---

Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
21

---

A Design of experiments
In the following, we provide further details on the design of experiments for the experiments in Section 5 as well as discussion on the baselines and model initialization and training.
A.1 Generating within-model objective functions
Within-model comparisons were introduced by Hennig and Schuler [26] to study the performance of BO methods functions that fullfil all model assumptions. With recent advances in latent space BO [52], optimizing such functions has gained relevance e.g., for drug discovery. To generate the within-model objective functions shown in Figure 3 and discussed in Section 5, we approximate prior samples ft with a M random Fourier features (RFFs) following Rahimi and Recht [50]. This yields parametric function
[())) Σ Wm Om (a) with (m))
V
2 M
cos O + T-
(21)
Here, 6m are sampled proportional to the kernel's spectral density and Tm ~ L(0, 2x). In all experiments, we use a equared-epponentaa kernel with lengthscales & = 0.1 and M = 1028 RFFs. For SAASBO, we report the out-of-model comparison results as the main mechanism of SAASBO is the way finds suitable hyperparameters for the given task (cf. Appendix A.6 and cf. Appendix D).
A.2 Generating constrained within-model objective functions
To generate the constrained within-model objective functions, we use the same approach as for the unconstrained case. Additionally to generating within-model objective function, we also generate
within-moded constraint function Ö(z). We then shift this function by one, i.e., c(x) = c(x) - 1 2 0 so that average only about 16% of the domain is feasible. 5 Note that multiple constraints possible, however, in the within-model setting with a shifted mean we do run the risk of generating an infeasible problem and therefore opted for only one constraint. In the other constrained benchmarks, we also consider multiple constraints.
also
A.3 Constrained versions of Ackley and Hartmann
For the constrained versions of Ackley and Hartmann, we use the pre-implemented benchmarks in BoTTrch [5] which are largely based on experiments in Letham et al. [37] and Eriksson and Poloczek [14]. For Ackley objectives, there are two inequality constraints and for Hartmann only one. Specifically, we use the following constraints for the Ackley and Hartmann function:
{
c1(z) = >i= c2(x)) = 5 =
Ackley:
Hartmann: { c1
- -1
(22)
For Ackley, we restrict the feasible region to [-5, 5, 10 [15] and use time horizon of T = 100 for = 5 and = 400 for d = 20. The feasible region for Hartmann is [0, 1]6 and we set T 100.
A.4 Speed Reducer
The Speed Reducer design problem aims to minimize the weight of a speed reducer mechanism subject to 11 mechanical constraints. The design variables are: face width (x1 [2.6, 3.6)), module of teeth (r2 E [0.7, 0.8)., number of teeth on pinion (I3 E [17, 28), integer which we treat as a continuous variable implemented in BoTorch and consistent with prior work [15]), length of shaft
1 (I4 E [7.3. 8.3), length of shaft 2 (x5 E [7-8, 8.3)). diameter shaft 1 (va E [2.9, 3.9). and diameter of shaft (x² E 5.5)) We run the benchmark for 200 iterations across 32 random seeds with 7' 200 and report the results in Table 2. For BayeSQP, we set ōy, 8 0.5 essentially reverting to expected value formulation. As discussed, including uncertainty will result in the final solution being robust in the sense of not directly laying on the boundary. We find that directly at the boundary, the objective for Speed Reducer significantly improves. A very practical approach could also be to schedule 5y. & over the number of iterations. For the full formulation of the optimization problem, we refer the BoTorch documentation (v.13.0) as well as the original paper [36].
'e have for output scale of that P(((z) 2 0) = P(Z 1) 0.16 where N/ 1).
22

---

A.5 Constrained Gramacy function
For the constrained Gramacy benchmark, we use the formulation introduced in Gramacy et al. [20]
and implemented in BoTorch. The problem is defined over the unit square x E [0, 1] with the objective of minimizing the sum of the two decision variables, f(x) 21 + z2. It imposes two nonlinear inequality constraints
(x) -1.55 - 31 2x² - 0.5 (2m(a) 2x²))),
(23)
c2(x)) = (x} + 2 1.5).
(24)
The problem is non-convex as shown in Figure 6. For these problems, local BO approaches are especially sensitive to the initialization; depending on the initialization, different local optima may be reached. To converge to global optimum over time, approaches such as restarting BayeSQP are promising. For a longer discussion, we refer to Appendix B.
A.6 Baselines
Expected improvement (EI) is a widely used acquisition function in BO. However, optimizing the EI acquisition function can be numerically unstable. To address this, Ament et al. [1] proposed a logarithmic transformation of El resulting a numerically more stable acquisition function even resulting in an increase in performance. We use the implementation of logEI from BoTorch [5] and adapt it for constrained optimization by using the same wrapping standard constrained EI [17, 19] resulting in C-logEI Ament et al. [1]. Neither logEI nor C-logEI require additional hyperparameters.
(C-)logEI
TuRBO and SCBO For the implementation of TuRBO and SCBO, we follow tutorials from BoTorch [5] which were provided by the authors of the respective methods. TuRBO (and SCBO) require various hyperparameters which specify when and by how much to shrink or expand the trust region. To set these hyperparameters, tutorial. With these recommendations. the initial length of the trust region is Linit = 0.8, the minimum and maximum length of the trust region are Lmin = 0.5² and Lmax = 1.6, respectively, the number of consecutive failures before the trust region is shrunk is That = [max (4, d)], the number of consecutive successes before the trust region is expanded is Tauee 3. The trust region is always centered around the best point found so far which in the context of constrained optimization follows (20) for SCBO. For posterior sampling, we evaluate the GP posterior within the trust region at 2000 points which are drawn from a Sobol sequence. For both methods, we use the recommended pertubation masking to cope with discrete sampling in high-dimensional spaces [54, 15), i.e., in order to not perturb all coordinates at once, we use the value in the Sobol sequence with probability min(1, 20/d /d for a given candidate and dimension, and the value of the center of the trust region otherwise which induces an exploitation bias [48]. While TuRBO (and SCBO) depend on all the above parameters for their trust region and sampling heuristies, we found that these suggested parameters work well across various tasks as also highlighted in previous work [40, 29].
follow the recommendations of Eriksson al. [15] in the mentioned
Sparse axis-aligned subspace Bayesian optimization (SAASBO) [13] is designed for high- dimensional BO by placing hierarchical sparsity priors on inverse lengthscales to identify and exploit -ow-dimensional structure. The method uses a global shrinkage parameter T HC(3) and pnmeesion--seciffi inverse lengthscales Pa ~ HC((T) for all d E 1d, where HC denotes the half-Cauchy distribution. This prior encourages small values while allowing heavy tails that enable relevant dimensions to escape shrinkage toward zero. We follow the BoTTorc tutorial which performs inference using Hamiltonian Monte Carlo (HMC) with the NUTS sampler from Pyro [7] and uses logEI [1] as acquisition function. SAASBO's computational cost scales cubically with the number of observations due to HMC resulting a significant computational scaling as shown in Figure 8.
SAASBO
MPD Local BO via maximizing probability of decent (MPD) [45] is a follow-up method to the discussed GIBO approach. Unlike GIBO, it defines a different acquisition function for the aub-sampling step that aims to maximize the probability that the posterior mean points in descent direction, rather than minimizing the uncertainty about the gradient. Additionally, reuses the estimated posterior gradient GP to iteratively move the current point along the most probable descent direction
Availlabee under the MIT license at https://boocch.orr//dossuttrraae
(BoTorch version v0.15.1).
23

---

TuRBO
logE1
SAASBO
MPD
BayeSQP
Best Value (+)
8
16
50
100
50
100
100
200
Best Value (-)
2
6
d
32
64
10
96
0
200
250
500
500
1000
Iteration
Iteration
Iteration
Figure 7: Unconstrained out-of-model comparison. As the dimensions grow, the benefit of local search increases, with BayeSQP outperforming the other baselines.
until the probability falls below a predefined threshold. We use the parameters implemented for the synthetic functions from the respective repository. We should note that we did not further tune these parameters to our specific problems. We hypothesize that some of the sub-optimal performance can be attributed the second stage of the algorithm, where the current GP estimate may be trusted for too many iterations before entering a sub-sampling step.
BayeSQp For BayeSQP, we set the hyperparameters 8y, K = and 0.05 in all experiments, unless stated otherwise. To solve B-SUB at each iteration, we use CVXOPT [3] with standard parameters for maximum iterations and tolerance. 8 Until we have reached feasible point, we set ō, 0.5 to focus on robust improvement of the constraints and switch to the specified value once we have observed feasible point. For the subsequent line search, we use constrained posterior sampling as in SCBO but without the perturbation masking as we only operate on one-dimensional line segment. On this line segment. 100 sample points from a Sobol sequence as candidate points and choose the next sample location following Friksson and Poloczek [14).
A.7 Model initialization and training
kernel. This kernel sufficiently smooth such that we can formulate the joint GP for BayeSQP as specified in (7). For the within-model comparisons, freeze the lengthscales of the kernel and do not perform any hyperparameter optimization. For the experiments classic benchmarks, kernel into scale kernel and initialize all lengthscales with vd. We furthermore as corse bounds on the lengthscales as l E [0.001, 2d] for all baselines and set the noise to a small values i.e., 2² 10-4. We subsequently optimize all hyperparameters by maximizing the marginal log- likelihood. In all experiments and for all baselines, we use standardization as a output transformation improve numerical stability of the GP and normalize the inputs to the unit hypercube [0, 1]4 [5].
For all algorithms, we use a
follow recent best-practice, wrap the exxonentxx
B Extensions to different settings and ideas for future work
BayeSQP provides a flexible framework that can be extended various settings. In this section, we briefly discuss some of these extensions as well as other interesting avenues for futuree work.
TThe repository is under the MIT license at htt.ggghhw...m///yy//////aaaaaaa $CVXPP is publicly available under modified GNU GENERAL PUBLIC LICENSE.
24

---

Termination and restarting The question of when to stop optimizing in BO has gained increasing attention in recent years [39, 63]. Addressing this question directly increases the practicality of an algorithm especially in the context of robotics where hardware experiments are often very expensive. For BayeSQP, we can draw inspiration from traditional SQP methods
develop an appropriate termination criterion such as stopping optimization once l l2 S Tool i.e., when the search direction becomes sufficiently small, indicating likely convergence to local optimum and little progress is to be expected in the subsequent line search. After termination, we can leverage ideas from TuRBO. Similar to how TuRBO restarts after trust region collapse, our algorithm can randomly reinitialize when optimization terminates, given that there remains sufficient computational budget.
Batch optimization Similar to TuRBO, implementing BayeSQP for batch BO is straightforward: We can utilize different initial conditions distinct starting points for local optimization. All resulting data points can be combined into the same GP model or as in TuRBO separated in different data sets. In general, scaling local BO methods to batch optimization is particularly promising as these algorithms inherently remain confined to the local region surrounding their initialization point likely generating high-diversity batches.
Localized GP model To combat model mismatch on real- problem, is possible to include sliding window on the training data of size Nmax as also proposed in [43, 45]. This effectively produces purely local model, which can better capture the local structure. However, some care has to be taken here to as this might result in unstable learning of kernel hyperparameters.
Active sub-sampling In its current version, BayeSQP relies a sub-sampling step to get good posterior estimates for the gradients and Hessians. While space-filling sampling using a Sobol sequence already yields good results (cf. Section 5), an active approach to the sub-sampling potentially with stopping criterion-is interesting. One approach would to build on ideas from Tang et al. [59]. However, an active sub-sampling likely will result in much slower overall runtime so it depends on the specific application this is desirable.
C Proof of Corollary
The proof of Corollary follows directly from the fact that - = O. With this, by and be. are zero and longer influence the constraints objective. Since by and be, are optimization variables in B-SUB, the cones can be trivially satisfied.
D Additional results
Out-of-model comparisons In addition the within-model comparisons from Section 5, we also perform out-of-model comparisons. In this setting, the objective still satisfies the assumption that the model is sample from GP, but instead of passing the correct lengthscales to the models, each baseline learns these lengthscales. For all baselines SAASBO, we initialize the lengthscales with the true length- scales of the objective. The results in Figure 7. We can observe the same trends as in the within-men mohell comparisons thought gap the final performance of TuRBO and BayeSQP is smaller, Look- ing the optimization time the out-of-model com- parisons in Figure 8, still is apparent with BayeSQP being two orders of magnitudes faster than TuRBO for the 96 dimensional problems. As stated in the main text, for dimensions larger than 32, SAASBO failed to solve the problem at hand within the 24h time cap of the server but Figure still shows the clear trend of the local approaches BayeSQP and TuRBO outperforming SAASBO in high dimensions.
TuRBO
SAASBO
BaycSQP
logEl
MPD
Optim. time [s]
10
103
10a
summarized
10
48
16
32
64
96
Dimension d
Figure 8: Optimization time of outeof-model comparison. BayeSQP shows significantly faster optimization compared to other base- lines also in the out-of setting.
25

---

10.0
Best Value (c-)
Param. K
Paran K
16
16
Dest Value
32
52
16
64
32
-1.0
128
128
206
11.5
3
5
9
0.4
Param.
Param. 6j
(a) Const. within-model.
(b) Const. out-of-model.
(c) Unconst. within-model.
Figure 9: Ablations for BayeSQP. Depending the specific problem, different parameter combina- tions may yield optimal performance. All results shown in the figures correspond to an objective (and constraint) function with 64 dimensions and each field reports the median over 20 seeds.
Ablation on the number of sub-samples K and line search samples M Figure 5c showed that for the specified K and M, i... number of sub-samples and number of line search samples, including some uncertainty in B-SUB will result in better performance. In Figure 9 we provide further ablations.
Figure 9a presents the results of the sensitivity analysis over the number of sub- uub-sampling steps (K) and line-search steps (M) for both the constrained within-mode model and out-of-mode settings. In the within- model case, increasing
beyond d tends degrade performance, whereas decreasing K which allows for more B-SUB solves under the same computational budget-con can lead to improvements. For out-of-model comparisons, the trends are less pronounced: increasing the number of line-search samples generally helps, while choosing M too small can be detrimental.
Overall, the results suggest that when strong priors for the surrogate models are available (e.g.. from domain knowledge), reducing K can enhance performance. In contrast, when such priors are absent, increasing M may offer better results. These findings further highlight the potential benefit of introducing suitable stopping criterion for the line search, enabling online adaptation of M. Finally, note that the performance BayeSQP in Figure 7 could likely be improved through hyperparameter tuning though similar improvements may be achievable for the other baselines as well.
Lastly, Figure 9c shows that, in the unconstrained within-model case, reducing K, as also observed in Figure 9a, can be benefi- cial. Moreover,
without accounting for uncertainty, the constrained case can achieve very good performance provided that the number of B-SUB solves is sufficiently large. Notably, for K = 64, incorporating uncertainty leads to improved re- sults, indicating that this configuration can more effectively handle scenarios with limited number of line-searches. The oft-oflmodel case in Figure 10 however highlights that in the absence of good prior knowledge, reducing the number of can be costly for small 5,. small óy and high uncertainty in the estimates in B-SUB will lead to small Ps- and with this only optimum. higher K results better hyperparameters and likely more confident estimates resulting in comparable process. A takeaway for practitioners in the unconstrained case is that reducing K below d can be advantageous. Furthermore, in the unconstrained case, employing the formulation in (9) (or equivalently setting 8, 0.5 for B-SUB) can be sufficient.
Param,
Dest Value
16
0.1
0.2
0.3
0.4
0.5
Param
by
Figure 10: Ablation for on uncon-
strained out-of-model functions.
limited progress towards the local
E
Numerical considerations
Ensuring positive definiteness The Hessian of the Lagrangian as described in (8) may become indefinite due to numerical issues, modeling inaccuracies, or nonconvexity in the surrogate models. To maintain numerical stability and ensure that eurvature information defines valid descent direction, we enforce positive definiteness through a simple eigenvalue modification. Concretely, let & RRxx denote the Hessian candidate. We form the spectral decomposition orthogonal (QQQ Q !) and A = diag(Ai..... An) contains the real eigenvalues ordered arbitrarily. The gigenvalue-clipping rule replaces each eigenvalue A by , max(A;, :) for a small threshold
QAQ where Q is
0 which we set 10 The modified matrix is then reconstructed as H QÅQ where
26

---

= diag(Xi..... A.)). By construction X; 2 >>> 0 for all i, hence H is symmetric positive definite. Note that more sophisticated modifications are also possible, but found that this simple approach already resulted in satisfactory performance.
Jitter on the joint covariance The joint covariance of the standard and derivative GP can be ill-conditioned. Here, we apply standard jitter to the diagonal if necessary, which is also standard for covariances in classic BO. This then assures that the Cholesky decomposition for B-SUB exists.
Ensure feasibility through slacked B-SUB formulation After every sub-sampling step, we aim to solve the B-SUB optimization problem. However, given the problem and the current linearization, the sub-problen composed of the surrogate model estimates may be infeasible. We therefore opt to solve the following cack-constrained version of B-SUB:
P:
arg min (ey, Asche
1 2P
H.P + PV, + py + 11-55 + P Σ 84
(25)
iElm
subject to
S by,
L
T C2
beg
E mm:
2
Mv + 6, bc.
- S Pep
Vi E Im:
by 2 0,
bc. 2 0,
8, 2 0,
Vi Im.
where by and Le, are Cholesky factorizations as
0 Evar
L, = EVil Evi
L D
E T
(26)
Eve
and p > 0 is the penalty parameter for slack variables. Here, we choose p = 100. By design, this subproblem is always feasible. With the search direction from this slacked optimization problem, we proceed as described in the main part the paper. It should further noted that the feasibility the problem will also depend the dy and &g. We therefore recommend that if the subproblem frequently fails to increase 5 and de and potentially use a form of scheduling for these hyperparameters.
27

---

Runtime breakdown of BayeSQP
To give an idea of the computational efficiency BayeSQP, we provide a runtime comparison against TuRBO across varying problem dimensions in Table 3. Overall, BayeSQP demonstrates a substantial reduction in total wall-clock time relative to TuRBO, particularly in higher-dimensional settings as also demonstrated in Figure 5b and Figure 8. The reported results are from the within-mon model setting, but training surrogate models incurs approximately the same computational cost per hyperparameter update for both methods since BayeSQP operates only with the marginal GP. Due to its sub-sampling strategy, BayeSQP requires fewer model training iterations within the same computational budget.
Table 3: Runtime comparison of BayeSQP to TuRBO across dimensions. The time of BayeSQP which unaccounted for is due logging overhead. A detailed per-step runtime breakdown is in Table 4.
Dimension
BayeSQP
SOCP
Hessian
Subsampling (6)
8.93+2.69
2.92+1.71
0.06:0.07 (2.2%)
0.13=0.06 (4.6%)
0.03+0.01 (1.0%)
1.4240.21 (48.5%)
8
13.74=2.49
1.431..75
0.06+0.07 (2.6%)
0.11=0.05 (4.4%)
0.02+0.01 (1.0%)
0.94+0.16 16 (38.7%)
16
27.39a5.60
2.74+1.75
0.08:0.07 (2.8%)
0.14 (5.3%)
0.0410.01 (1.4%)
1.1740.15 (42.7%)
32
39.2424.18
2.80#1.71
0.08:0.07 (2.9%)
0.14=0.07 (5.1%)
0.10:0.08 (3.6%)
1.10#0.16 (39.3%)
64
86.02=16.24
2,98+1.65
0.09=0.06 (3.1%)
0.22=0.08 (7.2%)
0.12+0.06 (4.0%)
1.21=0.16 (40.5%)
96
310.52486.77
6.7242.30
0.2640.06 (3.9%)
0.53 (7.8%)
0.14+0.08 (2.1%)
1.6940.23 (25.1%)
To provide deeper insight into the computational characteristics each core component of BayeSQP. analyze the per-step runtime costs in Table 4. This breakdown demonstrates how the computational burden shifts as problem dimensionality increases, with Thompson sampling rmaining the most expensive component but showing decreasing relative contribution in higher dimensions as the Hessian computation becomes more expensive. With an increased number of constraints, the contribution of evaluating Hessians will also further increase linearly in the number constraints.
Table 4: Runtime breakdown per BayeSQP step.
Dimension
SOCP (s/step)
Hessian (s/step)
Subsampling (s/step)
TS (s/step)
0.0043 (3.9%)
0.0090 (8.1%)
0.0019 (1.7%)
0.0954 (86.3%)
0.0065 (5.6%)
0.0110 (9.5%)
0.0025 (2.2%)
0.0959 (82.7%)
0.0072 (5.3%)
0.0137 10)
0.0035 (2.6%)
0.1111 (82.0%)
32
0.0087 (5.7%)
0.0154 (10.0%)
0.0108 (7.0%)
0.1187 (77.3%)
64
0.0100 (5.6%)
0.0236 (13.2%)
0.0130 (7.3%)
0.1323 (74.0%)
96
0.0258 (10.0%)
0.0521 (20.1%)
0.0142 (5.5%)
0.1671 (64.5%)
28


Reviews:
[{"review_id": "pTw0YmNc81", "reviewer": "Reviewer_Br7S", "summary": "The authors present BayesSQP, which combines sequential quadratic programming with Bayesian optimization to solve constrained optimization problems. BayesSQP capitalizes on the relationship between a Gaussian Process (GP) and its derivatives to estimate the gradient and curvature of functions using zeroth-order information (function evaluations). In doing so, the QP subproblem that needs to be solved is uncertainty-aware because it uses function approximations (GPs) that naturally express uncertainty. The method is compared with TuRBO and logEI in a simulation study and using variants of the Ackley and Hartmann test functions", "questions": "1. Is TuRBO incorrectly labeled as SCBO on line 96?\n\n2. There is substantial work on scalable BO which is not discussed and never compared with. For example [1] and [3]. While the focus is on constrained optimization, there is a section with unconstrained problems, so more baselines need to be included.\n\n3. I didn't see the assumptions on the objective and constraints in the main text. Is this included somewhere and if not can the authors please include this information?\n\n4.It would be nice to see more baselines included to know how much benefit one gets from incorporating this over existing, potentially simpler methods. For the unconstrained problems, why not compare with methods assuming an additive structure to the function [1], SAASBO [2], or random projections [3]?\n\n5. What is the computational complexity of the algorithm? This is not explicitly stated in the text as far as I can see. \n\nIf my questions 2-5 were addressed sufficiently, e.g. extra baselines included and complexity concerns thoroughly addressed, then I would be willing to increase my score. If I have missed anything or misunderstood any components please let me know and I will reevaluate. \n\n[1] Kandasamy, Kirthevasan, Jeff Schneider, and Barnabás Póczos. \"High dimensional Bayesian optimisation and bandits via additive models.\" International conference on machine learning. PMLR, 2015.\n\n[2] Eriksson, David, and Martin Jankowiak. \"High-dimensional Bayesian optimization with sparse axis-aligned subspaces.\" Uncertainty in Artificial Intelligence. PMLR, 2021.\n\n[3] Wang, Ziyu, et al. \"Bayesian optimization in a billion dimensions via random embeddings.\" Journal of Artificial Intelligence Research 55 (2016): 361-387.", "limitations": "I think the limitations of the method in terms of assumptions and scalability concerns should be fleshed out and I have included details on this in my questions.", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths:\n- Polished presentation. The figures and organization are good. Figures 1 and 2, in particular, are well-structured and convey their information efficiently.\n- The explanation of the method follows a clear structure, and the core idea is presented in a good format.\n- The simulation study (Figure 4) does an excellent job of sweeping over the problem dimensions, and the results are impressive in this setting.\n\nWeaknesses:\n- The exploration of high-dimensional Bayesian Optimization solutions is lacking and should be fleshed out. See question 2 for details.\n- Assumptions need to be stated more clearly. See question 3 for details.\n- Scalability concerns. See question 5 for details.", "was_revised": true, "final_justification": "The authors directly addressed my concerns 1-5 so I have raised my score to border line accept. I feel the paper could still use more comparisons—that's why my final score is not 5.", "rebuttal": "We thank the reviewer for the comments and time taken for the review. We address their questions and stated weaknesses below.\n\n> Q1: Is TuRBO incorrectly labeled as SCBO on line 96?\n\nYes, thank you very much for catching this! We have fixed this typo in the document.\n\n> Q2 (and Weakness 1): There is substantial work on scalable BO which is not discussed and never compared with. For example [1] and [3]. While the focus is on constrained optimization, there is a section with unconstrained problems, so more baselines need to be included.\n\nThank you for pointing out these additional references (including [2]). We agree that there is a large body of work on making BO more scalable. In our submitted version, we focused primarily on local methods—a subset of scalable BO approaches—since BayeSQP extends prior work that leverages first-order optimization techniques. However, your comment highlights that we should expand our related work discussion. We therefore include [1-3] and subsequent approaches in the related work section, adding a paragraph on scalable BO methods before discussing the local approaches. We also provide explicit experimental comparisons by incorporating an additional baseline from this set of methods to our experiments.\n\nHere, we decided to include SAASBO in our unconstrained experiments as a popular baseline in the sparse subspace setting. We find that it outperforms BayeSQP up to 8D, performs similarly on 16D, and afterwards BayeSQP also outperforms this baseline. Thank you for suggesting this baseline in Q4.\n\n> Q3 (and Weakness 2): I didn't see the assumptions on the objective and constraints in the main text. Is this included somewhere and if not can the authors please include this information?\n\nNo, you are correct that we did not explicitly formulate assumptions on the objective and constraints in our submitted version. While BayeSQP does not impose explicit structural assumptions (such as sparse subspaces) on the objective and constraints, our modeling approach does embed implicit assumptions about these functions. The sub-problem is based on the assumption that the modeling approach of the objective with a GP is reasonable. Since we rely also on derivatives of the GP, severe model misspecification could potentially harm performance more than other approaches which we in part see on the Ackley experiments. We will make sure to make this more explicit in the final version. Specifically, we will add the following to Section 2: “Here and in the following, we will assume that $f(x)$ and $c(x)$ are samples from Gaussian processes and we will base our algorithm as well as modelling approach on this assumption”. We will further add a small discussion to the existing one on limitations stating that a severe model mismatch could harm the performance of BayeSQP.\n\n> Q4: It would be nice to see more baselines included to know how much benefit one gets from incorporating this over existing, potentially simpler methods. For the unconstrained problems, why not compare with methods assuming an additive structure to the function [1], SAASBO [2], or random projections [3]?\n\nAs discussed above, we agree that a discussion on these methods was missing and we now have added it to the manuscript. We furthermore agree that a comparison to one of the methods mentioned above would benefit the current presentation. Therefore, we included SAASBO as a representative recent method from this category. SAASBO matched or exceeded BayeSQP's performance up to 16 dimensions, but BayeSQP significantly outperformed SAASBO at 32 dimensions. For dimensions 64 and 96, SAASBO could not complete runs within our 24-hour computational budget. This computational limitation has also been noted in prior evaluations of SAASBO [4].\n\n> Q5 (and Weakness 3): What is the computational complexity of the algorithm? This is not explicitly stated in the text as far as I can see.\n\nThe computational complexity with respect to the number of data points follows standard BO at $O(N^3)$ for Gram matrix inversion. Since this inverse is cached, it can be reused for evaluating both the gradient GP and the Hessian mean. The training time is identical to a standard GP as the log marginal likelihood objective only operates on the marginal standard GP. Evaluating the mean of the Hessian at a test point scales as $O(N^2 d^2)$. The $d^2$ factor arises from computing all d×d Hessian elements, while $N^2$ comes from the matrix operations with the cached Gram matrix inverse. The optimization subproblem involves solving a second-order cone program (SOCP) with $d+1+m$ decision variables (recall that m is the number of constraints), which has polynomial complexity independent of the number of data points. For our problem scale, CVXOPT provides adequate performance for solving the optimization subproblem, as demonstrated in the runtime results below.\n\nOn that note, we want to highlight that we found a more efficient implementation to evaluate the Hessian leveraging Einsum notation vs. iteratively building up the Hessian. With this, we were able to significantly reduce runtime by eliminating redundant computations and leveraging vectorized tensor operations. The optimized implementation achieves almost a $d^2$ speed up as it avoids nested loops and we can now report runtimes on, e.g., the 96D within-model functions setting of TuRBO (WM: 310sec, OHO: 12 206sec), logEI (WM: 186sec, OHO: 1743 sec), BayeSQP (WM: 7sec, OHO: 110sec) where WM is pure within-model, i.e., perfect model knowledge and OHO is with online hyperparameter optimization. As stated above, hyperparameter optimization is similar for all methods.\n\n---\n### References\n\n[1] Kandasamy, Kirthevasan, Jeff Schneider, and Barnabás Póczos. \"High dimensional Bayesian optimisation and bandits via additive models.\" International conference on machine learning. PMLR, 2015. \n\n[2] Eriksson, David, and Martin Jankowiak. \"High-dimensional Bayesian optimization with sparse axis-aligned subspaces.\" Uncertainty in Artificial Intelligence. PMLR, 2021.\n\n[3] Wang, Ziyu, et al. \"Bayesian optimization in a billion dimensions via random embeddings.\" Journal of Artificial Intelligence Research 55 (2016): 361-387.\n\n[4] Papenmeier, Leonard, Luigi Nardi, and Matthias Poloczek. \"Increasing the scope as you learn: Adaptive Bayesian optimization in nested subspaces.\" Advances in Neural Information Processing Systems 35 (2022): 11586-11601."}, {"review_id": "9kR26eoYs6", "reviewer": "Reviewer_r1jx", "summary": "They propose BayeSQP, a method that uses classical Sequential Quadratic Programming (SQP) ideas in Bayesian Optimization for high-dimensional problems, with or without constraints. To apply SQP in BO, they require a second-order Gaussian Process, i.e., a surrogate that returns the mean, variance, gradient, and Hessian, but uses only zero-order evaluations. This GP is used to approximate the gradient and Hessian at each iteration, which allows them to construct a robust quadratic subproblem. In this subproblem they apply Value-at-Risk to the objective and chance-constraints to the restrictions, and reformulate it as a second-order cone program (SOCP). They also implement a fallback strategy with a slack variable that guarantees initial feasibility of the subproblem. The SQP subproblem they solve indicates the direction to move in iteration p_t. Then, they determine how far to move in that direction by performing a constrained line search along the segment using Thompson Sampling with a budget of M=3 evaluations, and selecting the point with the highest feasibility. Finally, in several synthetic experiments, they demonstrate the performance of their method.", "questions": "Questions:\n\n - Could BayeSQP's strong performance simply come from its initial feasibility guarantee rather than the core algorithm itself?\n\n - Why does logEI converge quickly and perform poorly in high-dimensional problems?\n\n - How many initial points were used for logEI in these experiments?\n\nOther comments:\n\n - On page 3, there is a typo: \"Approaches such as SCBO [11] and SCBO [10]” (one of those should be a different method reference).\n\n - On page 5, Eq. 9: \"E [ ∇c_i (x_t )^⊤ p ≥ −c_i (x_t ) ]\" is incorrect because the expectation should not include the inequality.\n\n - On page 5, Eq. 11: \"µ^⊤_{∇c_i} (x_t)\" should be: \"µ_{∇c_i} (x_t)^⊤\".\n\n - On page 6: “on a step size α **which** which we can update”.\n\n - At the end of page 7: there is a “p_k” but everywhere else the iteration index is “p_t”.\n\n - Table 1: bold highlight the best value.", "limitations": "I have described the limitations I found in the Weaknesses section.", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strenghts:\n\n - Their method scales effectively to high-dimensional problems, with experiments demonstrated up to around 100 dimensions\n\n - Their method can handle both unconstrained and constrained problems\n\n - They formulate an uncertainty aware SQP subproblem and integrate it into a BO algorithm.\n\n - They implement a fallback mechanism that guarantees feasible evaluations in the first iteration.\n\nWeaknesses:\n\n - The performance of BayeSQP may come from guaranteeing feasible points at the start rather than from the rest of the pieces of its algorithm. It would be useful to see an ablation study without the fallback to isolate its effect.\n\n - I do not understand why logEI converges so quickly and performs poorly as dimensionality increases (its performance in d=16 and d=32 is almost the same). In the within-model experiments, fixing the lengthscale of all GPs at 0.1 might not be appropriate for logEI (even if it is the best for RFF). That small scale could limit its exploration. It would be better to follow logEI's original BoTorch configuration [4] (which is the one used in logEI [1]). Moreover, it would be better to use a Matern-5/2 ARD kernel, the one used in logEI [1].\n\n - The number of initial points may be affecting logEI. How many were used?\n\n- I think the number of repetitions, 20, is insufficient. They should use at least 30, and preferably 50 or 100, as is common in the literature.\n\n\n[1] Sebastian Ament, Samuel Daulton, David Eriksson, Maximilian Balandat, and Eytan Bakshy. Unexpected improvements to expected improvement for Bayesian optimization. In Advances in Neural Information Processing Systems (NeurIPS), 2023.\n\n[4] Maximilian Balandat, Brian Karrer, Daniel R. Jiang, Samuel Daulton, Benjamin Letham, Andrew Gordon Wilson, and Eytan Bakshy. BoTorch: a framework for efficient Monte-Carlo Bayesian optimization. In Advances in Neural Information Processing Systems (NeurIPS), 2020.", "was_revised": true, "final_justification": "I think that the authors have adequately addressed my comments (initial feasibility mechanism, logEI's high-dimensional performance, number of repetitions) and the important concerns the other reviewers pointed out (inclusion of SAASBO and MPD baselines, new ablation studies, new details about complexity and scalability, and a new real-world experiment). If the authors implement all of the changes mentioned, I consider it appropriate to raise the score I originally gave them.", "rebuttal": "We thank the reviewer for the comments and time taken for the review. We address their questions and stated weaknesses below.\n\n> Q1 (and W1): Could BayeSQP's strong performance simply come from its initial feasibility guarantee rather than the core algorithm itself?\n\nThank you for the question but we think there may be a misunderstanding. We think our explanation of feasibility may not have been sufficiently clear, as there are actually two distinct types of feasibility at play in BayeSQP that we should distinguish better. First, there is the feasibility of a sample point x_t with respect to the original problem constraints. Second, there is the feasibility of the SOCP subproblem that BayeSQP solves at each iteration. Upon reviewing our manuscript, we recognize that our terminology was not consistently precise in distinguishing these concepts, which we will address in a final version. To directly address the question: BayeSQP implements a slacked version of the SOCP subproblem to ensure computational feasibility of the subproblem itself, not to guarantee that the resulting sample points satisfy the original problem constraints. Therefore, BayeSQP's performance cannot be attributed to an initial feasibility guarantee for the sample points. In fact, we often have multiple iterations at the start, e.g., on the constrained Hartmann, where the sample points remain infeasible similar to the other algorithms.\n\n> Q2 (and W2): Why does logEI converge quickly and perform poorly in high-dimensional problems?\n\nLogEI's poor performance in high-dimensional settings stems primarily from its global search approach. Essentially, as nicely stated in [1], when the true objective function has small lengthscales relative to domain, global search becomes increasingly ineffective as dimensionality grows. Given a small length scale for the unit hypercube and a large number of dimensions, a global search is bound to fail to a certain extent. Therefore, as problem complexity increases, logEI's behavior increasingly resembles random search which becomes particularly evident in the \"pure\" within-model setting as also highlighted in prior work (e.g., [2]). \n\nTo investigate whether enforcing larger lengthscales could mitigate this issue, we conducted an experiment where we initialized the lengthscales as sqrt(d) [3] on the 64D within model objective. We can observe that this does improve performance slightly however, both TuRBO and BayeSQP are able to also outperform this baseline with BayeSQP showing the strongest performance.\n\n> Q3: How many initial points were used for logEI in these experiments?\n\nFor all experiments and all baselines, we use D initial points. Such a scaling of the number of initial points with the number of dimensions is quite common.\n\n> W1: The performance of BayeSQP may come from guaranteeing feasible points at the start rather than from the rest of the pieces of its algorithm. It would be useful to see an ablation study without the fallback to isolate its effect.\n\nAs discussed above, we believe there might have been misunderstanding about which feasibility was meant in the submitted version when discussing the slacked version of the SOCP. We adjust the text to make the distinction more precise. If we however misunderstood the point you were making, please let us know so that we can further clarify.\n\n> W2: I think the number of repetitions, 20, is insufficient. They should use at least 30, and preferably 50 or 100, as is common in the literature.\n\nWe appreciate this feedback and have increased the number of repetitions from 20 to 32, following the logEI paper [3]. We are happy to increase this number even further for the final version but opted for 32 due to the limited time frame for the rebuttal. Importantly, our results remain consistent with the results reported in the submitted version of the paper.\n\nAlso thank you for spotting these typos and the recommendation for Table 1! We have fixed the typos and included bold highlighting in the manuscript!\n\n---\n### References\n\n[1] Wüthrich, Manuel, Bernhard Schölkopf, and Andreas Krause. \"Regret bounds for Gaussian-process optimization in large domains.\" Advances in Neural Information Processing Systems 34 (2021): 7385-7396.\n\n[2] Müller, Sarah, Alexander von Rohr, and Sebastian Trimpe. \"Local policy search with Bayesian optimization.\" Advances in Neural Information Processing Systems 34 (2021): 20708-20720.\n\n[3] Ament, Sebastian, et al. \"Unexpected improvements to expected improvement for bayesian optimization.\" Advances in Neural Information Processing Systems 36 (2023): 20577-20612."}, {"review_id": "w0XbmR02Ig", "reviewer": "Reviewer_DWmx", "summary": "This paper introduces BayeSQP, a novel approach that blends ideas from Sequential Quadratic Programming (SQP) and Bayesian Optimization (BO) to tackle constrained black-box optimization problems. The core idea is to build joint Gaussian Process (GP) surrogates that model not just the function values, but also their gradients and Hessians, estimated purely from zero-order observations. At each iteration, BayeSQP constructs a local subproblem that incorporates model uncertainty and solves it as a second-order cone program (SOCP) to determine a promising search direction. A subsequent constrained Thompson sampling line search is used to select the next evaluation point. The approach is tested on synthetic problems, showing a competive performance.", "questions": "1. *Line search budget and efficiency:* The method uses constrained Thompson sampling with a fixed number of samples (M=3) for line search. Could the authors elaborate on how sensitive BayeSQP is to this choice? Would a larger budget significantly improve performance?\n\n2. *Alternative uses of uncertainty:* The paper leverages uncertainty using a Value-at-Risk (VaR) formulation, which is both intuitive and computationally tractable. However, one could imagine alternative approaches, such as using Conditional Value-at-Risk (CVaR) or distributionally robust optimization formulations. Have the authors considered or tested such alternatives?\n\n3. *Guidance on uncertainty parameters:* Although the authors include an illustrative analysis of the role of $\\delta_f$ and $\\delta_c$, a more systematic sensitivity study would be helpful. Providing practical recommendations for setting these parameters could make the method more accessible to non-expert users.\n\n4. *Clarification on problem formulation* In Equation (1), the use of the expectation $\\mathbb{E}[f(x)]$ may be misleading, especially since the rest of the paper assumes a deterministic function with noisy observations. Could the authors clarify the meaning of this expectation, or consider reformulating the objective if the method is not intended to handle true stochastic optimization problems?", "limitations": "Yes.", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 2, "originality": 3, "strengths_and_weaknesses": "**Strenghts:**\n\n1. *Technical depth and novelty:* The paper presents a well-thought-out integration of classical second-order optimization with Bayesian optimization. To my knowledge, this is the first work to take this approach in a probabilistic SQP framework.\n\n2. *Strong empirical performance:* BayeSQP outperforms strong baselines like TuRBO and logEI on several high-dimensional benchmarks, particularly when constraints are involved.\n\n**Weaknesses:**\n\n1. *Lack of real-world benchmarks:* The experiments are well-executed, but they’re limited to synthetic functions. Including one or two real-world applications would help demonstrate the practical impact and usability of the method.\n\n2. *Baseline selection should be broader:* The comparison focuses on established baselines like TuRBO and logEI, but doesn’t include some newer methods for local BO, such as the algorithm proposed by Nguyen et al. (2022). A broader comparison would help position BayeSQP more clearly in the current landscape.", "was_revised": true, "final_justification": "I believe this paper makes a significant contribution. Moreover, the authors were able to address most of my concerns during the discussion period, so I have decided to raise my score from 4 to 5.", "rebuttal": "We thank the reviewer for the comments and time taken for the review. We address their questions and stated weaknesses below.\n\n> Q1: The method uses constrained Thompson sampling with a fixed number of samples (M=3) for line search. Could the authors elaborate on how sensitive BayeSQP is to this choice? Would a larger budget significantly improve performance?\n\nThank you for the interesting question. We performed a sweep over different values with $M \\in \\\\{ 1, …, 10\\\\}$ to test the influence on the final performance on the 64D constrained within model objective functions. We find that for a within-model setting the median performance does not significantly change but by increasing $M$ the 90% upper quantile of the performance can be reduced. In general, we hypothesize that there is a deep connection of $M$ and the number of subsampling steps $K$. We therefore conducted a second sweep over combinations (similar to Figure 5c) with $K \\in \\\\{4, 8, 16, 32, 64, 128, 256\\\\}$ and found that one can reduce the number of subsampling steps when slightly increasing $M$ for a “pure” within-model (WM) setting, i.e., with perfect model knowledge. We performed the same sweep also with online hyperparameter optimization (OHO). Here, we find that increasing $M$ can improve performance but decreasing $K$ comes at a higher cost as a higher number of subsamples facilitates better learning of hyperparameters. Increasing $K$ beyond D also does not improve performance as it reduces the number of line searches given a fixed budget. That said, using D+1 subsampling points is a compelling heuristic as it offers a clean theoretical justification, as discussed in the paper, however, if more inductive biases are available one can probably reduce this number. We will include these additional sensitivity ablations in the main body of the paper. \n\n> Q2: However, one could imagine alternative approaches, such as using Conditional Value-at-Risk (CVaR) or distributionally robust optimization formulations. Have the authors considered or tested such alternatives?\nWe have not yet experimented with CVaR or distributionally robust formulations, but we agree that this is a promising direction. We see BayeSQP as a flexible framework that could benefit from various extensions and refinements, including improvements to the subproblem formulation. In particular, incorporating risk-sensitive objectives like CVaR or potentially leveraging approximate Hessian uncertainties for more informed subsampling are exciting directions which seem promising for future research. We hope that this initial version of BayeSQP can serve as a foundation for further research in these directions.\n\n> Q3: Although the authors include an illustrative analysis of the role of  and , a more systematic sensitivity study would be helpful. Providing practical recommendations for setting these parameters could make the method more accessible to non-expert users.\n\nThank you for that comment. With the illustrative example, we wanted to provide intuition on the behavior of BayeSQP for different confidence parameters and then highlight with Figure 5 c) that the exact choice of the delta does not seem to be critical as long as some delta_c < 0.5 is specified. Still, we think that we can further improve in explicit recommendations for practitioners as well as researchers. We believe that the above mentioned results will help with this and we will make sure to state suggestions more explicitly in a final version of the manuscript.\n\n> Q4: Could the authors clarify the meaning of this expectation, or consider reformulating the objective if the method is not intended to handle true stochastic optimization problems?\n\nThis was a mistake from our side and thank you for spotting this! We have fixed this in the manuscript. To clarify:  we assume a fixed f from which we can only obtain noisy zero-order observations.\n\n> W1: The experiments are well-executed, but they’re limited to synthetic functions. Including one or two real-world applications would help demonstrate the practical impact and usability of the method.\n\nWe agree that such an application would further strengthen the contributions of our paper. Unfortunately, we were not able to set up a new experiment in the limited time of the rebuttal and rather focused on adding baselines (MPD as requested in W2 and SAASBO requested by reviewer Br7S to which we also refer to for a longer discussion on this baseline) as well as conducting experiments that further provide insights into the behavior of BayeSQP as discussed in Q1. We also added additional out-of-model comparisons (similar to the mentioned experiments in Q1) to highlight BayeSQPs performance also in the OHO setting. In these experiments, we observe a slight drop in performance, which is to be expected, but the rankings among the algorithms remains consistent with BayeSQP showing the best performance in the high-dimensional setting.\n\n> W2: The comparison focuses on established baselines like TuRBO and logEI, but doesn’t include some newer methods for local BO, such as the algorithm proposed by Nguyen et al. (2022). A broader comparison would help position BayeSQP more clearly in the current landscape.\n\nThank you for the comment and suggestion. We agree that including more recent local BO methods would better position BayeSQP within the current landscape and can improve the empirical section of our paper. Adding the method from the paper above is especially interesting as, in a sense, BayeSQP is an extension of this class of first-order BO methods to a second order method. We therefore added MPD (method from the mentioned paper). We find that BayeSQP outperforms this baseline in all experiments with dimension greater than 8 and we furthermore observe that in the high-dimensional experiments, also TuRBO outperformed this baseline. We should also note that extending MPD (and similar methods like GIBO) to constrained optimization settings is not straightforward."}, {"review_id": "aGrZJnLmZB", "reviewer": "Reviewer_khbc", "summary": "This paper proposes BayeSQP, a method for constrained black-box optimization that unifies Sequential Quadratic Programming (SQP) with Bayesian Optimization (BO). BayeSQP models the objective and constraints using GPs that jointly encode the function value, gradient, and Hessian — all inferred from zeroth-order (function value only) observations. Each iteration constructs an uncertainty-aware subproblem, formulated as a second-order cone program (SOCP) using posterior means and covariances. A constrained Thompson sampling line search follows to select the next evaluation point. Empirical results show BayeSQP performs well in high-dimensional constrained settings (up to 96D), consistently returning feasible solutions when baselines fail.", "questions": "1. Can you provide runtime breakdowns (e.g., % time in SOCP solve vs. GP inference vs. sampling)?\n\n2. Would first-order-only BayeSQP (without Hessians) perform similarly? An ablation would clarify the benefit of the full model.\n\n3. Have you considered using initialization heuristics or multi-start strategies to mitigate local-only behavior?\n\n4. Why not include at least one CEC-style constrained benchmark, even if surrogate mismatch is expected, to demonstrate robustness?\n\n5. Could you expand on when the uncertainty-aware formulation (δ < 0.5) is most beneficial — e.g., high noise, narrow constraints, ill-conditioning?", "limitations": "Yes, the authors have explicitly discussed several limitations in Section 6. However, one key limitation not explicitly discussed is the lack of runtime profiling. While wall-clock performance is reported, the paper does not isolate the cost of SOCP solving, GP inference, or Hessian computation — which could become dominant at scale. Adding such profiling would better characterize the method’s practical bottlenecks.", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 4, "originality": 4, "strengths_and_weaknesses": "Strengths:\n\nS1. The proposed method tightly couples second-order classical optimization and probabilistic modeling. Unlike first-order local BO methods like GIBO, BayeSQP explicitly leverages the structure of the Hessian and models uncertainty in all components — a rare and nontrivial contribution in BO.\n\nS2. The subproblem is carefully reformulated to incorporate uncertainty using value-at-risk (VaR) and probabilistic constraint satisfaction, yielding an SOCP solvable via standard solvers. This derivation is mathematically clean and computationally viable.\n\nS3. Experiments show that BayeSQP outperforms TuRBO and C-logEI in synthetic constrained GP benchmarks, especially in high dimensions. Unlike C-logEI, it consistently finds feasible solutions across all runs.\n\nS4. The method includes fallback strategies (e.g., slack-variable variants for infeasibility), Sobol-based local resampling, and extensive consideration of hyperparameter robustness (e.g., δf, δc), showing good awareness of practical deployment.\n\nWeaknesses:\n\nW1. BayeSQP is a fully local method, relying on an initial point and never exploring globally. While this is common in high-dimensional BO, it risks getting trapped in poor local optima, especially in multimodal landscapes (e.g., Ackley5D). No hybrid exploration scheme (e.g., TuRBO-like restarts or multi-start) is employed.\n\nW2. While Figure 5b reports wall-clock times for constrained within-model benchmarks, the runtime data has limited generality:\n* It covers only a single problem class (synthetic GP-compatible functions with smooth constraints), and is not averaged across diverse benchmark types such as Ackley.\n* It does not break down runtime by component (e.g., GP posterior inference, SOCP solving, line search).\n\nW3. The use of within-model benchmarks is justified to study surrogate performance in ideal conditions, but limits generalizability. While standard BO functions like Ackley and Hartmann are included later (Table 1), the paper does not include more realistic constrained benchmarks, such as those from the CEC competitions. These benchmarks — although often non-smooth — are standard in constrained global optimization and would better demonstrate real-world feasibility and robustness, even if only as a stress test.\nThat said, the choice to avoid CEC is understandable: these benchmarks often violate the assumptions of GP-based modeling (smoothness, differentiability), and the authors state that second-order modeling would not be meaningful in such cases. Still, some acknowledgment of this limitation in the main text would help balance the paper’s scope.\n\nW4. The paper would benefit from additional empirical studies isolating the contribution of second-order information. How much does using the Hessian improve performance over gradient-only or function-only variants? Is the uncertainty-aware formulation (δ < 0.5) materially better than the expectation-only one (δ = 0.5)? These comparisons are essential to justify the extra complexity of modeling second-order structure.\n\nW5. While this is common in BO, the paper lacks any theoretical discussion of convergence, regret, or feasibility rates under noise. Given the probabilistic nature of the surrogate and the robustness of the SOCP formulation, it would be valuable to include even informal intuition about when the method is guaranteed to make progress.", "was_revised": true, "final_justification": "My concerns have all been resolved, so I’ll keep my score as it is.", "rebuttal": "We thank the reviewer for the comments and time taken for the review. We address their questions and stated weaknesses below.\n\n> Q1 (and W2): Can you provide runtime breakdowns (e.g., % time in SOCP solve vs. GP inference vs. sampling)?\n\nWe provide runtime breakdowns for BayeSQP on within-model comparisons across different dimensions below. Note that we identified a more efficient Hessian evaluation using Einsum notation, which significantly reduces runtime by avoiding nested loops. Accordingly, an updated Figure 5b reflects these improved runtimes.\n| Dimension | TuRBO (s) | BayeSQP (s) | SOCP (s) | Hessian eval (s) | Subsampling (s) | TS (s) |\n|-----------|----------|---------------------|----------|-------------|-----------------|--------|\n| 4 | 8.93±2.69 | 2.92±1.71 | 0.06±0.07 (2.2%) | 0.13±0.06 (4.6%) | 0.03±0.01 (1.0%) | 1.42±0.21 (48.5%) |\n| 8 | 13.74±2.49 | 2.43±1.75 | 0.06±0.07 (2.6%) | 0.11±0.05 (4.4%) | 0.02±0.01 (1.0%) | 0.94±0.16 (38.7%) |\n| 16 | 27.39±5.60 | 2.74±1.75 | 0.08±0.07 (2.8%) | 0.14±0.10 (5.3%) | 0.04±0.01 (1.4%) | 1.17±0.15 (42.7%) |\n| 32 | 39.24±4.18 | 2.80±1.71 | 0.08±0.07 (2.9%) | 0.14±0.07 (5.1%) | 0.10±0.08 (3.6%) | 1.10±0.16 (39.3%) |\n| 64 | 86.02±16.24 | 2.98±1.65 | 0.09±0.06 (3.1%) | 0.22±0.08 (7.2%) | 0.12±0.06 (4.0%) | 1.21±0.16 (40.5%) |\n| 96 | 310.52±86.77 | 6.72±2.30 | 0.26±0.06 (3.9%) | 0.53±0.33 (7.8%) | 0.14±0.08 (2.1%) | 1.69±0.23 (25.1%) |\n\nFrom the results above, we can observe that the active sampling-based line search (TS) is the most compute intensive part of BayeSQP, followed by the building the Hessian, and lastly solving the SOCP subproblem. The remaining runtime is attributed to logging overhead (note that with only 3sec of total runtime already 1.5sec of logging are quite significant) and initializing the GP models which we reinitialize in each iteration following various BoTorch tutorials. Maintaining the same model over iterations and appending the data should further decrease runtime. The time for initialization is however similar for all algorithms as it all is done through the same constructors from BoTorch. Also the time for logging is similar across methods. We will include this runtime breakdown in a similar form to Figure 5b in the paper. However, for this figure, we will report the time per BayeSQP steps. Note that above, for dimension 4, we have a longer runtime than for dimension 8 because the number of substeps is larger. The table below shows the runtime breakdown per substep. The percentages are based on the sum of the core components.\n\n\n| Dimension | SOCP (s/step) | Hessian eval (s/step) | Subsampling (s/step) | TS (s/step) |\n|-----------|---------------|------------------|---------------------|-------------|\n| 4 | 0.0043 (3.9%) | 0.0090 (8.1%) | 0.0019 (1.7%) | 0.0954 (86.3%) |\n| 8 | 0.0065 (5.6%) | 0.0110 (9.5%) | 0.0025 (2.2%) | 0.0959 (82.7%) |\n| 16 | 0.0072 (5.3%) | 0.0137 (10.1%) | 0.0035 (2.6%) | 0.1111 (82.0%) |\n| 32 | 0.0087 (5.7%) | 0.0154 (10.0%) | 0.0108 (7.0%) | 0.1187 (77.3%) |\n| 64 | 0.0100 (5.6%) | 0.0236 (13.2%) | 0.0130 (7.3%) | 0.1323 (74.0%) |\n| 96 | 0.0258 (10.0%) | 0.0521 (20.1%) | 0.0142 (5.5%) | 0.1671 (64.5%) |\n\nThe per-step analysis reveals expected dimensional scaling patterns that provide insights into algorithmic behavior. While Hessian evaluation time increases with dimension (from 8.1% to 20.1%), absolute computation times remain efficient even at higher dimensions (0.05 seconds at 96 dimensions). We appreciate this suggestion, as the breakdown provides valuable insights into our algorithm's computational characteristics of our core components.\n\nNote that learning hyperparameters adds the same runtime for all baselines as the log marginal likelihood objective only operates on the marginal standard GP. Upon request of reviewer Br7S we now also have added a discussion on the complexity of the different parts of the algorithm.\n\n> Q2: Would first-order-only BayeSQP (without Hessians) perform similarly? An ablation would clarify the benefit of the full model.\n\nThis is an interesting question and follows along the lines of the request by reviewer DWmx in W2. In the revised version, we have now added an additional first-order BO method, namely MPD [1], which has the substeps of learning the current gradient and then performing gradient descent using the predicted gradients. To an extent, BayeSQP is a second-order extension of this method. We can observe that after 8 dimensions, BayeSQP also outperforms this baseline leveraging the second-order approximation in a SQP style.\n\n> Q3: Have you considered using initialization heuristics or multi-start strategies to mitigate local-only behavior?\n\nWe have not explored initialization heuristics or multi-starts but see it as an exciting future direction of research. In general, we believe that BayeSQP can be also be applied in a batched setting using strategies similar to TuRBO [2], which uses for each entry in a batch a separate trust region.\n\n> Q4 (and W3): Why not include at least one CEC-style constrained benchmark, even if surrogate mismatch is expected, to demonstrate robustness?\n\nThank you for the comment. By including Ackley, we aimed to demonstrate our method's limitations on such “adversarial” objective functions for purely local approaches. For this revision, we prioritized implementing two additional baselines (SAASBO as requested by reviewer Br7S, and MPD discussed above) as well as performing ablations to further provide insights into the behavior of BayeSQP. As you note, the CEC-style functions are often stated in a way in which these local approaches tend to yield not meaningful results. We will follow your advice and further expand the discussion of BayeSQP's limitations in such settings within the main paper. \n\n> Q5: Could you expand on when the uncertainty-aware formulation (δ < 0.5) is most beneficial — e.g., high noise, narrow constraints, ill-conditioning?\n\nWe see the benefits of our uncertainty-aware formulation in two different mechanisms. First, incorporating uncertainty typically improves optimization behavior by providing more reliable search directions when gradient estimates are uncertain. Second, the choice of $\\delta$ will have a direct influence on the robustness of the final solution. The magnitude of $\\delta$ will likely depend on the application at hand and will mostly influence the final optimal solution. If, e.g., one aims to find a robust constraint satisfying solution one can directly encode this through the choice of very small $\\delta$. By design, it will yield a solution that has some safety margin. This touches also on the high-noise setting, where we naturally want some margin to the constraint. With Figure 3, we aimed to provide some qualitative insights on the influence of delta.\n\n---\n### References\n\n[1] Nguyen, Quan, et al. \"Local Bayesian optimization via maximizing probability of descent.\" Advances in neural information processing systems 35 (2022): 13190-13202.\n\n[2] Eriksson, David, et al. \"Scalable global optimization via local Bayesian optimization.\" Advances in neural information processing systems 32 (2019)."}]

Author rebuttal:

Meta-review: The paper presents a novel algorithm that combines Sequential Quadratic Programming (SQP) with Bayesian Optimization for sample-efficient black-box optimization problems. The key scientific contribution is a novel formulation of SQP in an uncertainty-aware fashion that leverages second-order Gaussian Process surrogates that jointly model function values, gradients, and Hessians using only zero-order observations. The method constructs uncertainty-aware subproblems formulated as second-order cone programs (SOCPs), incorporating Value-at-Risk for objectives and chance constraints for restrictions. The line search to select the next evaluation point is performed via Thompson sampling. Empirical results demonstrate improved performance against strong TuRBO and logEI baselines on synthetic benchmarks, with particular strength in constrained high-dimensional settings.

Strengths:
- Motivation and practical relevance
- Novelty of the methodological approach
- SOCP formulation is a strong technical contribution
- Natural inclusion of constraints
- Strong optimization performance (esp. In higher dimensions, and in constrained settings)
- Clarity of presentation (writing, exposition, polished figures)

Weaknesses
- Empirical evaluation largely limited to synthetic problems (one real-world problem added in rebuttal)
- Local optimization behavior makes method reliant on initialization
- Scalability to (very-)high-dimensional problems is limited

The primary concerns raised by the reviewers were the omission of other scalable Bayesian optimization baselines and real-world benchmarks, as well as the insufficient discussion of some assumptions and experimental details. In their rebuttal, authors comprehensively addressed these concerns by providing comparisons against SAASBO and MPD as additional baselines, additional ablation studies and runtime comparison, and discussion of complexity and scalability. All reviewers were satisfied with the authors’ responses and there were no major outstanding concerns. 

The paper makes a substantial contribution to the area of sample-efficient black-box optimization. The technical contribution is significant - successfully bridging classical second-order optimization with modern Bayesian optimization represents a meaningful advance for the constrained optimization community. The uncertainty-aware SOCP formulation is mathematically elegant, computationally tractable, and naturally allows for the inclusion of constraints. While the evaluation focuses on synthetic problems, the consistent strong performance across multiple high-dimensional settings and the addition of baselines like SAASBO and MPD during rebuttal strengthen the empirical evidence for the value of this method. Finally, the presentation is clear and polished and the key limitations are discussed explicitly. 

I believe this is an excellent paper - especially if the authors include the additional baselines and discussions produced during the rebuttal period.
Decision: Accept (spotlight)