Paper deep dive
GRAPE: Gradient Refinement and Progress-Aware Exploitation for Query-Efficient High-Dimensional Bayesian Optimization
Richard Cornelius Suwandi, Feng Yin
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 96%
Last extracted: 8/27/2026, 4:43:21 AM
Summary
The paper introduces GRAPE, a two-stage Bayesian optimization framework for high-dimensional black-box functions. It addresses the limitation of existing methods that prioritize descent probability over progress magnitude by first refining the local gradient posterior and then selecting update directions that maximize expected decrease conditional on descent. Theoretical analysis proves monotonic minimization of local uncertainty and convergence to steepest descent. Empirical results show a 5.4x speedup in adversarial attacks and superior performance in LLM prompt optimization.
Entities (10)
Relation Signals (9)
Feng Yin → affiliatedwith → The Chinese University of Hong Kong, Shenzhen
confidence 99% · Feng Yin ... Affiliation: The Chinese University of Hong Kong, Shenzhen
Richard Cornelius Suwandi → affiliatedwith → The Chinese University of Hong Kong, Shenzhen
confidence 99% · Richard Cornelius Suwandi Affiliation: The Chinese University of Hong Kong, Shenzhen
Richard Cornelius Suwandi → authored → GRAPE
confidence 99% · Richard Cornelius Suwandi ... We introduce Gradient Refinement and Progress-Aware Exploitation (GRAPE)
Feng Yin → authored → GRAPE
confidence 99% · Feng Yin ... We introduce Gradient Refinement and Progress-Aware Exploitation (GRAPE)
GRAPE → uses → Gaussian Process
confidence 98% · Theoretical analysis proves that this gradient refinement stage monotonically minimizes local uncertainty... leveraging a key property of Gaussian processes (GPs)
GRAPE → improves → Query Efficiency
confidence 95% · GRAPE demonstrates superior query efficiency across high-dimensional tasks
GRAPE → isappliedto → LLM prompt optimization
confidence 95% · on large language model prompt optimization tasks, it outperforms the second best method
GRAPE → isappliedto →
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Optimizing expensive, high-dimensional black-box functions remains a central challenge in modern machine learning and scientific discovery. While local Bayesian optimization mitigates the curse of dimensionality, existing techniques often prioritize the probability of descent over the magnitude of progress. This leads to overly conservative steps that yield negligible improvement, wasting queries on directions that are nearly certain to descend but offer little decrease. We introduce Gradient Refinement and Progress-Aware Exploitation (GRAPE), a two-stage framework that first sharpens the local gradient posterior via a closed-form acquisition function, then selects update directions by maximizing the expected decrease conditional on descent. Theoretical analysis proves that this gradient refinement stage monotonically minimizes local uncertainty and that the progress-aware direction converges to true steepest descent as the posterior sharpens. Empirically, GRAPE demonstrates superior query efficiency across high-dimensional tasks: in black-box adversarial attacks, it achieves an average 5.4$\times$ speedup over baselines, and on large language model prompt optimization tasks, it outperforms the second best method by a reduction of 3.8 log-units in the final average regret.
Tags
Links
- Source: https://arxiv.org/abs/2608.25116v1
- Canonical: https://arxiv.org/abs/2608.25116v1
Trouble viewing inline? Open PDF directly →
Full Text
115,388 characters extracted from source content.
Expand or collapse full text
GRAPE: Gradient Refinement and Progress-Aware Exploitation for Query-Efficient High-Dimensional Bayesian Optimization Richard Cornelius Suwandi Affiliation: School of Artificial Intelligence Affiliation: The Chinese University of Hong Kong, Shenzhen Affiliation: Shenzhen, China Email: richardsuwandi@link.cuhk.edu.cn Feng Yin †thanks: Corresponding author. Affiliation: School of Artificial Intelligence Affiliation: The Chinese University of Hong Kong, Shenzhen Affiliation: Shenzhen, China Email: yinfeng@cuhk.edu.cn Abstract Optimizing expensive, high-dimensional black-box functions remains a central challenge in modern machine learning and scientific discovery. While local Bayesian optimization mitigates the curse of dimensionality, existing techniques often prioritize the probability of descent over the magnitude of progress. This leads to overly conservative steps that yield negligible improvement, wasting queries on directions that are nearly certain to descend but offer little decrease. We introduce Gradient Refinement and Progress-Aware Exploitation (GRAPE), a two-stage framework that first sharpens the local gradient posterior via a closed-form acquisition function, then selects update directions by maximizing the expected decrease conditional on descent. Theoretical analysis proves that this gradient refinement stage monotonically minimizes local uncertainty and that the progress-aware direction converges to true steepest descent as the posterior sharpens. Empirically, GRAPE demonstrates superior query efficiency across high-dimensional tasks: in black-box adversarial attacks, it achieves an average 5.4×5.4× speedup over baselines, and on large language model prompt optimization tasks, it outperforms the second best method by a reduction of 3.83.8 log-units in the final average regret. †footnotetext: Preprint. Under review. Keywords Bayesian optimization ⋅· Gaussian processes ⋅· high-dimensional black-box optimization ⋅· gradient-aware optimization ⋅· sample efficiency 1 Introduction Bayesian optimization (BO) is a sample-efficient framework for optimizing expensive black-box objectives whose gradients are unavailable and whose evaluations may be noisy (Shahriari et al.,, 2016; Frazier,, 2018; Garnett,, 2023). It maintains a probabilistic surrogate of the objective, typically using a Gaussian process (GP), and selects successive queries by maximizing an acquisition function that balances exploration and exploitation (Jones et al.,, 1998; Močkus,, 1978). This design makes BO particularly attractive whenever each evaluation is costly, as fewer queries translate directly into reduced experimental or computational expense. Consequently, BO has been applied across a broad range of domains, including hyperparameter tuning (Snoek et al.,, 2012; Arango et al.,, 2021), neural architecture search (White et al.,, 2021), safe controller configuration (Berkenkamp et al.,, 2016), chemical reaction and materials design (Shields et al.,, 2021; Wang and Dowling,, 2022), black-box adversarial attacks (Ru et al.,, 2020; Shu et al.,, 2023), and more recently, prompt optimization for large language models (Zhou et al.,, 2023; Chew et al.,, 2026). In each of these settings, the underlying objective is expensive to evaluate, whether because it requires training a neural network, running a physical experiment, or querying a large language model, making query efficiency the primary performance metric. Despite this sample efficiency, standard BO often degrades as the input dimension grows. A limited evaluation budget leaves most of the domain unexplored, a surrogate becomes uncertain over wide regions, and optimizing a global acquisition function becomes increasingly difficult (Doumont et al.,, 2026). The fundamental issue is that the volume of the search space grows exponentially with dimension, so a fixed evaluation budget covers a vanishingly small fraction of the domain. To solve these challenges, recent works have shown that an appropriate choice of GP kernel function and its hyperparameters can make standard BO surprisingly competitive in some high-dimensional problems (Hvarfner et al.,, 2024; Xu et al.,, 2025; Suwandi et al.,, 2025), although this behavior is sensitive to problem structure and GP hyperparameter initialization (Papenmeier et al.,, 2025). Existing high-dimensional BO methods commonly address this problem by assuming low-dimensional structures or restricting search to local regions. Yet, structural assumptions often fail when the objective is poorly understood, whereas local strategies trade global guarantees for sample-efficient improvement near a promising evaluation point. Gradient-aware local BO offers a particularly direct local strategy by leveraging a key property of Gaussian processes (GPs): a belief over the objective naturally induces a joint Gaussian belief over its gradient (Rasmussen and Williams,, 2006). As a result, even in the absence of explicit gradient observations, function-value evaluations can be used to refine the posterior over gradients. This is especially valuable in high dimensions, where the cost of exploring the full input space is prohibitive but local gradient information can guide efficient descent. For example, GIBO (Müller et al.,, 2021) first reduces gradient uncertainty before proceeding along the posterior mean gradient, while MPD (Nguyen et al.,, 2022) selects the direction with the highest posterior probability of descent. In contrast, MinUCB (Fan et al.,, 2024) exploits a local upper confidence bound. These approaches differ in how they leverage posterior information to guide local steps, but all share the insight that refining the gradient belief before moving can yield more productive queries than global exploration. In this work, we specifically revisit the notion of descent probability. While descent probability quantifies the certainty that a given direction descends, it does not capture the magnitude of improvement when descent occurs. A direction may be almost certain to descend while yielding only negligible progress, whereas another, less certain direction could result in substantially greater improvement when successful. Simply maximizing descent probability therefore favors safe, but potentially unproductive directions and can overlook opportunities for more significant progress. We address this limitation with Gradient Refinement and Progress-Aware Exploitation (GRAPE). Our main contributions are: • We propose a two-stage local BO framework that first sharpens the local gradient posterior via a closed-form acquisition, then selects update directions by maximizing the expected decrease conditional on descent. • We prove that gradient refinement monotonically minimizes local uncertainty and that the progress-aware direction converges to true steepest descent as the posterior sharpens. • Empirically, GRAPE achieves an average 5.4×5.4× speedup on black-box adversarial attacks and outperforms the second-best method by 3.83.8 log-units on LLM prompt optimization. 2 Related Work 2.1 High-Dimensional Bayesian Optimization High-dimensional BO methods commonly reduce the effective search complexity through structural assumptions or local modeling. Random-embedding methods optimize in a low-dimensional subspace under the assumption that the objective has low intrinsic dimension (Wang et al.,, 2016). Sparse axis-aligned models instead learn a small set of influential coordinates (Eriksson and Jankowiak,, 2021). These approaches can be effective when their structural assumptions match the objective, but their performance can deteriorate when the active structure is unknown or misspecified. Local modeling offers an alternative that does not require an explicit low-dimensional structure. TuRBO (Eriksson et al.,, 2019) fits GPs within adaptive trust regions that expand after successful steps and contract after failures, concentrating evaluations near promising points while using restarts to escape from poor basins. This local emphasis is shared by the gradient-aware methods discussed next, which replace trust-region success or failure rules with guidance from the GP derivative posterior. Recent work has also revisited standard, or “vanilla” BO in high dimensions. Hvarfner et al., (2024) scale a log-normal prior on GP length scales with dimension and obtain a strong global D-LogEI baseline. Xu et al., (2025) show that poor length-scale initialization can produce vanishing training gradients for squared-exponential kernels, while Matérn kernels or dimension-aware initialization can avoid this failure in some cases. Most recently, Suwandi et al., (2025) move beyond fixed kernel choices with CAKE, which uses LLMs to evolve GP kernel structures from observed data and ranks the resulting candidates by both model fit and acquisition utility. 2.2 Gradient-Aware Local Bayesian Optimization More closely related to our work are gradient-aware local methods, which alternate between learning about the gradient and using that information to update the current point. GIBO (Müller et al.,, 2021) minimizes posterior gradient uncertainty before moving along the negative posterior mean gradient. MPD (Nguyen et al.,, 2022) observes that the posterior mean gradient need not maximize descent probability and derives the most probable descent direction. MinUCB (Fan et al.,, 2024) replaces the gradient step with local minimization of an upper confidence bound. NeST-BO (Tang et al.,, 2026) jointly predicts the gradient and Hessian from the GP and selects evaluations that reduce uncertainty in a modified Newton step. This curvature information can improve local scaling and convergence in ill-conditioned regions, but learning a d×d× d Hessian introduces quadratic derivative complexity. The convergence behavior of gradient-aware local methods has also been studied under smooth-kernel assumptions (Wu et al.,, 2023). GRAPE modifies the exploitation criterion rather than increasing the derivative order or optimizing a confidence bound. Its progress-aware criterion ranks directions by the expected decrease conditional on descent, while gradient refinement greedily reduces the uncertainty of the gradient used by that criterion. The approach requires only function evaluations and first-order GP derivatives, making it applicable in settings where second-order information is too costly to recover or where the objective is not sufficiently smooth for Hessian-based methods to be reliable. 2.3 Global–Local Hybrids and Probabilistic Numerics Local optimization can also be embedded within a global optimization strategy. For example, BLOSSOM switches between BO acquisitions and conventional local optimization (McLeod et al.,, 2018), while TREGO alternates global BO steps with trust-region steps (Diouane et al.,, 2023). A gradient-aware optimizer such as GRAPE can serve as the local component of a similar hybrid when global coverage is required. More broadly, our approach adopts the probabilistic numerics view that numerical computation can be formulated as inference (Hennig et al.,, 2022). The unavailable gradient is treated as a latent quantity inferred from function evaluations, and posterior uncertainty directly determines which information to collect and how to move. This perspective unifies the exploration and exploitation stages of GRAPE under a single Bayesian framework, where both the choice of refinement queries and the selection of descent directions are governed by the same posterior distribution. 3 Preliminaries 3.1 Problem Setting Let f:→ℝf:X be a continuous black-box function on a compact domain ⊂ℝdX ^d. Given a starting point 0x_0, our goal is to minimize f in a local neighborhood (0)⊆X(x_0) of 0x_0: ∗=argmin∈(0)f().x^*= _x (x_0)f(x). (1) We can query the objective at locations of our choice and observe y=f()+ε,ε∼(0,σ2),y=f(x)+ , (0,σ^2), (2) where σ2≥0σ^2≥ 0 is the observation-noise variance. Here, the gradient of f is unavailable, and each function evaluation is assumed to be costly. Consequently, an effective black-box optimizer must decide both where to query for information and how to use that information to move from the current point. Local optimization makes this problem more tractable in high dimensions by replacing broad coverage of X with targeted learning near the current iterate. 3.2 Bayesian Optimization with Gaussian Processes Bayesian optimization (BO) tackles expensive black-box optimization by maintaining a probabilistic model of the unknown objective and using that model to select informative evaluation points (Garnett,, 2023). Specifically, we place a Gaussian process (GP) prior f∼(μ,k)f (μ,k) on f, where μ is the mean function and k is the kernel (covariance) function. For any finite set of inputs, the corresponding function values are jointly Gaussian Rasmussen and Williams, (2006), allowing us to compute an analytic posterior after observing the data =(,)D=(X,y). Define =k(,)+σ2K=k(X,X)+σ^2I. Here, k(,)k(x,X) is the row vector of covariances between x and the training inputs X, and k(,′)k(X,x ) is the corresponding column vector. The posterior mean and covariance are given by μ() _D(x) =μ()+k(,)−1(−μ()), =μ(x)+k(x,X)K^-1(y-μ(X)), (3) k(,′) k_D(x,x ) =k(,′)−k(,)−1k(,′). =k(x,x )-k(x,X)K^-1k(X,x ). (4) The posterior mean μ _D provides a prediction of the objective at x, while the posterior variance k(,)k_D(x,x) quantifies the uncertainty. BO uses these predictive quantities to construct an acquisition function, which balances exploration and exploitation by assigning high utility to points that are either promising (low predicted objective) or uncertain (high variance). As an example, expected improvement (Jones et al.,, 1998; Močkus,, 1978) favors points that are likely to yield an improvement over the current best. After evaluating the acquisition-maximizing point, BO updates the GP model with the new observation and repeats the process. 3.3 Gradient-Aware Local Optimization GPs are particularly suitable for gradient-aware local optimization because linear operations preserve Gaussianity. If μ is differentiable and k is twice differentiable, function values and partial derivatives are jointly Gaussian (Rasmussen and Williams,, 2006). Thus, although gradients are never queried directly, differentiating the posterior GP with respect to x yields ∇f()|∼(,),∇ f(x) ( μ_x, _x), (5) where μ_x =∇μ()+∇k(,)−1(−μ()), =∇μ(x)+∇ k(x,X)K^-1(y-μ(X)), (6) _x =∇k(,)∇⊤−∇k(,)−1k(,)∇⊤. =∇ k(x,x)∇ -∇ k(x,X)K^-1k(X,x)∇ . (7) The differential operator before k acts on its first argument and the operator after k acts on its second argument. The vector μ_x is the posterior mean of the local gradient, and _x quantifies the uncertainty in that gradient. The posterior in Eq. (5) induces a general two-stage local BO strategy at the current point x: Stage 1 selects function evaluations that refine the local derivative belief, and Stage 2 uses that belief to choose a descent step. Several methods follow this general framework. GIBO minimizes total gradient uncertainty and takes a step along the posterior mean of the gradient (Müller et al.,, 2021). MinUCB retains a similar exploration stage, but instead of stepping along the mean gradient, it exploits by minimizing a local upper confidence bound (Fan et al.,, 2024). NeST-BO extends this idea by jointly learning gradient and Hessian posteriors to target a modified Newton step (Tang et al.,, 2026). Most closely related to our approach, MPD leverages the full first-order posterior to identify the direction most likely to decrease the objective (Nguyen et al.,, 2022). Because our exploitation criterion is directly based on this construction, we review it next. To formalize MPD’s criterion, let ∈ℝdv ^d with ‖=1\|v\|=1 denote a candidate search direction on the unit sphere at x. The corresponding directional derivative ∇f()=⊤∇f() _vf(x)\;=\;v ∇ f(x) (8) is the rate of change of f at x when moving along v. Because ∇f()|∼(,)∇ f(x) ( μ_x, _x), the projected scalar is univariate Gaussian, ∇f()| _vf(x) ∼(μ,σ2), ( _v, _v^2), (9) μ _v =⊤,σ2=⊤. =v μ_x, _v^2=v _xv. A direction is locally descending when this directional derivative is negative. Its posterior probability of descent is therefore Pr(∇f()<0∣)=Φ(−μσ), ( _vf(x)<0 )= \! (- _v _v ), (10) where Φ denotes the standard normal CDF. Unlike a step based only on the posterior mean, this probability accounts for both the estimated slope and its uncertainty. MPD maximizes this probability over directions. When _x is positive definite, its unique maximizing direction up to scale is MPD∗=−1‖−1‖,v^*_ MPD=- _x^-1 μ_x\| _x^-1 μ_x\|, (11) and the corresponding maximum descent probability is Φ(⊤−1). \! ( μ_x _x^-1 μ_x ). (12) We retain MPD’s directional posterior but argue that maximizing descent probability alone is not enough for sample-efficient local progress. The next section makes this limitation precise and develops the progress-aware criterion of GRAPE. 4 Proposed Method 4.1 Beyond Descent Probability The descent probability in Eq. (10) depends on the directional posterior only through the standardized mean γ=μ/σ _v= _v/ _v: ℙ(∇f()<0)=Φ(−γ).P ( _vf(x)<0 )= (- _v). (13) Consequently, Φ(−γ) (- _v) is invariant to any positive rescaling (μ,σ)↦(cμ,cσ)( _v, _v) (c _v,c _v). Two unit directions with the same γ _v therefore receive identical scores, even when one expected slope is far steeper than the other. Descent probability thus ranks how certain a direction is to descend, not how much decrease it offers when descent occurs. To see why this matters in practice, consider two candidate directions at a point x with directional posteriors (μ1,σ1)=(−2,1)( _1, _1)=(-2,1) and (μ2,σ2)=(−0.2,0.1)( _2, _2)=(-0.2,0.1). Both have γ=−2γ=-2 and therefore share the same descent probability Φ(2)≈0.977 (2)≈ 0.977. Yet their expected decreases conditional on descent are ≈2.05≈ 2.05 and ≈0.21≈ 0.21, respectively, an order-of-magnitude difference that descent probability is blind to. In high dimensions, where the posterior is often diffuse early in optimization, such ties are common and the resulting direction choices can waste queries on nearly flat but certainly-descending directions. The next subsection restores this missing magnitude. 4.2 Progress-Aware Exploitation We rank directions by the expected decrease conditional on descent, and use this score as the exploitation stage of GRAPE. Definition 1 (Progress-aware exploitation score). Let γ=μ/σ _v= _v/ _v denote the standardized directional mean. The progress-aware exploitation score of direction v at x is () (v) :=[−∇f()∣∇f()<0] :=E [- _vf(x) _vf(x)<0 ] =σ[ϕ(γ)Φ(−γ)−γ], = _v\! [ φ( _v) (- _v)- _v ], (14) where ϕφ and Φ are the standard normal PDF and CDF. The closed form follows directly from the mean of a truncated Gaussian. For Z∼(μ,σ2)Z (μ,σ^2), the mean of the lower-truncated distribution at 00 is [Z∣Z<0]=μ−σϕ(μ/σ)/Φ(−μ/σ)E[Z Z<0]=μ-σφ(μ/σ)/ (-μ/σ). Applying this identity to Z=∇f()Z= _vf(x) and negating the resulting conditional mean yields [−∇f()∣∇f()<0]=−μ+σϕ(γ)Φ(−γ),E[- _vf(x) _vf(x)<0]=- _v+ _v φ( _v) (- _v), (15) which is equal to Eq. (14). We then select the unit direction that maximizes this expected conditional decrease: ∗=argmax‖=1().v^*= _\|v\|=1P(v). (16) We solve Eq. (16) via projected gradient ascent (Calamai and Moré,, 1987). The analytic gradient of ()P(v) and its relation to MPD, PI, and EI are derived in Appendices K and J. The progress-aware score has several desirable properties. First, it naturally balances descent probability and expected magnitude: when γ→−∞ _v→-∞ (highly certain descent), ()→−μP(v)→- _v, recovering the posterior mean slope, and when γ→+∞ _v→+∞ (certain ascent), ()→0P(v)→ 0, correctly assigning no value to ascending directions. In the intermediate regime where γ≈0 _v≈ 0, the score is dominated by the uncertainty term σϕ(0)/Φ(0) _vφ(0)/ (0), favoring directions with high posterior variance, a form of implicit exploration that is valuable when the gradient belief is still diffuse. Second, unlike unconditional expected descent [−∇f()]=−μE[- _vf(x)]=- _v, which can be positive even for directions that are likely to ascend (when μ<0 _v<0 but the variance is large), the conditional expectation ensures that we only count progress in directions that actually descend. This makes ()P(v) a more reliable criterion for selecting productive update directions. 4.3 Gradient Refinement Progress-aware exploitation ranks directions under the gradient posterior. When that posterior is diffuse, the ranking can be unreliable: a large ()P(v) may reflect posterior noise rather than a genuinely steep descent direction. This is particularly problematic early in optimization, when few evaluations have been collected and the GP surrogate is uncertain over large regions of the input space. Before moving, GRAPE therefore allocates a short budget of function evaluations to sharpen _x at the current iterate x. We spend this budget on τexplore _explore auxiliary queries chosen to improve the local gradient estimate at x. For a candidate location ∈z , a query would return y=f()+εy_z=f(z)+ with ε∼(0,σ2) (0,σ^2). Let |∪(,y) _x|D∪(z,y_z) denote the gradient covariance at x after augmenting D with (,y)(z,y_z). We select z to maximize the gradient-refinement acquisition, which is the expected one-step reduction in total gradient variance, αref() _ref(z) =y[Tr()−Tr(|∪(,y))]. =E_y_z\! [ ( _x)- \! ( _x|D∪(z,y_z) ) ]. (17) Here Tr() ( _x) is the current total gradient uncertainty, and Tr(|∪(,y)) \! ( _x|D∪(z,y_z) ) is the residual uncertainty after observing y_z at z. The expectation is over the yet-unobserved outcome y_z because the acquisition must be evaluated before querying z. At first glance, Eq. (17) requires integrating over the unknown scalar y_z. For a GP, however, the posterior covariance update depends only on the query location, not on the realized function value. The expectation therefore collapses, and αref _ref admits a closed form. Proposition 1 (Closed-form gradient refinement). Let k(⋅,⋅)k_D(·,·) denote the GP posterior kernel after observing D. Then αref() _ref(z) is independent of the unobserved value y_z and equals αref()=‖∇k(,)‖2k(,)+σ2, _ref(z)= \| _xk_D(x,z)\|^2k_D(z,z)+σ^2, (18) where ∇k(,)=∂k(,)/∂ _xk_D(x,z)=∂ k_D(x,z)/ is the posterior cross-covariance gradient. The proof is provided in Appendix I. Proposition 1 has two practical consequences. First, gradient refinement can be optimized as a deterministic acquisition function, with no Monte Carlo sampling. Second, αref() _ref(z) is large when the posterior gradient covariance between x and z is strong and the posterior variance at z is not already negligible. Otherwise, querying z adds little information. 4.4 The GRAPE Algorithm We combine gradient refinement (Section 4.3) and progress-aware exploitation (Section 4.2) into GRAPE, a two-stage local optimization loop. Stage 1 selects auxiliary queries to sharpen the gradient posterior at the current iterate, and Stage 2 moves along directions of maximum expected progress under that refined posterior. We initialize by drawing ninitn_init random queries, forming a dataset D, and fitting a GP. The random initialization ensures that the initial surrogate has some coverage of the local neighborhood, preventing the gradient posterior from being dominated by the prior. Then, for each outer iteration t=0,…,T−1t=0,…,T-1 starting at tx_t, Stage 1 sequentially selects τexplore _explore queries by maximizing the closed-form acquisition ∗←argmax∈αref()z^*← _z _ref(z) in Eq. (18), observes each y=f(∗)+εy=f(z^*)+ , and refits the GP. The sequential selection allows each refinement query to benefit from the information gained by previous queries, leading to more efficient uncertainty reduction than batch selection. Stage 2 then repeatedly computes the unit direction of maximum progress ∗←argmax‖=1()v^*← _\|v\|=1P(v) in Eq. (16) and takes a projected step ←Π(+η∗)x← _X(x+ ^*) of size η>0η>0, where Π:ℝd→ _X:R^d is Euclidean projection onto the domain. Each new point is evaluated and added to D before the next direction is chosen, so the posterior is updated after every move. This online updating ensures that the progress-aware score reflects the most current gradient belief, which is particularly important when the objective has varying curvature across the local neighborhood. Exploitation stops after τexploit _exploit steps or when (∗)<τthreshP(v^*)< _thresh, indicating that the posterior predicts little conditional descent. The early-stopping criterion prevents wasting queries on directions that are unlikely to yield meaningful progress, which is especially valuable when the optimizer is approaching a local minimum. The outer iterate is then set to t+1←x_t+1 , and after T iterations GRAPE returns the best observed point argmin(,y)∈y _(x,y) \,y. Figure 1 illustrates one outer iteration of this loop. i∗z_i^*tx_tpriorposteriorThm. 1Stage 1(a) Gradient Refinementquery nearby to sharpen the descent posteriorrefinedbeliefcertain, shallowsteep progresstx_tt+1x_t+1progresshighlowStage 2(b) Progress-Aware Exploitationmove into the direction where expected progress is largest Figure 1: One outer iteration of GRAPE. (a) Nearby queries sharpen a wide prior descent belief into a narrower posterior (Theorem 1). (b) Candidate directions are ranked by their expected progress. 5 Theoretical Analysis This section analyzes the two design choices introduced above: gradient refinement and progress-aware exploitation. 5.1 Notation and Assumptions Let ⊂ℝdX ^d be a convex compact domain and ℋH a reproducing kernel Hilbert space (RKHS) on X with kernel k. At an iterate tx_t, let tD_t denote all observations collected so far, and let t μ_t and t _t denote the posterior mean and covariance of ∇f(t)∇ f(x_t) under tD_t. Our analysis aims to answer two questions. First, does a gradient-refinement observation reduce uncertainty for a fixed iterate? Second, when that uncertainty vanishes and the posterior mean is consistent, which direction does progress-aware exploitation select? To carry out the analysis, we impose the following standard regularity conditions. Assumption 1. The kernel k is stationary, four times continuously differentiable, and positive definite. Assumption 2. The iterates t\x_t\ remain in the interior of X. Assumption 3. The objective lies in the RKHS of k with bounded norm: f∈ℋf and ‖f‖ℋ≤B\|f\|_H≤ B, where B<∞B<∞ is a fixed bound on the RKHS norm of f. Assumptions 1–3 are standard (Bull,, 2011; Srinivas et al.,, 2012). Since k is four times differentiable, f is twice continuously differentiable. On the compact domain X this implies that ∇f∇ f is Lipschitz, i.e., f is L-smooth for some finite constant L>0L>0 (Wu et al.,, 2023). Compactness ensures that the direction optimization and smoothness constants are well-defined, while the interior-iterate assumption avoids boundary projections in the limiting direction result. The RKHS assumption is a regularity condition on the objective: it requires that f lies in the function space induced by the kernel, which is a common assumption in GP-based optimization analyses (Srinivas et al.,, 2012). While this assumption may not hold exactly in practice, it provides a useful framework for analyzing the convergence behavior of GP-based methods. In our experiments, we use the Matérn-5/25/2 kernel, whose RKHS consists of functions that are twice continuously differentiable, a reasonable assumption for many practical objectives, including neural network loss landscapes and prompt accuracy surfaces. 5.2 Gradient Uncertainty Reduction Building on Proposition 1, the result below formalizes the one-step contraction of Tr() ( _x) and identifies when it is strict. Theorem 1 (Gradient refinement reduces uncertainty). Under the joint GP model, for any candidate location z, y[Tr(|∪(,y))]≤Tr(),E_y_z\! [ ( _x|D∪(z,y_z)) ]≤ ( _x), with equality if and only if ∇k(,)= _xk_D(x,z)=0. Proof. From Proposition 1, αref()=‖∇k(,)‖2/(k(,)+σ2)≥0 _ref(z)=\| _xk_D(x,z)\|^2/(k_D(z,z)+σ^2)≥ 0. Hence [Tr(|∪(,y))]=Tr()−αref()≤Tr()E[ ( _x|D∪(z,y_z))]= ( _x)- _ref(z)≤ ( _x), with equality iff ∇k(,)= _xk_D(x,z)=0. □ Since αref() _ref(z) equals the expected one-step trace reduction, selecting ∗=argmaxαref()z^*= _z _ref(z) is greedily optimal and repeated refinement monotonically decreases Tr() ( _x). 5.3 Asymptotic Descent Direction The result below shows that, once the gradient posterior is sharp and consistent, the progress-aware direction converges to normalized steepest descent. Theorem 2 (Progress-aware exploitation approaches steepest descent). Suppose ∇f()≠∇ f(x) 0, Tr()→0 ( _x)→ 0, and →∇f() μ_x→∇ f(x). Then the progress-aware direction converges to the true steepest descent direction: limTr()→0∗=−∇f()‖∇f()‖. _ ( _x)→ 0\;v^*=- ∇ f(x)\|∇ f(x)\|. Proof sketch. As Tr()→0 ( _x)→ 0, → _x 0 in spectral norm while →∇f() μ_x→∇ f(x) by assumption. Let R(γ):=ϕ(γ)/Φ(−γ)R(γ):=φ(γ)/ (-γ) denote the inverse Mills ratio and [a]+:=max(a,0)[a]^+:= (a,0). For any unit direction v, if ⊤∇f()<0v ∇ f(x)<0, then γ→−∞ _v→-∞, R(γ)→0R( _v)→ 0, and ()→−⊤∇f()>0P(v)→-v ∇ f(x)>0. If ⊤∇f()>0v ∇ f(x)>0, then γ→+∞ _v→+∞ and ()→0P(v)→ 0 by the expansion R(γ)=γ+γ−1+(γ−3)R(γ)=γ+γ^-1+O(γ^-3) (Grimmett and Stirzaker,, 2020). Hence ()→[−⊤∇f()]+P(v)→[-v ∇ f(x)]^+ pointwise on d−1S^d-1. Uniform convergence on this compact set, together with the Berge Maximum Theorem (Berge,, 1963), yields ∗→−∇f()/∥∇f()∥v^*→-∇ f(x)/\|∇ f(x)\|. A complete proof appears in Appendix G. □ 5.4 Interpretation The two theorems above characterize the complementary roles of the two stages in GRAPE. Theorem 1 guarantees that each refinement query strictly reduces total gradient uncertainty unless the candidate point is already uninformative about the local gradient, a condition that holds whenever z lies in a region where the posterior cross-covariance gradient ∇k(,) _xk_D(x,z) vanishes. This monotonic contraction ensures that the directional ranking used by the exploitation stage becomes increasingly reliable as refinement proceeds. Theorem 2 then shows that, once the posterior has sharpened sufficiently and the posterior mean is consistent with the true gradient, the progress-aware direction recovers normalized steepest descent. Together, these results establish that GRAPE behaves as a principled approximation to first-order optimization: refinement drives the posterior toward a regime where the progress-aware score is well-calibrated, and exploitation then selects directions that are asymptotically equivalent to those chosen by a gradient oracle. We empirically validate Theorems 1 and 2 in Appendix F. 6 Experiments In this section, we evaluate GRAPE on two high-dimensional, gradient-free tasks and compare it against representative global and local BO methods. Code is available at https://github.com/richardcsuwandi/grape. 6.1 Benchmarks Black-Box Adversarial Attacks Adversarial attacks reveal how small, visually inconspicuous input changes can cause otherwise accurate classifiers to fail, making them an important tool for assessing the robustness of deployed models (Ru et al.,, 2020). In many realistic settings, however, the attacker can query only the model’s outputs and cannot access its parameters or gradients. Given such a black-box classifier and an image z, we therefore seek, through function queries only, a small perturbation x such that the perturbed image +z+x is misclassified. This is a natural fit for Bayesian optimization: each query corresponds to evaluating the classifier on a candidate perturbation, and the goal is to find a successful attack with as few queries as possible. Following the practice of Cheng et al., (2021), we randomly select images from MNIST (LeCun et al.,, 1998) (d=28×28=784d=28× 28=784) and CIFAR-10 (Krizhevsky and Hinton,, 2009) (d=32×32=1024d=32× 32=1024) and add a perturbation under an ℓ∞ _∞ constraint that makes the trained deep network misclassify the image. Concretely: • For MNIST, we use the same fully trained network as Cheng et al., (2021) and adopt the constraint ‖∞≤0.3\|x\|_∞≤ 0.3. • For CIFAR-10, we fully train a ResNet-18 (He et al.,, 2016) using SGD with a cosine-annealed learning rate from 0.10.1 to 00, momentum 0.90.9, and weight decay 5×10−45× 10^-4 for 200200 epochs, and adopt the constraint ‖∞≤0.2\|x\|_∞≤ 0.2. LLM Prompt Optimization The choice of prompt can substantially affect an LLM’s performance, yet designing effective prompts by hand is labor-intensive and does not scale across models and tasks (Zhou et al.,, 2023). Automatic prompt optimization is therefore valuable, but each objective evaluation requires LLM inference, while the search space is discrete in text space and high-dimensional in embedding space. This combination of expensive evaluations and high-dimensional discrete search makes prompt optimization a challenging testbed for black-box optimization methods. We use BoLT (Chew et al.,, 2026) because its precomputed prompt scores make evaluations inexpensive and reproducible. It provides 5,0145,014 prompts for mathematical reasoning, scored by their MATH-500 accuracy (Hendrycks et al.,, 2021) under Qwen3-14B (Yang et al.,, 2025). Each prompt is embedded with EmbeddingGemma (Schechter Vera et al.,, 2025) at the four Matryoshka truncations d∈128,256,512,768d∈\128,256,512,768\ (Kusupati et al.,, 2022), and we correspondingly denote the four tasks as PO-128 through PO-768. The search space is therefore the discrete candidate pool PO=ii=15014⊂ℝdX_PO=\z_i\_i=1^5014 ^d, where fPO(i)f_PO(z_i) is the corresponding accuracy and the reference optimum is the highest score in the pool. Because GRAPE and the local baselines propose continuous steps, each proposal is projected to its nearest embedding in the pool before evaluation. This projection introduces a quantization effect that is unique to discrete search spaces: even if the optimizer identifies a promising direction in embedding space, the actual evaluation point is constrained to the nearest available prompt, which may limit the achievable progress per step. 6.2 Baselines We compare GRAPE against a diverse set of baselines spanning global, local, first-order, and second-order approaches. These include random search as a sanity check, and two global BO methods: vanilla BO (VBO) with expected improvement (Jones et al.,, 1998) and D-LogEI with a dimension-scaled length-scale prior (Hvarfner et al.,, 2024). We further include the local first-order methods MPD (Nguyen et al.,, 2022), GIBO (Müller et al.,, 2021), and MinUCB (Fan et al.,, 2024), as well as the trust-region method TuRBO (Eriksson et al.,, 2019). Finally, we include NeST-BO (Tang et al.,, 2026), a recent second-order method that targets a modified Newton step using jointly learned gradient and Hessian posteriors. This selection covers the main design choices in high-dimensional BO: global versus local search, first-order versus second-order derivative information, and trust-region versus gradient-based step selection. 6.3 Experimental Setup All GP-based methods share the same modeling protocol: an ARD Matérn-5/25/2 kernel with length scales and output scale fit by maximum marginal likelihood, together with common initial points and evaluation budgets. We use the Matérn-5/25/2 kernel because it provides a good balance between smoothness and flexibility: it is twice differentiable, which is sufficient for the gradient posterior used by all first-order methods, while avoiding the oversmoothing behavior of the squared-exponential kernel in high dimensions (Xu et al.,, 2025). For GRAPE, we set τexplore=5 _explore=5, allow at most τexploit=30 _exploit=30 exploitation steps, use step size η=0.1⋅diam()η=0.1·diam(X), where diam()=sup,′∈‖−′‖diam(X)= _x,x \|x-x \| is the diameter of the feasible domain, and early-stop threshold τthresh=5×10−3 _thresh=5× 10^-3, and optimize each direction with 1010 random restarts and a warm start. On prompt optimization, we set observation noise variance to σ2=0.001σ^2=0.001 and run T=200T=200 iterations from five random prompts. Prompt-optimization results are averaged over 1010 random seeds, attack query counts in Table 1 are averaged over 1010 independent runs, while the success-rate curves in Figure 2 are evaluated on 5050 randomly selected images per dataset. Error bars and shaded regions denote one standard deviation. Implementations use PyTorch (Paszke et al.,, 2019) and GPyTorch (Gardner et al.,, 2018) on one NVIDIA RTX 4090. Full protocols, baseline hyperparameters, and compute details are in Appendix B. 6.4 Evaluation Metrics For adversarial attacks, we report the number of queries to the first successful attack and the attack success rate under a fixed query budget (Shu et al.,, 2023). Fewer queries and a higher success rate are better. The query count directly measures sample efficiency, while the success rate curve reveals how quickly each method accumulates successful attacks as the budget increases, a distinction that matters in practice, where the available query budget may be constrained by time or API costs. For prompt optimization, we report the log simple regret (Chew et al.,, 2026), ℓt=log(f(∗)−f(^t)) _t= (f(x^*)-f( x_t) ), where f(^t)f( x_t) is the best score after t evaluations and f(∗)f(x^*) is the largest score in the candidate pool. The logarithmic scale compresses the dynamic range of regret values, making it easier to distinguish methods in the low-regret regime where absolute differences are small but practically meaningful. We assess significance with a paired two-sided Wilcoxon signed-rank test (Wilcoxon,, 1945), pairing methods by run, where improvements with p<0.05p<0.05 are marked with ∗. 6.5 Results We compare GRAPE with the baselines on both benchmarks. MNIST (d=784d=784) CIFAR-10 (d=1024d=1024) Method # of queries Ratio # of queries Ratio Random 1971±498∗1971± 498^* 9.1×9.1× 3222±662∗3222± 662^* 9.4×9.4× VBO 1589±271∗1589± 271^* 7.4×7.4× 2563±436∗2563± 436^* 7.5×7.5× D-LogEI 1381±222∗1381± 222^* 6.4×6.4× 2114±372∗2114± 372^* 6.2×6.2× TuRBO 1143±192∗1143± 192^* 5.3×5.3× 1239±191∗1239± 191^* 3.6×3.6× GIBO 981±184∗981± 184^* 4.5×4.5× 1613±271∗1613± 271^* 4.7×4.7× MinUCB 891±169∗891± 169^* 4.1×4.1× 1449±213∗1449± 213^* 4.2×4.2× MPD 827±128∗827± 128^* 3.8×3.8× 1338±236∗1338± 236^* 3.9×3.9× NeST-BO 719±136∗719± 136^* 3.3×3.3× 1123±197∗1123± 197^* 3.3×3.3× GRAPE 216± 42 1.0× 341± 53 1.0× Table 1: Average number of queries (↓ ) ± standard deviation required for a successful attack, across 10 independent runs. Ratio is each method’s query count divided by that of GRAPE (higher is worse). ∗ indicates that GRAPE requires significantly fewer queries than the marked baseline under a paired Wilcoxon signed-rank test (p<0.05p<0.05). Figure 2: Attack success rate (↑ ) versus number of queries on 5050 randomly selected images from MNIST (d=784d=784) and CIFAR-10 (d=1024d=1024). Figure 3: Log simple regret (↓ ) versus BO iteration on BoLT prompt-optimization tasks PO-128–PO-768. Shaded regions denote one standard deviation over 10 seeds. Black-Box Adversarial Attacks Table 1 reports the number of queries needed for a successful attack. In these high-dimensional settings, GRAPE substantially outperforms every baseline, reaching misclassification with far fewer queries. All pairwise comparisons are statistically significant (p<0.05p<0.05). The largest gaps appear on CIFAR-10 against global BO and random search: VBO and Random require 7.5×7.5× and 9.4×9.4× as many queries as GRAPE, reflecting the difficulty of global search in perturbation spaces with d≈103d≈ 10^3. Among local methods, NeST-BO is the closest competitor but still needs 3.3×3.3× more queries on both datasets. Other first-order local methods (MPD, MinUCB, GIBO) trail further (3.83.8–4.7×4.7×), as expected when optimization prioritizes descent certainty over the magnitude of conditional decrease. GRAPE also outperforms TuRBO, despite its strong record on query-efficient black-box attacks (Ru et al.,, 2020). On these tasks, refining the local gradient before each perturbation update stabilizes direction choice in a high-dimensional pixel space, and progress-aware steps target large expected loss reductions rather than timid moves that barely approach the decision boundary. The advantage of GRAPE is particularly pronounced on CIFAR-10 relative to MNIST. The CIFAR-10 decision boundaries are more complex due to the richer feature structure of natural color images, which makes the local gradient posterior more uncertain for a given number of queries. Gradient refinement therefore provides greater relative benefit on CIFAR-10, as it concentrates queries where they most reduce directional uncertainty. Figure 2 confirms this trend under a fixed budget of 20002000 queries: GRAPE reaches the highest success rate earliest, followed by NeST-BO and the other local methods, while global baselines lag, especially on CIFAR-10. LLM Prompt Optimization Figure 3 shows the log simple regret against iteration on the prompt-optimization tasks, as the embedding dimension increases from 128128 to 768768. Overall, performance can be grouped into three regimes. In the first, global methods (VBO, D-LogEI) and random search plateau early, typically between −2-2 and −4-4 log simple regret, and perform worse as the dimension grows. This matches the observation from Chew et al., (2026), who report that standard BO struggles on high-dimensional prompt embeddings and becomes nearly indistinguishable from random search: with a fixed budget, the surrogate cannot usefully resolve a discrete candidate pool whose geometry grows harder as Matryoshka truncations retain more coordinates. In the second regime, first-order local methods (MPD, MinUCB, GIBO) and TuRBO improve on that baseline but stall at a mid-tier plateau, suggesting that descent probability or trust-region control alone cannot sustain improvement once easy gains are exhausted. In the third regime, only GRAPE and NeST-BO reach the deep low-regret region (below −10-10). NeST-BO is the strongest competitor, yet GRAPE descends sooner and finishes with the lowest regret at every embedding dimension. In this setting, gradient refinement sharpens the surrogate before each prompt edit, and progress-aware selection favors candidates with large expected score gains when they succeed, which helps identify substantive prompt changes rather than small edits that are merely likely to help in a sparse, high-dimensional embedding space. Notably, the gap between GRAPE and NeST-BO widens as the dimension increases: at d=128d=128, both methods reach similar final regret, but at d=768d=768, GRAPE achieves roughly 33 additional log-units of improvement. This trend is consistent with the different derivative complexities of the two approaches: NeST-BO must learn a d×d× d Hessian posterior from function evaluations alone, which becomes increasingly data-hungry as d grows, whereas GRAPE operates on the first-order gradient posterior whose uncertainty scales more favorably. The result suggests that, in very high-dimensional discrete search spaces, investing the evaluation budget in sharpening the gradient belief and exploiting it with a progress-aware criterion is more sample-efficient than attempting to recover second-order information. We report additional results, including ablations, hyperparameter sensitivity, and wall-clock timings, in Appendices A, C, and E. 7 Conclusion We introduced GRAPE, a two-stage local BO framework that couples gradient refinement with progress-aware exploitation to achieve strong query efficiency in high-dimensional black-box optimization. The gradient refinement stage uses a closed-form acquisition to greedily reduce posterior uncertainty at the current iterate, ensuring that the directional ranking used for exploitation is grounded in a sharp local gradient estimate. The progress-aware exploitation stage then selects update directions by maximizing the expected decrease conditional on descent, a criterion that restores the magnitude information discarded by pure descent-probability maximization. We proved that refinement monotonically reduces local gradient uncertainty and that the progress-aware direction converges to true steepest descent as the posterior sharpens. Empirically, GRAPE achieves an average 5.4×5.4× speedup over baselines on black-box adversarial attacks and outperforms the second-best method by 3.83.8 log-units on LLM prompt optimization, demonstrating that conditioning on descent while accounting for progress magnitude yields substantially more productive queries. Discussion. Despite these gains, several limitations point to natural extensions. First, as a purely local method, GRAPE inherits the well-known sensitivity to initialization: poor starting points or highly multimodal landscapes can trap the iterates in suboptimal basins. While trust-region restarts (Eriksson et al.,, 2019) and local–global switching schemes (McLeod et al.,, 2018; Diouane et al.,, 2023) provide principled escape mechanisms, integrating them with the two-stage loop requires careful budget allocation between exploration, exploitation, and restart phases. Second, the progress-aware score conditions on a Gaussian directional derivative, which assumes the local linearization is adequate over the chosen step size. In regions of high curvature, a truncated-Gaussian approximation may misestimate the true conditional decrease, suggesting that second-order corrections or adaptive step-size rules could improve robustness. Third, the cubic scaling of exact GP inference, (||3)O(|D|^3) per outer iteration, becomes a bottleneck as the evaluation budget grows. While rank-one and sparse updates (Csató and Opper,, 2002; Titsias,, 2009) alleviate this cost, they introduce additional approximation error whose interaction with the refinement–exploitation loop has not been characterized. We leave these directions to future work. References Arango et al., (2021) Arango, S. P., Jomaa, H. S., Wistuba, M., and Grabocka, J. (2021). HPO-B: A large-scale reproducible benchmark for black-box HPO based on OpenML. In Advances in Neural Information Processing Systems Datasets and Benchmarks Track. Berge, (1963) Berge, C. (1963). Topological Spaces. Oliver and Boyd. Berkenkamp et al., (2016) Berkenkamp, F., Schoellig, A. P., and Krause, A. (2016). Safe controller optimization for quadrotors with Gaussian processes. In IEEE International Conference on Robotics and Automation, pages 491–496. Bull, (2011) Bull, A. D. (2011). Convergence rates of efficient global optimization algorithms. Journal of Machine Learning Research, 12:2879–2904. Calamai and Moré, (1987) Calamai, P. H. and Moré, J. J. (1987). Projected gradient methods for linearly constrained problems. Mathematical Programming, 39(1):93–116. Carlini and Wagner, (2017) Carlini, N. and Wagner, D. (2017). Towards evaluating the robustness of neural networks. In IEEE Symposium on Security and Privacy, pages 39–57. Cheng et al., (2021) Cheng, S., Wu, G., and Zhu, J. (2021). On the convergence of prior-guided zeroth-order optimization algorithms. In Advances in Neural Information Processing Systems, volume 34, pages 14620–14631. Chew et al., (2026) Chew, R. W. T., Chen, Z., Hemachandra, A., and Low, B. K. H. (2026). BoLT: A benchmark to democratize black-box optimization research for expensive LLM tasks. arXiv preprint arXiv:2605.17000. Csató and Opper, (2002) Csató, L. and Opper, M. (2002). Sparse on-line Gaussian processes. Neural Computation, 14(3):641–668. Diouane et al., (2023) Diouane, Y., Picheny, V., Le Riche, R., and Scotto Di Perrotolo, A. (2023). TREGO: A trust-region framework for efficient global optimization. Journal of Global Optimization, 86(1):1–23. Doumont et al., (2026) Doumont, C., Fan, D., Maus, N., Gardner, J. R., Moss, H., and Pleiss, G. (2026). We still don’t understand high-dimensional Bayesian optimization. In International Conference on Artificial Intelligence and Statistics. Eriksson and Jankowiak, (2021) Eriksson, D. and Jankowiak, M. (2021). High-dimensional Bayesian optimization with sparse axis-aligned subspaces. In Conference on Uncertainty in Artificial Intelligence, volume 161 of Proceedings of Machine Learning Research, pages 493–503. PMLR. Eriksson et al., (2019) Eriksson, D., Pearce, M., Gardner, J., Turner, R. D., and Poloczek, M. (2019). Scalable global optimization via local Bayesian optimization. In Advances in Neural Information Processing Systems, volume 32. Fan et al., (2024) Fan, Z., Wang, W., Ng, S. H., and Hu, Q. (2024). Minimizing UCB: a better local search strategy in local Bayesian optimization. In Advances in Neural Information Processing Systems, volume 37, pages 130602–130634. Frazier, (2018) Frazier, P. I. (2018). Bayesian optimization. In Recent Advances in Optimization and Modeling of Contemporary Problems, INFORMS TutORials in Operations Research, pages 255–278. INFORMS. Gardner et al., (2018) Gardner, J. R., Pleiss, G., Bindel, D., Weinberger, K. Q., and Wilson, A. G. (2018). GPyTorch: Blackbox matrix-matrix Gaussian process inference with GPU acceleration. In Advances in Neural Information Processing Systems, volume 31. Garnett, (2023) Garnett, R. (2023). Bayesian Optimization. Cambridge University Press. Grimmett and Stirzaker, (2020) Grimmett, G. R. and Stirzaker, D. R. (2020). Probability and Random Processes. Oxford University Press, 4th edition. He et al., (2016) He, K., Zhang, X., Ren, S., and Sun, J. (2016). Deep residual learning for image recognition. In IEEE Conference on Computer Vision and Pattern Recognition, pages 770–778. Hendrycks et al., (2021) Hendrycks, D., Burns, C., Kadavath, S., Arora, A., Basart, S., Tang, E., Song, D., and Steinhardt, J. (2021). Measuring mathematical problem solving with the MATH dataset. In Advances in Neural Information Processing Systems Datasets and Benchmarks Track. Hennig et al., (2022) Hennig, P., Osborne, M. A., and Kersting, H. P. (2022). Probabilistic Numerics: Computation as Machine Learning. Cambridge University Press. Hvarfner et al., (2024) Hvarfner, C., Hellsten, E. O., and Nardi, L. (2024). Vanilla Bayesian optimization performs great in high dimensions. In International Conference on Machine Learning, volume 235 of Proceedings of Machine Learning Research, pages 20793–20817. PMLR. Jones et al., (1998) Jones, D. R., Schonlau, M., and Welch, W. J. (1998). Efficient global optimization of expensive black-box functions. Journal of Global Optimization, 13(4):455–492. Krizhevsky and Hinton, (2009) Krizhevsky, A. and Hinton, G. (2009). Learning multiple layers of features from tiny images. Technical report, University of Toronto. Kusupati et al., (2022) Kusupati, A., Bhatt, G., Rege, A., Wallingford, M., Sinha, A., Ramanujan, V., Howard-Snyder, W., Chen, K., Kakade, S., Jain, P., and Farhadi, A. (2022). Matryoshka representation learning. In Advances in Neural Information Processing Systems, volume 35, pages 30233–30249. LeCun et al., (1998) LeCun, Y., Bottou, L., Bengio, Y., and Haffner, P. (1998). Gradient-based learning applied to document recognition. Proceedings of the IEEE, 86(11):2278–2324. McLeod et al., (2018) McLeod, M., Roberts, S., and Osborne, M. A. (2018). Optimization, fast and slow: Optimally switching between local and Bayesian optimization. In International Conference on Machine Learning, volume 80 of Proceedings of Machine Learning Research, pages 3443–3452. PMLR. Močkus, (1978) Močkus, J. (1978). The application of Bayesian methods for seeking the extremum. In Dixon, L. C. W. and Szegö, G. P., editors, Towards Global Optimisation 2, pages 117–129. North-Holland. Müller et al., (2021) Müller, S., von Rohr, A., and Trimpe, S. (2021). Local policy search with Bayesian optimization. In Advances in Neural Information Processing Systems, volume 34, pages 20708–20720. Nguyen et al., (2022) Nguyen, Q., Wu, K., Gardner, J., and Garnett, R. (2022). Local Bayesian optimization via maximizing probability of descent. In Advances in Neural Information Processing Systems, volume 35, pages 13190–13202. Papenmeier et al., (2025) Papenmeier, L., Poloczek, M., and Nardi, L. (2025). Understanding high-dimensional Bayesian optimization. In International Conference on Machine Learning, volume 267 of Proceedings of Machine Learning Research, pages 47902–47923. PMLR. Paszke et al., (2019) Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., et al. (2019). PyTorch: An imperative style, high-performance deep learning library. In Advances in Neural Information Processing Systems, volume 32. Rasmussen and Williams, (2006) Rasmussen, C. E. and Williams, C. K. I. (2006). Gaussian Processes for Machine Learning. MIT Press. Ru et al., (2020) Ru, B., Cobb, A. D., Blaas, A., and Gal, Y. (2020). BayesOpt adversarial attack. In International Conference on Learning Representations. Schechter Vera et al., (2025) Schechter Vera, H. et al. (2025). EmbeddingGemma: Powerful and lightweight text representations. arXiv preprint arXiv:2509.20354. Shahriari et al., (2016) Shahriari, B., Swersky, K., Wang, Z., Adams, R. P., and de Freitas, N. (2016). Taking the human out of the loop: A review of Bayesian optimization. Proceedings of the IEEE, 104(1):148–175. Shields et al., (2021) Shields, B. J., Stevens, J., Li, J., Parasram, M., Damani, F., Martinez Alvarado, J. I., Janey, J. M., Adams, R. P., and Doyle, A. G. (2021). Bayesian reaction optimization as a tool for chemical synthesis. Nature, 590(7844):89–96. Shu et al., (2023) Shu, Y., Dai, Z., Sng, W., Verma, A., Jaillet, P., and Low, B. K. H. (2023). Zeroth-order optimization with trajectory-informed derivative estimation. In International Conference on Learning Representations. Snoek et al., (2012) Snoek, J., Larochelle, H., and Adams, R. P. (2012). Practical Bayesian optimization of machine learning algorithms. In Advances in Neural Information Processing Systems, volume 25. Srinivas et al., (2012) Srinivas, N., Krause, A., Kakade, S. M., and Seeger, M. W. (2012). Gaussian process optimization in the bandit setting: No regret and experimental design. Journal of Machine Learning Research, 13:3183–3212. Suwandi et al., (2025) Suwandi, R., Yin, F., Wang, J., Li, R., Chang, T.-H., and Theodoridis, S. (2025). Adaptive kernel design for Bayesian optimization is a piece of CAKE with LLMs. In Advances in Neural Information Processing Systems, volume 38, pages 132690–132723. Tang et al., (2026) Tang, W.-T., Kudva, A., and Paulson, J. A. (2026). NeST-BO: Fast local Bayesian optimization via Newton-step targeting of gradient and Hessian information. In International Conference on Artificial Intelligence and Statistics. Titsias, (2009) Titsias, M. K. (2009). Variational learning of inducing variables in sparse Gaussian processes. In Proceedings of the Twelfth International Conference on Artificial Intelligence and Statistics, volume 5 of Proceedings of Machine Learning Research, pages 567–574. PMLR. Wang and Dowling, (2022) Wang, K. and Dowling, A. W. (2022). Bayesian optimization for chemical products and functional materials. Current Opinion in Chemical Engineering, 36:100728. Wang et al., (2016) Wang, Z., Hutter, F., Zoghi, M., Matheson, D., and de Freitas, N. (2016). Bayesian optimization in a billion dimensions via random embeddings. Journal of Artificial Intelligence Research, 55:361–387. White et al., (2021) White, C., Neiswanger, W., and Savani, Y. (2021). BANANAS: Bayesian optimization with neural architectures for neural architecture search. In AAAI Conference on Artificial Intelligence, volume 35, pages 10293–10301. Wilcoxon, (1945) Wilcoxon, F. (1945). Individual comparisons by ranking methods. Biometrics Bulletin, 1(6):80–83. Wu et al., (2023) Wu, K., Kim, K., Garnett, R., and Gardner, J. R. (2023). The behavior and convergence of local Bayesian optimization. In Advances in Neural Information Processing Systems, volume 36, pages 73497–73523. Xu et al., (2025) Xu, Z., Wang, H., Phillips, J. M., and Zhe, S. (2025). Standard Gaussian process is all you need for high-dimensional Bayesian optimization. In International Conference on Learning Representations. Yang et al., (2025) Yang, A. et al. (2025). Qwen3 technical report. arXiv preprint arXiv:2505.09388. Zhou et al., (2023) Zhou, Y., Muresanu, A. I., Han, Z., Paster, K., Pitis, S., Chan, H., and Ba, J. (2023). Large language models are human-level prompt engineers. In International Conference on Learning Representations. Appendix This appendix contains supplementary material for the main text. We lead with the ablation that isolates each design choice, then provide protocols and robustness checks, followed by theoretical material and technical derivations. Component analysis. Section A reports an ablation experiment on the two main stages of our method: (1) the use of gradient refinement in the exploration stage, and (2) the use of progress-aware exploitation in the exploitation stage. This analysis isolates the contribution of each component. Protocols and robustness. Section B gives complete experimental details, including benchmark protocols, GP and method-specific hyperparameters, and compute information. Section C studies sensitivity to τexplore _explore, τexploit _exploit, η, and τthresh _thresh. Section D compares the robustness of GRAPE under Matérn and squared-exponential kernels. Section E reports wall-clock timings. Theory. Section F empirically validates Theorems 1 and 2. Section G provides the full proof of Theorem 2 (convergence of the progress-aware direction to steepest descent). Method details. Section H gives explicit componentwise expressions for the left- and right-acting kernel derivatives used in the gradient posterior. Section I proves the closed-form expression for the gradient-refinement acquisition (Proposition 1). Section J expands the connection between progress-aware exploitation and related criteria (MPD, unconditional expected descent, PI, and EI). Section K derives the analytic gradient of the progress-aware exploitation score used in the projected gradient ascent solver. Section L discusses the optimization procedure and computational cost. Appendix A Ablation Study To separate the two design choices in GRAPE, we compare the following four variants on the black-box adversarial attack and LLM prompt optimization tasks used in the main text: • MPD (baseline): original exploration and most-probable-descent exploitation (Nguyen et al.,, 2022). • MPD-Refine: Stage 1 uses gradient refinement; Stage 2 keeps MPD exploitation. • GRAPE-RandExp: Stage 1 draws uniform random queries (same budget τexplore _explore); Stage 2 uses progress-aware exploitation. • GRAPE (full): gradient refinement and progress-aware exploitation. All the experimental protocols including the hyperparameters, seeds, and budgets match the main text. Variant MNIST (#Q ↓ ) PO-128 (ℓT↓ _T ) PO-768 (ℓT↓ _T ) MPD (baseline) 827±128827± 128 −7.1±1.2-7.1± 1.2 −5.9±1.4-5.9± 1.4 MPD-Refine 741±119741± 119 −7.8±1.3-7.8± 1.3 −6.6±1.5-6.6± 1.5 GRAPE-RandExp 371±68371± 68 −11.8±1.6-11.8± 1.6 −12.6±1.7-12.6± 1.7 GRAPE (full) ±216± 42 −13.9±1.3-13.9± 1.3 −15.7±1.5-15.7± 1.5 Table 2: Ablation of gradient refinement vs. progress-aware exploitation (mean ± std over 1010 runs). Lower is better. MPD and full GRAPE match the main-paper settings; each hybrid changes one stage. MPD-Refine offers only a modest improvement over MPD, consistent with the findings of Nguyen et al., (2022), who observed that combining trace-based refinement with most-probable-descent moves did not reliably outperform GIBO. In contrast, GRAPE-RandExp closes much of the gap to full GRAPE across all three benchmarks, suggesting that the direction score is the more impactful of the two design choices. The improvement from GRAPE-RandExp to full GRAPE is smaller but consistent: αref _ref continues to provide a benefit when Stage 2 ranks directions by conditional progress, which leverages both μ_x and _x, instead of relying solely on descent probability. In summary, the two stages interact: refinement is most helpful when exploitation can utilize a sharper posterior, but provides minimal benefit to MPD alone. Appendix B Experimental Details B.1 Black-Box Adversarial Attacks Following Cheng et al., (2021), we attack correctly classified images from MNIST (LeCun et al.,, 1998) (d=28×28=784d=28× 28=784) and CIFAR-10 (Krizhevsky and Hinton,, 2009) (d=32×32=1024d=32× 32=1024) under an ℓ∞ _∞ perturbation constraint. • MNIST. We use the same fully trained CNN as Cheng et al., (2021) with ‖∞≤0.3\|x\|_∞≤ 0.3 and a query budget of 20002000 per attack. • CIFAR-10. We fully train a ResNet-18 (He et al.,, 2016) with SGD (cosine-annealed learning rate from 0.10.1 to 00, momentum 0.90.9, weight decay 5×10−45× 10^-4, 200200 epochs) and use ‖∞≤0.2\|x\|_∞≤ 0.2 with a query budget of 40004000. Observation model. All methods are strictly black-box: they may query only the classifier’s logits (pre-softmax scores). Hard predicted labels, true gradients, and internal network parameters are unavailable. Attack objective. Given a correctly classified image z with true label c, we search for a perturbation x that misclassifies +z+x. We minimize the Carlini–Wagner-style margin loss (Carlini and Wagner,, 2017; Cheng et al.,, 2021) L()=logitc(+)−maxj≠clogitj(+).L(x)=logit_c(z+x)- _j≠ clogit_j(z+x). (19) Intuitively, L()L(x) is the gap between the true-class logit and the strongest competing logit. Thus L()>0L(x)>0 means +z+x is still classified as c, while L()<0L(x)<0 means some other class scores higher and the attack has succeeded. The GP surrogate models L as a function of x over the feasible perturbation set =:‖∞≤ε,X= \x:\|x\|_∞≤ \, (20) with ε=0.3 =0.3 on MNIST and ε=0.2 =0.2 on CIFAR-10. Query counting and success. Each classifier forward pass that returns logits for a candidate x counts as one query. An attack is declared successful at the first queried point with L()<0L(x)<0; the reported query count is the number of evaluations up to and including that point (or the full budget if the attack fails). All methods minimize L over X under the same budgets above. Feasible-set projection. Exploitation steps that leave X are projected back by Euclidean projection onto the ℓ∞ _∞ ball. For GRAPE, this is the update ←Π(+η∗)x← _X(x+ ^*) described in Section 4.4 (with Π _X acting coordinate-wise by clipping to [−ε,ε]d[- , ]^d). Evaluation protocol. Reported query counts average over 1010 independent runs; success-rate curves are evaluated on 5050 randomly selected correctly classified images per dataset under a fixed budget of 20002000 queries. B.2 LLM Prompt Optimization We use the BoLT benchmark (Chew et al.,, 2026): 5,0145,014 prompts for mathematical reasoning, scored by MATH-500 accuracy (Hendrycks et al.,, 2021) under Qwen3-14B (Yang et al.,, 2025). Each prompt is embedded with EmbeddingGemma (Schechter Vera et al.,, 2025) at Matryoshka truncations d∈128,256,512,768d∈\128,256,512,768\ (Kusupati et al.,, 2022), yielding tasks PO-128 through PO-768. The search space is the discrete candidate pool PO=ii=15014⊂ℝdX_PO=\z_i\_i=1^5014 ^d, and fPO(i)f_PO(z_i) is the corresponding accuracy. Continuous proposals are projected to the nearest embedding in the pool. We run T=200T=200 iterations from five random initial prompts, with observation noise σN=0.001 _N=0.001, and report log simple regret over 1010 seeds. B.3 GP Hyperparameters All GP-based methods share the same modeling protocol unless noted otherwise. We use an ARD Matérn-5/25/2 kernel k(,′) k(x,x ) =s2(1+5r+53r2)exp(−5r), =s^2 (1+ 5\,r+ 53r^2 ) (- 5\,r), (21) r r =∑i=1d(xi−xi′)2ℓi2. = _i=1^d (x_i-x _i)^2 _i^2. (22) with per-dimension length scales ℓi _i and output scale s2s^2 estimated by maximum marginal likelihood (L-BFGS with 1010 restarts), matching BoTorch’s default Matérn surrogate with a gamma prior on the length scales. This choice is consistent with evidence that Matérn kernels are less prone to vanishing training gradients than squared-exponential (SE) kernels in high dimensions (Xu et al.,, 2025), and with the Matérn-5/25/2 option in the vanilla-BO definition of Hvarfner et al., (2024). For adversarial attacks the noise variance is fixed at σ2=10−4σ^2=10^-4; for prompt optimization we use σN=0.001 _N=0.001 as above. All methods start from the same initial points and share the same evaluation budget. An SE-kernel comparison appears in Section D. B.4 Hyperparameters for GRAPE The default hyperparameters for GRAPE are summarized in Table 3. Adaptive schedules (when enabled) can scale τexplore _explore, τthresh _thresh, and τexploit _exploit with dimension, GP length scale, dataset size, and recent improvement. All reported results use the fixed settings above (_= adaptive\_params=False). Parameter Value Refinement budget τexplore _explore 55 Exploitation budget τexploit _exploit 3030 (max) Step size η 0.1⋅diam()0.1·diam(X) Early-stop threshold τthresh _thresh 5×10−35× 10^-3 Direction restarts 1010 (with warm start) Refinement candidates NcandN_cand 100100 (uniform in X) Initial dataset size ninitn_init 1010 (attacks); 55 (prompt opt.) Outer iterations T until query budget (attacks); 200200 (prompt opt.) Table 3: Hyperparameters for GRAPE. B.5 Baseline Hyperparameters All baselines use the same Matérn-5/25/2 GP and marginal-likelihood training protocol. • MPD (Nguyen et al.,, 2022): 55 exploration samples per iterate; descent-probability threshold 0.650.65; step size η=0.01⋅diam()η=0.01·diam(X); 1010 random restarts for direction optimization. • GIBO (Müller et al.,, 2021): 55 gradient-information samples per iterate; local radius δ=0.1δ=0.1; learning rate 0.10.1; normalized mean-gradient steps. • MinUCB (Fan et al.,, 2024): local ball radius 0.1⋅diam()0.1·diam(X); exploration interval 55; UCB parameter β=3β=3; multi-start L-BFGS for local UCB minimization. • TuRBO (Eriksson et al.,, 2019): 55 trust regions; ARD length scales; success/failure thresholds and restart criteria as in the original paper. • NeST-BO (Tang et al.,, 2026): jointly learned gradient and Hessian posteriors targeting a modified Newton step; default settings from the authors’ implementation. • VBO: vanilla BO with Expected Improvement (Jones et al.,, 1998); multi-start L-BFGS acquisition optimization (1010 restarts, 512512 raw samples). • D-LogEI (Hvarfner et al.,, 2024): Expected Improvement with a dimension-scaled log-normal prior on length scales. • Random: uniform sampling within the feasible set (continuous ℓ∞ _∞ ball for attacks; uniform over the prompt pool for BoLT). B.6 Compute Resources All experiments run on an Ubuntu Linux workstation with an AMD Ryzen 9 7950X CPU (16 cores / 32 threads), 64 GB system RAM, and an NVIDIA GeForce RTX 4090 GPU (24 GB VRAM; driver 595.58.03595.58.03, CUDA 13.213.2). Software versions follow requirements.txt: PyTorch 2.13.02.13.0 (Paszke et al.,, 2019), GPyTorch 1.15.21.15.2 (Gardner et al.,, 2018), and BoTorch 0.18.10.18.1. Randomness is controlled by setting torch.manual_seed and numpy.random.seed to the trial seed s∈0,…,9s∈\0,…,9\ (or base_seed+trial with base_seed=42=42). Wall-clock time is dominated by GP refits and, for attacks, classifier forward passes. Per-method wall-clock measurements are reported in Section E. Appendix C Hyperparameter Sensitivity We sweep one hyperparameter at a time around the defaults in Table 3, holding the others fixed (_= adaptive\_params=False), on two representative tasks: MNIST and PO-128. Hyperparameter Value MNIST (#Q ↓ ) PO-128 (ℓT↓ _T ) τexplore _explore 11 287±61287± 61 −11.4±1.5-11.4± 1.5 33 238±49238± 49 −13.1±1.4-13.1± 1.4 55 (default) ±216± 42 −13.9±1.3-13.9± 1.3 1010 241±51241± 51 −14.0±1.4-14.0± 1.4 τexploit _exploit (max) 1010 268±58268± 58 −12.1±1.6-12.1± 1.6 2020 229±47229± 47 −13.4±1.4-13.4± 1.4 3030 (default) ±216± 42 −13.9±1.3-13.9± 1.3 5050 211±44211± 44 −14.0±1.3-14.0± 1.3 η/diam()η\,/\,diam(X) 0.010.01 392±74392± 74 −10.6±1.7-10.6± 1.7 0.050.05 251±53251± 53 −12.9±1.5-12.9± 1.5 0.10.1 (default) ±216± 42 −13.9±1.3-13.9± 1.3 0.20.2 274±62274± 62 −12.4±1.6-12.4± 1.6 τthresh _thresh 10−310^-3 209±45209± 45 −14.1±1.4-14.1± 1.4 5×10−35× 10^-3 (default) ±216± 42 −13.9±1.3-13.9± 1.3 10−210^-2 244±52244± 52 −13.0±1.5-13.0± 1.5 5×10−25× 10^-2 318±69318± 69 −10.9±1.8-10.9± 1.8 Table 4: One-at-a-time sensitivity of GRAPE on MNIST and PO-128 (mean ± std over 1010 runs). Default rows match Table 3 and are bolded within each block. Table 4 is consistent with defaults that lie in a reasonably flat region. Varying τexplore _explore over 3,5,10\3,5,10\ changes MNIST queries by only a few dozen evaluations and keeps PO-128 within about one log-unit of the default; τexplore=1 _explore=1 is the clear under-refinement failure case. The exploitation cap behaves similarly once it is large enough for early stopping to dominate: 3030 and 5050 are nearly indistinguishable under this protocol. The step size η is the most sensitive setting in the grid, where 0.010.01 is too conservative and 0.20.2 overshoots, but 0.050.05 and 0.10.1 remain close. For τthresh _thresh, values near 5×10−35× 10^-3 look similar, whereas 5×10−25× 10^-2 truncates exploitation early and hurts prompt-opt regret. We did not sweep NcandN_cand because Stage 1 scoring is cheap relative to black-box queries and the default matches TuRBO’s candidate budget; see Section L for the full protocol. Appendix D Kernel Robustness: Matérn vs. SE Main results use ARD Matérn-5/25/2. Following Xu et al., (2025), who show that SE kernels are more fragile under default length-scale initialization in high dimensions, we re-run GRAPE (and, for reference, MPD) with an ARD squared-exponential kernel under the same MLE fitting protocol, on MNIST and PO-128. Method Kernel MNIST (#Q ↓ ) PO-128 (ℓT↓ _T ) GRAPE Matérn-5/25/2 (default) ±216± 42 −13.9±1.3-13.9± 1.3 GRAPE SE / RBF 268±58268± 58 −12.2±1.6-12.2± 1.6 MPD Matérn-5/25/2 827±128827± 128 −7.1±1.2-7.1± 1.2 MPD SE / RBF 961±151961± 151 −5.9±1.5-5.9± 1.5 Table 5: Kernel robustness: Matérn-5/25/2 vs. squared-exponential (mean ± std over 1010 runs). Lower is better. Switching from Matérn-5/25/2 to SE worsens both methods on both metrics, which is consistent with the high-dimensional SE fragility discussed by Xu et al., (2025). Under our local protocol the degradation is limited rather than severe: GRAPE’s MNIST query count rises from 216216 to 268268, and PO-128 final regret worsens by roughly 1.71.7 log-units; MPD shifts in the same direction. The method ranking is unchanged on these two tasks, so the main-paper comparison is not overturned by the kernel swap, but we treat this only as a sanity check on two datasets rather than a general claim that GRAPE is kernel-agnostic. We keep Matérn-5/25/2 as the default to match the main experiments. Appendix E Wall-Clock Timing We report wall-clock cost on the same hardware as Section B (AMD Ryzen 9 7950X, NVIDIA RTX 4090). For each method we measure (i) total runtime per trial and (i) mean per-iteration overhead of surrogate updates and direction / acquisition optimization, excluding the black-box query itself when that query is a classifier forward pass or a precomputed prompt score. Attack timings use MNIST (d=784d=784); prompt timings use PO-128 (T=200T=200). On MNIST, a trial stops at attack success, so total time scales with the number of queries issued; on PO-128 every method runs a fixed iteration budget. MNIST PO-128 Method Total (s) Opt. / iter (s) Total (s) Opt. / iter (s) GRAPE 94±1794± 17 0.41±0.070.41± 0.07 88±1388± 13 0.37±0.060.37± 0.06 NeST-BO 1370±2401370± 240 1.91±0.241.91± 0.24 421±47421± 47 1.98±0.271.98± 0.27 MPD 251±41251± 41 0.30±0.050.30± 0.05 69±1169± 11 0.27±0.050.27± 0.05 MinUCB 420±70420± 70 0.47±0.080.47± 0.08 97±1497± 14 0.42±0.070.42± 0.07 GIBO 355±60355± 60 0.36±0.060.36± 0.06 81±1281± 12 0.33±0.060.33± 0.06 TuRBO 620±105620± 105 0.54±0.100.54± 0.10 119±17119± 17 0.51±0.090.51± 0.09 VBO 1810±3101810± 310 1.14±0.171.14± 0.17 251±29251± 29 1.09±0.161.09± 0.16 D-LogEI 1820±3001820± 300 1.32±0.201.32± 0.20 289±33289± 33 1.26±0.181.26± 0.18 Random 51±951± 9 — 7±17± 1 — Table 6: Wall-clock timing (mean ± std over 1010 runs). “Opt. / iter” is surrogate fit plus acquisition / direction optimization, excluding black-box queries. MNIST totals accumulate this overhead over the queries-to-success from the main text; PO-128 totals use a fixed budget T=200T=200. NeST-BO incurs additional Hessian-related cost; VBO / D-LogEI pay for global multi-start acquisition over the full domain. GRAPE’s overhead is dominated by GP refits and candidate scoring of αref _ref / projected gradient ascent for P. On a per-iteration basis, GRAPE sits with the other first-order local methods, close to GIBO / MPD and well below NeST-BO, whose joint gradient–Hessian updates dominate. Global acquisitions (VBO, D-LogEI) are slower still per step because each iteration runs multi-start L-BFGS over the full high-dimensional domain. On MNIST, total time additionally tracks query count: GRAPE finishes in under two minutes mainly because it succeeds with far fewer queries, while VBO / D-LogEI combine a high per-step cost with long runs. On PO-128 the iteration count is fixed, so totals largely mirror per-iteration cost; GRAPE remains close to MPD / GIBO and several times faster than NeST-BO. We do not claim wall-clock as a primary contribution, query efficiency is the main metric, only that Stage 1/Stage 2 does not make GRAPE an outlier among first-order local baselines on this hardware. Appendix F Empirical Validation of Theoretical Claims We complement the proofs of Theorems 1 and 2 with two controlled synthetic checks that isolate each claim. Both experiments use the smooth quadratic f()=12‖2+0.15z0z1f(z)= 12\|z\|^2+0.15\,z_0z_1 in d=2d=2 and report mean ± standard deviation over 1010 random seeds. F.1 Theorem 1: Monotone Gradient-Uncertainty Contraction Theorem 1 states that each gradient-refinement query weakly decreases the total gradient uncertainty Tr() ( _x), with a strict decrease whenever ∇k(,)≠ _xk_D(x,z) 0. At a fixed interior point x, we fit a squared-exponential GP (lengthscale 0.550.55, outputscale 11, noise variance 10−410^-4) on three random initialization points, then run eight greedy refinement steps that each select ⋆=argmaxαref()z = _z _ref(z) from Ncand=100N_cand=100 candidates sampled in a local neighborhood of x and append the noisy observation. Figure 4 plots Tr() ( _x) against the refinement step. The trace decreases at every step and shows diminishing returns as the local gradient posterior concentrates, matching Theorem 1. Figure 4: Empirical validation of Theorem 1. Total gradient uncertainty Tr() ( _x) at a fixed interior point decreases monotonically under greedy αref _ref refinement (mean ± std over 1010 seeds). F.2 Theorem 2: Convergence to Steepest Descent Theorem 2 states that if Tr()→0 ( _x)→ 0 and →∇f() μ_x→∇ f(x), then the progress-aware direction converges to true steepest descent, ⋆→sd:=−∇f()/∥∇f()∥v _sd:=-∇ f(x)/\|∇ f(x)\|. To probe this limit directly, we instantiate the theorem’s hypotheses at the same interior point x: set the gradient posterior mean to a noisy perturbation of the analytic ∇f()∇ f(x) and the covariance to a shrinking anisotropic matrix t=εtAt _t= _tA_t with εt→0 _t→ 0. At each step we compute ⋆=argmax‖=1()v = _\|v\|=1P(v) and record the angle ∠(⋆,sd) (v ,v_sd). Figure 5 shows that this angle decays toward zero as the posterior concentrates, in agreement with Theorem 2. Figure 5: Empirical validation of Theorem 2. Angle between the progress-aware direction ⋆v and true steepest descent sdv_sd decreases toward zero as Tr()→0 ( _x)→ 0 and →∇f() μ_x→∇ f(x) (mean ± std over 1010 seeds). Appendix G Full Proof of Theorem 2 We restate the theorem for convenience. Theorem (Progress-aware exploitation approaches steepest descent). Suppose ∇f()≠∇ f(x) 0, Tr()→0 ( _x)→ 0, and →∇f() μ_x→∇ f(x). Then the progress-aware direction converges to the true steepest descent direction: limTr()→0∗=−∇f()‖∇f()‖. _ ( _x)→ 0\;v^*\;=\;- ∇ f(x)\|∇ f(x)\|. Proof. As Tr()→0 ( _x)→ 0, we have → _x 0 in spectral norm. Together with the assumption →∇f() μ_x→∇ f(x), the directional moments converge for every unit v: μ→c:=⊤∇f() _v→ c_v:=v ∇ f(x) and σ→0 _v→ 0. Step 1: Pointwise convergence of J()J(v). Fix a unit vector v. Write J()=()=σ[R(γ)−γ]J(v)=P(v)= _v [R( _v)- _v ]. Case 1: c<0c_v<0. Then γ=μ/σ→−∞ _v= _v/ _v→-∞. Since ϕ(γ)→0φ(γ)→ 0 and Φ(−γ)→1 (-γ)→ 1 as γ→−∞γ→-∞, we have R(γ)→0R( _v)→ 0. Therefore J()=σR(γ)−μ→0−c=−c>0.J(v)= _vR( _v)- _v→ 0-c_v=-c_v>0. Case 2: c>0c_v>0. Then γ→+∞ _v→+∞. Using the asymptotic expansion R(γ)=γ+γ−1+(γ−3)R(γ)=γ+γ^-1+O(γ^-3) for large γ Grimmett and Stirzaker, (2020): J()=σ[γ−1+(γ−3)]=σ2μ+(σ4μ3)→0.J(v)= _v [ _v^-1+O( _v^-3) ]= _v^2 _v+O\! ( _v^4 _v^3 )→ 0. Case 3: c=0c_v=0. We have J()≤σR(0)=σ2/π→0J(v)≤ _vR(0)= _v 2/π→ 0. Combining all three cases, J() J(v) →Tr()→0[−c]+ \; ( _x)→ 0\;[-c_v]^+ =[−⊤∇f()]+pointwise on d−1. =[-v ∇ f(x)]^+ on S^d-1. Step 2: Uniform convergence. J()J(v) is jointly Lipschitz in v and in (,)( μ_x, _x), since R is smooth on any compact set away from Φ(−γ)=0 (-γ)=0, and Φ(−γ)≥δ>0 (- _v)≥δ>0 uniformly over the sphere once σ _v is small enough relative to min|c| _v|c_v| (away from the boundary of the descent hemisphere). Because (,)( μ_x, _x) converges under the stated assumptions on compact X, the convergence of J()J(v) is uniform over d−1S^d-1. Step 3: Convergence of the argmax. The limit function g():=[−⊤∇f()]+g(v):=[-v ∇ f(x)]^+ is continuous on the compact sphere d−1S^d-1. By uniform convergence of J()→g()J(v)→ g(v) and the Berge Maximum Theorem Berge, (1963), the argmax correspondence is upper hemicontinuous: any cluster point of ∗\v^*\ lies in argmaxg() _vg(v). Step 4: Uniqueness of the argmax of g. Since ∇f()≠∇ f(x) 0, the function g()=−⊤∇f()g(v)=-v ∇ f(x) is strictly positive on the open descent hemisphere :⊤∇f()<0\v:v ∇ f(x)<0\ and zero elsewhere. It is uniquely maximized at ∗=−∇f()/∥∇f()∥v^*=-∇ f(x)/\|∇ f(x)\|, because g()=−⊤∇f()≤∥∇f()∥g(v)=-v ∇ f(x)≤\|∇ f(x)\| with equality iff =−∇f()/∥∇f()∥v=-∇ f(x)/\|∇ f(x)\|. Therefore the argmax is a singleton and ∗→−∇f()/∥∇f()∥v^*→-∇ f(x)/\|∇ f(x)\|. ∎ Appendix H Explicit Kernel Derivatives for the Gradient Posterior The main text writes the gradient posterior mean and covariance using the compact operators ∇k∇ k and k∇⊤k∇ . This section derives the corresponding explicit componentwise definitions. Let k be a twice-differentiable kernel on ⊆ℝdX ^d, and write =(x1,…,xd)⊤x=(x_1,…,x_d) and ′=(x1′,…,xd′)⊤x =(x _1,…,x _d) . A differential operator placed before k acts on its first argument, while an operator placed after k acts on its second: [∇k(,′)]i [∇ k(x,x ) ]_i =∂xik(,′),i=1,…,d, = ∂ x_i\,k(x,x ), i=1,…,d, (23) [k(,′)∇⊤]j [k(x,x )∇ ]_j =∂xj′k(,′),j=1,…,d. = ∂ x _j\,k(x,x ), j=1,…,d. (24) Thus ∇k(,′)∈ℝd∇ k(x,x ) ^d is a column vector and k(,′)∇⊤∈ℝ1×dk(x,x )∇ ^1× d is a row vector. The mixed second-derivative matrix appearing in the prior gradient covariance is [∇k(,′)∇⊤]i,j=∂2∂xi∂xj′k(,′),i,j=1,…,d, [∇ k(x,x )∇ ]_i,j= ∂^2∂ x_i\,∂ x _j\,k(x,x ), i,j=1,…,d, (25) so ∇k(,)∇⊤∈ℝd×d∇ k(x,x)∇ ^d× d. When the second (resp. first) argument is a training set =(1,…,n)⊤X=(x_1,…,x_n) , the same convention yields matrices. Writing k(,)∈ℝ1×nk(x,X) ^1× n for the row vector of covariances and k(,)∈ℝnk(X,x) ^n for the corresponding column vector, [∇k(,)]i,ℓ [∇ k(x,X) ]_i, =∂xik(,ℓ),∇k(,)∈ℝd×n, = ∂ x_i\,k(x,x_ ), ∇ k(x,X) ^d× n, (26) [k(,)∇⊤]ℓ,j [k(X,x)∇ ]_ ,j =∂xjk(ℓ,),k(,)∇⊤∈ℝn×d. = ∂ x_j\,k(x_ ,x), k(X,x)∇ ^n× d. (27) By symmetry of k, k(,)∇⊤=(∇k(,))⊤k(X,x)∇ = (∇ k(x,X) ) . Substituting Eqs. (23)–(27) into the conditioning formulas of the main text recovers the usual joint GP over function values and gradients (Rasmussen and Williams,, 2006): μ_x =∇μ()+∇k(,)−1(−μ()), =∇μ(x)+∇ k(x,X)\,K^-1 (y-μ(X) ), (28) _x =∇k(,)∇⊤−∇k(,)−1k(,)∇⊤, =∇ k(x,x)∇ -∇ k(x,X)\,K^-1\,k(X,x)∇ , (29) with =k(,)+σ2K=k(X,X)+σ^2I. Appendix I Proof of Proposition 1 (Closed-Form Gradient Refinement) We restate the claim for convenience. Recall that the gradient-refinement acquisition of the main text is the expected one-step reduction in total gradient variance at the current iterate x, αref()=y[Tr()−Tr(|∪(,y))], _ref(z)=E_y_z\! [ ( _x)- \! ( _x|D∪(z,y_z) ) ], (30) where y=f()+εy_z=f(z)+ with ε∼(0,σ2) (0,σ^2) is the yet-unobserved noisy evaluation at a candidate ∈z , and |∪(,y) _x|D∪(z,y_z) denotes the gradient covariance at x after augmenting D with (,y)(z,y_z). Proposition (Closed-form gradient refinement). Let k(⋅,⋅)k_D(·,·) denote the GP posterior covariance after observing D. Then αref() _ref(z) is independent of the unobserved value y_z and equals αref()=‖∇k(,)‖2k(,)+σ2, _ref(z)= \| _xk_D(x,z)\|^2k_D(z,z)+σ^2, (31) where ∇k(,)=∂k(,)/∂ _xk_D(x,z)=∂ k_D(x,z)/ . Proof. Write k_D for the posterior covariance after conditioning on D. The gradient posterior covariance at x is the mixed Hessian of this kernel, =∇′⊤k(,′)|′=. _x= _x _x k_D(x,x ) |_x =x. (32) Now augment D by a single noisy observation at z. The standard rank-one GP covariance update (Rasmussen and Williams,, 2006) gives the updated posterior covariance k′(,′)=k(,′)−k(,)k(,′)k(,)+σ2,k_D (x,x )=k_D(x,x )- k_D(x,z)\,k_D(z,x )k_D(z,z)+σ^2, (33) which depends on the location z but not on the realized value y_z. Applying Eq. (32) to k′k_D therefore yields |∪(,y) _x|D∪(z,y_z) =∇′⊤k′(,′)|′= = _x _x k_D (x,x ) |_x =x =−(∇k(,))(∇′k(,′))⊤|′=k(,)+σ2. = _x- ( _xk_D(x,z) ) ( _x k_D(z,x ) ) |_x =xk_D(z,z)+σ^2. (34) By symmetry of k_D, ∇′k(,′)|′==∇k(,) _x k_D(z,x ) |_x =x= _xk_D(x,z). Writing :=∇k(,)∈ℝdg_z:= _xk_D(x,z) ^d for this common vector, Eq. (34) simplifies to the rank-one downdate |∪(,y)=−⊤k(,)+σ2. _x|D∪(z,y_z)= _x- g_zg_z k_D(z,z)+σ^2. (35) In particular, the right-hand side does not involve y_z, so the expectation in Eq. (30) is vacuous: αref()=Tr()−Tr(|∪(,y)). _ref(z)= ( _x)- \! ( _x|D∪(z,y_z) ). (36) Substituting Eq. (35) and using Tr(⊤)=‖2 (g_zg_z )=\|g_z\|^2 gives αref()=‖2k(,)+σ2=‖∇k(,)‖2k(,)+σ2, _ref(z)= \|g_z\|^2k_D(z,z)+σ^2= \| _xk_D(x,z)\|^2k_D(z,z)+σ^2, (37) which is Eq. (31). ∎ The closed form makes the information-theoretic content of αref _ref explicit: the numerator rewards candidates whose posterior cross-covariance gradient with x is large (so observing f()f(z) is informative about ∇f()∇ f(x)), while the denominator down-weights locations that are already well determined by D. Appendix J Connections to Related Acquisition Criteria This section expands the relation between the progress-aware score ()P(v) and three closely related quantities: the most-probable descent criterion of MPD (Nguyen et al.,, 2022), the unconditional expected descent of a directional derivative, and the classical Probability of Improvement (PI) / Expected Improvement (EI) pair (Močkus,, 1978; Jones et al.,, 1998). The goal is to make precise which factor each criterion maximizes, and why the resulting directions generally disagree. J.1 Setup and Notation Recall that under a twice-differentiable GP prior, the directional derivative along a unit vector v is univariate Gaussian, ∇f()=⊤∇f()∼(μ,σ2), _vf(x)\;=\;v ∇ f(x)\; \;N( _v, _v^2), (38) with μ=⊤ _v=v μ_x, σ=⊤ _v= v _xv, and standardized mean γ=μ/σ _v= _v/ _v. Write ϕφ and Φ for the standard normal PDF and CDF, and define the inverse Mills ratio R(γ):=ϕ(γ)Φ(−γ).R(γ)\;:=\; φ(γ) (-γ). (39) Throughout, (a)+=max(a,0)(a)_+= (a,0). J.2 Derivation of the Progress-Aware Score For Z∼(μ,σ2)Z (μ,σ^2), the mean of the lower-truncated law at zero is the standard truncated-Gaussian identity [Z∣Z<0]=μ−σϕ(μ/σ)Φ(−μ/σ)=μ−σR(μ/σ).E[Z Z<0]=μ-σ\, φ(μ/σ) (-μ/σ)=μ-σ\,R(μ/σ). (40) Taking Z=∇f()Z= _vf(x) and negating both sides yields () (v) :=[−∇f()∣∇f()<0] :=E [- _vf(x) _vf(x)<0 ] =−μ+σR(γ)=σ[R(γ)−γ], =- _v+ _v\,R( _v)= _v [R( _v)- _v ], (41) which is Definition 1. Two elementary properties follow immediately: 1. ()>0P(v)>0 whenever σ>0 _v>0, because the Mills ratio satisfies R(γ)>γR(γ)>γ for all finite γ when the conditioning event has positive probability under a nondegenerate Gaussian. 2. ()P(v) is homogeneous of degree one in the directional scale (μ,σ)( _v, _v): if (μ,σ)↦(cμ,cσ)( _v, _v) (c _v,c _v) for c>0c>0, then P scales by c. Thus the score measures a magnitude of conditional decrease, not a dimensionless probability. J.3 Decomposition into Probability and Conditional Magnitude The unconditional expected descent of a directional derivative is the first truncated moment without conditioning: [(−∇f())+] \! [(- _vf(x))_+ ] (42) =[−∇f()⋅∇f()<0]. =E [- _vf(x)·1_\ _vf(x)<0\ ]. By the law of total expectation, we obtain [(−∇f())+] \! [(- _vf(x))_+ ] =ℙ(∇f()<0)[−∇f()∣∇f()<0] =P( _vf(x)<0)\,E [- _vf(x) _vf(x)<0 ] =Φ(−γ)(). = (- _v)\,P(v). (43) This factorization is the key structural observation. It separates three natural criteria on the same Gaussian directional posterior: MPD maximizes the descent probability Φ(−γ) (- _v); progress-aware exploitation maximizes the conditional magnitude ()P(v); and unconditional expected descent maximizes their product Φ(−γ)() (- _v)\,P(v). Because the first factor is a probability in [0,1][0,1] and the second is a scale-dependent magnitude, maximizing either factor alone is not equivalent to maximizing their product. In particular: • MPD is invariant to positive rescaling of (μ,σ)( _v, _v) and therefore ignores how large a decrease would be if descent occurs. • Progress-aware exploitation is deliberately not probability-weighted: a direction with modest descent probability can still win if its negative tail is heavy enough that ()P(v) is large. • Unconditional expected descent recovers a probability-weighted magnitude, analogous in spirit to EI (see below), but still acts on directional derivatives rather than on objective values. J.4 Worked Example: Scale Blindness and Ranking Disagreement We first make the scale-invariance of descent probability concrete. Consider two unit directions with identical standardized means γ=γ=−5 _u= _w=-5, ∇f() _uf(x) ∼(−0.01,0.0022), (-0.01,0.002^2), (44) ∇f() _wf(x) ∼(−5,1). (-5,1). Both receive the same descent probability Φ(5)≈1 (5)≈ 1, so MPD is indifferent, yet their conditional magnitudes differ by orders of magnitude: ()≈0.01P(u)≈ 0.01 while ()=[R(−5)+5]≈5P(w)= [R(-5)+5 ]≈ 5. This is exactly the scale discarded by Φ(−γ) (- _v). The three criteria can also disagree in ranking. Take ∇1f() _v_1f(x) ∼(−0.01,0.0022), (-0.01,0.002^2), (45) ∇2f() _v_2f(x) ∼(0.5,1). (0.5,1). For 1v_1, γ1=−5 _1=-5, so Φ(−γ1)=Φ(5)≈1 (- _1)= (5)≈ 1 and (1)=0.002[R(−5)+5]≈0.002⋅5=0.01,P(v_1)=0.002 [R(-5)+5 ]≈ 0.002· 5=0.01, (46) since R(−5)=ϕ(−5)/Φ(5)R(-5)=φ(-5)/ (5) is negligible. Thus MPD strongly prefers 1v_1, while the conditional magnitude is only about 0.010.01. For 2v_2, γ2=0.5 _2=0.5, so Φ(−γ2)=Φ(−0.5)≈0.309 (- _2)= (-0.5)≈ 0.309 and (2)=1⋅[R(0.5)−0.5]=ϕ(0.5)Φ(−0.5)−0.5≈0.64.P(v_2)=1· [R(0.5)-0.5 ]= φ(0.5) (-0.5)-0.5≈ 0.64. (47) Progress-aware exploitation therefore prefers 2v_2. The unconditional expected descents are [(−∇1f)+] [(- _v_1f)_+] ≈1⋅0.01=0.01, ≈ 1· 0.01=0.01, (48) [(−∇2f)+] [(- _v_2f)_+] ≈0.309⋅0.64≈0.20, ≈ 0.309· 0.64≈ 0.20, (49) so the product criterion also prefers 2v_2, but by a smaller margin than P alone. Probability-only ranking is the most conservative of the three: it is blind to slope scale when γ is fixed, and can prefer a nearly certain but tiny direction over an uncertain but large-tailed one. J.5 Analogy with Probability of Improvement and Expected Improvement Classical BO acquisitions act on the objective-value posterior f()∼(μ(),σ()2)f(z) (μ(z),σ(z)^2) at a candidate z, relative to the best observed value f∗f_*. With γimp()=(f∗−μ())/σ() _imp(z)=(f_*-μ(z))/σ(z) for minimization (or the sign-flipped analogue for maximization), PI() (z) =Φ(γimp()), = ( _imp(z) ), (50) EI() (z) =σ()[γimp()Φ(γimp())+ϕ(γimp())] =σ(z)\, [ _imp(z)\, ( _imp(z) )+φ ( _imp(z) ) ] (51) =ℙ(f()<f∗) =P(f(z)<f_*) ⋅[f∗−f()∣f()<f∗]. ·E [f_*-f(z) f(z)<f_* ]. Eq. (51) is again an unconditional truncated-Gaussian moment: EI factors as “probability of improvement × conditional magnitude of improvement,” exactly as Eq. (43) factors unconditional expected descent (Jones et al.,, 1998; Močkus,, 1978). The directional analogues are therefore PI ↔ MPD (probability only), conditional improvement ↔ ()P(v) (magnitude given the event), and EI ↔ unconditional expected descent (probability × magnitude). Progress-aware exploitation is thus closer to the conditional factor of EI than to EI itself: it deliberately discards the probability weight Φ(−γ) (- _v) and retains only the conditional magnitude. Two distinctions remain important: 1. Domain of the random variable. PI/EI truncate the posterior of f()f(z) about a baseline f∗f_*; P truncates the posterior of a directional derivative about zero. The shared algebra is truncated-Gaussian, but the optimized quantity is different. 2. Role in the algorithm. PI/EI are typically used as global acquisition functions over X. Progress-aware exploitation is a local direction-selection criterion applied after gradient refinement has already concentrated (,)( μ_x, _x). Theorem 2 then shows that, as Tr()→0 ( _x)→ 0 with consistent mean, P recovers normalized steepest descent, a limit that has no direct PI/EI analogue. In short, MPD stands to progress-aware exploitation as PI stands to the conditional half of EI, while unconditional expected descent stands to P as EI stands to that same conditional half. GRAPE optimizes the conditional magnitude on directional derivatives, preceded by closed-form gradient refinement. Appendix K Gradient of the Progress-Aware Score To solve ∗=argmax‖=1()v^*= _\|v\|=1P(v) (52) via projected gradient ascent (Calamai and Moré,, 1987), we derive the analytic gradient of J()=()=σ[R(γ)−γ]J(v)=P(v)= _v [R( _v)- _v ] with respect to v, where R(γ):=ϕ(γ)/Φ(−γ)R(γ):=φ(γ)/ (-γ) is the inverse Mills ratio and γ=μ/σ _v= _v/ _v. Recall μ=⊤ _v=v μ_x, σ=⊤ _v= v _xv, and γ=μ/σ _v= _v/ _v. Their partial derivatives are ∂σ∂ ∂ _v =σ, = _xv _v, (53) ∂γ∂ ∂ _v =σ−μ(/σ)σ2=σ−γσ2. = μ_x _v- _v( _xv/ _v) _v^2= _v μ_x- _v _xv _v^2. (54) Using the identity R′(γ)=−R(γ)[γ+R(γ)]R (γ)=-R(γ) [γ+R(γ) ] Grimmett and Stirzaker, (2020), the chain rule gives ∇J() _vJ(v) =∂σ∂[R(γ)−γ]+σ∂[R(γ)−γ] = ∂ _v [R( _v)- _v ]+ _v ∂ [R( _v)- _v ] =σ[R(γ)−γ] = _xv _v [R( _v)- _v ] +σR′(γ)−1σ−γσ2 + _v \R ( _v)-1 \ _v μ_x- _v _xv _v^2 =σ[R(γ)−γ] = _xv _v [R( _v)- _v ] −σ−γσ2 - _v μ_x- _v _xv _v^2 ×R(γ)[γ+R(γ)]+1. × \R( _v) [ _v+R( _v) ]+1 \. (55) Projected gradient ascent. At each iteration update ′←+α∇J()v +α _vJ(v) and project back to the unit sphere: ←′/‖′‖v /\|v \|, where α>0α>0 is the step size. We use 1010 random restarts with warm-starting from the previous iteration’s best direction, and return the unit vector achieving the highest P value. The step size is set by Armijo backtracking line search with initial value α0=0.1 _0=0.1 and reduction factor 0.50.5. Appendix L Optimization and Computational Cost The sphere-constrained maximization of the progress-aware score Eq. (52) has no general closed form. We optimize ()P(v) with projected gradient ascent, 1010 random restarts, and warm-starting from the previous iterate (Section K). Gradient refinement is deterministic once a query location is chosen: αref()=‖∇k(,)‖2/(k(,)+σ2) _ref(z)=\| _xk_D(x,z)\|^2/(k_D(z,z)+σ^2) requires no Monte Carlo integration. Maximizing αref _ref in the reported experiments. In all reported runs, we maximize αref() _ref(z) by evaluating the closed-form expression on Ncand=100N_cand=100 candidates sampled uniformly from the feasible set X, and selecting the highest-scoring candidate. Each αref _ref evaluation is inexpensive relative to the subsequent GP hyperparameter refit and black-box query, so we use this lightweight random scoring rather than continuous multi-start optimization over z. Continuous or local multi-start maximization of αref _ref, for example, L-BFGS initialized from the top random candidates in a neighborhood of x, remains a drop-in upgrade and was not used for the main-paper numbers. Each outer iteration requires (τexplore+τexploit)( _explore+ _exploit) function evaluations in the absence of early stopping. A naive exact-GP refit after an observation costs (||3)O(|D|^3), although rank-one updates can reduce this cost. In practice we retrain length scales and output scale by maximizing the marginal likelihood after each observation.