Stochastic Subspace Descent Accelerated via Bi-fidelity Line Search
Authors: Nuojin Cheng
Alireza Doostan
Stephen Becker
Venue: Accepted by TMLR ()

Abstract: Efficient optimization remains a fundamental challenge across numerous scientific and engineering domains, particularly when objective function evaluations are computationally expensive and gradient information is inaccessible. While zeroth-order optimization methods address the lack of gradients, their performance often suffers due to the high cost of repeated function queries. This work introduces a bi-fidelity line search scheme tailored for zeroth-order optimization. Our method constructs a temporary surrogate model by strategically combining inexpensive low-fidelity (LF) evaluations with accurate high-fidelity (HF) evaluations of the objective function. This surrogate enables an efficient backtracking line search for step size selection, significantly reducing the required number of costly HF queries. We provide theoretical convergence guarantees for this scheme under standard assumptions. Furthermore, we integrate this bi-fidelity strategy into the stochastic subspace descent framework, proposing the bi-fidelity stochastic subspace descent (BF-SSD) algorithm. A comprehensive empirical evaluation of BF-SSD is conducted across four distinct problems: a synthetic optimization benchmark, dual-form kernel ridge regression, black-box adversarial attacks on machine learning models, and transformer-based black-box language model fine-tuning. The numerical results consistently demonstrate that BF-SSD achieves superior optimization performance, particularly in terms of solution quality obtained per HF function evaluation, when compared against relevant baseline methods. This study highlights the efficacy of integrating bi-fidelity strategies within zeroth-order optimization frameworks, positioning BF-SSD as a promising and computationally efficient approach for tackling large-scale, high-dimensional optimization problems encountered in various real-world applications.

Published in Transactions on Machine Learning Research (09/2025)
Stochastic Subspace Descent Accelerated via Bi-fidelity Line Search
Nuojin Cheng
nuojin.cheng@eollmddo.ee
Deportment of Applied Mathematics University of Colorado Bouider
Alireza Doostan
srnrzaadoosttan@@ooorrddoo edu
Smead Acrospace Engineering Sciences Department University of Colorado Boulder
Stephen Becker
Ceepheen.bcckrr@oolrree edu
Department of Applied Mathematics University of Colorado Boulder
Reviewed on OpenReview: https: // openreview. net/ net/feruu forum? =u0QUs5mmm
Abstract
Efficient optimization remains a fundamental challenge across numerous scientific and en- gincering domains, particularly when objective function evaluations are computationally expensive and gradient information inaccessible. While zeroth-order optimization meth- ods address the lack of gradients, their performance often suffers due to the high cost of repeated function queries. This work introduces bi-fidelity line search scheme tailored for zeroth-order optimization. Our method constructs a temporary surrogate model by strategically combining inexpensive low-fidelity (LF) evaluations with accurate high-fidelity (HF) evaluations the objective function. This surrogate enables an efficient backtracking line search for step size selection, significantly reducing the number of costly HF queries required. We provide theoretical convergence guarantees for this scheme under standard sumptions. Furthermore, we integrate this bi-fidelity strategy into the stochastic subspace descent framework, proposing the bi-fidelity stochastic subspace descent (BF-SSD) algo- rithm. A comprehensive empirical evaluation of BF-SSD is conducted across four distinct problems: a synthetic optimization benehmark, dual-form kernel ridge regression, black-box adversarial attacks on machine learning models, and transformer-based black-box language model fine-tuning. The numerical results consistently demonstrate that BF-SSD achieves superior optimization performance, particularly in terms of solution quality obtained per HF function evaluation, when compared against relevant. baseline methods. This study highlights the efficacy of integrating bi-fidelity strategies within zeroth-order optimization Trameworks, positioning BF-SSD as promising and computationally efficient approach for tackling large-scale, hign-dimensional optimization problems encountered in various real- world applications.
1 Introduction
In this work, we interested in the unconstrained optimization problem
x:" E arg min f(x),
(1.1)
øERD
where the objective function f: R is assumed to be L-smooth (i.e., its gradient V L-Lipschitz con- tinuous). Crucially, operate in black-box setting where direct access to the gradient V is unavailable

---

