1 Introduction↩︎

In recent years, Bayesian optimization (BO) has emerged as a powerful framework for black-box optimization ranging from applications in robotics [1]–[3] to hyperparameter tuning [4], [5] and drug discovery [6]–[8]. To address high-dimensional problems emerging in these fields, a variety of high-dimensional BO (HDBO) approaches have been proposed, including the use of local BO (LBO) methods [9], [10] or methods that exploit specific structure in the objective [11]. Recently, a 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 [12], [13]. However, as shown by [14], 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 a matter of choosing either HDBO approaches or standard approaches, but rather of leveraging recent advances in 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 [10], [15]–[18]—which can be interpreted as combining BO with first-order optimization methods—to LBO with second-order methods. We introduce BayeSQP, a 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 1). At each iteration, a local subproblem is constructed using the GP posterior estimates and solved to 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:

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

2. An uncertainty-aware subproblem for BayeSQP that accounts for the variance and covariance in function and gradient estimates, resulting in a tractable second-order cone program.

3. Empirical experiments demonstrating that BayeSQP outperforms state-of-the-art BO methods in specific high-dimensional constrained settings.

2 Problem formulation↩︎

We consider the problem of finding an optimizer to the general non-convex optimization problem \[\begin{align} \label{eq:opt95problem} {\boldsymbol{x}}^* = \mathop{\mathrm{arg\,min}}_{{\boldsymbol{x}}\in {\mathbb{X}}} f({\boldsymbol{x}}) \quad \text{subject to} \quad c_i({\boldsymbol{x}}) \geq 0, \quad \forall i \in {\mathbb{I}}_m \mathrel{\vcenter{:}}= \{1, \dots, m\} \end{align}\tag{1}\] where \(f: {\mathbb{X}}\to {\mathbb{R}}\) and constraints \(c_i : {\mathbb{X}}\to {\mathbb{R}}\) for all \(i \in {\mathbb{I}}_m\) are black-box functions defined over the compact set \({\mathbb{X}}\subset {\mathbb{R}}^d\). At each iteration \(t \in {\mathbb{I}}_T\) where \(T\) is the total budget for the optimization, an algorithm selects a query point \({\boldsymbol{x}}_t \in {\mathbb{X}}\) and receives noisy zeroth-order feedback following the standard observation model in BO as \(f_t = f({\boldsymbol{x}}_t) + \varepsilon_f\) for the objective, and \(c_{i,t} = c_i({\boldsymbol{x}}_t) + \varepsilon_{c_i}\) for all \(i \in {\mathbb{I}}_m\) for the constraints where \(\varepsilon_f\) and \(\varepsilon_{c_i}\) are independent realizations from a zero-mean Gaussian distribution with possibly different noise variances. From these observations, we construct independent datasets for the objective function \({\mathcal{D}}_{f}^t = \{({\boldsymbol{x}}_j, f_j)\}_{j=1}^{t}\) and for each constraint \({\mathcal{D}}_{c_i}^t = \{({\boldsymbol{x}}_j, c_{i,j})\}_{j=1}^{t}\) for all \(i \in {\mathbb{I}}_m\), which any zero-order method can leverage to solve 1 .

3 Preliminaries↩︎

3.1 Sequential quadratic programming↩︎

SQP represents a 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 a wide range of optimization problems. The foundation of constrained optimization rests on the Lagrangian function, defined as \(\mathcal{L}({\boldsymbol{x}}, \boldsymbol{\xi}) = f({\boldsymbol{x}}) - \sum_{i=1}^m \xi_i c_i({\boldsymbol{x}})\). It combines the objective with the constraints, where each constraint is weighted by its Lagrange multiplier \(\xi_i\). 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 \({\boldsymbol{H}}_t = \nabla_{{\boldsymbol{x}}{\boldsymbol{x}}}^2\mathcal{L}({\boldsymbol{x}}_t, \boldsymbol{\xi})\) (or an appropriate approximation thereof) and linearizes the constraints around the current point \({\boldsymbol{x}}_t\). This generates the following subproblem: \[\label{eq:standard95sqp} \begin{align} {\boldsymbol{p}}_{t} = \mathop{\mathrm{arg\,min}}_{{\boldsymbol{p}}\in \mathbb{R}^d} \quad & \frac{1}{2}{\boldsymbol{p}}^{\top}{\boldsymbol{H}}_t {\boldsymbol{p}}+ \nabla f({\boldsymbol{x}}_t)^{\top}{\boldsymbol{p}}+ f({\boldsymbol{x}}_t) \\ \text{subject to} \quad & \nabla c_i({\boldsymbol{x}}_t)^{\top}{\boldsymbol{p}}\geq -c_i({\boldsymbol{x}}_t), \,\, \forall i \in {\mathbb{I}}_m \end{align}\tag{2}\] The solution \({\boldsymbol{p}}_t\) provides a search direction, and the next iterate is typically computed as \({\boldsymbol{x}}_{t+1} = {\boldsymbol{x}}_t + \alpha_t {\boldsymbol{p}}_t\), where \(\alpha_t\) is 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 [19].

Under standard assumptions, SQP exhibits local superlinear convergence when using exact Hessian information, and various quasi-Newton approximation schemes (such as BFGS or SR1 updates, cf. [19]) 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-parametric way. A GP is defined as a collection of random variables, any finite number of which have a joint Gaussian distribution. Formally, a GP is defined by its mean function \(m({\boldsymbol{x}}) \mathrel{\vcenter{:}}= \mathbb{E}[f({\boldsymbol{x}})]\) and kernel \(k({\boldsymbol{x}}, {\boldsymbol{x}}') \mathrel{\vcenter{:}}= \mathrm{Cov}[f({\boldsymbol{x}}), f({\boldsymbol{x}}')]\) [20]. In BayeSQP, we use GPs to model the objective function \(f\) and the constraints \(c_i\) as standard in BO [21]. Contrary to standard BO, we aim to leverage the following property of GPs: They are closed under linear operations, i.e., the derivative of a GP is again a GP given that the kernel is sufficiently smooth [20]. This enables us to derive a distribution for the gradient and Hessian. We can formulate the following joint prior distribution: \[\begin{align} \begin{bmatrix} {\boldsymbol{y}}\\ f \\ \nabla f \\ \nabla^2 f \end{bmatrix} & \sim \mathcal{N}\left( \begin{bmatrix} m({\boldsymbol{X}}) \\ m({\boldsymbol{x}}) \\ \nabla m({\boldsymbol{x}}) \\ \nabla^2 m({\boldsymbol{x}}) \end{bmatrix}, \begin{bmatrix} k({\boldsymbol{X}}, {\boldsymbol{X}}) + \sigma^2 {\boldsymbol{I}}& \bullet & \bullet & \bullet \\ k({\boldsymbol{x}}, {\boldsymbol{X}}) & k({\boldsymbol{x}}, {\boldsymbol{x}}) & \bullet & \bullet \\ \nabla k({\boldsymbol{x}}, {\boldsymbol{X}}) & \nabla k({\boldsymbol{x}}, {\boldsymbol{x}}) & \nabla^2 k({\boldsymbol{x}}, {\boldsymbol{x}}) & \bullet \\ \nabla^2 k({\boldsymbol{x}}, {\boldsymbol{X}}) & \nabla^2 k({\boldsymbol{x}}, {\boldsymbol{x}}) & \nabla^3 k({\boldsymbol{x}}, {\boldsymbol{x}}) & \nabla^4 k({\boldsymbol{x}}, {\boldsymbol{x}}) \end{bmatrix} \right) \end{align}\] where \({\boldsymbol{y}}\in {\mathbb{R}}^n\) are the \(n\) function observations, \({\boldsymbol{X}}= [{\boldsymbol{x}}_1, \ldots, {\boldsymbol{x}}_n] \in {\mathbb{R}}^{d \times n}\) is the matrix of all training inputs with \({\boldsymbol{x}}_i \in {\mathbb{R}}^d\), and \({\boldsymbol{x}}\in {\mathbb{R}}^d\) is the test point.¹ Here and in the following, we use \(\bullet\) to denote symmetric entries for improved readability. The joint conditional distribution is then:\[\begin{bmatrix} f \\ \nabla f \\ \nabla^2 f \end{bmatrix} \ifcsname fenced@\string|\endcsname \@firstoftwo \else \@secondoftwo \fi {\amsmath@bigm\csname fenced@\string|\endcsname}{\bigm|} {\boldsymbol{x}}, {\boldsymbol{X}}, {\boldsymbol{y}} \sim \mathcal{N}\left( \begin{bmatrix} \mu_f({\boldsymbol{x}}) \\ {\boldsymbol{\mu}}_{\nabla f}({\boldsymbol{x}}) \\ {\boldsymbol{\mu}}_{\nabla^2 f}({\boldsymbol{x}}) \end{bmatrix}, \begin{bmatrix} \sigma_{f}^2({\boldsymbol{x}}) & \bullet & \bullet \\ {\boldsymbol{\Sigma}}_{\nabla f, f}({\boldsymbol{x}}) & {\boldsymbol{\Sigma}}_{\nabla f}({\boldsymbol{x}}) & \bullet \\ {\boldsymbol{\Sigma}}_{\nabla^2 f, f}({\boldsymbol{x}}) & {\boldsymbol{\Sigma}}_{\nabla^2 f, \nabla f}({\boldsymbol{x}}) & {\boldsymbol{\Sigma}}_{\nabla^2 f}({\boldsymbol{x}}) \end{bmatrix} \right).\] Following standard conditioning of multivariate normal distributions, we can directly compute the mean and covariance functions of the marginals of the posterior as \[\begin{align} \text{\textit{Marginal GP of }} f: & \,\, \left\{ \begin{array}{ll} \mu_f({\boldsymbol{x}}) = m({\boldsymbol{x}}) + k({\boldsymbol{x}}, {\boldsymbol{X}}){\boldsymbol{K}}^{-1}({\boldsymbol{y}}- {\boldsymbol{m}}({\boldsymbol{X}})) & \in \mathbb{R}, \\ [4pt] \sigma_f^2({\boldsymbol{x}}) = k({\boldsymbol{x}}, {\boldsymbol{x}}) - k({\boldsymbol{x}}, {\boldsymbol{X}}){\boldsymbol{K}}^{-1}k({\boldsymbol{X}}, {\boldsymbol{x}}) & \in \mathbb{R} \end{array} \right. \\[4pt] \text{\textit{Marginal GP of }} \nabla f: & \,\, \left\{ \begin{array}{ll} {\boldsymbol{\mu}}_{\nabla f}({\boldsymbol{x}}) = \nabla m({\boldsymbol{x}}) + \nabla k({\boldsymbol{x}}, {\boldsymbol{X}}){\boldsymbol{K}}^{-1}({\boldsymbol{y}}- {\boldsymbol{m}}({\boldsymbol{X}})) & \in \mathbb{R}^d, \\ [4pt] {\boldsymbol{\Sigma}}_{\nabla f}({\boldsymbol{x}}) = \nabla^2 k({\boldsymbol{x}}, {\boldsymbol{x}}) - \nabla k({\boldsymbol{x}}, {\boldsymbol{X}}){\boldsymbol{K}}^{-1}\nabla k({\boldsymbol{X}}, {\boldsymbol{x}}) & \in \mathbb{R}^{d \times d} \end{array} \right. \end{align}\] Here, we defined the Gram matrix as \({\boldsymbol{K}}\mathrel{\vcenter{:}}= k({\boldsymbol{X}}, {\boldsymbol{X}}) + \sigma^2 {\boldsymbol{I}}\) with entries \([k({\boldsymbol{X}}, {\boldsymbol{X}})]_{ij} = k({\boldsymbol{x}}_i, {\boldsymbol{x}}_j)\) for \(i,j \in {\mathbb{I}}_n\), and use the notation that \(k({\boldsymbol{X}}, {\boldsymbol{x}}) \in {\mathbb{R}}^{n \times 1}\) is the vector of kernel evaluations between each training point and the test point, with entries \([k({\boldsymbol{X}}, {\boldsymbol{x}})]_i = k({\boldsymbol{x}}_i, {\boldsymbol{x}})\). Similarly, we can obtain the covariance term \({\boldsymbol{\Sigma}}_{\nabla f, f}({\boldsymbol{x}})\) as well as the mean estimate of the Hessian.² As noted in [10], we must perform the inversion of the Gram matrix \({\boldsymbol{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 set size.

3.3 Related work↩︎

3.3.0.1 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 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 [11], [22]–[24]. 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 [12]–[14]. 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 [25]–[27] or sparse GPs [28], [29] have been discussed; addressing these scalability issues, however, is not the focus of this work.

3.3.0.2 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 [9] and SCBO [30] 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, [10] 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 [15], [17], [18], [31], theoretically investigated [16], and extended with crash constraints [32]. 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 a GP. Similar ideas have been leveraged in a quasi-Newton methods [33]. However, by incorporating ideas from SQP, BayeSQP is directly applicable to both unconstrained and constrained optimization problems—something which is not possible with GIBO.

3.3.0.3 Bayesian optimization and Gaussian processes in classical optimization

There have been various papers integrating BO with first-order optimization, e.g., for line search [34], [35]. GPs have been successfully applied and leveraged in optimization—both for local optimization [36], [37] and global optimization (essentially BO) [21], [38]. All of these can be classified as a subfield of probabilistic numerics [39], [40]. Similar to our approach, [41] merged classical methods with BO by lifting the constraints into the objective using an augmented Lagrangian approach which later got extended to a slacked [42] and recently a relaxed version [43]. These approaches are based on expected improvement (EI) and, crucially, [30] showed that these approaches do not scale well to high-dimensional problems. BayeSQP 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.

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 a single model. BayeSQP then leverages this model at each iteration to construct a quadratic uncertainty-aware subproblem for a search direction that yields improvement with high probability. In the following, we will first discuss our modeling approach. Based on this, we will construct the subproblem, followed by a discussion on line search. In the end, we touch on further practical extensions and give intuition on the optimization behavior.

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: \[\label{eq:joint95gp} \begin{bmatrix} f \\ \nabla f \\ \mathrm{vec}(\nabla^2 f) \end{bmatrix} \ifcsname fenced@\string|\endcsname \@firstoftwo \else \@secondoftwo \fi {\amsmath@bigm\csname fenced@\string|\endcsname}{\bigm|} {\boldsymbol{x}}, {\boldsymbol{X}}, {\boldsymbol{y}} \sim \mathcal{N}\left( \begin{bmatrix} \mu_f({\boldsymbol{x}}) \\ {\boldsymbol{\mu}}_{\nabla f}({\boldsymbol{x}}) \\ \mathrm{vec}({\boldsymbol{\mu}}_{\nabla^2 f}({\boldsymbol{x}})) \end{bmatrix}, \begin{bmatrix} \sigma_{f}^2({\boldsymbol{x}}) & \bullet & \times \\ {\boldsymbol{\Sigma}}_{\nabla f, f}({\boldsymbol{x}}) & {\boldsymbol{\Sigma}}_{\nabla f}({\boldsymbol{x}}) & \times \\ \times & \times & \times \end{bmatrix} \right)\tag{3}\] We use surrogate models of the same form for each constraint \(c_i({\boldsymbol{x}})\). We do not compute the covariance of the Hessian (\(\times\)) 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-order information.

Crucially, it is not required to always evaluate the full posterior distribution for each test point. In a SQP framework, we can approximate the Hessian of the Lagrangian once at our current iterate as \[\label{eq:lag} {\boldsymbol{H}}_t = {\boldsymbol{\mu}}_{\nabla^2 f}({\boldsymbol{x}}_t) - \textstyle\sum_{i=1}^m \xi_i^{(t-1)} {\boldsymbol{\mu}}_{\nabla^2 c_i}({\boldsymbol{x}}_t),\tag{4}\] where \(\xi_i^{(t-1)}\) are the Lagrange multipliers from the solution of the last subproblem, but for the subsequent line search, we can directly work with the cheap marginal GP \(f\sim \mathcal{N} (\mu_f({\boldsymbol{x}}), \sigma_{f}^2({\boldsymbol{x}}))\).

4.2 Deriving the subproblem for BayeSQP↩︎

Standard SQP approaches typically require exact knowledge of the objective function, constraints, and their respective gradients. In our case, we only have access to zero-order feedback and the question arises how to formulate a suitable subproblem given our choice of surrogate model.

4.2.0.1 Expected value SQP subproblem

A straightforward approach is to simply formulate a subproblem using expectations, leading to the following expected value subproblem: \[\begin{align} {\boldsymbol{p}}_t \in \mathop{\mathrm{arg\,min}}_{{\boldsymbol{p}}\in \mathbb{R}^d} \quad & \mathbb{E}\left[\frac{1}{2}{\boldsymbol{p}}^{\top}{\boldsymbol{H}}_t{\boldsymbol{p}}+ \nabla f({\boldsymbol{x}}_t)^{\top}{\boldsymbol{p}}+ f({\boldsymbol{x}}_t)\right] \\ \text{subject to} \quad & \mathbb{E}\left[c_i({\boldsymbol{x}}_t) + \nabla c_i({\boldsymbol{x}}_t)^{\top}{\boldsymbol{p}}\right]\geq 0, \,\, \forall i \in {\mathbb{I}}_m \, . \end{align} \label{eq:eqp}\tag{5}\] While intuitive, this formulation fails to account for the inherent uncertainty in the estimates. As discussed by [15] and [18], taking into account the uncertainty of, e.g., the gradient, can be crucial for improving with high probability for LBO approaches.

4.2.0.2 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: \[\label{eq:robust95sqp} \begin{align} {\boldsymbol{p}}_t \in \mathop{\mathrm{arg\,min}}_{{\boldsymbol{p}}\in \mathbb{R}^d} \quad & \underbracket{\mathrm{VaR}_{1-\delta_f}\left[\frac{1}{2}{\boldsymbol{p}}^{\top}{\boldsymbol{H}}_t{\boldsymbol{p}}+ \nabla f({\boldsymbol{x}}_t)^{\top}{\boldsymbol{p}}+f({\boldsymbol{x}}_t)\right]}_{\text{Objective value-at-risk with confidence level }1-\delta_f} \\ \text{subject to} \quad & \underbracket{\mathbb{P}\left(c_i({\boldsymbol{x}}_t) + \nabla c_i({\boldsymbol{x}}_t)^{\top}{\boldsymbol{p}}\geq 0\right) \geq 1-\delta_c}_{\text{Constraint satisfaction with confidence }1-\delta_c}, \,\, \forall i \in {\mathbb{I}}_m \, . \end{align}\tag{6}\] This formulation accounts for uncertainty through two mechanisms: employing value-at-risk (VaR) for the objective function and enforcing probabilistic feasibility for the constraints. The resulting search direction minimizes the worst-case objective value while ensuring the constraints are satisfied with high probability.

4.2.0.3 Tractability through joint Gaussian process

The robust formulation in 6 remains intractable without distributional assumptions. By modeling the objective and constraints as jointly Gaussian with their gradients, we can transform 6 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 \(\mathbb{P}\left( c_i({\boldsymbol{x}}_t) + \nabla c_i({\boldsymbol{x}}_t)^{\top}{\boldsymbol{p}}\geq 0 \right) \geq 1-\delta_c\). Since we have \({\boldsymbol{z}}^\top {\boldsymbol{v}}\sim \mathcal{N}({\boldsymbol{\mu}}_z^\top {\boldsymbol{v}}, {\boldsymbol{v}}^\top {\boldsymbol{\Sigma}}_z {\boldsymbol{v}})\) for a multivariate Gaussian random variable \({\boldsymbol{z}}\sim \mathcal{N}({\boldsymbol{\mu}}_z, {\boldsymbol{\Sigma}}_z)\) and \({\boldsymbol{v}}\) is a deterministic vector, we know that \(c_i({\boldsymbol{x}}_t) + \nabla c_i({\boldsymbol{x}}_t)^{\top}{\boldsymbol{p}}\) is also normal distributed with moments \[\begin{align} \mathbb{E}\left[c_i({\boldsymbol{x}}_t) + \nabla c_i({\boldsymbol{x}}_t)^{\top}{\boldsymbol{p}}\right] & = \mu_{c_i}({\boldsymbol{x}}_t) + \boldsymbol{\mu}_{\nabla c_i}^{\top}({\boldsymbol{x}}_t) \,{\boldsymbol{p}}\\ \text{Var}\left[c_i({\boldsymbol{x}}_t) + \nabla c_i({\boldsymbol{x}}_t)^{\top}{\boldsymbol{p}}\right] & = \sigma_{c_i}^2({\boldsymbol{x}}_t) \, + {\boldsymbol{p}}^{\top}\boldsymbol{{\boldsymbol{\Sigma}}}_{\nabla c_i}({\boldsymbol{x}}_t) \,{\boldsymbol{p}}+ 2{\boldsymbol{p}}^{\top}\boldsymbol{{\boldsymbol{\Sigma}}}_{c_i,\nabla c_i}({\boldsymbol{x}}_t) \, \end{align}\] where the last term accounts for the covariance between the function and its gradient. In the following, we drop the explicit evaluation at \({\boldsymbol{x}}_t\) for all moments for notational convenience, i.e., \(\mu_{c_i} = \mu_{c_i}({\boldsymbol{x}}_t)\).

For a Gaussian random variable to remain non-negative with probability at least \(1 - \delta\), we require its mean to exceed its standard deviation multiplied by the corresponding quantile. This yields: \[\begin{align} \mu_{c_i} + \boldsymbol{\mu}_{\nabla c_i}^{\top} \,{\boldsymbol{p}}& \geq q_{1-\delta} \sqrt{\sigma_{c_i}^2 + {\boldsymbol{p}}^{\top} \boldsymbol{{\boldsymbol{\Sigma}}}_{\nabla c_i}\, {\boldsymbol{p}}+ 2{\boldsymbol{p}}^{\top} \boldsymbol{{\boldsymbol{\Sigma}}}_{c_i,\nabla c_i}} \end{align}\] where \(q_{1-\delta} = \Phi^{-1}(1 - \delta)\) denotes the \((1 - \delta)\)-quantile of the standard normal distribution. Rearranging the terms and introducing an auxiliary variable \(t_{c_i}\) to upper-bound the square root term allows the constraint to be reformulated as a set of two inequalities: \[\begin{align} -\boldsymbol{\mu}_{\nabla c_i}^{\top} \,{\boldsymbol{p}}+ q_{1-\delta}b_{c_i} \leq \mu_{c_i} \quad \text{and} \quad \sqrt{\sigma_{c_i}^2 + {\boldsymbol{p}}^{\top}\boldsymbol{{\boldsymbol{\Sigma}}}_{\nabla c_i} \,{\boldsymbol{p}}+ 2{\boldsymbol{p}}^{\top}\boldsymbol{{\boldsymbol{\Sigma}}}_{c_i,\nabla c_i} \,} \leq b_{c_i} \end{align}\] To express the square root term more compactly, we consider the full covariance matrix associated with the joint Gaussian distribution of \(c_i\) and its gradient \(\nabla c_i\). Specifically, we can state \[\begin{align} \sigma_{c_i}^2 + {\boldsymbol{p}}^{\top}\boldsymbol{{\boldsymbol{\Sigma}}}_{\nabla c_i} \,{\boldsymbol{p}}+ 2{\boldsymbol{p}}^{\top}\boldsymbol{{\boldsymbol{\Sigma}}}_{c_i,\nabla c_i} \, & = \begin{bmatrix} 1 \\ {\boldsymbol{p}}\end{bmatrix}^{\top} \begin{bmatrix} \sigma_{c_i}^2 & \bullet \\ \boldsymbol{{\boldsymbol{\Sigma}}}_{c_i,\nabla c_i} & \boldsymbol{{\boldsymbol{\Sigma}}}_{\nabla c_i} \end{bmatrix} \begin{bmatrix} 1 \\ {\boldsymbol{p}}\end{bmatrix} \end{align}\] By Cholesky decomposition of the covariance matrix, we can express the square root term as a second-order cone constraint: \[\begin{align} \sqrt{\sigma_{c_i}^2 + {\boldsymbol{p}}^{\top}\boldsymbol{{\boldsymbol{\Sigma}}}_{\nabla c_i} \,{\boldsymbol{p}}+ 2{\boldsymbol{p}}^{\top}\boldsymbol{{\boldsymbol{\Sigma}}}_{c_i,\nabla c_i} \,} & = \sqrt{\begin{bmatrix} 1 \\ {\boldsymbol{p}}\end{bmatrix}^{\top} {\boldsymbol{L}}_{c_i}{\boldsymbol{L}}_{c_i}^{\top} \begin{bmatrix} 1 \\ {\boldsymbol{p}}\end{bmatrix}} = \left\|{\boldsymbol{L}}_{c_i}^{\top}\begin{bmatrix} 1 \\ {\boldsymbol{p}}\end{bmatrix}\right\|_2 \leq b_{c_i} \end{align}\] Using the same reasoning, we can reformulate the objective function by introducing the auxiliary variable \(b_f\). In the end, we obtain the following formulation that we refer to as \(\boldsymbol{B-SUB}\).

This formulation also naturally incorporates the subproblem formulation in 5 .

4.3 Line search through constrained posterior sampling↩︎

With the search direction given as the solution of the \(\boldsymbol{B-SUB}\) subproblem, the next step is to decide on a step size \(\alpha\) which which we can update the current iterate as \({\boldsymbol{x}}_{t+1} = {\boldsymbol{x}}_t + \alpha_t {\boldsymbol{p}}_t\). To implicitly decide on the step size, we perform constrained posterior sampling [30] on the one-dimensional line segment spanned by \({\boldsymbol{p}}_t\). Specifically, we aim to solve \[\mathop{\mathrm{arg\,min}}_{\left\{ {\boldsymbol{x}}_t + \alpha {\boldsymbol{p}}_t \;\middle|\; \alpha \in [0, 1] \right\}} f({\boldsymbol{x}}) \quad \text{subject to} \quad c_i({\boldsymbol{x}}) \geq 0, \,\, \forall i \in {\mathbb{I}}_m.\] This is similar to LineBO [44] 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 \(\alpha_t\) that yields progress given a limited evaluation budget \(M\) for the line search which we set to 3 in all experiments. Similar to [30], we either choose the next point to be the index of the best feasible point, or, if none of the points are feasible, as the point with the least amount of constraint violations as \[\label{eq:best95point} {\boldsymbol{x}}_{k+1} \gets \begin{cases} \mathop{\mathrm{arg\,min}}_{{\boldsymbol{x}}_t^{(j)} \in \mathcal{F}} f({\boldsymbol{x}}_t^{(j)}), & \text{if } \mathcal{F} \neq \emptyset, \\ \mathop{\mathrm{arg\,min}}_{1 \leq j \leq M} \sum_{i \in {\mathbb{I}}_m} \max\left(0, -c_i\left({\boldsymbol{x}}_t^{(j)}\right)\right), & \text{otherwise}, \end{cases}\tag{7}\] where \(\mathcal{F} = \left\{ {\boldsymbol{x}}_t^{(j)} \middle| \, c_i\left({\boldsymbol{x}}_t^{(j)}\right) \geq 0, \,\forall i \in {\mathbb{I}}_m \right\}\) denotes the set of feasible points among the \(M\) samples.

4.4 Practical considerations and intuition on optimization behavior↩︎

4.4.0.1 Local sub-sampling

Unlike GIBO-style methods [10], [15], [18], 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 a \(d\)-dimensional ball of radius \(\varepsilon\) centered at \({\boldsymbol{x}}_t \in \mathbb{R}^d\). For this, we first draw a Sobol sequence from the hypercube \([0,1]^{d+1}\). Each Sobol point \((\tilde{{\boldsymbol{x}}}, u) \in [0,1]^{d} \times [0,1]\) is then transformed such that \(\tilde{{\boldsymbol{x}}}\) approximates a standard normal vector to yield a unit direction \(\bar{\boldsymbol{x}}\), and \(u\) determines the individual radius as \(r = \varepsilon \cdot u^{1/d}\). The final sample is then \({\boldsymbol{x}}= {\boldsymbol{x}}_t + r \cdot \bar{\boldsymbol{x}}\).

4.4.0.2 Slack variable fallback strategy

The subproblem \(\boldsymbol{B-SUB}\) may become infeasible due to constraint linearization or high uncertainty in gradient estimates. To address this, we implement a slack variable version of \(\boldsymbol{B-SUB}\) as a fallback, which guarantees feasibility by design (cf. Appendix 12). This approach aligns with established practices in classical SQP methods [19]. 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.

4.4.0.3 Intuition on optimization behavior

To gain intuition about the parameters \(\delta_f\) and \(\delta_c\) and their influence on the optimization process, we study BayeSQP on a small toy example. We generate a two-dimensional within-model objective function (cf. Appendix 8) with a quadratic constraint, resulting in only a 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 (\(\delta_f, \delta_c = 0.5\)), the linearization of the quadratic constraint results in tangential directions \(p_k\), leading to limited or no improvement. We observe that incorporating uncertainty into the subproblem pushes the search direction toward the feasible set. Additionally, selecting a very low value for \(\delta_c\) effectively robustifies the constraints, as shown by the resulting directions \(p_k\).

5 Empirical evaluations↩︎

We next quantitatively evaluate our proposed method BayeSQP. Our evaluation first considers unconstrained and then constrained optimization problems using BoTorch [45]. We benchmark against four baselines: logarithmic EI (logEI) [38], [46], TuRBO [9], SAASBO [11], and MPD [15]. These baselines are widely used [7], [47]–[49] and represent complementary approaches—logEI employs a classic global optimization strategy, TuRBO implements a pseudo-local approach, SAASBO aims to automatically identify and exploit low-dimensional structure within high-dimensional search spaces through a hierarchical sparsity prior, and MPD is a fully local BO approach. Additionally, logEI and TuRBO can be readily adapted for constrained optimization through their respective variants: C-logEI [46], [50], [51] and SCBO [30] to which we compare on the constrained optimization problems.

In all subsequent plots, we present the median alongside the \(5^\text{th}\) to \(95^\text{th}\) percentile range (\(90\%\) inner quantiles) computed across 32 independent random seeds. For BayeSQP, we set the hyperparameters \(\delta_f, \delta_c = 0.2\) (unless stated otherwise) and \(K = d + 1\), following [16].

5.0.0.1 Unconstrained optimization

We first consider unconstrained within-model problems [52] for which we adapt \(\boldsymbol{B-SUB}\) accordingly. We generate the functions using random Fourier features following [53], [54] (cf. Appendix 8 for all details). Optimizing such functions has gained relevance with recent advances in latent space BO [7], [55], [56], where GP priors are enforced in the latent space [57]. Figure 4 summarizes the results. BayeSQP outperforms the other baselines from dimension 16 onward. Furthermore, we can observe the step-like behavior of BayeSQP resulting from the subsampling followed by solving \(\boldsymbol{B-SUB}\) and the subsequent line search which yields the improvement.

5.0.0.2 Constrained optimization

Similarly, we can perform within-model comparisons for the constrained case. Here, also the constraint function \(c({\boldsymbol{x}})\) is a sample from an GP. Again, all details are provided in Appendix 8. Figure 5 summarizes the constrained within-model results. As in the unconstrained case, BayeSQP outperforms the baselines at high dimensions (Figure 5 (a)), while remaining orders of magnitude faster than SCBO and C-logEI despite computing full Hessians per \(\boldsymbol{B-SUB}\)(Figure 5 (b)). However, as we keep increasing dimensions, computing the Hessians of size \({d \times d}\) will results in a 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 13. Lastly, in Figure 5 (c) we can observe the influence of the parameters \(\delta_f\) and \(\delta_c\) of \(\boldsymbol{B-SUB}\) on the performance for \(d=64\). We can observe as visualized in Figure 3, not considering uncertainty especially in the constraints (\(\delta_c=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.

5.0.0.3 Performance on standard benchmarks

Lastly, we also evaluate BayeSQP on standard BO benchmarks. Here, we follow recent best practices and initialize lengthscales with \(\sqrt{d}\) for all baselines [12]–[14]. The results are summarized in Table 1. We can clearly observe that BayeSQP is sensitive to initialization highlighted by the large 90% quantile especially for Ackley. This is to be expected as the algorithm is 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.

r0.55

To demonstrate the real-world applicability of BayeSQP, we compare constrained optimization baselines on the 7-dimensional Speed Reducer benchmark [59], which minimizes the weight of a speed reducer subject to 11 mechanical design non-linear constraints (more details in Appendix 8.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, BayeSQP demonstrates a clear runtime advantage even in the presence of 11 constraints—each requiring separate Hessian evaluations—and a substantially larger \(\boldsymbol{B-SUB}\).

6 Discussion on limitations↩︎

While BayeSQP provides a novel framework combining classic optimization methods with BO, there are several limitations and addressing them will be interesting future research.

6.0.0.1 Initialization matters

As with any local approach, the initialization of BayeSQP will directly influence its performance (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 [41] in Figure 6 (details in Appendix 8.5). Depending on the initialization, the algorithm converges to a different local optimum of the constrained problem.

Although global approaches can also exhibit sensitivity to initialization, this sensitivity is amplified in LBO approaches, particularly in constrained optimization. However, this sensitivity provides practitioners with the option to incorporate some expert knowledge into the optimization 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.

6.0.0.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 \(\boldsymbol{B-SUB}\), for example following efficient schemes such as [33], [60]; whether this would lead to empirical performance improvements remains an open question. Future work could focus on 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 a BFGS-type update scheme.

6.0.0.3 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 \(\boldsymbol{B-SUB}\) directly relies on the accuracy of both gradient and Hessian estimates, a poorly chosen kernel can lead to unreliable curvature information and ultimately to suboptimal search directions. While standard kernels such as the squared-exponential 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., [12], [13]) are expected to directly improve the robustness and effectiveness of BayeSQP.

In Appendix 9, we list possible extensions of 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.

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 Forschungsgemeinschaft (DFG, German Research Foundation)–RTG 2236/2 (UnRAVeL). Simulations were performed in part with computing resources granted by RWTH Aachen University under projects rwth1579, p0022034, and p0021919.

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

8.1 Generating within-model objective functions↩︎

Within-model comparisons were introduced by [52] to study the performance of BO methods on functions that fullfil all model assumptions. With recent advances in latent space BO [57], 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 \(f_i\) with a \(M\) random Fourier features (RFFs) following [53]. This yields a parametric function \[f_i({\boldsymbol{x}}) = \sum_{m=1}^M w_m \phi_m({\boldsymbol{x}}) \quad \text{with} \quad \phi_m({\boldsymbol{x}}) = \sqrt{\frac{2}{M}} \cos(\boldsymbol{\theta}_m^\top {\boldsymbol{x}}+ \tau_m).\] Here, \(\boldsymbol{\theta}_m\) are sampled proportional to the kernel’s spectral density and \(\tau_m \sim \mathcal{U}(0, 2\pi)\). In all experiments, we use a squared-exponential kernel with lengthscales \(\ell_i = 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 it finds suitable hyperparameters for the given task (cf. Appendix 8.6 and cf. Appendix 11).

8.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 a within-model objective function, we also generate a within-model constraint function \(\hat{c}({\boldsymbol{x}})\). We then shift this function by one, i.e., \(c({\boldsymbol{x}}) = \hat{c}({\boldsymbol{x}}) -1 \geq 0\) so that on average only about 16% of the domain is feasible.³ Note that multiple constraints are also 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.

8.3 Constrained versions of Ackley and Hartmann↩︎

For the constrained versions of Ackley and Hartmann, we use the pre-implemented benchmarks in BoTorch [45] which are largely based on experiments in [58] and [30]. 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: \[\begin{align} \text{Ackley: } \left\{ \begin{array}{ll} c_1({\boldsymbol{x}}) & = -\sum_{i=1}^d {\boldsymbol{x}}_i , \\ [4pt] c_2({\boldsymbol{x}}) & = 5 - \| {\boldsymbol{x}}\|_2 \end{array} \right. \qquad \text{Hartmann: } \left\{ \begin{array}{ll} c_1({\boldsymbol{x}}) & = 1 - \| {\boldsymbol{x}}\|_2^2. \end{array} \right. \end{align}\] For Ackley, we restrict the feasible region to \([-5, 10]^d\) [9] and use a time horizon of \(T=100\) for \(d=5\) and \(T=400\) for \(d=20\). The feasible region for Hartmann is \([0, 1]^6\) and we set \(T=100\).

8.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 (\(x_1 \in [2.6, 3.6]\)), module of teeth (\(x_2 \in [0.7, 0.8]\)), number of teeth on pinion (\(x_3 \in [17, 28]\), integer which we treat as a continuous variable as implemented in BoTorch and consistent with prior work [9]), length of shaft 1 (\(x_4 \in [7.3, 8.3]\)), length of shaft 2 (\(x_5 \in [7.8, 8.3]\)), diameter of shaft 1 (\(x_6 \in [2.9, 3.9]\)), and diameter of shaft 2 (\(x_7 \in [5.0, 5.5]\)). We run the benchmark for 200 iterations across 32 random seeds with \(T=200\) and report the results in Table 2. For BayeSQP, we set \(\delta_f, \delta_c = 0.5\) essentially reverting to a 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 \(\delta_f, \delta_c\) over the number of iterations. For the full formulation of the optimization problem, we refer to the BoTorch documentation (v.13.0) as well as the original paper [59].

8.5 Constrained Gramacy function↩︎

For the constrained Gramacy benchmark, we use the formulation introduced in [41] and implemented in BoTorch. The problem is defined over the unit square \({\boldsymbol{x}}\in [0, 1]^2\) with the objective of minimizing the sum of the two decision variables, \(f({\boldsymbol{x}}) = x_1 + x_2\). It imposes two nonlinear inequality constraints \[\begin{align} c_1({\boldsymbol{x}}) &= -\big(1.5 - x_1 - 2x_2 - 0.5 \sin\big(2\pi(x_1^2 - 2x_2)\big)\big), \\ c_2({\boldsymbol{x}}) &= -(x_1^2 + x_2^2 - 1.5). \end{align}\] 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 a global optimum over time, approaches such as restarting BayeSQP are promising. For a longer discussion, we refer to Appendix 9.

8.6 Baselines↩︎

8.6.0.1 (C-)logEI

Expected improvement (EI) is a widely used acquisition function in BO. However, optimizing the EI acquisition function can be numerically unstable. To address this, [46] proposed a logarithmic transformation of EI resulting a numerically more stable acquisition function even resulting in an increase in performance. We use the implementation of logEI from BoTorch [45] and adapt it for constrained optimization by using the same wrapping on standard constrained EI [50], [51] resulting in C-logEI [46]. Neither logEI nor C-logEI require additional hyperparameters.

8.6.0.2 TuRBO and SCBO

For the implementation of TuRBO and SCBO, we follow tutorials from BoTorch [45] 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, we follow the recommendations of [9] in the mentioned tutorial. With these recommendations, the initial length of the trust region is \(L_{\text{init}} = 0.8\), the minimum and maximum length of the trust region are \(L_{\text{min}} = 0.5^7\) and \(L_{\text{max}} = 1.6\), respectively, the number of consecutive failures before the trust region is shrunk is \(\tau_{\text{fail}} = \lceil \max\{4, d\}\rceil\), the number of consecutive successes before the trust region is expanded is \(\tau_{\text{succ}} = 3\). The trust region is always centered around the best point found so far which in the context of constrained optimization follows 7 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 [9], [61], 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\}\) for a given candidate and dimension, and the value of the center of the trust region otherwise which induces an exploitation bias [14]. While TuRBO (and SCBO) depend on all the above parameters for their trust region and sampling heuristics, we found that these suggested parameters work well across various tasks as also highlighted in previous work [7], [47].

8.6.0.3 SAASBO

Sparse axis-aligned subspace Bayesian optimization (SAASBO) [11] is designed for high-dimensional BO by placing hierarchical sparsity priors on inverse lengthscales to identify and exploit low-dimensional structure. The method uses a global shrinkage parameter \(\tau \sim \mathcal{HC}(\beta)\) and dimension-specific inverse lengthscales \(\rho_d \sim \mathcal{HC}(\tau)\) for all \(d \in \mathbb{I}_d\), where \(\mathcal{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 BoTorch tutorial⁴ which performs inference using Hamiltonian Monte Carlo (HMC) with the NUTS sampler from Pyro [62] and uses logEI [46] 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.

8.6.0.4 MPD

Local BO via maximizing probability of decent (MPD) [15] is a follow-up method to the discussed GIBO approach. Unlike GIBO, it defines a different acquisition function for the sub-sampling step that aims to maximize the probability that the posterior mean points in a descent direction, rather than minimizing the uncertainty about the gradient. Additionally, it reuses the estimated posterior gradient GP to iteratively move the current point along the most probable descent direction 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 to the second stage of the algorithm, where the current GP estimate may be trusted for too many iterations before entering a sub-sampling step.

8.6.0.5 BayeSQP

For BayeSQP, we set the hyperparameters \(\delta_f, \delta_c= 0.2\), \(K=d+1\), \(M=3\), and \(\varepsilon=0.05\) in all experiments, unless stated otherwise. To solve \(\boldsymbol{B-SUB}\) at each iteration, we use CVXOPT [63] with standard parameters for maximum iterations and tolerance.⁶ Until we have not reached a feasible point, we set \(\delta_f=0.5\) to focus on robust improvement of the constraints and switch to the specified value once we have observed a 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 a one-dimensional line segment. On this line segment, we \(100\) sample points from a Sobol sequence as candidate points and choose the next sample location following [30]. The implementation is provided under https://github.com/brunzema/bayesqp and easily accessible via PyPI, i.e., using pip install bayesqp.

8.7 Model initialization and training↩︎

For all algorithms, we use a squared-exponential kernel. This kernel is sufficiently smooth such that we can formulate the joint GP for BayeSQP as specified in 3 . For the within-model comparisons, we freeze the lengthscales of the kernel and do not perform any hyperparameter optimization. For the experiments on classic benchmarks, we follow recent best-practice, wrap the squared-exponential kernel into a scale kernel and initialize all lengthscales with \(\sqrt{d}\). We furthermore as corse bounds on the lengthscales as \(\ell_i \in [0.001, 2 d]\) for all baselines and set the noise to a small values i.e., , \(\sigma^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 to improve numerical stability of the GP and normalize the inputs to the unit hypercube \([0, 1]^d\) [45].

9 Extensions to different settings and ideas for future work↩︎

BayeSQP provides a flexible framework that can be extended to various settings. In this section, we briefly discuss some of these extensions as well as other interesting avenues for futuree work.

9.0.0.1 Termination and restarting

The question of when to stop optimizing in BO has gained increasing attention in recent years [64]–[66]. 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 to develop an appropriate termination criterion such as stopping optimization once \(\|{\boldsymbol{p}}_t \|_2 \leq \tau_{\text{tol}}\), i.e., when the search direction becomes sufficiently small, indicating likely convergence to a 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.

9.0.0.2 Batch optimization

Similar to TuRBO, implementing BayeSQP for batch BO is straightforward: We can utilize different initial conditions as 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.

9.0.0.3 Localized GP model

To combat model mismatch on real-world problem, it is possible to include a sliding window on the training data of size \(N_{\text{max}}\) as also proposed in [10], [15]. This effectively produces a 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.

9.0.0.4 Active sub-sampling

In its current version, BayeSQP relies on a sub-sampling step to get good posterior estimates for the gradients and Hessians. While a space-filling sampling using a Sobol sequence already yields good results (cf. Section 5), an active approach to the sub-sampling—potentially with a stopping criterion—is interesting. One approach would be to build on ideas from [67]. However, an active sub-sampling likely will result in a much slower overall runtime so it depends on the specific application if this is desirable.

10 Proof of Corollary 1↩︎

The proof of Corollary 1 follows directly from the fact that \(q_{1-\delta}(0.5) = \Phi^{-1}(0.5)=0\). With this, \(b_f\) and \(b_{c_i}\) are zero and no longer influence the constraints or objective. Since \(b_f\) and \(b_{c_i}\) are optimization variables in \(\boldsymbol{B-SUB}\), the cones can be trivially satisfied.

11 Additional results↩︎

11.0.0.1 Out-of-model comparisons

In addition to 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 a sample from a 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 lengthscales of the objective. The results are summarized in Figure 7. We can observe the same trends as in the within-model comparisons thought gap in the final performance of TuRBO and BayeSQP is smaller. Looking at the optimization time of the out-of-model comparisons 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 7 still shows the clear trend of the local approaches BayeSQP and TuRBO outperforming SAASBO in high dimensions.

11.0.0.2 Ablation on the number of sub-samples \(K\) and line search samples \(M\)

Figure 5 (c) showed that for the specified \(K\) and \(M\), i.e., number of sub-samples and number of line search samples, including some uncertainty in \(\boldsymbol{B-SUB}\) will result in better performance. In Figure 9 we provide further ablations.

Figure 9 (a) presents the results of the sensitivity analysis over the number of sub-sampling steps (\(K\)) and line-search steps (\(M\)) for both the constrained within-model and out-of-model settings. In the within-model case, increasing \(K\) beyond \(d\) tends to degrade performance, whereas decreasing \(K\)—which allows for more \(\boldsymbol{B-SUB}\) solves under the same computational budget—can lead to improvements. For the 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 a suitable stopping criterion for the line search, enabling online adaptation of \(M\). Finally, note that the performance of 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 9 (c) shows that, in the unconstrained within-model case, reducing \(K\), as also observed in Figure 9 (a), can be beneficial. Moreover, even without accounting for uncertainty, the unconstrained case can achieve very good performance provided that the number of \(\boldsymbol{B-SUB}\) solves is sufficiently large. Notably, for \(K=64\), incorporating uncertainty leads to improved results, indicating that this configuration can more effectively handle scenarios with a limited number of line-searches. The out-of-model case in Figure 10 however highlights that in the absence of good prior knowledge, reducing the number of can be costly for small \(\delta_f\). A small \(\delta_f\) and high uncertainty in the estimates in \(\boldsymbol{B-SUB}\) will lead to small \({\boldsymbol{p}}_k\) and with this only limited progress towards the local optimum. A higher \(K\) results in 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 5 (or equivalently setting \(\delta_f=0.5\) for \(\boldsymbol{B-SUB}\)) can be sufficient.

12 Numerical considerations↩︎

12.0.0.1 Ensuring positive definiteness

The Hessian of the Lagrangian as described in 4 may become indefinite due to numerical issues, modeling inaccuracies, or nonconvexity in the surrogate models. To maintain numerical stability and ensure that curvature information defines a valid descent direction, we enforce positive definiteness through a simple eigenvalue modification. Concretely, let \(H\in\mathbb{R}^{n\times n}\) denote the Hessian candidate. We form the spectral decomposition \(H=Q\Lambda Q^\top\), where \(Q\) is orthogonal (\(Q^\top Q=I\)) and \(\Lambda=\operatorname{diag}(\lambda_1,\dots,\lambda_n)\) contains the real eigenvalues ordered arbitrarily. The eigenvalue-clipping rule replaces each eigenvalue \(\lambda_i\) by \(\tilde{\lambda}_i=\max(\lambda_i,\varepsilon)\) for a small threshold \(\varepsilon>0\) which we set as \(\varepsilon=10^{-5}\). The modified matrix is then reconstructed as \(\tilde{H}=Q\tilde{\Lambda} Q^\top\), where \(\tilde{\Lambda}=\operatorname{diag}(\tilde{\lambda}_1,\dots,\tilde{\lambda}_n)\). By construction \(\tilde{\lambda}_i\ge\varepsilon>0\) for all \(i\), hence \(\tilde{H}\) is symmetric positive definite. Note that more sophisticated modifications are also possible, but we found that this simple approach already resulted in satisfactory performance.

12.0.0.2 Jitter on the joint covariance

The joint covariance of the standard and derivative GP can be ill-conditioned. Here, we apply a standard jitter to the diagonal if necessary, which is also standard for covariances in classic BO. This then assures that the Cholesky decomposition for \(\boldsymbol{B-SUB}\) exists.

12.0.0.3 Ensure feasibility through slacked \(\boldsymbol{B-SUB}\) formulation

After every sub-sampling step, we aim to solve the \(\boldsymbol{B-SUB}\) optimization problem. However, given the problem and the current linearization, the sub-problem composed of the surrogate model estimates may be infeasible. We therefore opt to solve the following slack-constrained version of \(\boldsymbol{B-SUB}\): \[\begin{align} {\boldsymbol{p}}_t \in \mathop{\mathrm{arg\,min}}_{{\boldsymbol{p}}, b_f, \{b_{c_i}\}_{i \in {\mathbb{I}}_m}, \{s_i\}_{i \in {\mathbb{I}}_m}} \quad & \frac{1}{2}{\boldsymbol{p}}^{\top}{\boldsymbol{H}}_t{\boldsymbol{p}}+ {\boldsymbol{\mu}}_{\nabla f}^{\top}{\boldsymbol{p}}+ \mu_f + q_{1-\delta_f} b_f + \rho\sum_{i \in {\mathbb{I}}_m} s_i \\ \text{subject to} \quad & \left\|{\boldsymbol{L}}_f^{\top}\begin{bmatrix} 1 \\ {\boldsymbol{p}}\end{bmatrix}\right\|_2 \leq b_f, \quad \left\|{\boldsymbol{L}}_{c_i}^{\top}\begin{bmatrix} 1 \\ {\boldsymbol{p}}\end{bmatrix}\right\|_2 \leq b_{c_i}, \quad \forall i \in {\mathbb{I}}_m, \notag \\ & -{\boldsymbol{\mu}}_{\nabla c_i}^{\top}{\boldsymbol{p}}+ q_{1-\delta_c}b_{c_i} - s_i \leq \mu_{c_i}, \quad \forall i \in {\mathbb{I}}_m, \notag \\ & b_f \geq 0, \quad b_{c_i} \geq 0, \quad s_i \geq 0, \quad \forall i \in {\mathbb{I}}_m. \notag \end{align}\] where \({\boldsymbol{L}}_f\) and \({\boldsymbol{L}}_{c_i}\) are Cholesky factorizations as \[\begin{align} {\boldsymbol{L}}_f{\boldsymbol{L}}_f^{\top} & = \begin{bmatrix} \sigma_f^2 & \bullet \\ {\boldsymbol{\Sigma}}_{\nabla f, f} & {\boldsymbol{\Sigma}}_{\nabla f} \end{bmatrix}, \quad {\boldsymbol{L}}_{c_i}{\boldsymbol{L}}_{c_i}^{\top} = \begin{bmatrix} \sigma_{c_i}^2 & \bullet \\ {\boldsymbol{\Sigma}}_{\nabla c_i, c_i} & {\boldsymbol{\Sigma}}_{\nabla c_i} \end{bmatrix}, \quad \forall i \in {\mathbb{I}}_m, \end{align}\] and \(\rho > 0\) is the penalty parameter for slack variables. Here, we choose \(\rho=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 of the paper. It should further be noted that the feasibility of the problem will also depend on the \(\delta_f\) and \(\delta_c\). We therefore recommend that if the subproblem frequently fails to increase \(\delta_f\) and \(\delta_c\) and potentially use a form of scheduling for these hyperparameters.

13 Runtime breakdown of BayeSQP↩︎

To give an idea of the computational efficiency of 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 5 (b) and Figure 8. The reported results are from the within-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.

To provide deeper insight into the computational characteristics of each core component of BayeSQP, we 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 on constraints.

References↩︎