Published in Transactions on Machine Learning Research (09/2025)
computationally infeasible to obtain with low complexity (e.g., via closed-form expressions automatic dif- ferentiation), thereby precluding the direct application of standard first-order or higher-order optimization techniques. Furthermore, we focus on high-dimensional scenarios where D is large (e.g., D 2 100), pos- ing significant challenges related to computational cost and scalability for many traditional derivativo-free optimization methods.
The primary focus of this work is on selecting an appropriate step size (a.k.a learning rate) os 0 for iterative descent schemes of the form:
Th+11 = 2z. - Ug:
where 2g E IR represents estimate of the true gradient (.( or another suitable descent direction at iteration k. Selecting appropriate step size dynamically can significantly improve the convergence performance of the optimization process. This is illustrated in Figure 1, where an example function is optimized using different methods with and without a step size tuming scheme. However, common practices in machine learning often involve either using a fixed step size throughout the optimization or employing predefined adaptive step size scheduling, e-g., Duchi et al. (2011). While convenient implement, these approaches often neglect the intrinsic local geometry and characteristics of the objective function.
In contrast, classical line search methods, including exact line search and backtracking algorithms, typically yield better step sizes by incorporating information from additional objective function evaluations within each iteration. However, the cost of these additional evaluations can render such methods impractical, particularly when the computational budget is limited - a common satiotion in black-box optimization or when evaluating f (or its gradients) expensive.
To address this limitation, we propose a novel approach to tune the step size by leveraging multi-fidelity evaluations of the objective function. Here, multi-fidelity modeling involves utilizing two more levels of objective function representations: the high-fidelity (HF) objective f, which provides accurate but expensive evaluations, and one or more low-fidelity (LF) objectives that serve as cheaper, albeit less accurate, approx- imations of f. We emphasize that this multi-fidelity setup (often bi-fidelity, involving one HF and one LF model) does not necessarily need to originate from the problem's inherent structure. Even for a problem initially defined with single fidelity, a practitioner can potentially construct a multi-fidelity extension (e.g., by creating simplified surrogate models) to accelerate the optimization procedure, particularly through more informed step size selection.
HF Function value
GD-LS
10-1
55D
CD-L5
HF-550
BF-550
5000
10000
15000
20000
25000
30000
Equivalent HF Function Calls
Figure/ Gradient descent (GD), coordinate descent (CD), and stochastic subspace descent (SSD), along with their respective back- tracking line search (LS) variants for step size tuning, as well the proposed Bi-fidelity SSD (BF-SSD), are evaluated on the "worst func- tion in the world" example, detailed in Section 4.1.
2

---

Published in Transactions on Machine Learning Research (09/2025)
For simplicity, we focus the bi-fidelity case, where only two fidelity levels are considered. The HF objective, JHP, is treated as the ground-truth objective function, so we treat JHl and the ʃ from Equation (1.1) synony- mously. We construct simple bi-fidelity surrogates after obtaining the gradient estimation vn. Specifically, given the LF objective ftF IRD + R, the current position 2k. D): an initial step size amax: and a budget ny for HF evaluations at this step, the local 1D surrogate of the HF objective (p((a) = fHF (ax - -00k)) - R is constructed as
(x( pf.F (æa avvs) + (a; a E [0, [max]
(1.2)
Here, P is a scalar, and ( nx)= : RR IR is a piecewise linear function constructed using ng HF evalu- ations. Once the surrogate @k is constructed, the step size is selected using backtracking line search by (approximately) solving
arg min (a) at aimel
Assuming the scalar p is properly chosen so that the difference
d(x) HHP (x) pfLF (x)
is Lipschitz continuous, we show that the convergence of this descent method is guaranteed, and K, O(L/e) iterations are needed to ensure that min V f (22 is e-small. Moreover, when the HF and LF functions are well-aligned, i.e., the Lipschitz constant W of d((x) is small, the required number of HF function evaluations N.- O(W L2/e + DLfe) is not large.
For implementation, we focus on high- nigh-dimensional zeroth-order optimization problems, using the stochastic subspace descent (SSD) method (Kozak et al., 2021) combined with the proposed step size tuning strategy, and call the resulting method bi-fidelity stochastic subspace descent (BF-SSD). BF-SSD demonstrates strong empirical performance across various tasks and holds great potential for future applications.
1.1 Related Work
Line Search for Optimization Line search is a widely used method for determining step sizes in opti- mization algorithms. Line searches can be either exact, meaning that chosen to exactly or almost exactly minimize ful (xx avu), or inexact. Exact line searches are computationally expensive, other than in special cases, they are rarely used in practice. Common inexact line search methods include backtracking line scarch (Nocedal & Wright, 1999), the Polyak step size (Polyak, 1987), spectral methods such as (Barzilai & Borwein, 1988), and learning rate scheduling (Duchi et al., 2011). Among these, backtracking line search is particularly popular due to its simplicity and explainable design, often employing stopping criteria like the Armijo and Wolfe conditions (Nocedal Wright, 1999). However, backtracking line search increases the overall computational costs considerably due the numerous function evaluations required at cach iteration. One way lo mitigate this issue is by constructing surrogate models to guide step size selection. For example, Yue & Meerbergen (2013) and Grundvig (2023) used reduced-order models to approximate the objective function during line scarch, while Mahscreci & Hennig (2017) and Cartis & Scheinberg (2018) employed a probabilistic Gaussian model for step size selection. Paquette & Scheinberg (2020) provided a theoretical analysis of line search in stochastic optimization. A recent work by Nguyen et al. (2025) extends the backtracking line search framework to the stochastic ISTA/FISTA method. These approaches do not account for the multi multi-fidelity structure of objective functions, which the focus of this work.
Perivative-Free and Zeroth-Order Optimization Derivative-free optimization refers to a family of optimization techniques that rely solely on function evaluations, without requiring gradient information, to find the optimum of objective function. This category includes methods such Bayesian optimiza- tion (Shahriari et al., 2015), direct scarch (Kolda et al., 2003), trust region methods (Conn et al., 2000), genetic a algorithms (Srinivas & Palnaik, 1994), and zeroth-order optimization (Liu et al., 2020). Among these, zeroth-order methods stand out for their scalability to high-dimensiona, problems and reliable convergence properties. Following Liu al. (2020), we refer to zeroth-order algorithms as the type of algorithms that approximate gradients using finite difference techniques and subsequently apply strategies similar to first- order methods. These methods have shown great promise in various machine learning applications where
3

---

Published in Transactions on Machine Learning Research (09/2025)
objective functions smooth but lack accessible easy-to-compute derivatives. Recent advances include their use in solving black-box adversarial attacks (Chen et al., 2017; 2023) and fine-tuning large models, such as MeZO, S-MeZOO Sub2O, etc. (Sun et al., 2022a;b; Malladi et al., 2023; Liu et al., 2024; Yu et al., 2024; Zhang et al., 2024), with minimal memory overhead. A recent work (Brilli et al., 2024) discloses the worst-case bound for derivalive-Iree optimization with a line search-style method.
Randomized Zeroth-Order Optimization for High-Dimensional Problems In high-dimensional zeroth-order optimization problems, estimating gradients via finite differences can be computationally pro- hibitive. To address this, randomized algorithms have been proposed to reduce the cost of gradient estima- tion. The simultancous perturbation stochastic approximation (SPSA) (Spall, 1992; 1998) uses Rademacher random vectors for gradient estimation, while Gaussian smoothing methods (Nesterov & Spokoiny, 2017) employ Gaussian random vectors. These algorithms typically provide gradient estimates projected onto one- dimensional subspaces. However, for certain problems, it is worth the increased number of function calls improve gradient estimates. SSD (Kozak al., 2021) explores this idea by projecting the gradient onto a random subspace of dimension & for any f E S D, providing a more generalized framework for randomized zeroth-order optimization.
Modeling and Optimization Multi-fidelity is a well-established concept in engineering and scientific computing for reducing computational costs. It has been widely applied across various domains, including aerodynamic design (Zhang et al., 2021), structural optimization (Ng Willcox, 2014; De et al., 2020), data sampling (Cheng et al., 2024b); and uncertainty quantification (Peherstorfer et al., 2018; Cheng et al., 2024a; De & Doostan, 2022; De et al., 2023; Cheng & Doostan, 2025; Cheng, 2025). As one sub- sub-branch, multi-fidelity optimization has been employed in hyperparameter tuning (Wu ct al., 2020), accelerating Bayesian oplimization ion (Kandasany el al., 2016; Takeno el al., 2020), and reinforcement learning (Cutler et al., 2014) within machine learning. However, despite its relevance in settings where function evaluations are costly, its application in zeroth-order optimization remains largely unexplored (de Montbrun & Gerchinovitz, 2024) and has not been applied to any randomized zeroth-order method.
1.2 Contributions
In this work, we propose a multi-fidelity line-search scheme. Unlike previous approaches that utilize static surrogate (e.g., reduced-order) models, which remain fixed throughout the optimization process (Yue & Mecrbergen, 2013; Grundvig, 2023), our method constructs a temporary surrogate model in cach iteration, specifically after the gradient (or search direction) has been estimated. This allows us to focus on building a one- one-dimensional surrogate along the search direction, a significantly simpler task compared to constructing expensive D-dimensional surrogates for the objective function f. By leveraging computationally cheaper LF model, we construct simple, local, linear surrogate using only a small number, n4. of HF evaluations per iteration. Assuming certain conditions between the LF and HF models hold, this 1D surrogate is then used to efficiently identify a suitable step size.
Specifically, this work makes the following contributions:
We develop the BF-SSD algorithm, a stochastic zeroth-order optimization method with a bi-fidelity line search that allows for choosing the approximation quality of the gradient by tuning & (reducing to deterministic gradient descent when é D).
When the error of the gradient estimate negligible (e.g., is sufficiently large), we give specific conditions on the relation between the HF and LF functions that will guarantee convergence to a stationary point (or a global minimizer when f is convex).
We highlight that many machine learning problems naturally have a corresponding LF model that can be used to construct surrogate model, improving optimization efficiency. Despite its high potential, this strategy has not received sufficient attention in prior work.
4. We compare BF-SSD with other zeroth-order optimization methods on one synthetic function and the following three real-world applications:
4

---

Published in Transactions on Machine Learning Research (09/2025)
Kernel ridge regression with byssedm-based LF approximation;
Black-box image-based adversarial attacks with an LF model trained via knowledge distillation;
Soft prompting of language models using smaller training set construct the bi-fidelity line search.
The rest of the paper is organized as follows. Scction introduces the proposed bi-fidelity line scarch method and provides convergence results. Section 3 delails the implementation of the proposed method with SSD. Section 4 presents the experimental results, and Section concludes the paper.
2 Line Search on Bi-fidelity Surrogate
In this section, we discuss the proposed algorithm and present the main theoretical results derived in this work. Unless specified otherwise, ||-I| denotes the Fuclidean norm for vectors and the spectral norm (i.e., the induced 2-norm) for matrices. For simplicity and to maintain foeus on our primary contributions during the theoretical analysis, assume in the proofs that v& provides an accurate estimate of the high-fidelity gradient, specifically U2 f This assumption is employed solely for the theoretical development presented herein and does not hold for the practical implementation discussed in Section 3 or the numerical results presented in Section 4.
2.1 Algorithm
First, we define the algorithm, which consists of three steps for each iteration k:
Given the current position æ& E RD, gradient UR RD, and initial step size Qmax E E, sample ny equi-spaced HF ovaluations in |0, amax and build the surrogate Øx R IR following Equation (1.2) (see Algo. 1 for details);
Given Armijo condition parameters c E (0, 8 1/2, and initial step size amax 2 c/(L cL), conduct bi-fidelity adjusted Armijo backtracking that
ag max cd Amax mEN
(2.1)
S.L. (2(c" Ommxx nx) < H (axx) 3c Om |0x/1
See Algo. for details.
3. Evaluate PHF at the new point and continue the iterations.
2.2 Convergence Results
For convergence, we make the following assumptions:
Assumption 2.1. The objective function fHF, RD - R attains its minimum f and f is L-Lipschitz continuous; i.e., there exists L E R such that
NI V
II (1))
All vll, Va, y E R.
Note that Assumption 2.1 is standard for analysis zeroth- and first-order methods. The constant L must be known to the algorithm since is used to set. Omax-
Assumption 2.2. The difference between fHT and f&. is assumed to be smooth with bounded Lipschitz constant. Specifically, we assume there exists W, p E R such that
H (x: OLLF ()))(())) (JHF (v) pLLF (v)) W| vll. Va, y R..
5

---

Published in Transactions on Machine Learning Research (09/2025)
The Assumption 2.2 allows for fl. to be uncalibrated, meaning that we do not require fl. (x) (*) since discrepancies can be reduced by building the surrogate PP..
Our next assumption, Assumption 2.3, is a sufficient condition that will be used in Lemma 2.7 to show that the surrogate accurate.
Assumption 2.3. For each iteration k, we assume the number of HF evaluations, nk, for building the surrogate Φ is sufficiently large such that
WL(1 c)0max
WT.
nkk
cB11ux111
i.c., ny 5
||2 12
where 22+) denotes lower bound up to a constant difference.
Using 2s fxxx) and with the above assumptions satisfied and sufficiently large initial step size Omax c/(cL L), the designed bi-fidelity line search leads the following result:
Thcorem 2.4. Given on initial point To, assuminn actual gradients are accurately estimated, the algorithm
in Section 2.1 generales sequence (xx) such that
| (xx))1 2L(1 + ( (20o) - -*)
min KO
ke
(K + 1)e3
That is to say, K, = OLLe)) ilerations are required to oblain minkg V H 13
Remark 2.5. Theorem 2.4 holds when = Vf(=x). The error in approximating Vf(=x) using finite difference methods with O(D) samples typically negligible in comparison to the optimization error (see Kozak et al. (2023) for a procise quantitative statement for the case of SSD). Hence, assuming we accurately estimate er = (a(a)) with O(D) samples per step, a bound for the total number of HF evaluations for é-convergence of the algorithm in Section 2.1 is
W 2²
DL
N.
(na + +OD)) 0
+
(2.2)
E
Remark 2.6. When using zeroth- gradient descent, with the same assumption that we accurately estimate 2% fa(z) with O(D) samples per step, bound for the total number of HF evaluations for t-convergence is
N. = Σ (log.- (Amax O(D)) = [LlogLL)
K.
+ (
DL
(2.3)
The proof of Remark 2.6 follows the convergence proof of gradient descent using backtracking line search. Comparing the results in Equations (2.2) and (2.3), we observe that the advantage of using our bi-fidelity surrogate depends on the value of W. Notice that if W is sufficiently small that WL S clog L, then the worst-case bound of our method is better than that of the zeroth-order gradient descent. We emphasize that our convergence result in Equation (2.2) is loose, due to the global nature of the Assumption 2.2 and difficulty in precisely describing the quality of the LF function relative to its HF counterpart. Hence, we view our convergence analysis as a reassurance that the method does converge, and rely on mumerical experiments elucidate when the method improves over baseline methods.
2.3 Proof of Theorem 2.4
Before the proof, we first introduce the following lemma:
Lemma 2.7. With Assumption 2.2 and Assumption 2.3 satisfied, for any a E [0, 0mxx), the 1D surrogate (a(a)) satisfies the following bound,
|
cB
cB c
(a((;; - s(aa) 5
min
(1+ L 80
(2.4)
2
(1+c)
2(1 cc

---

Published in Transactions on Machine Learning Research (09/2025)
The proof of Lemma 2.7 is in Appendix A. Following this lemma, the proof of Theorem 2.4 is given below.
Proof. Given Lemmn 2.7, we have
v [ S { au Soll/oxl1 112 + abL |
cB
(x2+1) 5x (ax;m))) = (p((n)) - (x(()))a) S
min c1+ (1 c)
30max
2
and, using the standard descent lemma for L-smooth functions (guaranteed by Assumption 2.1),
Therefore, using the triangle inequality, the surrogate 6b is bounded as
) & fHF (IK+1) 1 +
+ (xx+1) -
s | (mx) + (aa++ + 0f
2(1 c) ||0x11?
When the step size satisfies og E [cc(LL 1/(L+ ([[++co) cb)),, the quadratic inequality -ox + L/2 e/ L) v -ax/2 holds, along with the fact that S 1/2, which implies the following bi-fidelity-adjusted Armijo condition
agL
Px (axinR)
H -okk +
2(1 + cj2T. ||0x
H (xx) {H))
||2x*1 11
(2.5)
f (#x) - -Ba 11
The last line in Equation (2.5) satisfics the bi-fideltt-ddjuuseed Armijo condition in Equation (2.1). Therefore, the bi-fidelity backtracking either terminates immediately with aj O'ma or else Or 2 c/(L cl), and implies
11 min {
Bc
6 (ag;me)= S H (axx) 0||00
(1+ Omax
,HF (Ir*) (HF((*))
(1 c) ||
(2.6)
where the last equality comes from Qmax 2 c)) Using c) (from Lemma 2.7) and Equation (2.6), have
Pp (ax: nx)) S 1cllex/11/(21)
|HP (xx+1) I 6a (ax;mE) + p (xx++)1)) 56 (ous
+ Bcl|xx
I Px (axi
(2.7)
2(1 + c)
Bel|voll Mu
HP (mx)
2(1 + c)
Equation (2.7) leads to the telescopic series
Bc
K
K
< Σ (fHl (xx) - f ((x++))
2(1
L
==
k=0
HHF (xx) - -HHP (#0 I JHF (xo0 f"
Hence,
Be
(K + 1)
min
(f11)( -
(0.......K
2(1 + c)T.)
2(1 + C) L)L (fHIf -f*).
Bc

---

Published in Transactions on Machine Learning Research (09/2025)
To guarantee m | f (#x) S e, the value of Ke should be
K2 2 2( (zo) (f + r)ه O (
Bee
Remark 2.8. Even if JHF is non-convex, Equation (2.7) implies that the method is descent method, meaning fH (I++1) I (xs). Hence, after K iterations, it is natural to use ææ as the output. This descent properly is not enjoyed by other methods, such as subgradient descent, stochastic gradient descent, or Polyak step size gradient descent.
Remark 2.9. If fHP is convex, then Theorem 2.4 implies convergence to global minimizer. Or, if fH satisfies the Polyak-Tojasiewiey inequality with parameter l (which includes some non-convex functions, well as all strongly convex functions), then Theorem 2.4 in conjunction with the descent property implies {HF (æx) f" S (K+1)ue3 (JHF (20) - f"), cf. Karimi et al. (2016).
L(l+e)
2.4 Examples of Possible Low-Fidelity Functions
In practice, the LF function fLP can be constructed various ways. The most straightforward approach is when a multi-fidelity structure is intrinsically present in the problem. For example, in Cheng et al. (2024a, Section 5.1), the LF model is the exact solution to simplified physical model that can be simulated with negligible cost, while the intended HF objective relies on relatively expensive finite element simulations. In most machine learning problems, the LF model is not explicitly given, making its construction necessary. In this section, we discuss multiple approaches for building the LF model and the associated upper bound on W.
Affine Bi-Fidelity Relationship The ideal case occurs when the HF model is an affine transformation of the LE model, i.e,, f (a) (x) c. In this case, the Lipschitz constant W, as defined in As- sumption 2.2, is zero, and the number of function evaluations required for convergence is proportional to the number of iterations, as n2 = 1 is sufficient.
Quadratic Objective with Raw-Rank LF Approximation Consider the case where the objective is quadratic with positive semi-definite matrix A E DDD and denote its rank-r approximation A E DxD and assume that rank(A) >> rank(A). The HF objective is f (z) = (xx Ax) + (ax a) for E X' where ) denotes the Euclidean inner-product, and the LF objective fLF (a) = (x, Ax) + (ax a). Assuming the input space X is bounded by a unit ball with radius R, the Lipschitz constant W is upper bounded as
W < V (x) - (Vf"((x)) = supl| A - ||||||=|| < Art1R,
where 2y41 is the (r 1)-th largest eigenvalue of A. The empirical problems in Section 4.1 and Section 4.2.1 fall into this category.
This use-case satisfies the assumptions mentioned in Section 2.2. Assumption 2.1 is satisfied since the minimum is achieved (fuFf is continuous and coercive) and the Lipschitz constant of the gradient is L = ||]]] Assumption 2.2 is satisfied using the value of W above, and Assumption 2.3 can be satisfied since it just parameler choice.
Full-Batch HF and Mini-Batch LF Objectives In many machine learning settings, the objective function is expressed as a sum over a large number of terms, each corresponding to the evaluation of a loss function on an individual data sample. In this case, a natural choice for the LF objective is the summation over a smaller subset of the data. Specifically, assuming that the HF objective sums over datapoints i 1,...,n and (without loss of generality, i.e., by relabeling) the LF objective sums over datapoints i = l....,r r for & n, the HF objective is f(1((z) (x) Enm f((æ), and the LF objective is fLF (x) = 1 f((x). Using the triangle inequality, the Lipschitz constant W is upper bounded as
W& < sup
Σ
V (x))
n n-r max <
[V+(=)] + max [V((()))
n i=1
13i5n
8

---

Published in Transactions on Machine Learning Research (09/2025)
The terms [[)]] are bounded if each f. is Lipschitz, or equivalently, if f. is continuous and x: is constrained to a compact set. An empirical problem with this setting is presented in Section 4.2.3. Our analysis is deterministic, so W is a worst-case bound, but if r is large and the LF subsamples are chosen uniformly at random, it would be reasonable to expect that, due to the law of large numbers, the average case behavior is significantly better than our bound.
This use-case also satisfies the assumptions mentioned in Section 2.2 under reasonable conditions. If each ff, is Lipschitz continuous with constant L; (i.e., this is always true if fi is continuous and is constrained to compact set), then ff is L Lipschitz continuous with L = 1 - L; via the triangle inequality, so under the mild assumption that the minimum is achieved, Assumption 2.1 is satisfied. Furthermore, Assumption 2.2 is satisfied using the value of W above, and Assumption 2.3 can again be automatically satisfied since it is just a parameter choice.
Generic Case Finally, we consider the most general case, without assuming specific relationships between the HF and LF objectives. By assuming the Lipschitz contimuity of both the HF and LF objectives, W can be bounded as
W = 1I (x) - pfLP (x)Il< S Il (z) - lpl (z))ll))
for any choice of p, where ||-Ilz denotes the Lipschitz constant. The proportionality should not be chosen to minimize this bound (since that leads to 0) but can instead be chosen by any heuristic, such as the one used in control variate techniques (Gorodetsky et al., 2020) where p p= -ê/0 where è an estimate of the covariance between jHP and fL.F, and è is an estimate of the variance fHP.
3 Bi-Fidelity Line Search with Stochastic Subspace Descent
In this manuscript, focus zeroth-order optimization, utilizing stochastic subspace descent (SSD) the implementation method. Following the algorithmic steps introduced in Section 2.1, combined with SSD, the entire process is divided into three main (iteratively implemented) components: gradient estimation to construct D2, bi-fidelity surrogate construction, and Armijo backtracking on the surrogate.
Gradient Estimation SSD employs a random projection matrix PA E Dxt with E D. The random matrix Ps satisfies the properties P In and P Pa (D/E)Tz A common choice for Pe is based the Haar measure, where P& is derived from the Gram-Schmit orthogonalization of random Gaussian matrix. The gradient estination is given by ta- = Pagg, where Os is the finile difference estimator of the gradient:
HF (F APi) - HF (ark) H (ætg Ap2) H (DL)
HHF (00k + Ap:)
HHF (xrs)
gk
(3.1)
A
A
where A R is small step size and p. is the i-th column of Pa. Estimating Uy using Equation (3.1) requires
l function evaluations - a more accurate ((A?) approximation is also possible at the cost of 2 function evaluations if more than 8 digits of precision are needed. Up to the finite-difference error, 9k ~ P V P (æx) so that U2 ~ P.P. V PHF f (æ*), hence E(vy] zs f (xx). It is also possible to construct the same estimator without reference to the Haar measure by rewriting 22 as o% = projool(q. (V)H(( where Q2 E Dxt is any random matrix with indopendent columns from an isotropic probability distribution (such as the standard normal).
Surrogate Construction Given the estimated gradient us and the current position æR, the goal of sur- rogate construction is to build 6x, denoted as
(3.2) The analysis in Scction 2 assumes that p is known and fixed. However, in practice, P tunod for better performance. In our application, we set Pk: ,HF (ax))/ LLF (œx). We model th as a piecewise linear function using nı additional HF evaluations at equispaced points {0, "o, an. amax]. Specifically,
(a) = pfLF (xx + avv) + (a).
PLFF
a (a(aa = h "((āg-1)+ v(āz), a aālāy-r,ūj) j j= j=1..... -.....,24 R.k
(3.3)
h
9

---

Published in Transactions on Machine Learning Research (09/2025)
where h Omax//ng and ((a) = p(a) - -f. fxx
avk). -av))) This piecewise linear interpolation is simple, yet effective, approach for interpolating şp in 1D and satisfies the bounds in Lemma 2.7 given sufficient The detailed algorithm is presented in Algorithm 1.
Figure illustrates the bi-fidelity backtracking line search process using the example problem in Section 4.2.1. The blue curve represents the bi-fidelity surrogate model f approximating the HF function jp (red curve). Rather than performing the line search directly on the computationally expensive HF function (red dots), the method utilizes the surrogate Ga to estimate optimal step size. While this surrogate is approximation and may require more surrogate function evaluations during the search itself, it substantially reduces the computational cost of line search. In this example, the cost is decreased from 4 HF function calls (for direct search) to only HP call (to build the surrogate) combined with 6 LF function calls.
Algorithm 1: Surrogate Construction
Input: fLV, fHF. Lk, Uk n& N, O'max Output: 1D surrogate 6
1: Define {(āj, ([(ū)))] as equispaced points between 0 and O'max (including endpoints), and compute
HF ovaluations ()ā;) +- jHP (xx + āzbk);
Ph f 1HF (xx)/ ( (RR);
3: v(āz) + p(ā,) Paf (xx + a,Uk): j= L...., Re:
4: Construct piecewise linear function the using Equation (3.3);
5: Return 6 using Equation (3.2).
HF gv(a)
5600
dula)
Armijo Condition
5400
Objective Value
5200
5000
4800
0.4
0.6
0.8
1.0
Step Size
Figure 2: Illustration of the bi-fidelity backtracking line search pro- cess using the example problem in Section 4.2.1. The blue curve represents the bi-fidelity surrogate model &. approximating the HF function (red curve). It significantly lowers computational cost (e.g,, reducing 4 HF calls to HF + LF calls).
Armijo Backtracking on the Surrogate Based the criteria in Equation (2.1), we set the maximum number of iterations for testing the Armijo condition to M N. The detailed procedure is presented in Algorithm 2.
Convergence Analysis of SSD with Line Scarch The convergence results of SSD with line search (on the exact ((a)) are presented in Appendix B, under three separate scenarios: strongly convex, convex, and non-convex. The proof shows that, in the SSD with line search setting, the value of 9 can be set L/2D.
10

---

Published in Transactions on Machine Learning Research (09/2025)
Algorithm 2: BF-Backtracking
Input: Ek, 8 > O (0, 1), amax > 0, UR. M EN N
// typical value of c 0.9
Output: Stop size Ock
1: Initialize a₁ f- Qmaxi
for m =0 M do
if (x(as) - pHl (xx) - a1Blloul1 then
Break;
5: else
6:
f- coo:
7: end if
8: end for
Return ox:
The proposed bi-fidelity line search algorithm, combined with SSD, will be referred to bi-fidelity SSD (BF-SSD), and is summarized in Algorithm 3. Our theory covers either / = D with bi-fidelity line search (Thm. 2.4) or 1 < S D with HF line search (Appendix B) under strongly convex, convex non-convex settings. Combining the two analyses is fairly complicated, and we defer it to future study. For reference, the bi-fidelity result is in Theorem 2.4 and below is one of the results from Appendix B.3 for the 1 1 D with HF line search case, particular, the result covering the non-convex case:
Theorem 3.1 (proof in Appendix B.3). With Assumption B.5 holding and backtracking implemented for line search, we have
(J(z)) DL(S(m))-(')) B[][vf(==]111+ max { (KK 1)30max
min
(K + 1)Ec8 }
kE (O....,))
That is, k = O(L/((((maa)) + (L/(el))) iterations are required to achieve EE [(xx)) S.
For practical purposes, since the parameters in Assumption 2.3 are often unknown, we set Pk as described above and choose n, 1 to minimize the cost of generating (a)). This choice also demonstrates excellent empirical performance across all HF and LF pairs have examined. We leave it as future work for adaptively choosing ng, but suggest one possible scheme inspired by similar methods used in trust-region algorithms: starting with a small n4, we find the stepsize via line search on the surrogale, as usual, and then compare the actual function value of the high-fidelity function with the predicted value from the surrogate. the predicted decrease from the surrogate was overly optimistic, we then rebuild the surrogate with a larger ni (re-using existing samples) and repeat.
Algorithm 3: Bi- Bi-Fideeity Line Search SSD Algorithm
Input: JHF fLF,C,c, e, M, O Output: HF minimum value
// by default, 1 and B f/2D
1: Initialize No and set of HF values D = Hf (xn)}
2: B E/2d;
for == K do
Sample random matrix Ps; Appraximate
P.P. V (xxx Normalize Uj uz/ /
using finite difference (& HF evaluations):
7: Construct G + suuroogate-consrrutiinn FLFF HHF 2k, Uk ) (n HF evaluations);
8: O% + BF-backtracking(77 8, Omax, Uk M);
9: Update Øk+1 + Øx 04
10: Evaluate FHF (Ix+1) and update D:
11: end for
12: Return min D;
11

---

Published in Transactions on Machine Learning Research (09/2025)
4 Empirical Experiments
In this section, we evaluate the proposed BF-SSD Algorithm 3 on four distinct: one synthetic optimization problem discussed in Section 4.1 and three machine learning-related problems across diverse scenarios pre- sented in Section 4.2. These include dual-form kernel ridge regression (Section 4.2.1), black-box adversarial attacks (Section 4.2.2), and transformer- based black-box language model fine-tuning (soft prompting) Section 4.2.3. We demonstrate that the BF-SSD algorithm consistently outperforms competing methods. To illustrate these advantages, we compare BF-SSD against the following baseline algorithms:
Gradient descent (GD): A zeroth-order gradient descent method, where the ball-batch gradient estimated using forward differences, and a fixed step size is used.
Nesterov accelerated gradient descent (NAG): An accelerated zeroth-order gradient descent method by Nesterov that incorporates a decaying momentum term, which often leads to faster convergence. The gradient is estimated using forward differences.
Coordinate descent (CD): Iteratively and individually optimizes each coordinate using finite- difference estimated coordinate gradients.
Stochastic subspace descent with fixed step size (FS-SSD): The standard stochastic subspace descent method, which samples subspaces from the Haar measure and uses a fixed step size.
Simultaneous perturbation stochastic approximation (SPSA): A randomized optimization method using& a Hadamard random variable to estimate the gradient, as proposed by Spall (1992) and with step sizes as described in Spall (1998); this is a time-tested, well-established zeroth-order method.
Gaussian smoothing (GS): A method popularized by Nesterov & Spokoiny (2017), which is nearly equivalent to SSD with L ==1, and uses a fixed step size.
High-Fidelity stochastic subspace descent (HF-SSD): A single-iidelity SSD method that utilizes high-fidelity function for backtracking line search, with its convergence analysis detailed in Appendix B.
Variance-reduced stochastic subspace descent (VR-SSD): A varianco-reduced version of the SSD method inspired by SVRG (Johnson & Zhang, 2013), as described in the technical report Kozak et al. (2019, Section 2.2) of the SSD authors.
Bi-fidelity stochastic subspace descent (BF-SSD): The proposed method detailed in Section 3.
The performance of the optimizers is assessed based on the number of HF objective function evaluations required, accounting for LF calls (in terms of fractional equivalent HF function calls) as appropriate. While we also measured wall-clock time, the results strongly aligned with the equivalent HF call data. Therefore, for clarity, we present performance primarily in terms of equivalent HF calls.
4.1 Synthetic Problem: Worst Function in the World
In this section, investigate the performance of our proposed BF-SSD algorithm the "worst function in the world" (Nesterov, 2013). With fixed Lipschitz constant L > 0, the function is
[(;;,, (
xr - - ) x I1
Lr
8
4
8(r + +1
where I; denotes the tth entry of the input æ and r < D a constant integer defining the intrinsic dimension of the problem. The function convex with the global minimum value 0. The Lipschitz constant of the gradient of this function is L. Nesterov (2013) has shown that wide range of iterative first-order methods perform poorly when minimizing f(x L.) with initial point æt = 0.
12

---

Published in Transactions on Machine Learning Research (09/2025)
Worst Function the World (ell=20)
Worst Function in the World (ell=50)
15,00
5754
5754
Best Function
10-500
Best Function
01-350
V1.53
V1555
10m
10
19800
1966e
20ien
1et0
Isien
20000
Ireco
Equivalenn HF Function Calls
Equivalinn HF Function Calls
(a)
Worst Function World (ell=100)
Worst Function World (ell=200)
Best Tunction
Punction
60
69
03
10
V0.
5000
13686
15600
20000
25300
30000
5000
18660
15660
20000
25000
30000
Equivalent HF Function Calls
Equivalent HF Function Calls
(c)
(d)
Figure 3: Convergence performance for different optimizers. The x- axis is the equivalent nunber of HF function evaluations, including the number of LF function evaluations based on the cost ratio FI/TH. The y-axis is the HF function evaluation value at the current stage. We investigate the results when l = 20, 50, 100, 200 with 9y = 2, ти = 100. The corresponding results are presented with their ti- tles indicating the specific choices. The shadow regions are the area between the best and the worst behavior over 10 trials.
We set the dimension D = 1,000, f= 20, and = 20. The intrinsic dimensions of the LF and HF functions are ти and rz, respectively. We choose < ти and assume the computational cost ratio between HF and LF evaluations is rи rz- For Gradient Descent, we choose the standard step size of 1/L = 0.05, and for the GS and SSD-based methods, the step size is t/(T.)). The backtracking parameter hyperparameter study conducted according different values of E (0.8, 0.9, 0.99) and f E (5, 10, 20). All the experiments repeated 10 times, with shaded regions denoting the worst and the best performance over 10 trials.
= f/(2D). The
Figure 3 illustrates the performance of various optimizers across different values of t. Detailed results for = 20 and 0.99 at N from 500 to 8,000 are presented in Table 1, while additional comparisons across different and configurations are included outperforms the other optimizers in most scenarios. For different. SSD methods, the effect of & on the final performance varies. Large values of improve the optimization results for FS-SSD and VR-SSD, while HF-SSD and BF-SSD prefer relatively smaller f, highlighted in Table 2.
Table 4. These results show that BF-SSD consistently
13

---

Published in Transactions on Machine Learning Research (09/2025)
Equivalent HF function evaluations N
Method
N 100
N 1000
N 10000
N 20000
N 30000
GD
2.48 + 0.00
2.48 = 0.00
0.62 = 0.00
0.43 = 0.00
0.34 + 0.00
CD
1.48 + 0.00
1.48 0.00
0.70 0.00
0.49 + 0.00
0.40 0.00
NAG
2.48 0.00
2.48
0.33
0.16 = 0.00
0.11
FS-SSD
2
2.45
2.20
2.08 = 0.00
1.93
SPSA
1.93
0.00
0.00
0.15 = 0.00
0.12
GS
+.00
2.00
2.36 + 0.00
2.25 = 0.00
2.15
HF-SSD
1.99 0.05
0.01
0.40 = 0.21
0.37 0.21
0.20
BF-SSD
2.00 + 0.08
0.02
0.17 = 0.00
0.11 0.00
0.09
VR-SSD
2.47 + 0.00
2.09 0.01
0.62 0.00
0.43 + 0.00
0.35 0.00
Table Performance values (mean std over 10 runs) showing the objective function for different optimization methods at various HF function evaluations N with & = 20 and 0.9. The minimum values in each column are highlighted in bold.
Table 2: Comparison of SSD methods for different values of & (Mean + Std at N = 20, 000). Bold values indicate the minimum mean for each SSD method, i.e., across each row.
Method
{ = 20
{ = 50
{ = 100
L = 200
FS-SSD
2.0787 0.0013
1.6621 0.0003
1.2942 + 0.0022
0.9473 0.0027
HF-SSD
0.3696 + 0.2109
0.1482 0.0098
0.1574 + 0.0157
0.3696 = 0.2368
BF-SSD
0.1143 + 0.0005
0.1226 + 0.0015
0.1260 + 0.0032
0.1298 + 0.0013
VR-SSD
0.4309 + 0.0025
0.4281 + 0.0013
0.4238 + 0.0009
0.4226 = 0.0013
4.2
Zeroth-Ordar Optimization for Machine Learning Problems
Next, we contrast the BF-SSD optimization results against other competing zeroth-order methods. Besides showing the advantages of the BF-SSD, we also show that it is often convenient design cheap LF model in many machine learning problems that can be leveraged to accelerate the convergence.
4.2.1 Dual Form of Kernel Ridge Regression
Consider a kernel ridge regression problem as follows. By the representer theorem, given data {(0,, W:) and a kernel function x R x RRth + R, the goal is to find the coefficients â such that
D Spredict (x) - a.))
(4.1)
One way to compute the coefficients is to solve the dual form of the kernel ridge regression,
à" arg min at Ka 2(a, y) + Al|ae|??
(4.2)
where K is the kernel matrix with K k(æ,, x,), [y); Vt and À is positive scalar regularization parameter. The solution of Equation (4.2) can be explicitly represented as
(K -
However, solving the explicit solution involves inverting the matrix K + AI, which takes O(D3) and can be extremely expensive when D is large. When D is sufficiently large, evaluating the function in Equation (4.1)
14

---

Published in Transactions on Machine Learning Research (09/2025)
takes O(D²) and becomes expensive. Therefore, alternative approach to solve this problem is to build low-rank approximation (surrogate) of the kernel matrix K, which we adopt for evaluating the LF objective. To do this, we employ the Nyströem method, which finds a subset S C [I...., with size / < D and builds the kernel surrogate K :| S (() Kis, :]. By implementing the Nyströem method, the complexity evaluating the objective function reduced to O(D).. Therefore, the ratio computational cost between HF and LFF function evaluation is D/l.
3.0
2.5
Egenvalues
0.5
0.0
200
600
1000
Eigenvalue Order
Figure 4: The eigenvalues of the kernel matrix implemented in Equa- tion (4.2).
We pretend the problem in a black-box format, where access to the HF function is only available through an API, thus hiding y (and/or K) in Equation (4.2) and making derivative information inaccessible. In this case, assume the values of y are unavailable for privacy reasons. For the regression data, we select the first D 1,000 samples from the California housing dataset provided in the scikit-learn library (Pedregosa et al., 2011). A subset S with size = 10 is randomll selected for the Nyströem approximation. We use a Gaussian (RBF) kernel with hscale 1.0 to generate the corresponding kernel matrix K. Figure 4 shows the decay of eigenvalues for K, with a rapid drop, especially within the first 100 eigenvalues, due to the Gaussian kernel's propertics. This fast decay motivates our focus on cases where the values are below 100. The starting point æn is set at the origin 0.
Black-box Kernel Regression (c=0.9)
Black box Kernel Regression (c=0.99)
Best MF Functionvall
20
200
60
GB 85
HF
FS-550
$7.4
10-
V9.500
20000
47528
5000O
1e06
136500
23000
40000
estde
89933
10905
Equivalent HF Function Calls
Equivalent HF Function Calls
(a)
(b)
Figure 5: Similar to Figure compare the optimizer performances with varying parameters e 100 and E {0.9, 0.99]. The corre- sponding results are presented with their titles indicating the specific choices. The shadow regions are the area between the best and the worst behavior over 10 trials.
The results of kernel ridge regression are shown in Figure 5, with values chosen from 0.9 and 0.99, and the value of is fixed as 100. According these results, BF-SSD shows advantages over other methods except NAG and HF-SSD in Figure 5a. When the backtracking factor decreases, the step sizes determined by the
15

---

Published in Transactions on Machine Learning Research (09/2025)
backtracking method become more conservative, leading to suboptimal results, especially for BF-SSD. We also implement different combinations of and & and collect the SSD performances in Table 3. The results suggest that BF-SSD outperforms other SSD methods for most cases, and larger values of & and c improve the performance of BF-SSD.
FS-SSD
HF-SSD
VR-SSD
BF-SSD
10
3500.91 + 5.74
9.49 + 0.21
23.31 + 0.83
8.18 + 0.60
0.9
50
1015.82 + 5.90
4.01 + 0.13
20.87 + 0.55
6.13 + 0.22
100
270.95 + 7.07
5.49 + 0.18
21.56 = 0.69
6.22 + 0.25
10
3493.88 + 11.27
18.39 0.28
23.20 0.30
2.94 + 0.49
0.95
50
1013.22 8.76
5.33 0.35
21.44 1.37
2.22 0.21
100
270.36 + 2.62
6.15 0.17
21.02 0.98
2.28 0.10
10
3503.99 + 3.77
29.55 + 0.05
23.35 = 0.38
1.40 + 0.29
0.99
50
1010.92 + 7.17
6.95 + 0.19
20.90 = 0.31
0.75 + 0.12
100
273.64 3.41
6.26 + 0.08
21.38 = 0.46 0.79 + 0.24
Table 3: Black-box kernel ridge regression HF function values (mean + std) for FS-SSD, HF-SSD, VR-SSD, and BF-SSD at various com- binations of l and at N 50, 000. Considering uncertainties, the minimum values in each row are highlighted in bold.
4.2.2 Black-box Adversarial Attack on MNIST
In practice, especially in explainable AI (XAI), researchers have found that many deep lcarning models are not robust to dala noise. Specifically, if test data is contaminated by a small perturbation imperceptible to humans, many previously well-performing deep learning models fail to produce reasonable results (Goodfellow et al., 2014). Generating such biased noise to confuse a trained neural network model is usually referred to as "attack" in adversarially robust training. This need not be "black hat" activity, as it can be used as part of hardening a system to prevent these attacks in the future.
There are primarily two types attacks: one a white box attack, in which we have knowledge of the model and inject the adversarial noise to confuse the given model. The standard approach under this scenario is to generate the shift. in pixel space based on the gradient of the objective function to maximize the loss. The other type of attack is called black-box attack, in which one does not have knowledge of the trained model and would like to generate adversarial data from it. The black-box scenario is more difficult due to the missing knowledge, and one way to solve it is to treat this problem as black box optimization. To generate an adversarial sample for given data instance æt E RD, with D reprosenting the number of pixels in the given image, common formulation of the adversarial attack is to find noise sample x' solving
T²
= argmax L (g(ar + +x ))) S
(4.3)
where yf is the correct label of æt, L is the attack loss function, and g represents the model for attack. Following the adversarial attack paradigm of Carlini & Wagner (2017) and its black-box extension (Chen et al., 2017), we use soft version of the problem (4.3) as follows:
arg min CE (g(x + at =l),"') + 7||0||?
where the cross entropy loss CE is assigned as the attack loss and g(-) outputs the probabilities of different classes, usually using a softmax function for normalization. is variable balancing the attack loss function
16

---

Published in Transactions on Machine Learning Research (09/2025)
CE and the attack norm. The optimization goal is to find small shift in pixel space that the output results are greatly changed.
In this study, we utilized two convolutional neural network (CNN) architectures to model the HF and LF representations for classification tasks on the MNIST dataset with 60,000 training data and 10,000 testing data. The HF model was deeper CNN consisting of two convolutional layers, the first with 32 filters and the second with 64 filters, both using 5 x 5 kernels, followed by ReLU activations and 2 max-pooling, The flattened output from the convolutional layers (7 x 7 x 64) was connected to fully connected layer with 1024 neurons, followed by 10-class output layer. In contrast, the LF model employed simplified architecture with single convolutional layer containing 2 filters and 3 x 3 kernel, followed by ReLU activation and 2 x 2 max-pooling. The output (13 x 13 x 2) was flattened and passed through a fully connected layer with 16 neurons, leading to a 10-class output layer with log-softmax activation. The HF model was designed to provide high-capacity representations, while the LF model served as lightweight alternative for computational efficiency. The LF model was trained using knowledge distillation (Hinton et al., 2014), leveraging only 1000 training samples and 1000 evaluations of the HF function. Knowledge distillation is technique where smaller, simpler model (the student) learns to replicate the outputs of a larger, more complex model (the teacher), effectively transferring knowledge while reducing computational costs. The classification accuracy for the HF and LF CNNs are 99.02% and 82 82.21%, respectively. There are 27,562 parameters for the LF CNN and 3,274,634 parameters for the HF CNN. We estimate the ratio between HF and LF computational costs as 3274634/27562 ss 118.8. The images in MNIST dataset are 28 x 28 with a single channel, hence the dimension is D 784. The starting points are initialized as the origin point for all experiments.
Black-box Adversarial Attack (Test Case 1)
Black- box Adversarial Attack (Test Case 2)
Best Function
Best Function
23-1
10-
6AD
640
RF-550
HF 620
2000
2000
3030
2030
2000
5000
5300
Equivalent
Equivalent
Calls
(a)
(b)
Figure 6: Optimization performances according to different attack targets. The images and their attack noises are presented in Figure 7.
Figure illustrates the convergence of various zeroth-order methods on two test images. For the SSD methods (including the line search version), the parameters are to l = 50 and 0max 2.0. Since BF-SSD uses 1, ,,000 HF evaluations for knowledge distillation training, it begins at N = 1,000. The convergence results demonstrate that HF-SSD) outperforms other methods in this task. Additionally, HF-SSD, BF-SSD, and SPSA exhibit clear advantages over other methods, underscoring the importance of tuning suitable step sizes for the optimization process.
In Figure 7, we present the adversarially attacked test images generated by different optimization approaches for N 2,000, 5,000, and 7,000. For the first test image (a-f), only HF-SSD, BF-SSD, and SPSA successfully flip the output of the HF model under limited HF evaluations (N < 7,000). Similarly, for the second test image (g-1), HF-SSD, BF-SSD, SPSA, and VR-SSD succeed in flipping the HF model output. However, in both cases, we observe that HF-SSD (Figure 7c and Figure 7i) and BS-SSD (Figure 7d and Figure 7j) tend to blur the images more than SPSA (Figure 7b and Figure 7h). This behavior may result from differences in sampling stratogics, such as Haar measure sampling versus Hadamard sampling. From an adversarial attack perspective, a successful attack should flip the model's output without excessively blurring the image. In this regard, SPSA outperforms HF-SSD and BF-SSD, although its loss function remains higher than that of
17

---

Published in Transactions on Machine Learning Research (09/2025)
CD
SPSA
HF-S5D
S
S
S
S
S
S
S
S
(a) GD (Label=5, Predict=5)
(b) SPSA (Label=5, Predict=6)
(c) HF-SSD (Label=5, Predict=6)
S
S
5
(d) BF-SSD (Label=5, Predict=6)
(e) CD (Label=5, Predict=5)
(f) VR-SSD (Label=5, Predict=5)
GD
SPSA
HF-S5D
e
(g) GD (Label=3, Predict=3)
(h) SPSA (Label=3, Predict=8)
(i) HF-SSD (Label=3, Predict=8)
BF-SSD
VR.SSD
3
(i) BF-SSD (Label-3, Predict-8)
(k) CD (Label-3, Predict-3)
(1) VR-SSD (Label=3, Predict-8)
Figure Adversarial examples for MNIST sample #8 (top two rows)
and MNIST sample #18 (bottom two rows) using different methods at N 2,000, 5,000, and 7,000.
the other methods. The reason for this observation requires further investigation, but it is out of the scope of this work.
18

---

Published in Transactions on Machine Learning Research (09/2025)
4.2.3 Soft Prompting Black-bor Language Model
Fine-tuning pre-trained models like BERT or GPT has become cornerstone of modern natural language processing (NLP). These models, trained massive corpora, achieve state-of-the-art performance across wide range of downstream tasks when adapted using task-specific fine-tuning. However, traditional fine- luning involves updating millions or even billions of parameters, making it computationally expensive and prone to overfitting, especially in low-resource settings. To address these challenges, soft prompting has emerged as a lightweight and efficient alternative. Instead of modifying the modcl's internal paramcters, soft prompting introduces learnable embeddings (soft prompts) that are prepended to the input sequence, enabling task adaptation with minimal computational cost. This approach is particularly appealing for tasks requiring minimal intervention in the model's architecture while leveraging its pre-trained knowledge.
Despite the efficiency of soft prompting, its practical applicability faces challenges when dealing with black- box models where gradients with respect to the model parameters are inaccessible. For instance, many com- mercial APIs or proprietary models only provide access to predictions or loss values, making gradient- based optimization infeasible. In such scenarios, zeroth-order optimization becomes a crucial tool. Specifically, in this section, we consider a black-box, pre-trained language classifier fc R - [0, 1], pre-trained tokenizer ft : str - IR L where str is any string of arbitrary length, and the sequence length L is positive integer up to 512 representing the length of the embedding. The goal is to find a soft prompt xx E R67 such that
x= = arg min ((z,)) (C
(4.4)
x
where -E(-,)) is the cross-entropy loss function, and the dataset (z, E str x 10. This particular formulation addresses a binary sentiment analysis task where a given string is classified as expressing a positive (1) or negative (0) sentiment. Since the classifier fo is pre-trained and treated as a black box, gradient information for the loss function is unavailable, necessitating the use of zeroth-order optimization to solve the problem.
For the pre-trained classifier and tokenizer, we employed the BERT model, a state-of-the-aar transformer- basod architecture trained on large corpora. Specifically, we focused on a simplified version of BERT, named DistilBERT. Sentiment analysis on the aclImdb dataset was used as a soft prompting task. This dataset comprises movie reviews categorized into positive and negative sentiments, forming a binary classification problem. A D = 768 dimension soft prompt æ is considered as the input. The transformer's 's parameters were kept frozen to focus optimization on the soft prompt, reducing the degrees of freedom and computational overhead. For computational convenience, we focus on small scale in this example. The HF model, as described in Equation (4.4), was evaluated using 10 samples from aclImdb to approximate the expectation, while the LF model leveraged only 2 samples that were randomly selected from them. Consequently, the evaluation cost ratio between HF and LF was 5:1.
We set the initial starting point at the origin. We let E 50, 0.99, and for the methods without line search we chose fixed step size of 1 x 10-2 The y-axis in the following figure represents the average cross-entropy loss. Figure 8 illustrates the performances of different competing methods with the BF-SSD demonstrating notable advantages.
5 Conclusion
In this work, we proposed and analyzed a bi-fidelity line search scheme designed to accelerate the convergence of zeroth-order optimization algorithms. By constructing surrogate model using both computationally expensive high-fidelity (HF) and inexpensive low-fidelity (LF) objective function evaluations, our approach enables an efficient backtracking line search on the surrogate to determine suitable step sizes. This method significantly reduces the reliance on costly HF evaluations, a common bottleneck in zeroth-order methods. We established thcoretical convergence guarantees for this scheme under standard assumptions and subsequently integrated it into the stochastic subspace descent framework, yielding the Bi-Fidelity Stochastic Subspace Descent (BF-SSD) algorithm.
19

---

Published in Transactions on Machine Learning Research (09/2025)
Soft Prompting DistiliBERT on Sentiment Analysis
6.825 -
Best HF Function value
6.823 x 10
6.822 10
GD
NAG
6.821
CD
FS-SSD
6.82 10
SPSA
GS
10-1
HF-SSD
BF-SSD
VR-SSD
6,818 x
200
400
600
800
1000
1200
1400
1600
Equivalent HF Function Calls
Figure 8: Average cross entropy loss for various zeroth-order optimiz- ers. The BF-SSD method achieves competitive performance while requiring substantially fewer expensive HF function evaluations.
The practical efficacy of BF-SSD was demonstrated through comprehensive experiments across diverse ap- plications, including synthetic optimization benchmark, dual-form kernel ridge regression, black-box ad- versarial attacks, and the fine-tuning of transformer-basedd language models. Across these tasks, BF-SSD consistently outperformed relevant baseline methods, including standard gradient descent, coordinate de- scent, SPSA, and a purely high-fidelity SSD variant, particularly when comparing solution quality achieved for given budget of HF evaluations.
These findings underscore the potential of leveraging bi-fidelity information within stochastic subspace meth- ods to effectively address large-scale, high-dimensional optimization problems where function evaluations are expensive. While the current analysis focuses on deterministic subspace steps, future work could extend the theoretical guarantees to encompass the stochastic nature of gradient estimation within the SSD frame- work. Overall, BF-SSD presents a promising and computationally efficient optimization tool for a variety of challenging real-world applications.
References
Jonathan Barzilai and Jonathan Borwein. Two-point step size gradient methods. IMA Journal of Numerical Analysis, 811):141-148, 01 1988.
Andrea Brilli, Morteza Kimiaei, Giampaolo Liuzzi, and Stefano Lucidi. Worst case complexity bounds for linesearch-type type derivative-free algorithms. Journal of Optimization Theory and Applications, 203(1): 419-454, 2024.
Nicholas Carlini and David Wagner. Towards evaluating the robustness of neural networks. In 2017 IEEE Symposium on Security and Privacy, pp. 39-57, 2017.
Coralia Cartis and Katya Scheinberg. Global convergence rate analysis of unconstrained optimization meth- ods based probabilistic models. Mathematical Programming, 169:337-375, 2018.
Aochuan Chen, Yimeng Zhang, Jinghan Jia, James Diffenderfer, Jiancheng Liu, Konstantinos Parasyris, Yihua Zhang, Zheng Zhang, Bhavya Kailkhura, and Sijia Liu. Deepzero: Scaling up zeroth-order opti- mization for deep model training. Xi preprint NXi:2330.02022 2023.
20

---

Published in Transactions on Machine Learning Research (09/2025)
Pin-Yu Chen, Huan Zhang, Yash Sharma, Jinfeng Yi, and Cho-Jui Hsich. Zoo: Zeroth order optimization based black-box attacks to deep neural networks without training substitute models. In Proceedings of the 10th ACM workshop on artificial intelligence and security, pp. 15 26, 2017.
Nuojin Cheng. Fulli-Fidelity Uncertainty Quantification and Oplimization. PhD thesis, University of Col- orado at Boulder, 2025.
Nuojin Cheng and Alireza Doostan. Langevin bi-fidelity importance sampling for failure probability estima- tion. arXiv preprint arXiv:2503.17796 2025.
Nuojin Cheng, Osman Asif Malik, Subhayan De, Stephen Becker, and Alireza Doostan. Bi-fidelity variational aulo-encoder for uncertainty quantification. Compuler Melhods in Applied Mechanics and Engineering, 421:116798, 2024a.
Nuojin Cheng, Osman Asif Malik, Yiming Xu, Stephen Becker, Alireza Doostan, and Akil Narayan. Subsam- pling of parametric models with bifidelity boosting. STAM/ASA Journal on Uncertainty Quantification. 12(2):213 241, 2024b.
Andrew R. Conn, Nicholas I. M. Gould, and Philippe L. Toint. Trust Region Methods. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2000. doi: 1011137/.97800997718855
Mark Cutler, Thomas J Walsh, and Jonathan How. Reinforcement learning with multi-fidelity simulators. In 2014 IEEE International Conference on Robotics and Automation (ICRA), pp.
IEEE, 2014.
Subhayan De and Alireza Doostan. Neural network training using 11-regularization and bi-fidelity data.
Journal of Computational Physics, 458:111010, 2022.
Subhayan De, Kurt Maute, and Alireza Doostan. Bi-fidelity stochastic gradient descent for structural opti-
mization under uncertainty. Computational Mechanics, 66:745-77 7777777,1111 2020.
Subhayan De, Matthew Reynolds, Malik Hassanaly, Ryan N King, and Alireza Doostan. Bi-fidelity modeling of uncertain and partially unknown systems using deeponets. Computational Mechanics, 71(6):1251 1267, 2023.
Étienne de Montbrun and Sébastion Gerchinovitz. Certified multifidelity zeroth-order optimization. STA M/ASA Journal on Uncerlainly Quantification, 12(4):1135-11)) December 2024. doi: 10.1137/ 22M1178833
John Duchi, Elad Hazan, and Yoram Singer. Adaptive subgradient methods for online learning and stochastic optimization. Journal of machine learning research, 12(7), 2011.
lan J. Goodfellow, Jonathon Shlens, and Christian Szegedy, Explaining and harnessing adversarial examples. arviv preprint ar Xiv: 14 12. 6572, 2014.
Alex A. Gorodetsky, Gianluca Gcraci, Michael S. Eldred, and John Jakeman. A generalized approximate control variate framework for multifidelity uncertainty quantification. Journal of Compulational Physics, 408:109257, 2020.
Dane Grundvig. Line search based optimization using function approximations with tunable accuracy. Master's thesis, Rice University, 2023.
Geoffrey Hinton, Oriol Vinyals, and Jeff Dean. Distilling the knowledge in neural network. In NIPS 2014 Deep Learning Workshop, 2014. arXiv preprint arXiv:1503.02591
Rie Johnson and Tong Zhang- Accelerating stochastic gradient descent using predictive variance reduction. Advances in Neural Information Processing Systems, 26, 2013.
Kirthevasan Kandasany, Gautam Dasarathy, Junier Oliva, Jeff Schneider, and Barnabás Póczos. Gaussian process bandit optimisation with multi- fideliyy evaluations. Advances in Neural Information Processing Systems, 29, 2016.
21

---

Published in Transactions on Machine Learning Research (09/2025)
Hamed Karimi, Julie Nutini, and Mark Schmidt. Linear convergence of gradient and moxximal-gradiin methods under the Polyak-Tajassewier condition. In Joinl European conference on machine learning and knowledge discovery in databases, pp. 795 811. Springer, 2016.
Tamara G. Kolda, Robert Michael Lewis, and Virginia. Torczon. Optimization by direct search: New per- spectives some classical and modern methods. SIAM Review, 45(3):385-455 2003. doi: 10.1137/
5003614450242889.
David Kozak, Stephen Becker, Alireza Doostan, and Luis Tenorio. Stochastic subspace descent. preprinl 1rXie:1904.01145, 2019.
arXiv
David Kozak, Stephen Becker, Alireza Doostan, and Luis Tenorio. A stochastic subspace approach to gradient-free optimization in high dimensions. Computational Optimization and Applications, :9(2):339- 368, 2021.
David Kozak, Cesare Molinari, Lorenzo Rosasco, Luis Tenorio, and Silvia Villa. Zeroth-order optimization with orthogonal random directions. Mathematical Programming. 199(1):1179 1219, 2023.
Sijia Liu, Pin-Yu Chen, Bhavya Kailkhura, Gaoyuan Zhang, Alfred O Hero III, and Pramod K Varshney. A primer zeroth-order optimization in signal processing and machine learning: Principals, recent advances, and applications. TEEE Signal Processing Magazine, 37(5):43-55) -4,, 2020.
Yong Liu, Zirui Zhu, Chaoyu Gong, Minhao Cheng, Cho-Jui Hsieh, and Yang You. Sparse MeZO: Less parameters for better performance in zeroth-order LLM fine-tuning. arXiv preprint arXiv:240.. 15751, 2024.
Maren Mahsereci and Philipp Hennig. Probabilistic line searches for stochastic optimization. Journal of Machine Learning Research, 18119)):1-59, 2017.
Sadhika Malladi, Tianyu Gao, Eshaan Nichani, Alex Damian, Jason D Lee, Danqi Chen, and Sanjeev Arora. Fine-tuning language models with just forward passes. Advances in Neural Information Processing Systems, 36:53038- 53075, 2023.
Y Nesterov. A method of solving convex programming problem with convergence rate o(1/k²). Proceedings of the USSR Academy of Sciences, 269:3.
Yurii Nesterov. Introductory lectures on conver optimization: basic course, volume 87. Springer Science & Business Media, 2013.
Yurii Nesterov and Vladimir Spokoiny. Random gradicnt-free minimization of convex functions. Foundations of Compulational Mathemalics, 17(2):527 566, 2017.
Leo W. T. Ng and Karen E. Willcox. Multifidelity approaches for optimization under uncertainty. Interna- tional Journal for numerical methods in Engineering, 104410)7746-77, 2014.
Lam M. Nguyen, Katya Scheinberg, and Trang H. Tran. Stochastic ISTA /TTA/FISTA adaptive step search algorithms for convex composite optimization. Journal of Optimization Theory and Applications, 205(1): 10, 2025.
Jorge Nocedal and Stephen Wright. Numerical optimization. Springer, 1999.
Courtney Paquette and Katya Scheinberg. A stochastic line search method with expected complexity anal- ysis. SIAM Journal on Optimization, 30(1):349-33) 376, 2020. doi: 10.1137/18M1221
Fabian Pedregosa, Gaöl Varoquaux, Alexandre Gramfort, Vincent Michel, Bertrand Thirion, Olivier Grisel, Mathieu Blondel, Peter Prettenhofer, Ron Weiss, Vincent Dubourg, et al. Scikit-learn: Machine learning in python. Journal of machine learning research, 12(0ct):2825-2888 2830, 2011.
Benjamin Peherstorfer, Karen Willcox, and Max Gunzburger. Survey of multifidelity methods in uncertainty propagation, inference, and optimization. STAM Review, 60(3):550-55) -915000--1110 2018.
22

---

Published in Transactions on Machine Learning Research (09/2025)
Boris T. Polyak. Introduction to optimization. New York, Optimization Software, 1987.
Bobak Shahriari, Kevin Swersky, Ziyu Wang, Ryan P. Adams, and Nando Freitas. Taking the human out of the loop: A review of Bayesian optimization. Proceedings of the IEEE, 104(():148- 11711,,,, January 2015. doi: 10.1 TT.100//RRCC.:25.2222200
James C. Spall. Multivariate stochastic approximation using a simultaneous perturbation gradient approxi- mation. IEEE transactions on automatic control, (7(3)3332 341, 1992.
James C. Spall. Implementation of the simultaneous perturbation algorithm for stochastic optimization. IEEE Transactions on aerospace and electromic systems, 34(3):817 823, 1998.
Mandavilli Srinivas and Lalit M. Patnaik. Genetic algorithms: A survey, Computer, 27(6))17 26, June 1994.
doi: 10.1109/2.209 2.299489..
Tianxiang Sun, Zhengfu He, Hong Qian, Yunhua Zhou, Xuan-Jing Huang, and Xipeng Qiu. BBTv2: Towards a gradient-free future with large language models. In Proceedings of the 2022 Conference on Empiricul Methods in Natural Language Processing, pp. 3916 39300, 2022a.
Tianxiang Sun, Yunfan Shao, Hong Qian, Xuanjing Huang, and Xipeng Qiu. Black-box tuning for language- aoseeies----ervice. International Conference on Machine Learning, pp. 20841 -0855.. PMLR, 2022b.
Shion Takeno, Hitoshi Fukuoka, Yuhki Tsukada, Toshiyuki Koyama, Motoki Shiga, Ichiro Takeuchi, and Masayuki Karasuyama. Multi-fidelity Bayesian optimization with max-value entropy search and its par- allelization. In International Conference on Machine Learning, pp. 9334 9345. PMLR, 2020.
Jian Wu, Saul Jescano-Palmernn Peter I. Frazier, and Andrew Gordon Wilson. Practical multi-fidelity
Bayesian optimization for hyperparameter tuning. In Uncertainty in Artificial Intelligence, pp. 788-78777 188-798. PMLR, 2020.
Ziming Yu, Pan Zhou, Sike Wang, Jia Li, and Hua Huang. Subzero: Random subspace zeroth-order opti- mization for memory-effecinnt Ilm finc-tuning. arXiv preprint arXiv:2410.0890 2024.
Yao Yue and Karl Meerbergen. Accelcrating optimization of parametric linear systems by model order reduction. SIAM Journal on Optimization, 23(2):1344 1370, 2013.
Xinshuai Zhang, Fangfang Xie, Tingwei Ji, Zaoxu Zhu, and Yao Zheng. Multi-fidelityy deep neural net- work surrogate model for aerodynamic shape optimization. Computer Methods in Applied Mechanics and Engineering, 373: 11345, 2021.
Yihua Zhang, Pingzhi Li, Junyuan Hong, Jiaxiang Li, Yimeng Zhang, Wenqing Zheng, Pin- Fin-Yu Chen, Jason D. Lee, Wotao Yin, Mingyi Hong, Zhangyang Wang, Sijia Liu, and Tianlong Chen. Revisiting zeroth-order optimization for mefory-effciiett LLM fine-tuning: A benchmark, 2024.
23

---

Published in Transactions on Machine Learning Research (09/2025)
A Proof of Lemma 2.7
Proof. We let ny evaluations positioned at equispaced points between 2g and 2x + UmaxXa, each sub-interval has length Omax/Tx: /R We also define (x)(a) = (fa) fl + a2%). For each sub-interval, we define the surrogate (); mx) as a linear function connecting these values.
WLOG, we prove the bound in Equation (2.4) holds the interval a E [0, hj with = Omax n and this result can be extended to other sub-intervals, The surrogate v(a) is defined
h h-a (0) + h (1), [),
ns)
and similar definitions of vh (x) hold when a in other sub-intervals. Such linear approximation va((a) satisfics
lpp(a) Pula; (x)) (((() v))(a h h
(v(00 vlaa) + h (w(h) x(a))
(A.1)
for any E [0, Since the Lipschitz constant of vp(a) is strictly controlled by W and the fact that va is unit vector, Equation (A.1) satisfies
Omnax ma - a
(a 0) + W a
W Womax
(p(a)) [x(a; & W
(Amax/ n - aa) S
Va E h].
amax /iaxx/Tkk n2k
Omaxx mmnxx/nx
2n
Since we have
WL(1 + c)a
n 2
cB 11 1
the sup-norm is bounded
(x(a) (a(( < cB11
|22ll
cB
min
30mmax
2(1 + e)T.
2
(1 C)t.
where the last equality stems from the fact that 1/2 and amax 2 c/(cLt L).
B Single-fidelity SSD with Line Search
B.1 Assuming Strong convexxyy
Assumption B.1. Assume the objective function fH and algorithm satisfies the following conditions
1. Pa E DxL are independent random matrices such that E[PxP]] = La and P Ps (D/E) with
D>
2. Objective function fu .R.. - R attains its minimum fr and V f is L-Lipschitz continuous;
Objective function fH JR.. > is y-strongly convex; note y I L.
Theorem B.2. (Single fidelity) With the assumptions of B.1, SSD with line search (either exact line search or backtracking) converges in the sense that f(xx) 4 f and f(xx) 4 f
In particular,
Elf(x+41)] wht (f(zo) f"),
for w E (0,1) where
7f
exact line search
DL
w
min {
2t DL
backtracking.
24

---

Published in Transactions on Machine Learning Research (09/2025)
Proof. Define the filtration F& o(P1..... Pk-1) and Fi {0, with Up the sample space. By Lipschitz continuity, we have
f(xx+1) & f + Vf(())) " (xk+1 - -x,) + xx+1 - xxIl?
(B.1)
2
By defining fe(x)) = f(x) - f" and plugging 2k+1 = 22 - -n P P (ax), Equation (B.1) yields
fe(0x+1) fo(wx)) -as -a.(V Pa P (xxx)) (LP) + LL (P& (
=-a V. + Da
P.PP
(B.2)
24
(-ox +
D
P.P Vf(x))
2l
where the fact P P.P.P.PPP = (D/0)P,,P is applied. We have two line search approaches to determine the step size Ok:
Exact line search:
Ok arg mig f(xe VPPPP Backtracking: for some fixed Omax > 0, 8 (0,8/2d), and E (0, 1),
(B.3)
ax max
max cm Cmmax mEN
(B.4)
S.t. [(xx 0 P. PAP V ((x*)) Omaal Px \((x)))
Ommx
We will prove the convergence for two line scarch methods separately. All the following analyses hold for any Pa satisfying Assumption B.1.
Exact line search According to Equation (B.3), the exact line search method can find the optimal og such that the quadratic term in Equation (B.2) yields -ak + Daf -L/(2dL) for any Pk, thereby
fe( fx())) S 2dL P.P. VP..
With condition on the current filtration Fk, the conditional expectation both sides turn to
[[f((
E (V
+ fa(x)
2dL
(B.5)
V ff(xx)]] (xxx)
2dL
where the equality is from the fact E.P.P. [i = I By invoking the Polyak- inequality,
78
fc (mx) + fe(xx) (1 DL fc(xx))
78
(B.6)
DL
Recursive application yields
Blfc((x+ < (1- 7P DL
fe(xx) = (1-
7l
Blfa (xx)) [ (1 (1- 7l K+1
E[fe(ao)].
(B.7)
DL
DL
Since LSD D and S L, the torm D-LL/DL is loss than 1. Equation (B.7) implics
E B[f(x++1) - ---7
7p
A+1
[((((z)) - f³]] (1- DL)
k+1 (f(xo) -f²),
DL
which proves f(xx) 33 fr and f(xx) 4 s
25

---

Published in Transactions on Machine Learning Research (09/2025)
Backtracking (Showing the existence of feasible set such that the Armijo condition is satisfied.) Follow- ing Equation (B.4), the backtracking method selects the maximal possible step size value that satisfies the Armijo condition with specified parameter B < e/2d and shrinking parameter с < 1. When 0 S ag S L/DL, -Ok + L/2 I -a holds with Haar measure probability one, which implics the following Armijo stopping condition
DajL 28
f(xx+1) - f(xx) an P.
V.
S f(ak) ak PxPP
agt
f(xx) 2d
$ f(xs)) - 3ax P.PP V fxx))l
Therefore, the backtracking terminates when as = amax or on 2 te/DL/ which implies
fa(xx+1) S Se((n)) - (0m (Vf(o)). P; P ((())))
(B.8)
or
tcB
fe(xx+1) S fe(xx)
DL P.Pk
(B.9)
Similar with Equation (B.5), by combining Equations (B.8) and (B.9), and taking expectations conditioned on the filtration Fk, we have
LeB
min
(VT(
+ fe(x))
DL
(B.10)
LcB
-im
3a
[ fa(x))
DL
Similar to Equation (B.6), by invoking the Polyak-Lojjsiowicz incquality,
200788 fe(ax)) + fc(w) (1 minn 27 DT.
26c7
}) Sa(wx).
() < mmin 27,00
DT.
(B.11)
By recursivoly implementing Equation (B.11), we have
[[f(( 1 (1-min min 278amaxx
200y3 DL
}) k+1
fa(wo).
Since 0 < <c<,, convergences f(xx) - f and f(xx)
and < 8 < 0.5, the term 0 < (1 min ( < 1, the
L) f guaranteed.
B.2 Assuming Convexity
Assumption B.3. For the non-strongly convex objective function THF and the algorithm, we make the following assumptions
Ps Dxt are independent random matrices such that E[P,P(] La and P. (D/E) with
D
Objective function PHF DRD - is convex and HF L-Lipschitz continuous;
3. The function fH attains its minimum(s) fo at x= that there exists known R satisfying max {|[ x=I| HE (x)) < (xo)) < R:
Given the above conwex-bat--ootsstoonnyyooooooo assumption, we have the following I convergence result:
26

---

Published in Transactions on Machine Learning Research (09/2025)
Theorem B.4. With backtracking implemented for line search, we have
[[f(x+)] - f I max {
2R
2dLR
k38
kleB
Proof. We use the same notation as in the proof of Thm. B.2. Starting from Equation (B.10), we have
Lep
B(fc(ex+1))Fll min 3
[
fe(x)).
DL
By convexity and the Cauchy-Schwartz inequality, [IV 2 fe(m))/R R. With the fact that (.f(((441) S fe (2x), we have
tcB
Ef2
(xx+1)) - fe (xx) S min
30maxx
DL
2
LcB
E²
fa(x))
S min
30maxx
DL
2R²
tcB
Effe (xx++ (eff((
< min
30max:
DT.
2R2
which further implies
LeB
+ 2R min
30
(B.12)
E. (1k+1)
DL
Applying (B.12) recursively, obtain
teB
2R
2R
2dLR
Bf(())41) < min Bomax:
= max
DL
k:Smmax
klc3
B.3 No convexity assumptions
Assumption B.5. We make the following assumptions:
Ps E Dxt are independent random matrices such that E[P;P[]] La and P Ps (D/0) with
D>6
2. The objective function fur RD - R (or f"1)) attains its minimum f and V f is L-Lipschitz continuous;
When Assumption B.5 holds, have the following T.2 convergence of the gradient result for SSD with line search:
Theorem B.6 (Same Theorem 3.1). With Assumption B.5 holding and backtracking implemented for line search, we have
< max (f(zv) DL(f(zo)-
min
[= (K 1)30max
kE (E.......K"
(K 1)Ec8
That is, l = ((1/(e50maa)) + (L/(clc))) iterations are required to achieve (f(x))] Sc.
Proof. Following Equation (B.10)
Lc3
min 30
| [V)())) s Se(w)) (
DL
27

---

Published in Transactions on Machine Learning Research (09/2025)
which leads to the telescope series
te3
K Σ
min
La 3mmax:
v=xx))]
& Σ (fe( -(fe((x))))
DL
km
km0
[(o)) [x((k)+)) (2000) -p"
Therefore,
PcB3
)) -1
(K + 1) min
B / <
min
30 DL
(f(x0) (²)
ke K)
(f(xo) - f*) DL(f(xo) f4)
max
3amax
tc
A sufficient condition let V (+=))l be e-smaal to let
(f(zo)-f9
D
h 2 max
cEe8
C
Worst Function in the World: Additional Data
0.8
0.9
c 0.99
Method
t= 5
= 10
{ 20
= 10
t 20
8 = 5
{ = 10
{ 20
GD
0.9026
0.9026
0.9026
0.9026
0.9026
0.9026
0.9026
0.9026
0.9026
CD
0.8984
0.8984
0.8984
0.8984
0.8984
0.8984
0.8984
0.8984
0.8984
FS-SSD
2.4495
2.4194
2.3611
2.4495
2.4196
2.3619
2.4497
2.4194
2.3622
SPSA
0.7713
0.7756
0.7502
0.6245
0.4623
0.5549
0.6046
0.6721
0.7006
GS
2.4598
2.4442
2.4129
2.4597
2.4447
2.4144
2.4596
2.4445
2.4150
HF-SSD
0.3620
0.2511
0.2194
0.8667
0.5109
0.3337
3.4582
1.7191
1.0331
BF-SSD
0.3177
0.2947
0.2932
0.2686
0.2526
0.2497
0.2316
0.2104
0.1984
VR-SSD
0.9885
0.9374
0.9178
0.9881
0.9464
0.9154
0.9925
0.9395
0.9121
Table 4: Performance values for different optimization methods across various e and & combinations at N = 5,000. The minimum value in each row is highlighted in bold.
28


Reviews:
[{"review_id": "GHWaJTpAWH", "reviewer": "Reviewer_iGJk", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "**Strengths**\n1. The proposed framework allows to efficiently use low-fidelity oracles in a clever way: it is used to construct an approximate surrogate loss, which is in turn used to tune the step-size using backtracking line search.\n2. A theoretical convergence result is given, showing faster convergence of the proposed method.\n3. The backtracing line search method is then applied to the stochastic subspace descent algorithm, showing promising results.\n4. An important effort has been put into (i) describing what could be used as low-fidelity oracles, and (ii) studying the efficiency of the method on a wide range of problems. Numerical experiments are quite strong and show that the proposed method is quite promising for more practical use.\n\n**Weaknesses**\n1. The study of the number of oracles calls presented in Remark 2.5 may be incorrect: it seems to rely on the bound on $n_k$ presented in Assumption 2.3, which is $n_k \\ge \\frac{WL(1+c)\\alpha_{max}}{c\\beta\\|\\|v_l\\|\\|^2}$, claimed to be $n_k = O(WL)$, which is not necessarily true as $v_k$ goes to zero (which is expected as it is typically the gradient).\n2. The claim in Remark 2.5 does not seem to say that the number of oracle calls can be reduced, since the second term is already the number of oracle calls of black-box GD.\n3. Theoretical results could be stated with more rigor: the assumptions that are used in Theorem 2.4 should be explicitly mentioned in the theorem statement, and Remark 2.5 should be given with a proper proof. In particular, the theorems only hold when $v_k$ is chosen to be the gradient, which is quite limiting for studying a black-box method.\n4. It is not very clear from the current result how the precision of the low-fidelity oracle affects the results: one could expect to see that if a very good low-fidelity oracle exists, there is a speed-up that depends on the ratio between the complexity of the high-fidelity and the low-fidelity oracles.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper proposes novel methods for black box optimization in the bi-fidelity framework, where the function to be optimized can be accessed through high-fidelity but expensive and low-fidelity but cheap oracles.\nThese new methods are based on line-search, which aim to find the best possible step size using backtracking, on a surrogate loss constructed using low-fidelity oracles together with a linear approximation of the function based on a few high-fidelity oracles.\nA convergence result is given in the case where the gradient estimate is the gradient itself, using an Armijo condition.\nAn additional method based on stochastic subspace descent (SSD) is then provided and studied empirically.\nExamples of low-fidelity oracles are provided, together with a numerical study on a synthetic function, on kernel ridge regression, adversarial attacks in machine learning, and soft-prompting for language models.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "The paper already presents nice ideas, and the empirical evaluation is already quite complete and relevant to the discussion.\nMost of my concerns lie in the theoretical part, with some results that may be incorrect:\n1. Assumption 2.3 seems to assume that Equation (2.4) holds, but Equation (2.4) is itself a result of a Lemma that assumes Assumption 2.3. In this sense, what does this assumption mean? Could the author clarify Assumption 2.3 and make it self-contained without reference to the result from a further lemma?\n2. Remark 2.5 looks erroneous, as it seems that it relies on Assumption 2.3, together with the assumption that $\\|\\|v_k\\|\\| = \\Omega(1)$, which is not true as $k \\rightarrow \\infty$. Could the authors provide a proof of this remark and clarify what this result means? More specifically, it seems there is kind of a \"homogeneity issue\", since $WL^2$ is a product of a Lipschitz and two smoothness constants, and $L$ is a smoothness constant, making the sum $WL^2 + L$ look disturbing.\n3. In Remark 2.5, it is claimed that the method allows for reducing the number of high-fidelity oracle calls, but the second term $DL/\\epsilon$ seems to be the complexity of vanilla black-box GD (with good estimators of the gradients). How is that an improvement over this algorithm in terms of oracle calls?\n4. In numerical experiments, all results are presented as \"equivalent number of HF oracle calls\". While this metric makes sense, would the same observations still hold if changing the x-axis, e.g., for wall-clock time? One could expect that the proposed method can increase the time of the training itself.\n5. In practice, how does one choose $n_k$?", "broader_impact_concerns": "No broader impact concerns.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "AUUwuFeyR4", "reviewer": "Reviewer_6TiG", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "This submission is written with clear language and structure. It has the following strengths:\n\n1. The authors proposed a multi-fidelity line-search scheme for zeroth order optimization methods. The authors also proved a sub-linear convergence rate of $1/k$ under standard assumptions, which is the same convergence order of the classical gradient descent method with standard backtracking line search scheme. \n\n2. The authors developed the BF-SSD algorithm, a stochastic zeroth-order optimization method with a bi-fidelity line search that allows for choosing the approximation quality of the gradient by tuning $l$. The authors proved the convergence results of SSD with line search under three separate scenarios: strongly convex, convex, and non-convex.\n\n3. Besides those theoretical convergence rates, the authors also compared the performance of BF-SSD with other zeroth-order optimization methods on one synthetic function and the three real-world applications objective function. The empirical results from numerical experiments are consistent with the theoretical analysis.\n\nThis submission has the following weaknesses:\n\n1. Assumption 2.3 is a bit strong. The authors need to explain more about this assumption 2.3. The authors need to discuss this Assumption 2.3 in the examples of section 2.4.\n\n2. The convergence results of SSD with line search under three separate scenarios: strongly convex, convex, and non-convex should be moved from Appendix B to the main part of the submission as important theoretical contributions.\n\n3. The authors should also add the Nesterov accelerated gradient descent method in the numerical experiments.\n\n4. The authors should change the presentation of Theorem B.2. to add the linear convergence rates.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "In this submission, the authors introduced a bi-fidelity line search scheme tailored for zeroth-order optimization. This method constructed a temporary surrogate model by strategically combining inexpensive low-fidelity (LF) evaluations with accurate high-fidelity (HF) evaluations of the objective function. The authors proved theoretical convergence rates of this bi-fidelity line search scheme under standard assumptions. Moreover, the authors integrated this bi-fidelity strategy into the stochastic subspace descent framework and proposed the bi-fidelity stochastic subspace descent (BF-SSD) algorithm. The authors proved the convergence rates of this novel algorithm under non-convex condition, convex condition and strongly convex condition. The authors also conducted comprehensive numerical experiments on the one synthetic function and three practical applications.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "As summarized in the weakness section, this submission need the following changes before accepted as publication:\n\n1. The authors need to explain more details about Assumption 2.3. In section 2.4, the authors should explain more about how those concrete function examples satisfy the Assumption 2.3.\n\n2. The authors need to move the convergence rates of SSD with line search from Appendix B to the main part of the submission. The authors also need to change the order: first introduce the results of non-convexity, then the results of convexity and finally the results of strong convexity.\n\n3. Please supplement the performance of accelerated gradient descent and compare it with all other zeroth-order optimization methods in the numerical experiments.\n\n4. The authors should add the linear convergence rates in the strongly convexity condition in Theorem B.2. instead of the current asymptotic convergence results. Theorem B.4 and Theorem B.6. include the non-asymptotic convergence rates. It is natural to supplement the non-asymptotic linear convergence rates under the strongly convex condition.", "broader_impact_concerns": "No Broader Impact Concerns", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "ncdm9qqVV4", "reviewer": "Reviewer_PMqd", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strength:\nExpensive black-box optimisation is indeed a pain-point in ML and engineering. A method that claims to reduce costly HF queries is automatically interesting.\n\nThe SSD samples a random low-dimensional subspace (dimension \\ell<<D), then it estimates the gradient projected onto that subspace using finite differences. The author argues that along such a descent direction v_k it costs less computational resource, since the 1D surrogate is a piecewise linear function with few points and a scalar blending with LF; its main difference from traditional surrogate-based 0-th order methods is that rather than modeling the full objective BF-SSD just models the function along the line only; this shares the spirit of arxiv:1301.1942 and is not new in the community. And as the authors acknowledge, SSD that approximates the gradient using finite differences in a random low-dimensional subspace. The integration of bi-fidelity modeling with stochastic subspace descent for zeroth-order optimization is a relatively fresh idea.\n\nDespite length, the narrative is mostly fluent and figures are legible (some fonts like Figure 3 legend could be larger, and use the same legend).\n\nWeaknesses:\nHowever, the idea of using a LF along a low-dimensional space is questionable; In addition, line search on a crude 1D surrogate may not generalize well across problem types, especially with nonlinear LF-HF mismatch; yet the paper provides convergence bounds under technical assumptions (smoothness of the HF objective and LF–HF alignment). It is not clear how this bound helps us to understand the advantage of bi-fidelity modeling better. What is a bit un-usual is to assume that we have access to almost exact gradient v_k's, I understand this is a technical assumption (Assumption 2.3), yet I have to say that this looks overly optimistic for *most* 0-th order applications. To give an example, it is hard to believe that LLM example would satisfy this..or it is a reasonable assumption to make for it.\n\nIn fact, the complexity bound depends on an unknown Lipschitz difference W between HF and scaled LF functions (Assumptions 2.2–2.3) . No guidance is given on estimating W or choosing n_k​; in experiments they just set n_k=1 “for simplicity” , undermining the theoretical premises. \n\nWhile the idea of using bi-fidelity models for step size tuning is novel in the context of SSD, the surrogate model is rudimentary and the analysis applies only under strong assumptions. The empirical results are encouraging but limited in scale and scope.\nI am not sure how practical this proposed method is based on the current manuscript.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper proposes a new algorithm, Bi-Fidelity Stochastic Subspace Descent (BF-SSD), for high-dimensional 0-th-order optimization where evaluating the objective function is expensive and gradients are unavailable. \n\nThe paper introduces a line search method that leverages a combination of cheap low-fidelity (LF) evaluations and a small number of expensive high-fidelity (HF) queries to build a simple 1-dimensional surrogate along the search direction. This surrogate is used to adaptively choose the step size, reducing HF query cost per iteration. The method is combined with Stochastic Subspace Descent (SSD) — a zeroth-order technique that estimates descent directions in randomly sampled low-dimensional subspaces — to form the full BF-SSD algorithm.\n\nThe paper proves theoretical convergence guarantees under the assumption that the descent direction is accurate (i.e., equivalent to the true gradient), and that the difference between HF and LF models is Lipschitz continuous.\n\nThe algorithm is tested on four tasks to support improved performance in terms of HF query efficiency compared to several zeroth-order baselines.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "All below changes are needed for sufficient novelty for the derivative-free community:\n(1) Please extend the experiment to higher dimensions, it is not clear to the reader how well the LF surrogate extends to higher dimensional input space. The mismatch/mis-align between LF and HF models will be more serious in high dimensions, and it will be helpful to see if BF-SSD can actually works for standard high-dimensional optimization (around 100-1000 dimensional functions).\n\n(2)It is interesting to link the method in the context of building suitable LF models in section 2.4; but in practice, defining an appropriate low-fidelity function LF​ is non-trivial. The paper assumes you can just “construct” an LF model that is aligned (Assumption 2.2.): but doesn’t provide a metric for LF-HF alignment; I think at least an assessment of failure modes when LF is poorly correlated with HF should be provided.\n\n(3)Table 3 reports the comparison between HF-SSD and BF-SSD, however, HF-SSD does not have similar performance as BF-SSD, which is counter-intuitive. I have also tried to understand the npz file provided in their supplementary files, yet I cannot reproduce the results reported in this table. Could the authors explain why we should expect HF (high-fidelity) is much worse than bi-fidelity? Should we have more information from HF case? \n\n(4)The abstract calls the approach “computationally efficient for large-scale problems”, yet no experiment demonstrates scale nor reports runtime.  LF models themselves are not free: training the distillation network uses 1000 HF calls upfront , but this is simply tacked onto the first iteration and then ignored in plots. Likewise, Nyström preprocessing for kernel ridge regression is not counted. Wall-clock time or GPU hours are never reported. Also it may be more convincing to have a larger scale of ML tasks, in transformer fine-tuning, they only used DistilBERT + 10 samples, which is hardly representative of a real black-box setting.", "broader_impact_concerns": "Not applicable", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Decision: Accept (poster)