arXiv CSopen access

From Cursed to Competitive: Closing the ZO-FO Gap via Input-to-State Stability

arXiv:2604.25372v1 [math.OC] 28 Apr 2026

From Cursed to Competitive: Closing the ZO–FO Gap via Input-to-State Stability

Amir Ali Farzin School of Engineering Australian National University Canberra, ACT 2601 [email protected]

Philipp Braun School of Engineering Australian National University Canberra, ACT 2601 [email protected]

Iman Shames Department of Electrical and Electronic Engineering University of Melbourne Melbourne, VIC 3010 [email protected]

Abstract While it is generally understood that zeroth-order (ZO) algorithms have an extra dependency on their number of iterations for any choice of parameters, compared to their first-order (FO) counterparts, in this work, we show that under several conditions, in expectation, ZO methods do not suffer from extra dimension dependencies in their convergence rates with respect to their FO counterparts. We look at optimisation algorithms from the dynamical systems perspective and analyse the conditions under which one can formulate the average of a ZO algorithm as the average of its FO counterpart with bounded perturbations with values dependent on design parameters. Then, using input-to-state stability properties, we show ZO methods follow the same decay rate as their FO counterparts and converge to a neighbourhood of the fixed point of FO methods, where its radius depends on the bound of the norm of the perturbations, which can be made arbitrarily small. The theoretical findings are illustrated via numerical examples.

1

Introduction

In this paper, we study a minimisation problem of the form min

x∈Rn

f (x)

(1)

where f : Rn → R is continuously differentiable, integrable, and bounded below. Such problems arise routinely in machine learning, signal processing, and control. While previously, in the literature, it was commonly understood that convergence complexity bounds for zeroth-order (ZO) methods have extra dimension dependencies compared to their first-order (FO) counterparts, in this work, we analyse conditions under which it is guaranteed that ZO methods, in expectation, have no extra dimension dependencies in their convergence rate compared to their FO counterparts. Here, we study iterative algorithms and their convergence properties from the dynamical system perspective and analyse them by studying input-to-state stability (ISS) properties. Motivation: In many modern applications, gradient information is either unavailable, prohibitively expensive to compute, or simply does not exist in closed form. Examples include black-box adversarial attacks [CZS+ 17, YHF+ 18, MMS+ 17], reinforcement learning with opaque simuPreprint.

lators [SHC+ 17, CRS+ 18], hyperparameter tuning [SLA12], simulation-based engineering design [CSV09], and memory-efficient fine-tuning of large language models [MGN+ 23]. In these settings, one must resort to ZO methods that construct gradient surrogates from function evaluations alone. While attractive for their simplicity and broad applicability, ZO methods have traditionally been seen as having fundamentally weaker convergence guarantees than their FO counterparts, independent of the parameter selection. The dimension-dependence gap: The seminal work of [NS17] established the modern framework for ZO optimisation based on Gaussian smoothing and showed that, for smooth convex problems, ZO gradient descent requires at most n times more iterations than standard gradient descent, where n is the dimension. This O(n) factor has become the standard penalty attributed to working without gradient information. Subsequent works have refined and extended these bounds across convex [DJWW15, BG18], nonconvex stochastic [GL13], and structured non-convex settings such as quasar-convex, submodular, and Polyak-Łojasiewicz functions [FPBS25a, FPB+ 25b, FS24]. In each case, the convergence bounds carry an explicit dependence on n beyond what FO methods require. Efforts to reduce this overhead via sparsity [WDBS18, BG22] or effective dimension [YYFL23] typically require additional structural assumptions (e.g., gradient sparsity or a fast-decaying Hessian spectrum) that may not hold in general. Related work: On the ZO optimisation side, building on [NS17, GL13], ZO variants of classical FO algorithms have been developed, including accelerated methods [DGG21], variance-reduced schemes [LKC+ 18], and minimax methods [LLC+ 20, WBMR20, FPB+ 25a, FPB+ 26]. Notably, the MeZO framework [MGN+ 23] demonstrated that ZO-SGD can fine-tune LLMs with billions of parameters at competitive performance, attributed to the favourable loss landscape induced by pre-training; subsequent work [ZCY+ 24] benchmarked broader ZO optimiser families in this setting. On the dynamical systems side, interpreting optimisation algorithms as discrete-time dynamical systems is well established [LRP16, WRJ21, FPBS26, DP18, FPBS25b]. Input-to-state stability (ISS), introduced in [Son89] and extended to discrete-time systems in [JW01], provides a natural framework for analysing robustness of stable equilibria to bounded perturbations and has been widely applied in nonlinear control [Son08, KB23]. Its application to bridge ZO and FO convergence theory is, to the best of our knowledge, novel. A detailed comparison of our results with those of prior ZO methods is provided in Table 2 (Appendix A). Contributions: The main contributions of this paper are as follows. (i) We show that, under suitable conditions, the averaged dynamics of a Gaussian ZO algorithm approximate its FO counterpart, where the approximation error is bounded and can be made arbitrarily small by tuning the smoothing parameter µ, the step size h, and the number of sampled directions t. (ii) Leveraging the ISS framework, we prove that ZO methods inherit the same transient convergence rate as their FO counterparts, i.e., the same exponential/asymptotic decay, and converge to a neighbourhood of the FO fixed point whose radius depends on the perturbation bound but carries no additional dependence on the number of iterations (Theorems 1 and 2). This provides a fundamentally different explanation of ZO convergence behaviour. The dimension dependence does not appear in the iteration complexity but rather in the size of the convergence neighbourhood. Moreover, the framework serves as a general-purpose tool. Whenever the FO algorithm is known to converge, the ISS machinery automatically yields convergence of its ZO counterpart to a neighbourhood of the FO fixed point, without requiring a separate, algorithm-specific proof. (iii) We apply this framework to GD, the heavy ball, and Nesterov’s accelerated gradient. Under strong convexity, we derive explicit perturbation bounds (Lemmas 2, 3, 4) and dimension-free convergence rate results (Theorems 3, 4, 5). For non-strongly-convex objectives, we show that L2 regularisation provides sufficient contraction, enabling the same dimension-free analysis with a controlled trade-off between regularisation strength and neighbourhood size (Lemma 5 and Theorem 6). (iv) We validate the theoretical findings through numerical experiments, confirming that ZO trajectories closely track their FO counterparts with no extra dimension dependence, and demonstrating how µ, h, and t control this proximity. Outline: Preliminaries are introduced in Section 2. In Section 3, the main results on the dimensional independency of ZO methods are given. In Section 4, some illustrative examples are studied. The conclusions and future research directions are discussed in Section 5. Additional details and the proof of the technical results can be found in Appendices A- E.

2

Notation: For x, y ∈ Rd , we write ⟨x, y⟩ = x⊤ y, and let ∥ · ∥ denote the Euclidean norm of a vector and the corresponding spectral norm of a matrix. For δ > 0 and z ∗ ∈ Rd , we define Bδ (z ∗ ) := {z ∈ Rd | ∥z − z ∗ ∥ ≤ δ}. For f : Rn → R, we denote the gradient and Hessian by ∇f and ∇2 f , respectively. The identity matrix of proper size is denoted by I, and Id denotes the identity function Id(s) = s. A continuous function α : [0, a) → [0, ∞) (with a > 0 or a = ∞) is of class K if it is strictly increasing and α(0) = 0. A function α is of class K∞ if α ∈ K and α(s) → ∞ as s → ∞. A continuous function β : [0, ∞) × [0, ∞) → [0, ∞) is of class KL if, for each fixed t ≥ 0, the map s 7→ β(s, t) belongs to K, and for each fixed s ≥ 0, the map t 7→ β(s, t) is decreasing with β(s, t) → 0 as t → ∞.

2

Preliminaries

In this section, we present preliminaries on dynamical systems, ISS, and ZO Gaussian oracles. 2.1

Dynamical systems and perturbations

We study optimisation problem (1) through the lens of discrete-time dynamical systems, interpreting an algorithm as (or z + = w(z)),

zk+1 = w(zk )

(2)

with state z ∈ Rd and update map w : Rd → Rd . A point z e ∈ Rd is an equilibrium (fixed point) of w if w(z e ) = z e . To analyse stability of equilibria, we use standard notions from discrete-time systems theory; see, e.g., [KB23, Ch. 2, Ch. 5]. Definition 1. An equilibrium z e ∈ Rd of (2) is stable if there exist α ∈ K∞ and δ > 0 such that ∥zk − ze ∥ ≤ α(∥z0 − ze ∥) for all k ∈ N and z0 ∈ Bδ (z e ). If z e is not stable, it is called unstable. The equilibrium is asymptotically stable if there exist β ∈ KL and δ > 0 such that ∥zk − ze ∥ ≤ β(∥z0 − ze ∥, k)

∀ k ∈ N,

∀ z0 ∈ Bδ (z e ).

(3)

If β in (3) has the form β(r, k) = M rγ k for M > 0 and γ ∈ (0, 1), then z e is exponentially stable. For an optimisation algorithm, KL-stability or asymptotic stability of an equilibrium of (2) is equivalent to local convergence to the associated critical point. We now consider a perturbed version of (2): zk+1 = w̄(zk , qk ) = w(zk ) + qk ,

(4)

where {qk }k∈N models perturbations (e.g., inexact gradients or noise). When these perturbations are bounded, robustness can be studied via ISS [SW96, JW01]. Definition 2 (Input-to-state stability). The system (4) with zk , qk ∈ Rd is (locally) input-to-state stable (ISS) if there exist δ > 0, β ∈ KL, and γ ∈ K such that, for all sequences {qk }k≥0 ,   ∥zk − z e ∥ ≤ β(∥z0 − z e ∥, k) + γ sup ∥qj ∥ ∀ k ∈ N, ∀ z0 ∈ Bδ (z e ). (5) 0≤j≤k

The ISS estimate (5) ensures that the state remains bounded for bounded inputs and recovers asymptotic stability when qk ≡ 0 [SW96, KB23]. We adopt the additive perturbation structure in (4), which is sufficient for our analysis and simplifies the ISS arguments. A convenient way to verify ISS is via Lyapunov functions. The following theorems show that (4) inherits the transient decay of (2) and converges to a neighbourhood of its equilibria, whenever the perturbations are uniformly bounded and the unperturbed system is exponentially or asymptotically stable. Theorem 1. Consider the system (4) with equilibrium z e ∈ Rn for qk = 0 for k ∈ N, and suppose there exist a locally Lipschitz function V : Rn → R≥0 , constants c ∈ (0, 1) and R > 0, and constants c1 , c2 ∈ R>0 , such that c1 ∥z − z e ∥2 ≤ V (z) ≤ c2 ∥z − z e ∥2 , V (w̄(z, 0)) − V (z) ≤ −c V (z), 3

∀ z ∈ BR (z e ) ∀ z ∈ BR (z e )

(6) (7)

whenever qk = 0 for all k ∈ N. Then there exist a class-K function γ and radii r̄, q̄ > 0 such that, for the perturbed system (4) with supk ∥qk ∥ ≤ q̄ and initial condition satisfying z0 ∈ Br̄ (z e ), the corresponding solution satisfies r   k c2 e ∥zk − z ∥ ≤ k ≥ 0. (8) ∥z0 − z e ∥(1 − c) 2 + γ sup ∥qj ∥ , c1 0≤j≤k In particular, z e is locally exponentially stable for qk = 0, k ∈ N, and locally ISS with respect to perturbations qk ∈ Bq̄ (0), k ∈ N. The function V in Theorem 1 is called Lyapunov function. A proof of Theorem 1 is given in Appendix B. The above theorem shows that the perturbed system exhibits the same local exponential decay rate as the unperturbed system and converges to a neighbourhood of the equilibrium of the unperturbed system, where the size of the neighbourhood depends on the supremum of the norm of the perturbations. For extension to the setting where the unperturbed system is asymptotically stable but not exponentially, we first present the following lemma. Lemma 1. Let φ : [0, ∞) → [0, ∞) be continuous, strictly increasing, with φ(0) = 0, and Id − φ a class-K function. Let {vk }k≥0 be a non-negative sequence satisfying vk+1 ≤ φ(vk ) + d, ∀ k ≥ 0, for some constant d ∈ [0, b], where b = lim r − φ(r). Then the following properties are satisfied. r→∞

(i) There exists a class-KL function β̂ such that, if d = 0, then vk ≤ β̂(v0 , k) for all k ≥ 0. ¯ > 0 such that, if d ≤ d(ε), ¯ (ii) For every ε > 0, there exists d(ε) then lim supk→∞ vk ≤ ε. (iii) For every R > 0, there exists a class-K function γ̂R such that, with β̂ as in (i), vk ≤ β̂(v0 , k) + γ̂R (d), ∀ k ≥ 0, ∀ v0 ∈ [0, R], ∀ d ≥ 0. Theorem 2. Consider the system (4) with equilibrium z e ∈ Rn for qk = 0 for k ∈ N, and suppose there exist a locally Lipschitz function V : Rn → R≥0 , class-K∞ functions α1 , α2 , α3 , and a constant R > 0 such that α1 (∥z − z e ∥) ≤ V (z) ≤ α2 (∥z − z e ∥), ∀ z ∈ BR (z e ), (9) e e V (w̄(z, 0)) − V (z) ≤ −α3 (∥z − z ∥), ∀ z ∈ BR (z ). (10) Then there exist a class-KL function β, a class-K function γ, and radii r̄, q̄ > 0 such that, for the perturbed system (4) with supk ∥qk ∥ ≤ q̄ and z0 ∈ Br̄ (z e ), the solution satisfies   ∥zk − z e ∥ ≤ β(∥z0 − z e ∥, k) + γ sup ∥qj ∥ , k ≥ 0. (11) 0≤j≤k

Proofs of Lemma 1 and Theorem 2 are given in Appendix B. Together with Theorem 1, they show that the perturbed system converges to a neighbourhood of the equilibrium of the unperturbed system with the same transient decay profile, without any extra dependence on the iteration index k, where the neighbourhood size is determined by the perturbation bound. This provides a general mechanism for transferring convergence guarantees from FO algorithms to their ZO counterparts. Whenever an FO algorithm is known to converge, one only needs to verify that the averaged ZO dynamics can be written in the form (4) with bounded perturbations, and the ISS machinery automatically guarantees convergence to a neighbourhood of the same fixed point, inheriting the same decay profile. In Section 3, we analyse conditions under which such a formulation is possible. Crucially, since the perturbations originate from the ZO approximation rather than a model mismatch, disturbances or sensor noise, their magnitude can be manipulated through appropriate parameter selections. 2.2

Gaussian smoothing

Following [NS17], we define the Gaussian-smoothed version of f with parameter µ > 0 as fµ (x) = Eu [f (x + µu)], where u ∼ N (0, of whether f is differentiable, fµ is always h I). Regardless i differentiable with ∇fµ (x) = Eu f (x+µu) u [NS17]. We use the two-point random oracle µ gµ (x) = µ1 (f (x + µu) − f (x)) u,

(12)

which satisfies Eu [gµ (x)] = ∇fµ (x). We call an algorithm a Gaussian ZO method if it uses (12) in place of gradients. 4

Assumption 1. The variance of the oracle gµ defined in (12) is uniformly upper bounded by σ 2 ≥ 0; i.e., Eu [∥gµ (x) − ∇fµ (x)∥2 ] ≤ σ 2 . This is a standard assumption; see, e.g., [MMSN21, LLC+ 20, FPB+ 25a]. The variance t independent oracle evaluations per iteration, yielding  can be reduced by2 averaging Eu ∥gµ (xk ) − ∇fµ (xk )∥ ≤ σ 2 /t while preserving unbiasedness [BG22]. The details of relaxation of Assumption 1 to a more relaxed notion are given in Remark 4.

3

ZO methods as perturbed FO methods

Without further assumptions, there is no guarantee that the expectation of a Gaussian ZO method can be formulated as a perturbed FO method with characterisable bounded perturbations. Here, we identify conditions under which such a formulation is possible. Then, Theorems 1 and 2 imply that the convergence of the Gaussian ZO method has no extra dependence on the problem dimension and follows the same transient decay as the FO method. The proofs of all the lemmas and theorems of this section are given in Appendix C. 3.1

Gradient descent and strong convexity

We focus on the GD algorithm and assume that f is strongly convex. The GD update rule is zk+1 = zk − h∇f (zk ),

(13)

where h > 0 is the step size. Suppose f has Lipschitz gradients with constant L1 > 0. The ZO counterpart of GD is z̃k+1 = z̃k − hgµ (z̃k ).

(14)

From [NS17], we know that E[gµ ] = ∇fµ . Let Uk = {u1 , . . . , uk } and define the averaged (deterministic) iterate of ZO-GD and its update rule as z̄k = EUk [z̃k ],

z̄k+1 = z̄k − hEUk [∇fµ (z̃k )].

(15)

This can be reformulated as z̄k+1 = z̄k − h∇f (z̄k ) + qk ,

(16)

where


qk = h ∇f (z̄k )−EUk [∇fµ (z̃k )] = h ∇f (z̄k )−∇fµ (z̄k ) +h ∇fµ (z̄k )−EUk [∇fµ (z̃k )] . (17) Lemma 2. Consider (13), (16), and
(17). Let f be β-strongly convex with Lipschitz gradients with β constant L1 > 0, and let h ∈ 0, 2L 2 . Then 1 s µ 2hσ 2 sup ∥qk ∥ ≤ h L1 (n + 3)3/2 + hL1 . (18) 2 t(β − 2L21 h) k Lemma 2 shows that the averaged ZO-GD iterates satisfy a perturbed version of the GD dynamics, where the perturbations are bounded and their magnitude can be manipulated through appropriate parameter selection. It is known that, for a smooth strongly convex objective, GD converges to the minimiser exponentially fast [N+ 18]. In [NS17], ZO-GD was analysed under the same setting, and it was shown that the iteration complexity carries an extra dimension-dependent factor for any choice of parameters. The following theorem shows that, instead, the dimension affects only the size of the convergence neighbourhood, and there are parameter choices such that ZO-GD closely follows GD. Theorem 3. Consider (13) and
(16). Let f be β-strongly convex with Lipschitz gradients with β constant L1 > 0, let h ∈ 0, 2L and for r > 0 let Br (z e ) denote a convergence neighbourhood 2 1 around the fixed point z e of GD. Then there exist constants µ̄, q̄ > 0 and t̄ ∈ N (depending on h, β, L1 , and the Lyapunov function of the GD dynamics) such that, for any µ ∈ (0, µ̄) and t ≥ t̄ satisfying supk ∥qk ∥ ≤ q̄ via the bound in Lemma 2, ZO-GD has the same exponential decay rate as GD converges to Br (z e ), and there exists a class-K function α such that s ! 2 µ 2hσ r = α h L1 (n + 3)3/2 + hL1 . (19) 2 t(β − 2L21 h) 5

3.2

Heavy ball method and Nesterov’s accelerated gradient under strong convexity

In this section, we consider both the heavy ball (HB) method and Nesterov’s accelerated gradient (NAG) under strong convexity. The two algorithms share the same augmented state-space structure and differ only in the point at which the gradient is evaluated, so we treat them in a unified framework. The HB and NAG update rules are, respectively, xk+1 = xk + h2 (xk − xk−1 ) − h1 ∇f (xk ),

(20)

xk+1 = xk + h2 (xk − xk−1 ) − h1 ∇f xk + h2 (xk − xk−1 ) , (21) and their ZO counterparts (ZO-HB, ZO-NAG) replace ∇f with gµ defined in (12). Define the  ⊤ x⊤ augmented state vector zk = x⊤ ∈ R2n , the state transition matrix k−1 k   0 In A= , −h2 In (1 + h2 )In and the output map T = [0 In ] ∈ Rn×2n , so that xk = T zk . Note that T Azk = xk + h2 (xk − xk−1 ) is the lookahead point used by NAG. Both algorithms can be written compactly as   0 zk+1 = Azk + , (22) −h1 ∇f (P zk ) where P = T for HB and P = T A for NAG. The ZO counterpart replaces ∇f (P zk ) with gµ (P zk ). Following the same steps as in Section 3.1, conditioning on zk , taking expectations, and adding and subtracting the gradient evaluated at the averaged state, we obtain the averaged dynamics   0 z̄k+1 = Az̄k + + qk , (23) −h1 ∇f (P z̄k ) where the perturbation term is " # 0   . qk = −h1 EUk [∇fµ (P zk )] − EUk [∇f (P zk )] + EUk [∇f (P zk )] − ∇f (P z̄k )

(24)

Equation (23) shows that the averaged trajectory follows the FO dynamical system (HB or NAG, depending on the choice of P ) perturbed by qk . Lemma 3. Consider (22)–(24) with P = T (HB). Let f be β-strongly convex with L1 -Lipschitz gradients, and let h1 , h2 > 0 satisfy ρHB := (1 + 2h2 )2 − h1 (1 + h2 )β + 4h1 h2 L1 + 2h21 L21 < 1. (25) Then r h1 µL1 (n + 3)3/2 h21 L1 σ 2 sup E[∥qk ∥] ≤ + √ . (26) 2 1 − ρHB t k

β−4c A sufficient condition for (25) is h2 ∈ (0, ch1 ) and h1 ∈ 0, 4c2 +4cL with c ∈ 0, β4 . 2 1 +2L 1

Lemma 4. Consider (22)–(24) with P = T A (NAG). Let f be β-strongly convex with L1 -Lipschitz gradients, and let h1 , h2 > 0 satisfy ρNAG := (1 − h1 β + 2h21 L21 )(1 + 2h2 )2 < 1. (27) Then r h1 µL1 (n + 3)3/2 (1 + 2h2 )h21 L1 σ 2 √ sup E[∥qk ∥] ≤ + . (28) 2 1 − ρNAG t k


β 1 A sufficient condition for (27) is h1 ∈ 0, 2L and h2 ∈ 0, 21 √ −1 . 2 2 2 1−h1 β+2h1 L1

1

Both lemmas show that the averaged ZO iterates satisfy a perturbed version of their FO dynamics with bounded, controllable perturbations. The key structural difference is that, for NAG, the strong convexity contraction acts directly on the lookahead error, eliminating the momentum cross term that appears in the HB analysis. For HB, it is well known that the method locally converges exponentially fast [Pol64]; for NAG, global convergence with an accelerated rate is established in [N+ 18]. The following theorems provide a unified perspective that for both methods, the dimension dependence does not appear in the iteration complexity, instead, the convergence is to a neighbourhood of FO dynamics, where the size of the neighbourhood can be made arbitrarily small. 6

Theorem 4. Consider (22) and (23) with P = T . Let f be β-strongly convex with L1 -Lipschitz gradients, let the hypotheses of Lemma 3 be satisfied, and for r > 0 let Br (z e ) denote a convergence neighbourhood around the fixed point z e of HB. Then there exist positive scalars µ̄ and q̄, and t̄ ∈ N such that, for any µ ∈ (0, µ̄) and t ≥ t̄ satisfying supk E[∥qk ∥] ≤ q̄ via the bound in Lemma 3, ZO-HB has the same exponential decay rate as HB and converges to Br (z e ), and there exists a class-K function α such that r  h µL (n + 3)3/2  h21 L1 σ 2 1 1 r=α + √ . (29) 2 1 − ρHB t Theorem 5. Consider (22) and (23) with P = T A. Let f be β-strongly convex with L1 -Lipschitz gradients, let the hypotheses of Lemma 4 be satisfied, and for r > 0 let Br (z e ) denote a convergence neighbourhood around the fixed point z e of NAG. Then there exist µ̄, q̄ > 0 and t̄ ∈ N such that, for any µ ∈ (0, µ̄) and t ≥ t̄ satisfying supk E[∥qk ∥] ≤ q̄ via the bound in Lemma 4, ZO-NAG has the same decay rate as NAG and converges to Br (z e ), and there exists a class-K function α such that r   h µL (n + 3)3/2 (1 + 2h2 )h21 L1 σ 2 1 1 √ + . (30) r=α 2 1 − ρNAG t 3.3

L2 regularisation

In the preceding sections, we showed that strong convexity is sufficient to guarantee that the averaged ZO dynamics can be written as a bounded perturbation of the FO dynamics. In this section, we show that adding an L2 regularisation term to the problem (1) provides sufficient contraction to recover the same perturbation framework, even when f itself is not strongly convex. We therefore consider, in place of (1), the regularised problem minn f (x) + λ∥x∥2 , x∈R

(31)

where λ > 0 is the regularisation parameter. We focus on the GD case; analogous results for HB and NAG are discussed in Remark 11. Setting λ = 1−c 2h , where h > 0 is the step size and c ∈ [0, 1) (thus 1 λ can be arbitrary small but not larger than 2h ), the GD update (13) applied to (31) becomes zk+1 = czk − h∇f (zk ),

(32)

z̃k+1 = cz̃k − hgµ (z̃k ),

(33)

z̄k+1 = cz̄k − h∇f (z̄k ) + qk ,

(34)

the ZO-GD update (14) becomes and the averaged ZO update is where


qk = h ∇f (z̄k ) − EUk [∇fµ (z̃k )] = h ∇f (z̄k ) − ∇fµ (z̄k ) + h ∇fµ (z̄k ) − EUk [∇fµ (z̃k )] . (35) Lemma 5. Consider (32), (34), and (35). Let f be a possibly non-strongly-convex function with √ 2 −c
Lipschitz gradients with constant L1 > 0, and let h ∈ 0, 2−c . Then 4L1 s µ 2h2 σ 2 sup ∥qk ∥ ≤ h L1 (n + 3)3/2 + hL1 . (36) 2 2 t(1 − c − 4cL1 h − 8L21 h2 ) k Remark 1. Since f has L1 -Lipschitz gradients, regularising with λ > L21 makes the overall objective strongly convex, and one could directly apply the results of Section 3.1. With the parametrisation √ √
2L1 (1−c) 2−c2 −c λ = 1−c , where c ∈ [0, 1) and h ∈ 0, , we have λ > √ ≥ 2L1 (1 − c). 2 2h 4L1 2−c −c 1 1 Consequently, for c ≤ 1 − 2√2 the regularised objective is strongly convex, whereas for c > 1 − 2√ 2 it need not be. This reveals a trade-off governed by c. A small value of c yields strong convexity, a larger admissible range for h, and more control over the magnitude of the perturbations, but it also introduces a larger shift of the minimiser from that of the original problem. Conversely, a large value of c preserves proximity to the original minimiser and still provides sufficient contraction for the perturbation to remain bounded, at the expense of a smaller admissible h. Thus, if fidelity to the original minimiser is the priority, one should choose c large and control the perturbation through t and h; if stability properties are more important, one can choose c smaller and benefit from smaller bounds on the perturbation. 7

Lemma 5 shows that the contraction induced by the regularisation term ensures the averaged ZO-GD iterates satisfy a perturbed version of the GD dynamics, where the perturbations are bounded and controllable through the algorithm parameters. The following theorem establishes that, whenever GD converges to a fixed point, ZO-GD inherits the same convergence rate with no extra dimension dependence in the iteration complexity; the dimension affects only the size of the convergence neighbourhood. Theorem 6. Consider (32) and (34). Let f be a possibly non-strongly-convex function with Lipschitz √ 2 −c
gradients with constant L1 > 0, and let h ∈ 0, 2−c , and for r > 0 let Br (z e ) denote 4L1 e a convergence neighbourhood around the fixed point z of GD. Suppose that GD applied to the regularised problem (31) converges to a fixed point z e . Then there exist constants µ̄, q̄ > 0 and t̄ ∈ N (depending on h, c, L1 , and the Lyapunov function of the GD dynamics) such that, for any µ ∈ (0, µ̄) and t ≥ t̄ satisfying supk ∥qk ∥ ≤ q̄ via the bound in Lemma 5, ZO-GD has the same decay rate as GD and converges to Br (z e ) and there exists a class-K function α such that s ! µ 2h2 σ 2 3/2 . (37) r = α h L1 (n + 3) + hL1 2 t(1 − c2 − 4cL1 h − 8L21 h2 ) In the next section, we illustrate the theoretical findings through numerical examples.

4

Numerical examples

We illustrate the main results through two settings: a strongly convex quadratic (Section 4.1) and a non-convex neural network with L2 regularisation (Section 4.2). Extended parameter studies are provided in Appendices E.1 and E.2. 4.1

Quadratic objective

We consider minx∈Rn x⊤ Ax, where n = 1000 and A is positive definite with condition number 100 and maximum eigenvalue 100, comparing GD and ZO-GD (averaged over 5 runs). Figure 1 shows the effect of varying h ∈ {10−4 , 10−5 , 10−7 } with µ = 10−6 and t = 1 fixed. When h is too large, the perturbation exceeds the admissible bound of Theorem 1 and ZO-GD diverges. As h decreases, ZO-GD converges to a progressively smaller neighbourhood of the GD trajectory, with no extra dimension dependence. Analogous sweeps over µ and t, as well as parameter compensation experiments, are presented in Appendix E.1. In all cases, the observations are consistent with the bound in Lemma 2: h and µ control the smoothing bias, while h and t control the variance term. GD ZO-GD h = 1e 04

1013 106

3.2 × 100

GD ZO-GD h = 1e 05

3 × 100

2.8 × 100

10 15

2 × 100

x x

x x

10 8

x x

GD ZO-GD h = 1e 07

3 × 100

10 1

10 22

2.6 × 100 2.4 × 100

10 29 2.2 × 100

10 36 0

100

200

Iteration

300

400

500

0

GD ZO-GD h = 1e 04

1013 106

250

500

750

1000

Iteration

1250

1500

1750

2000

0

3 × 100

20000

30000

40000

Iteration

50000

GD ZO-GD h = 1e 07

3 × 100

10 1

2.8 × 100

10 15

2 × 100

x x

x x

10 8

x x

10000

3.2 × 100

GD ZO-GD h = 1e 05

10 22

2.6 × 100 2.4 × 100

10 29 2.2 × 100

10 36 0.0

0.2

0.4

Wall-clock time (s)

0.6

0.8

0.0

0.5

1.0

1.5

Wall-clock time (s)

2.0

2.5

0

20

40

Wall-clock time (s)

60

80

Figure 1: Quadratic objective (n = 1000): effect of the step size h on GD (solid) and ZO-GD (dashed). Reducing h shrinks the perturbation and brings ZO-GD closer to GD. 4.2

Binary classification with a neural network

We consider binary classification on MNIST (digits 0 vs 1) using a two-layer fully connected network with ReLU activations and sigmoid output (d = 784, H = 128 (hidden layer neurons), n = 100,609 8

parameters). We minimise the L2 -regularised logistic loss with λ = 10−2 and M = 2000 training samples. Since the logistic loss composed with a neural network is non-convex, this example falls in the framework of Section 3.3. Convergence comparison. Figure 2 shows the training loss for all six methods (GD, HB, NAG and their ZO counterparts). We set h = 1 × 10−3 for GD and ZO-GD, h1 = 1 × 10−3 , h2 = 5 × 10−4 for HB/NAG and their ZO counterparts with smoothing parameter µ = 10−7 and number of sampled directions t = 7. The ZO trajectories closely track their FO counterparts, and all methods reach 99.55% training accuracy and 99.40% test accuracy. It takes approximately 4 times longer to run the ZO methods compared to their FO counterpart. Parameter sensitivity experiments confirming the same phenomena as in the quadratic case are provided in Appendix E.2. Gradient Descent

Heavy Ball

8 × 10 1

GD (FO) ZO-GD ( =1e-07, t=7)

7 × 10 1

6 × 10 1

6 × 10 1

6 × 10 1

Training Loss

7 × 10 1

5 × 10 1

5 × 10 1

0

200

400

600

Iteration k

800

1000

NAG (FO) ZO-NAG ( =1e-07, t=7)

5 × 10 1

4 × 10 1

4 × 10 1

4 × 10 1

Nesterov's Accelerated Gradient

8 × 10 1

HB (FO) ZO-HB ( =1e-07, t=7)

7 × 10 1

Training Loss

Training Loss

8 × 10 1

0

200

400

600

Iteration k

800

1000

0

200

400

600

Iteration k

800

1000

Figure 2: Training loss for GD, HB, NAG and their ZO counterparts on a two-layer neural network (n = 100,609). Solid: FO; dashed: ZO. All methods achieve 99.55% train / 99.40% test accuracy. Dimension scaling. To directly validate the absence of extra dimension dependence, we vary H ∈ {16, 32, 64, 128, 256} (n from ≈12,700 to ≈201,000) and run GD and ZO-GD with identical parameters for 1000 iterations. We set h = 5 × 10−4 , t = 10, and µ = 10−8 . Figure 3 (left) shows that ZO-GD tracks GD at the same rate across all dimensions. Figure 3 (right) confirms the periteration time ratio remains approximately constant. Table 1 reports iteration counts to reach within 20% of the FO final loss: the ratio ZO/FO stays close to 1, confirming no extra dimension-dependent factor in the iteration complexity. 9 × 10 1

5.8

ZO-GD (n=50,305) GD (n=100,609) ZO-GD (n=100,609) GD (n=201,217) ZO-GD (n=201,217)

5.6

Time ratio (ZO / FO) per run

GD (n=12,577) ZO-GD (n=12,577) GD (n=25,153) ZO-GD (n=25,153) GD (n=50,305)

Training Loss

8 × 10 1

7 × 10 1

5.4 5.2 5.0 4.8 4.6

0

200

400

Iteration k

600

800

1000

Number of parameters n

Figure 3: Dimension scaling: Training loss for GD (solid) and ZO-GD (dashed) across five network sizes.

5

Conclusion and future directions

In this work, we revisited the commonly held belief that ZO methods inherently suffer from additional dimension dependence in their iteration complexity compared to their FO counterparts. By adopting a dynamical systems perspective, we showed that, under suitable conditions, the averaged dynamics of ZO algorithms can be interpreted as perturbed versions of their FO counterparts with bounded, controllable perturbations. Leveraging ISS arguments, we proved that ZO methods inherit the same convergence rates as FO methods in expectation, converging to a neighbourhood of the FO fixed point 9

105

Table 1: Dimension scaling: iterations to reach within 20% of FO final loss reduction. H

n

FO final

ZO final

Iter ZO/FO

Time ZO/FO

16 32 64 128 256

12,577 25,153 50,305 100,609 201,217

0.6922 0.6720 0.6761 0.6529 0.6866

0.6920 0.6721 0.6798 0.6584 0.7017

0.99 1 1.03 1.03 1.05

5.5 4.9 5.4 5.7 4.5

whose size can be made arbitrarily small by tuning the algorithm parameters. These findings suggest that the practical efficiency of ZO methods may be closer to that of FO methods than previously understood, particularly when gradient information is unavailable or expensive. Promising directions for future work include extending the analysis to broader classes of non-convex and non-smooth problems, establishing high-probability guarantees rather than results in expectation, and exploring adaptive schemes for tuning the smoothing and step-size parameters to optimise the trade-off between convergence speed and neighbourhood size.

References [BG18] Krishnakumar Balasubramanian and Saeed Ghadimi. Zeroth-order (non)-convex stochastic optimization via conditional gradient and gradient updates. In Advances in Neural Information Processing Systems, volume 31, pages 3455–3464, 2018. [BG22] Krishnakumar Balasubramanian and Saeed Ghadimi. Zeroth-order nonconvex stochastic optimization: Handling constraints, high dimensionality, and saddle points. Foundations of Computational Mathematics, 22(1):35–76, 2022. [CRS+ 18] Krzysztof Choromanski, Mark Rowland, Vikas Sindhwani, Richard Turner, and Adrian Weller. Structured evolution with compact architectures for scalable policy optimization. In International Conference on Machine Learning, pages 970–978. PMLR, 2018. [CSV09] Andrew R Conn, Katya Scheinberg, and Luis N Vicente. Introduction to derivative-free optimization. SIAM, 2009. [CZS+ 17] Pin-Yu Chen, Huan Zhang, Yash Sharma, Jinfeng Yi, and Cho-Jui Hsieh. 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, pages 15–26, 2017. [DGG21] Pavel Dvurechensky, Eduard Gorbunov, and Alexander Gasnikov. An accelerated directional derivative method for smooth stochastic convex optimization. European Journal of Operational Research, 288(2):601–621, 2021. [DJWW15] John C. Duchi, Michael I. Jordan, Martin J. Wainwright, and Andre Wibisono. Optimal rates for zero-order convex optimization: The power of two function evaluations. IEEE Transactions on Information Theory, 61(5):2788–2806, 2015. [DP18] Constantinos Daskalakis and Ioannis Panageas. The limit points of (optimistic) gradient descent in min-max optimization. Advances in neural information processing systems, 31, 2018. [FPB+ 25a] Amir Ali Farzin, Yuen-Man Pun, Philipp Braun, Antoine Lesage-Landry, Youssef Diouane, and Iman Shames. Min-max optimisation for nonconvex-nonconcave functions using a random zeroth-order extragradient algorithm. Transactions on Machine Learning Research, 2025. [FPB+ 25b] Amir Ali Farzin, Yuen-Man Pun, Philipp Braun, Tyler Summers, and Iman Shames. Minimisation of submodular functions using gaussian zeroth-order random oracles. arXiv preprint arXiv:2510.15257, 2025. 10

[FPB+ 26] Amir Ali Farzin, Yuen-Man Pun, Philipp Braun, Tyler Summers, and Iman Shames. Solving the offline and online min-max problem of non-smooth submodular-concave functions: A zeroth-order approach. arXiv preprint arXiv:2601.21243, 2026. [FPBS25a] Amir Ali Farzin, Yuen-Man Pun, Philipp Braun, and Iman Shames. Minimisation of quasar-convex functions using random zeroth-order oracles. arXiv preprint arXiv:2505.02281, 2025. [FPBS25b] Amir Ali Farzin, Yuen-Man Pun, Philipp Braun, and Iman Shames. Properties of fixed points of generalised extra gradient methods applied to min-max problems. IEEE Control Systems Letters, 2025. [FPBS26] Amir Ali Farzin, Yuen-Man Pun, Philipp Braun, and Iman Shames. On the stability connection between discrete-time algorithms and their resolution odes: Applications to min-max optimisation. arXiv preprint arXiv:2603.01430, 2026. [FS24] Amir Ali Farzin and Iman Shames. Minimisation of polyak-łojasewicz functions using random zeroth-order oracles. In 2024 European Control Conference (ECC), pages 3207–3212. IEEE, 2024. [GL13] Saeed Ghadimi and Guanghui Lan. Stochastic first- and zeroth-order methods for nonconvex stochastic programming. SIAM Journal on Optimization, 23(4):2341–2368, 2013. [HJ12] Roger A Horn and Charles R Johnson. Matrix analysis. Cambridge university press, 2012. [JW01] Zhong-Ping Jiang and Yuan Wang. Input-to-state stability for discrete-time nonlinear systems. Automatica, 37(6):857–869, 2001. [KB23] Christopher M Kellett and Philipp Braun. Introduction to Nonlinear Control: Stability, Control Design, and Estimation. Princeton University Press, 2023. [Kel14] Christopher M Kellett. A compendium of comparison function results. Mathematics of Control, Signals, and Systems, 26(3):339–374, 2014. [LKC+ 18] Sijia Liu, Bhavya Kailkhura, Pin-Yu Chen, Paishun Ting, Shiyu Chang, and Lisa Amini. Zeroth-order stochastic variance reduction for nonconvex optimization. Advances in Neural Information Processing Systems, 31, 2018. [LLC+ 20] Sijia Liu, Songtao Lu, Xiangyi Chen, Yao Feng, Kaidi Xu, Abdullah Al-Dujaili, Mingyi Hong, and Una-May O’Reilly. Min-max optimization without gradients: Convergence and applications to black-box evasion and poisoning attacks. In International Conference on Machine Learning, pages 6282–6293. PMLR, 2020. [LRP16] Laurent Lessard, Benjamin Recht, and Andrew Packard. Analysis and design of optimization algorithms via integral quadratic constraints. In SIAM Journal on Optimization, volume 26, pages 57–95. SIAM, 2016. [MGN+ 23] 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. [MMS+ 17] Aleksander Madry, Aleksandar Makelov, Ludwig Schmidt, Dimitris Tsipras, and Adrian Vladu. Towards deep learning models resistant to adversarial attacks. arXiv preprint arXiv:1706.06083, 2017. [MMSN21] Alejandro I Maass, Chris Manzie, Iman Shames, and Hayato Nakada. Zeroth-order optimization on subsets of symmetric matrices with application to mpc tuning. IEEE Transactions on Control Systems Technology, 30(4):1654–1667, 2021. [N+ 18] Yurii Nesterov et al. Lectures on convex optimization, volume 137. Springer, 2018. [NS17] Yurii Nesterov and Vladimir Spokoiny. Random gradient-free minimization of convex functions. Foundations of Computational Mathematics, 17:527–566, 2017. 11

[Pol64] Boris T Polyak. Some methods of speeding up the convergence of iteration methods. Ussr computational mathematics and mathematical physics, 4(5):1–17, 1964. [SHC+ 17] Tim Salimans, Jonathan Ho, Xi Chen, Szymon Sidor, and Ilya Sutskever. Evolution strategies as a scalable alternative to reinforcement learning. arXiv preprint arXiv:1703.03864, 2017. [SLA12] Jasper Snoek, Hugo Larochelle, and Ryan P Adams. Practical bayesian optimization of machine learning algorithms. Advances in Neural Information Processing Systems, 25, 2012. [Son89] Eduardo D. Sontag. Smooth stabilization implies coprime factorization. IEEE Transactions on Automatic Control, 34(4):435–443, 1989. [Son08] Eduardo D Sontag. Input to state stability: Basic concepts and results. In Nonlinear and optimal control theory: lectures given at the CIME summer school held in Cetraro, Italy June 19–29, 2004, pages 163–220. Springer, 2008. [SW96] Eduardo D Sontag and Yuan Wang. New characterizations of input-to-state stability. IEEE transactions on automatic control, 41(9):1283–1294, 1996. [WBMR20] Zhongruo Wang, Krishnakumar Balasubramanian, Shiqian Ma, and Meisam Razaviyayn. Zeroth-order algorithms for nonconvex minimax problems with improved complexities. arXiv preprint arXiv:2001.07819, 2020. [WDBS18] Yining Wang, Simon S. Du, Sivaraman Balakrishnan, and Aarti Singh. Stochastic zeroth-order optimization in high dimensions. In International Conference on Artificial Intelligence and Statistics, pages 1356–1365. PMLR, 2018. [WRJ21] Ashia C. Wilson, Benjamin Recht, and Michael I. Jordan. A lyapunov analysis of accelerated methods in optimization. Journal of Machine Learning Research, 22(113):1– 34, 2021. [YHF+ 18] Haishan Ye, Zhichao Huang, Cong Fang, Chris Junchi Li, and Tong Zhang. Hessianaware zeroth-order optimization for black-box adversarial attack. arXiv preprint arXiv:1812.11377, 2018. [YYFL23] Pengyun Yue, Long Yang, Cong Fang, and Zhouchen Lin. Zeroth-order optimization with weak dimension dependency. In The Thirty Sixth Annual Conference on Learning Theory, pages 4429–4472. PMLR, 2023. [ZCY+ 24] Yihua Zhang, Pingzhi Chen, Haowen Yuan, Kaidi Jin, Mingyi Hong, Sijia Liu, et al. Revisiting zeroth-order optimization for memory-efficient LLM fine-tuning: A benchmark. arXiv preprint arXiv:2402.11592, 2024.

A

Related works and Basic definitions

In this section, we provide a detailed comparison of our results with prior work on zeroth-order optimisation, followed by the basic definitions used throughout the paper. A.1

Comparison with prior work

Table 2 summarises the convergence guarantees of representative ZO methods and highlights how our framework differs from existing analyses. In all prior works, the iteration complexity carries an explicit factor depending on the problem dimension n (or an effective dimension neff ), and convergence is to the exact optimiser. By contrast, our ISS-based framework eliminates the dimension factor from the iteration complexity entirely; the dimension dependence instead appears in the radius of the convergence neighbourhood, which can be made arbitrarily small by tuning the algorithm parameters µ, h, and t. Moreover, our framework applies uniformly to ZO-GD, ZO-HB, and ZO-NAG, whereas each prior result requires a separate, algorithm-specific analysis. 12

Table 2: Comparison of convergence guarantees for ZO methods. Reference

Assumptions on f

Iteration Dim. dependence

Convergence target

[NS17] [DJWW15] [GL13]

convex / s.c. convex nonconvex convex / s.c; Hessian regularity convex; sparse gradients

O(n) factor O(n) factor O(n) factor

Minimiser Minimiser ∥∇f ∥2 ≤ ε

O(neff ), neff ≤ n

Minimiser

O(s log n), s-sparse

Minimiser

No explicit dependence

Neighbourhood of FO fixed point (can be made arbitrarily small)

[YYFL23] [WDBS18] This work

A.2

s.c. (Thms. 3–5); nonconvex + L2 reg. (Thm. 6)

Basic definitions

We now recall the standard regularity and convexity conditions assumed throughout the paper. Definition 3 (Lipschitz continuity and gradients). A continuous function f : Rn → R is globally Lipschitz with constant L0 > 0 if ∥f (x) − f (y)∥ ≤ L0 ∥x − y∥ for all x, y ∈ Rn . If f is C 1 , the gradient of f is globally Lipschitz with constant L1 > 0 if ∥∇f (x) − ∇f (y)∥ ≤ L1 ∥x − y∥,

∀ x, y ∈ Rn .

(38)

Definition 4 (Strong convexity). Let f : Rn → R be a C 1 function. Then f is said to be β-strongly convex with constant β > 0 if there exists a constant β > 0 such that f (y) ≥ f (x) + ∇f (x)⊤ (y − x) +

B

β ∥x − y∥2 , 2

∀ x, y ∈ Rn .

(39)

Additional details and proofs of Section 2

In this section, we present additional explanations and proofs of the results in Section 2. Lemma 6 ([Kel14, Lem. 10]). Let α : R≥0 → R≥0 be a class-K function. Then for any a, b ≥ 0 and any ε > 0,
(40) α(a + b) ≤ α((1 + ε)a) + α (1 + 1ε )b . Proof of Theorem 1. Let R > 0 be such that (6)–(7) hold on BR (z e ). Since V is locally Lipschitz on BR (z e ), there exist constants LV > 0 satisfying |V (x) − V (y)| ≤ LV ∥x − y∥ e

for all x, y ∈ BR (z ). Write λ := 1−c ∈ (0, 1) and δ := supj≥0 ∥qj ∥. Let r̄, q̄ > 0, let z0 ∈ Br̄ (z e ), let {zk }k∈N be defined through zk+1 = w(zk ) + qk for k ∈ N and {qk }k∈N ⊂ Bq̄ (z e ), and let r̄, q̄ be defined such that {zk }k∈N ⊂ Br̄ (z e ). (The existence of r̄, q̄ > 0 with this property will be shown later.) Using the Lipschitz property of V and (7), for the perturbed successor zk+1 = w(zk ) + qk , the estimate
V (zk+1 ) = V w(zk ) + qk ≤ V (w(zk )) + LV ∥qk ∥ ≤ λ V (zk ) + LV ∥qk ∥ (41) is obtained. Iterating (41) from k = 0 yields V (zk ) ≤ λk V (z0 ) + LV

k−1 X

λk−1−j ∥qj ∥ ≤ λk V (z0 ) + LV δ

j=0

k−1 X

λj ≤ λk V (z0 ) +

j=0

P∞

j

LV δ, c

where the geometric series was bounded by j=0 λ = 1/c. Applying c1 ∥zk − z e ∥2 ≤ V (zk ) and the sandwich bound V (z0 ) ≤ c2 ∥z0 − z e ∥2 , we obtain e 2 LV ∥zk − z e ∥2 ≤ λk c−1 1 c2 ∥z0 − z ∥ + cc1 δ

13

or e

r

∥zk − z ∥ ≤ where γ(δ) =

q

Lv cc1

√

  k c2 ∥z0 − z e ∥λ 2 + γ sup ∥qj ∥ , c1 0≤j≤k

k ≥ 0,

δ, which is a K function.

It remains to choose r̄ > 0 and q̄ > 0 so that {zk }k∈N ⊂ BR (z e ). For z ∈ Rn such that V (z) ≤ 2 ℓ = Rc2 , (6) implies that ∥z − z e ∥2 ≤ R2 , i.e., x ∈ BR (z e ), and we can conclude {z ∈ Rn | V (z) ≤ R2 R2 e n c2 } ⊂ BR (z ). From (41), the set {z ∈ R | V (z) ≤ c2 } is forward-invariant for the perturbed

system whenever λ ℓ + LV q̄ ≤ ℓ, i.e. whenever q̄ ≤ c ℓ/LV . Setting r̄ :=

ℓ , c1

q̄ :=

cℓ LV

ensures that ∥z0 − z e ∥ ≤ r̄ implies V (z0 ) ≤ ℓ, and the sublevel set {z ∈ Rn | V (z) ≤ ℓ} is forward-invariant under the perturbed dynamics, so the Lipschitz bounds and the decrease condition remain valid for all k ≥ 0. This completes the proof. Proof of Lemma 1. Item (i): Since φ is continuous with φ(0) = 0 and φ(s) < s for s > 0, the iterates φ(k) (s) := φ ◦ · · · ◦ φ(s) (k times) satisfy φ(k) (s) ↓ 0 as k → ∞ for each fixed s > 0. Define β̂(s, k) := φ(k) (s). This is continuous, non-decreasing in s, decreasing to 0 in k, and β̂(0, k) = 0; so it can be upper-bounded by a class-KL function, which we denote β̂. This gives (i). ¯ := µ/2. If d ≤ d(ε) ¯ and Item (ii): Fix ε > 0. Since φ(ε) < ε, set µ := ε − φ(ε) > 0 and take d(ε) vk ≤ ε, then vk+1 ≤ φ(ε) + d ≤ φ(ε) + µ/2 = ε − µ/2 < ε. This in particular implies that φ(s) ∈ [0, ε] for all s ∈ [0, ε]. Meanwhile, if vk > ε we have vk+1 ≤ φ(vk ) + d < vk − µ + d ≤ vk − µ/2 (using µ ≤ s − φ(s) for s ≥ ε by continuity and compactness on [ε, v0 ]). Hence, for vk > ε, vk decreases by at least µ/2 per step and there exists K ∈ N such that vK ∈ [0, ε]. Item (iii): Fix R > 0 throughout and let β̂ ∈ KL be the function from part (i). Since Id − φ : [0, ∞) → [0, b] is class-K and d ≤ b, its inverse is strictly increasing and is defined at d. Define

σ(d) := (Id − φ)−1 (d),

so that φ(σ(d)) + d = σ(d). We record two facts, both following directly from Id − φ being strictly increasing. We will analyse vk+1 , for the cases where vk ∈ [0, σ(d)] or vk ∈ (σ(d), ∞). Forward invariance of [0, σ(d)] under the mapping ϕ(·)+d: if vk ≤ σ(d), then vk+1 ≤ φ(σ(d))+d = σ(d), so [0, σ(d)] is forward-invariant. Strict decrease of vk on (σ(d), ∞): if vk > σ(d), then vk − φ(vk ) > σ(d) − φ(σ(d)) = d, hence ϕ(vk ) + d < vk and vk+1 ≤ φ(vk ) + d < vk . Combining these two properties, the sequence {vk } is bounded above by max(v0 , σ(d)) for all k ≥ 0: it cannot increase above v0 once it starts decreasing (which happens immediately if v0 > σ(d)), and it is trapped once it enters [0, σ(d)]. In particular, for v0 ∈ [0, R]:
vk ≤ max R, σ(d) for all k ≥ 0. (42) We claim that for all v0 ∈ [0, R], d ≥ 0, and k ≥ 0, vk ≤ β̂(v0 , k) + ηR (d),

(43)

for some function ηR with ηR (0) = 0 that depends only on d (and R), not on v0 or k. To see this, write vk = β̂(v0 , k) + ek where ek := vk − β̂(v0 , k) is the excess of the perturbed trajectory over the unperturbed bound. We need to show that ek is bounded by a function of d (uniformly over v0 ∈ [0, R] and k ≥ 0). Consider the unperturbed comparison sequence u0 = v0 , uk+1 = φ(uk ), so uk = β̂(v0 , k) by part (i). Both sequences satisfy: vk+1 ≤ φ(vk ) + d, uk+1 = φ(uk ), 14

v0 = s ∈ [0, R], u0 = s.

Subtracting and using φ(vk ) − φ(uk ) ≤ vk − uk (which follows from Id − φ being non-decreasing, so vk − φ(vk ) ≥ uk − φ(uk ) when vk ≥ uk ): vk+1 − uk+1 ≤ (vk − uk ) + d

whenever vk ≥ uk .

Since v0 − u0 = 0, iterating this inequality gives vk − uk ≤ kd whenever vk ≥ uk persists for all steps up to k. At first glance, kd grows without bound. But recall that vk eventually enters and remains in [0, σ(d)]. Let K ∗ be the first time vk ≤ σ(d); from that point on, vk ≤ σ(d) for all subsequent times. The strict decrease above σ(d) gives a minimum decrease per step: for vk ∈ [σ(d), max(R, σ(d))], we have vk − vk+1 ≥ vk − φ(vk ) − d ≥ (Id − φ)(vk ) − d. On the compact set [σ(d), max(R, σ(d))], the continuous function (Id − φ)(s) − d achieves its minimum, which is strictly positive (being zero only at s = σ(d), which is the left endpoint). Hence the entry time satisfies   R − σ(d)  K∗ ≤   min (Id − φ)(s) − d
   s∈[σ(d),R]

∗

when R > σ(d) (and K = 0 otherwise). In particular, K ∗ depends only on R and d (not on v0 beyond the constraint v0 ≤ R). Now we bound the excess for all k: • For k ≤ K ∗ : using vk ≤ max(R, σ(d)) from (42) and β̂(v0 , k) ≥ 0, we get ek ≤ max(R, σ(d)). But more precisely, ek = vk − uk ≤ K ∗ · d by the iterated bound above. • For k > K ∗ : vk ≤ σ(d) and β̂(v0 , k) ≥ 0, so ek ≤ σ(d). Therefore, for all k ≥ 0:

ek ≤ max K ∗ · d, σ(d) =: ηR (d).

Since K ∗ depends only on R and d, so does ηR (d). Moreover, ηR (0) = 0 because K ∗ = 0 and σ(0) = 0 when d = 0. The function ηR is non-negative, finite, nondecreasing in d, and satisfies ηR (0) = 0. By regularisation, there exists γ̂R ∈ K with ηR (d) ≤ γ̂R (d) for all d ≥ 0. Substituting into (43) and using uk = β̂(v0 , k): vk ≤ β̂(v0 , k) + γ̂R (d),

k ≥ 0,

for all v0 ∈ [0, R] and d ≥ 0. Lemma 7 ([JW01, Lemma B.1]). For any K∞ function α, there exists a K∞ function α̂ such that the following holds: • α̂(r) ≤ α(r) for all r ≥ 0; and, • Id − α̂ ∈ K. Proof of the Theorem 2. Let R > 0 be chosen so that (6)–(7) hold on BR (z e ), and let LV > 0 be a Lipschitz constant for V on BR (z e ). Write δ := supj≥0 ∥qj ∥. Provided the trajectory stays in BR (z e ), the Lipschitz property of V gives V (zk+1 ) = V (w(zk ) + qk ) ≤ V (w(zk )) + LV ∥qk ∥. Combining (44) with the decrease condition (7) and the upper sandwich bound (6) yields
V (zk+1 ) ≤ V (zk ) − α3 (∥zk − z e ∥) + LV δ ≤ V (zk ) − α3 α2−1 (V (zk )) + LV δ,

(44)

(45)

where the second inequality uses ∥zk − z e ∥ ≥ α2−1 (V (zk )) from the right-hand side of (6). Define ρ := α3 ◦ α2−1 , which is class-K∞ , and set vk := V (zk ). Then (45) reads vk+1 ≤ vk − ρ(vk ) + LV δ. 15

Define φ(s) := s − ρ(s). Since ρ is class-K∞ , using Lemma 7, without loss of generality we can say φ ∈ K, thus φ is continuous, φ(0) = 0, and Id − φ ∈ K∞ . Thus, we have vk+1 ≤ φ(vk ) + LV δ, which is exactly the setting of Lemma 1 with d = LV δ. Since Theorem 2 is concerned with properties on a compact set, let v0 ∈ [0, R̄], where R̄ > 0 such that {z : V (z) ≤ R̄} ⊂ BR (z e ). Applying part (iii) of the lemma gives a class-KL function β̂ and a class-K function γ̂ such that V (zk ) ≤ β̂(V (z0 ), k) + γ̂R̄ (LV δ). According to the sandwich bounds (9) the inequalities V (z0 ) ≤ α2 (∥z0 − z e ∥) and ∥zk − z e ∥ ≤ α1−1 (V (zk )) are satisfied. Using Lemma 6 with ε > 0 fixed, we obtain

∥zk − z e ∥ ≤ α1−1 (1 + ϵ) β̂(α2 (∥z0 − z e ∥), k) + α1−1 (1 + 1ϵ )γ̂R̄ (LV δ) . With the definitions
β(s, k) := α1−1 (1 + ϵ)β̂(α2 (s), k) ,

γ(s) := α1−1 (1 + 1ϵ )γ̂R̄ (LV s) ,

we see that β is class-KL (as a composition of class-K∞ and class-KL functions) and γ is class-K, giving   ∥zk − z e ∥ ≤ β(∥z0 − z e ∥, k) + γ sup ∥qj ∥ . 0≤j≤k

To complete the proof, we need trajectories to remain in {z ∈ Rn : V (z) ≤ R̄} ⊂ BR (z e ). From the proof of Lemma 1(iii), we know that vk ≤ max(v0 , σ(LV δ)), so forward invariance of {z ∈ Rn : V (z) ≤ R̄} holds provided σ(LV q̄) ≤ R̄, i.e., q̄ ≤ ρ(R̄)/LV . Setting r̄ := α2−1 (R̄),

q̄ :=

ρ(R̄) , LV

we have that ∥z0 − z e ∥ ≤ r̄ implies V (z0 ) ≤ R̄, and supk ∥qk ∥ ≤ q̄ guarantees σ(LV q̄) = R̄, so the sublevel set {z ∈ Rn : V (z) ≤ R̄} is forward-invariant under the perturbed dynamics. Hence, the proof is complete. Lemma 8. Let f be β-strongly convex. Then fµ is β-strongly convex. Proof. From the strong convexity of f , we have f (y) ≥ f (x) + ⟨∇f (x), y − x⟩ + β2 ∥y − x∥2 . Replacing x and y with x + µu and y + µu, respectively, and taking expectation with respect to u, we obtain fµ (y) ≥ fµ (x) + ⟨∇fµ (x), y − x⟩ + β2 ∥y − x∥2 , which completes the proof.

C

Proof of lemmas and theorems of Section 3

Proof of Lemma 2. Consider qk given in (17), which satisfies ∥qk ∥ ≤ h∥∇f (z̄k ) − ∇fµ (z̄k )∥ + h∥∇fµ (z̄k ) − EUk [∇fµ (z̃k )]∥.

(46)

Using [NS17, Lemma 3], we know that ∥∇f (z̄k ) − ∇fµ (z̄k )∥ ≤ µ2 L1 (n + 3)3/2 . Thus, the first term on the right-hand side (46) is bounded and can be manipulated through the choice of h and µ. For the second term on the right-hand side of (46), we have ∥EUk [∇fµ (z̃k )] − ∇fµ (z̄k )∥ ≤ ∥EUk [∇fµ (z̃k ) − ∇fµ (z̄k )]∥ ≤ EUk [∥∇fµ (z̃k ) − ∇fµ (z̄k )∥] ≤ L1 EUk [∥z̃k − z̄k ∥] (47) Now let ek = z̃k − z̄k , i.e., we need to find an upper bound on E[∥ek ∥] or E[∥ek ∥2 ] to prove the assertion. Combining (14), (16) and (17), we know that ek+1 = ek − hgµ (z̃k ) + hEUk [∇fµ (z̃k )] = ek − h(gµ (z̃k ) − ∇fµ (z̃k )) − h(∇fµ (z̃k ) − EUk [∇fµ (z̃k )]). 16

To proceed, we denote vk = gµ (z̃k ) − ∇fµ (z̃k ) and lk = ∇fµ (z̃k ) − EUk [∇fµ (z̃k )] and thus ek+1 = ek − h(vk + lk ). Considering Assumption 1, the variance of the random oracle is bounded, 2 and we have E[∥vk ∥2 ] ≤ σt . Moreover, for lk we have EUk [∥lk ∥] = EUk [∥∇fµ (z̃k ) − ∇fµ (z̄k ) + ∇fµ (z̄k ) − EUk [∇fµ (z̃k )]∥] ≤ EUk [∥∇fµ (z̃k ) − ∇fµ (z̄k )∥] + EUk [∥∇fµ (z̄k ) − ∇fµ (z̃k )∥] ≤ 2L1 EUk [∥ek ∥]

(48)

Thus we get σ (49) E[∥ek+1 ∥] ≤ (1 + 2Lh)E[∥ek ∥] + h √ , t which does not provide a uniform bound on EUk [∥z̃k − z̄k ∥] yet. To proceed, we use the strong convexity properties of f . Focusing on EUk [∥ek ∥2 ] and using the fact that E[vk ] = 0, we get EUk [∥ek+1 ∥2 ] = EUk [∥ek ∥2 ] + h2 EUk [∥vk + lk ∥2 ] − 2hEUk [⟨ek , lk ⟩].

(50)

Moreover, we have EUk [⟨ek , lk ⟩] = EUk [⟨ek , ∇fµ (z̃k )⟩] since EUk [⟨ek , EUk [∇fµ (z̃k )]⟩] = 0. Now, the strong convexity of f and consequently strong convexity of fµ (Lemma 8) implies that fµ (z̄k ) ≥ fµ (z̃k ) + ⟨∇fµ (z̃k ), z̄k − z̃k ⟩ +

β ∥z̄k − z̃k ∥2 2

or

β EU [∥ek ∥2 ]. 2 k Moreover, combining strong convexity of fµ and Jensen’s inequality, leads to the estimates EUk [⟨ek , ∇fµ (z̃k )⟩] ≥ EUk [fµ (z̃k )] − fµ (z̄k ) +

EUk [fµ (z̃k )] ≥ fµ (EUk [z̃k ]) ≥ fµ (z̄k ), and hence EUk [⟨ek , ∇fµ (z̃k )⟩] ≥

β EU [∥ek ∥2 ]. 2 k

(51)

Thus, substituting (51) in (50), we have EUk [∥ek+1 ∥2 ] ≤ (1 − hβ)EUk [∥ek ∥2 ] + h2 EUk [∥vk + lk ∥2 ].

(52)

2

We know that ∥vk + lk ∥2 ≤ 2∥vk ∥2 + 2∥lk ∥2 and E[∥vk ∥2 ] ≤ σt . Similar to the estimates in (48), we can hence conclude the following chain of equalities and inequalities: E[∥lk ∥2 ] = EUk [∥∇fµ (z̃k ) − ∇fµ (z̄k ) + ∇fµ (z̄k ) − EUk [∇fµ (z̃k )]∥2 ] = E[∥∇fµ (z̃k ) − ∇fµ (z̄k )∥2 + ∥∇fµ (z̄k ) − EUk [∇fµ (z̃k )]∥2 ] − 2E[⟨∇fµ (z̃k ) − ∇fµ (z̄k ), EUk [∇fµ (z̃k )]∇fµ (z̄k )⟩]

(53)

2

≤ E[∥∇fµ (z̃k ) − ∇fµ (z̄k )∥ ] ≤ L21 E[∥ek ∥2 ]. Thus 2h2 σ 2 (54) t β and where ρ is defined as ρ = 1 − hβ + 2L21 h2 . To have ρ < 1, we require h < 2L 2 . Since e0 = 0, 1 by induction we get EUk [∥ek+1 ∥2 ] ≤ ρEUk [∥ek ∥2 ] +

2h2 σ 2 2h2 σ 2 2hσ 2 = = . 2 2 t(1 − ρ) t(hβ − 2L1 h ) t(β − 2L21 h) Hence, combining all the estimates above, we get s µ 2hσ 2 3/2 sup ∥qk ∥ ≤ h L1 (n + 3) + hL1 2 t(β − 2L21 h) k EUk [∥ek ∥2 ] ≤

which completes the proof. 17

(55)

Proof of Theorem 3. By Lemma 2, the averaged ZO-GD dynamics (16) take the form of the perturbed system (4) with perturbation satisfying the bound (18). It remains to verify that the unperturbed GD dynamics satisfy the hypotheses of Theorem 1.
2 For β-strongly convex f with L1 -Lipschitz gradients and h ∈ 0, β+L , from [N+ 18, Thm. 2.1.15] 1 we have  2hβL1  ∥zk+1 − z ∗ ∥2 ≤ 1 − ∥zk − z ∗ ∥2 , β + L1 where z ∗ = arg minx∈Rn f (x). Consider the Lyapunov function V (z) = ∥z − z ∗ ∥2 . Then (6) holds 1 with c1 = c2 = 1, and the above contraction gives (7) with c = 2hβL β+L1 ∈ (0, 1). All hypotheses of Theorem 1 are therefore satisfied. The ISS bound (8) then guarantees that ZO-GD converges to a neighbourhood of z ∗ with the same
exponential rate (1 − c)k/2 as GD. The radius of this neighbourhood is governed by γ supk ∥qk ∥ , which, combined with the perturbation bound from Lemma 2, yields (19). Since the bound (18) can be made smaller than q̄ by choosing µ sufficiently small and t sufficiently large, such parameter choices always exist. Proof of Lemma 3. Consider qk given in (24) for HB, which satisfies ∥qk ∥ ≤ h1 EUk [∇fµ (xk ) − ∇f (xk )] + h1 EUk [∇f (xk )] − ∇f (x̄k ) .

(56)

Using [NS17, Lem. 3], we know that ∥∇fµ (x) − ∇f (x)∥ ≤ µ2 L1 (n + 3)3/2 uniformly in x. Thus, the first term in the right-hand side (56) is bounded and can be ensured to be arbitrarily small by selecting h1 and µ appropriately. For the second term in the right-hand side of (56), by L1 -Lipschitz continuity of ∇f and Jensen’s inequality, we have EUk [∇f (xk )] − ∇f (x̄k )

≤ L1 EUk [∥xk − x̄k ∥].

(57)

Now let ek = zk − z̄k so that xk − x̄k = T ek where T = [0 In ]. Thus we are interested in bounding E[∥T ek ∥] or E[∥T ek ∥2 ]. Subtracting the averaged dynamics (23) (for HB) from the ZO-HB update gives   0 ek+1 = A ek − h1 , (58) vk + l k where vk = gµ (xk ) − ∇fµ (xk ) with E[vk ] = 0, and lk = ∇fµ (xk ) − EUk [∇fµ (xk )]. As a next step, let xek = T ek . From the lower block of (58) we have xek+1 = (1 + h2 ) xek − h2 xek−1 − h1 (vk + lk ).

(59)

Let ak = E[∥xek ∥2 ]. Squaring (59), taking expectations, and using E[⟨y, vk ⟩] = 0 for any y independent of uk , leads to the expression     ak+1 = E ∥(1 + h2 )xek − h2 xek−1 ∥2 + h21 E[∥vk + lk ∥2 ] − 2h1 E ⟨(1 + h2 )xek − h2 xek−1 , lk ⟩ . (60) We derive bounds for each term in (60) separately. For the deterministic quadratic term, expanding and applying Young’s inequality provides the upper bound ∥(1 + h2 )xek − h2 xek−1 ∥2 = (1 + h2 )2 ∥xek ∥2 + h22 ∥xek−1 ∥2 − 2h2 (1 + h2 )⟨xek−1 , xek ⟩ ≤ (1 + h2 )(1 + 2h2 ) ∥xek ∥2 + h2 (1 + 2h2 ) ∥xek−1 ∥2 .

(61)

For the inner product, since E[xek ] = 0 and EUk [∇fµ (xk )] is deterministic, we get E[⟨xek , lk ⟩] = E[⟨xk − x̄k , ∇fµ (xk )⟩]. By β-strong convexity of fµ (from Lemma 8), the inequality β ∥xk − x̄k ∥2 2 is thus satisfied. Taking expectations and applying Jensen’s inequality to E[fµ (xk )] ≥ fµ (x̄k ) then yields E[⟨xek , lk ⟩] ≥ β2 ak . (62) e e For the cross term E[⟨xk−1 , lk ⟩], using E[∥lk ∥] ≤ 2L1 E[∥xk ∥] (from L1 -smoothness of fµ and the triangle inequality) and Young’s inequality, we know that
|E[⟨xek−1 , lk ⟩]| ≤ L1 ak−1 + ak . (63) ⟨∇fµ (xk ), xk − x̄k ⟩ ≥ fµ (xk ) − fµ (x̄k ) +

18

Since we assume that the variance of the random oracle is bounded by E[∥vk ∥2 ] ≤ σ 2 /t, for lk , 2 similar to (53), it holds that E[∥lk ∥2 ] ≤ L21 ak and thus E[∥vk + lk ∥2 ] ≤ 2σt + 8L21 ak . Substituting all bounds into (60) and letting ϕ = (1 + h2 )(1 + 2h2 ) − h1 (1 + h2 )β + 2h1 h2 L1 + 2h21 L21 , c0 = h2 (1 + 2h2 ) + 2h1 h2 L1 ,

(64) (65)

we get the two-step recursion 2h21 σ 2 . (66) t As a next step, we define mk = max{ak , ak−1 }. Since ak ≤ mk and ak−1 ≤ mk , (66) gives 2h2 σ 2 ak+1 ≤ (ϕ + c0 ) mk + 1t . Since additionally ak ≤ mk holds, mk+1 can be upper bounded by   2h2 σ 2 (67) mk+1 = max{ak+1 , ak } ≤ max (ϕ + c0 ) mk + 1 , mk , t ak+1 ≤ ϕ ak + c0 ak−1 +

From m0 = 0 (as e0 = 0), we get 2h21 σ 2 , (68) t which is a standard one-step contraction provided ρHB := ϕ + c0 < 1. Computing ρHB explicitly, we get ρHB = ϕ + c0 = (1 + 2h2 )2 − h1 (1 + h2 )β + 4h1 h2 L1 + 2h21 L21 . (69) Thus ρHB < 1 requires mk+1 = max{ak+1 , ak } ≤ (ϕ + c0 ) mk +

h1 (1 + h2 )β > 4h2 + 4h22 + 4h1 h2 L1 + 2h21 L21

(70)

to hold. To find a sufficient condition for (70) to be satisfied let h2 ∈ (0, ch1 ). Then it is sufficient to have h1 β > 2h21 L21 + 4ch1 + 4c2 h21 + 4ch21 L1

or equivalently

β > 2h1 L21 + 4c + 4c2 h1 + 4ch1 L1 .

β−4c Thus we need c ∈ (0, β4 ) and h1 ∈ (0, 4c2 +4cL 2 ) to guarantee that (70) is satisfied, providing 1 +2L1 the bounds stated in Lemma 3.

Since m0 = 0 (as e0 = 0), iterating (68) leads to sup ak ≤ sup mk ≤ k

k

2h21 σ 2 t(1 − ρHB )

and therefore, by Jensen’s inequality we have √

h1 σ sup E[∥xk − x̄k ∥] ≤ sup ak ≤ √ t k k

r

2 . 1 − ρHB

Combining (56), the smoothing bias bound, (57), and (71) we can thus conclude that r h1 µ L1 (n + 3)3/2 h1 σ 2 sup E[∥qk ∥] ≤ + h 1 L1 · √ 2 1 − ρHB t k r h1 µ L1 (n + 3)3/2 h21 L1 σ 2 = + √ , 2 1 − ρHB t

(71)

(72)

which completes the proof. Proof of Theorem 4. By Lemma 3, the averaged ZO-HB dynamics (23) take the form of the perturbed system (4) with bounded perturbation. For β-strongly convex f with L1 -Lipschitz gradients and parameters satisfying the hypotheses of Lemma 3, HB locally converges to the minimiser z ∗ exponentially fast [Pol64]. By the converse Lyapunov theorem for exponentially stable discrete-time systems (see, e.g., [KB23, Thm. 5.6]), there exists a Lyapunov function V satisfying (6)–(7) on some ball BR (z e ), R > 0, in the augmented state space R2n . All hypotheses of Theorem 1 are therefore satisfied, and the conclusion follows by the same argument as in the proof of Theorem 3. 19

Proof of Lemma 4. Consider qk given in (24) for NAG, which satisfies ∥qk ∥ ≤ h1 EUk [∇fµ (T Azk ) − ∇f (T Azk )] + h1 EUk [∇f (T Azk )] − ∇f (T Az̄k ) .

(73)

Using [NS17, Lem. 3], we know that ∥∇fµ (x) − ∇f (x)∥ ≤ µ2 L1 (n + 3)3/2 uniformly in x. Thus, the first term in the right-hand side (73) is bounded and can be controlled by the choice of h1 and µ. For the second term in the right-hand side of (73), by L1 -Lipschitz continuity of ∇f and Jensen’s inequality, we have EUk [∇f (T Azk )] − ∇f (T Az̄k ) ≤ L1 EUk [∥T Azk − T Az̄k ∥] = L1 EUk [∥T Aek ∥],

(74)

where ek = zk − z̄k . Thus, we need to derive an upper bound on E[∥T Aek ∥] or E[∥T Aek ∥2 ]. Let xek = T ek and define the lookahead error yke = T Aek = (1 + h2 )xek − h2 xek−1 . Subtracting the averaged dynamics (23) (for NAG) from the ZO-NAG update gives   0 ek+1 = A ek − h1 , (75) vk + l k where vk = gµ (T Azk ) − ∇fµ (T Azk ) with E[vk ] = 0, and lk = ∇fµ (T Azk ) − EUk [∇fµ (T Azk )]. From the lower block of (75) we obtain xek+1 = yke − h1 (vk + lk ).

(76)

Let ak = E[∥xek ∥2 ] and bk = E[∥yke ∥2 ]. Squaring (76), taking expectations, and using the fact that E[⟨yke , vk ⟩] = 0 leads to the equation ak+1 = bk + h21 E[∥vk + lk ∥2 ] − 2h1 E[⟨yke , lk ⟩].

(77)

To proceed, we derive upper bounds for each term in the right-hand side. For the inner product, since E[yke ] = 0 and EUk [∇fµ (T Azk )] is deterministic, we get E[⟨yke , lk ⟩] = E[⟨T Aek , ∇fµ (T Azk )⟩]. By β-strong convexity of fµ (Lemma 8), it additionally holds that ⟨∇fµ (T Azk ), T Azk − T Az̄k ⟩ ≥ fµ (T Azk ) − fµ (T Az̄k ) +

β ∥T Aek ∥2 . 2

Taking expectations and applying Jensen’s inequality to E[fµ (T Azk )] ≥ fµ (E[T Azk ]) = fµ (T Az̄k ) allows us to write β E[⟨yke , lk ⟩] ≥ bk . (78) 2 This is the key advantage of the NAG structure compared to the HB structure. The strong convexity contraction acts directly on bk = E[∥yke ∥2 ], which is the same quantity that appears with a positive sign in (77). There is no separate momentum cross term to bound. For the noise, since lk = ∇fµ (T Azk ) − EUk [∇fµ (T Azk )], using L1 -smoothness and similar to the process in (48) and (53), E[∥lk ∥] ≤ 2L1 E[∥yke ∥], so E[∥lk ∥2 ] ≤ L21 bk and with E[∥vk ∥2 ] ≤ σ 2 /t, 2 we get E[∥vk + lk ∥2 ] ≤ 2σt + 2L21 bk . Substituting these estimates into (77) provides the upper bound 2h2 σ 2 ak+1 ≤ (1 − h1 β + 2h21 L21 ) bk + 1 . (79) t Let ψ = 1 − h1 β + 2h21 L21 . For the quadratic term bk = E[∥yke ∥2 ], expanding and applying Young’s inequality exactly as in (61) yields bk ≤ (1 + h2 )(1 + 2h2 ) ak + h2 (1 + 2h2 ) ak−1 .

(80)

Substituting (80) into (79) and letting ϕ = ψ (1 + h2 )(1 + 2h2 ), c0 = ψ h2 (1 + 2h2 ), we get the two-step recursion ak+1 ≤ ϕ ak + c0 ak−1 + 20

2h21 σ 2 . t

(81)

Next, we define mk = max{ak , ak−1 }. Since ak ≤ mk and ak−1 ≤ mk , (81) gives ak+1 ≤ 2h2 σ 2 (ϕ + c0 ) mk + 1t and from ak ≤ mk we have  2h2 σ 2 (82) mk+1 = max{ak+1 , ak } ≤ max (ϕ + c0 ) mk + 1t , mk . Since m0 = 0 (as e0 = 0), we get 2h21 σ 2 , (83) t which is a standard one-step contraction provided ρNAG := ϕ + c0 < 1 and computing ρNAG explicitly yields mk+1 ≤ (ϕ + c0 ) mk +

ρNAG = ψ (1 + 2h2 )2 = (1 − h1 β + 2h21 L21 )(1 + 2h2 )2 .

(84)

2

Thus ρNAG < 1 requires ψ(1 + 2h2 ) < 1, i.e., (1 + 2h2 )2 <

1 . 1 − h1 β + 2h21 L21

(85)

β 1 −1/2 Choosing h1 ∈ (0, 8L − 1)) guarantees (85). 2 ) ensures ψ ∈ (0, 1), and then any h2 ∈ (0, 4 (ψ 1

Since m0 = 0, iterating (83), we then have sup ak ≤ sup mk ≤ k

k

2h21 σ 2 . t(1 − ρNAG )

2

To bound E[∥T Aek ∥], note that ∥T Aek ∥ = ∥yke ∥2 = bk , and from (80) we obtain sup E[∥yke ∥2 ] ≤ (1 + 2h2 )2 sup mk ≤ k

k

2(1 + 2h2 )2 h21 σ 2 . t(1 − ρNAG )

Finally, by Jensen’s inequality, r (1 + 2h2 )h1 σ 2 √ , 1 − ρNAG t k and combining (73), the smoothing bias bound, (74), and (86) it holds that r (1 + 2h2 )h1 σ h1 µ L1 (n + 3)3/2 2 √ + h1 L1 · sup E[∥qk ∥] ≤ 2 1 − ρNAG t k r h1 µ L1 (n + 3)3/2 (1 + 2h2 )h21 L1 σ 2 √ = + , 2 1 − ρNAG t sup E[∥T Aek ∥] ≤

(86)

(87)

which completes the proof. Proof of Theorem 5. By Lemma 4, the averaged ZO-NAG dynamics (23) take the form of the perturbed system (4) with bounded perturbation. For β-strongly convex f with L1 -Lipschitz gradients and parameters satisfying the hypotheses of Lemma 4, NAG converges to the minimiser z e exponentially fast [N+ 18]. By the converse Lyapunov theorem for exponentially stable discrete-time systems (see, e.g., [KB23, Thm. 5.6]), there exists a Lyapunov function V satisfying (6), (7) on some ball BR (z e ), R > 0, in the augmented state space R2n . All hypotheses of Theorem 1 are therefore satisfied, and the conclusion follows by the same argument as in the proof of Theorem 3. Proof of Lemma 5. Consider qk given in (35), which satisfies ∥qk ∥ ≤ h∥∇f (z̄k ) − ∇fµ (z̄k )∥ + h∥∇fµ (z̄k ) − EUk [∇fµ (z̃k )]∥

(88)

Using [NS17, Lem. 3], we know that ∥∇f (z̄k ) − ∇fµ (z̄k )∥ ≤ µ2 L1 (n + 3)3/2 . Thus, the first term

in the right-hand side of (88) is bounded and can be manipulated by choosing h and µ appropriately. For the second term in the right-hand side of (88), we have ∥EUk [∇fµ (z̃k )] − ∇fµ (z̄k )∥ ≤ ∥EUk [∇fµ (z̃k ) − ∇fµ (z̄k )]∥ ≤ EUk [∥∇fµ (z̃k ) − ∇fµ (z̄k )∥] ≤ L1 EUk [∥z̃k − z̄k ∥] 21

(89)

Now, let ek = z̃k − z̄k . As a next step, we derive an upper bound on E[∥ek ∥2 ]. We know that ek+1 = cek − hgµ (z̃k ) + hEUk [∇fµ (z̃k )] = ek − h(gµ (z̃k ) − ∇fµ (z̃k )) − h(∇fµ (z̃k ) − EUk [∇fµ (z̃k )]) To proceed, we use the notation vk = gµ (z̃k ) − ∇fµ (z̃k ) and lk = ∇fµ (z̃k ) − EUk [∇fµ (z̃k )] and thus ek+1 = cek − h(vk + lk ). We assume the variance of the random oracle is bounded and get 2 E[∥vk ∥2 ] ≤ σt . Moreover, for lk the following chain of inequalities is satisfied: E[∥lk ∥] ≤ E[∥∇fµ (z̃k ) − ∇fµ (z̄k )∥] + E[∥∇fµ (z̄k ) − ∇fµ (z̃k )∥] ≤ 2L1 E[∥ek ∥].

(90)

2

From the fact that E[vk ] = 0, we get the following bound on EUk [∥ek+1 ∥ ], EUk [∥ek+1 ∥2 ] = c2 EUk [∥ek ∥2 ] + h2 EUk [∥vk + lk ∥2 ] − 2chEUk [⟨ek , lk ⟩] 2

2

2

2

≤ c EUk [∥ek ∥ ] + h EUk [∥vk + lk ∥ ] + 2chEUk [∥ek ∥∥lk ∥] 2

2

2

(91) (92)

2

≤ (c + 4cL1 h)EUk [∥ek ∥ ] + h EUk [∥vk + lk ∥ ], (93) where the last inequality uses Jensen’s inequality and steps similar to same as (90). Moreover, we 2 know that ∥vk + lk ∥2 ≤ 2∥vk ∥2 + 2∥lk ∥2 and E[∥vk ∥2 ] ≤ σt . Similar to the process in (90), we have E[∥lk ∥2 ] ≤ 4L21 E[∥ek ∥2 ] and thus EUk [∥ek+1 ∥2 ] ≤ ρEUk [∥ek ∥2 ] +

2h2 σ 2 , t

where ρ = c2 + 4cL1 h + 8L21 h2 . For c ∈ [0, 1), to have ρ < 1, we require 0 < h < e0 = 0, by induction we get 2h2 σ 2 2h2 σ 2 = . t(1 − ρ) t(1 − c2 − 4cL1 h − 8L21 h2 ) Hence, combining the derivations above, we get s µ 2h2 σ 2 3/2 sup ∥qk ∥ ≤ h L1 (n + 3) + hL1 2 2 t(1 − c − 4cL1 h − 8L21 h2 ) k EUk [∥ek ∥2 ] ≤

(94) √

2−c2 −c . Since 4L1

(95)

(96)

which completes the proof. Proof of Theorem 6. By Lemma 5, the averaged ZO-GD dynamics (34) take the form of the perturbed system (4) with perturbation satisfying the bound (36). It remains to verify that the unperturbed regularised GD dynamics (32) satisfy the hypotheses of Theorem 2. By assumption, GD applied to the regularised problem converges to a fixed point z e . Since we have convergence for GD, the converse Lyapunov theorem for asymptotically stable discrete-time systems (see, e.g., [KB23, Thm. 5.5 and Thm. 5.6]) guarantees the existence of a Lyapunov function V satisfying (9)–(10) on some ball BR (z e ), R > 0. All hypotheses of Theorem 2 are therefore satisfied. The ISS bound (11) then guarantees that ZO-GD converges to a neighbourhood of z e with the
same transient decay profile as GD. The radius of this neighbourhood is governed by γ supk ∥qk ∥ , which, combined with the perturbation bound from Lemma 5, yields (37). Since the bound (36) can be made smaller than q̄ by choosing µ sufficiently small and t sufficiently large, such parameter choices always exist.

D

Complementary remarks on the results of Section 3

The following remarks provide additional technical details on the results presented in Section 3. Remark 2. While the first term in (19) can be manipulated by choosing µ small, the second term is governed by h and t. In the regime of large variance σ, for ZO-GD to operate with the largest admissible step size of GD, one must increase t, resulting in more function evaluations per iteration. When the optimal step size of GD is of order σ −2/3 , a small value of t suffices for ZO-GD
to track β GD closely with the same computational cost. Moreover, Theorem 3 requires h ∈ 0, 2L2 , whereas 1
+ 2 for strongly convex functions GD requires h ∈ 0, β+L [N 18]. Thus, if β ∼ L1 , the two step-size 1 bounds are of the same order. 22

Remark 3. We note that the expected performance of ZO-GD for strongly convex objective functions in minimisation problems has been analysed in [NS17, Thm. 8]. Compared to Theorem 3, the result in [NS17, Thm. 8] exhibits an explicit dependence of the iteration complexity on the dimension, and does not account for the variance bound σ or the number of sampled directions t. Moreover, the prescribed step size scales inversely with n, which is not consistent with standard GD. In contrast, Theorem 3 shows that ZO-GD, in expectation, achieves the same decay rate as GD and converges to a neighbourhood of the GD fixed point, whose size depends on σ and can be manipulated via h, µ, and t. Furthermore, for any step size consistent with GD, one can select µ and t such that the perturbation is arbitrarily small and closely track GD. Remark 4 (Relaxation of Assumption 1). Assumption 1 requires the variance of the oracle gµ defined in (12) to be uniformly bounded by a constant σ 2 . In some scenarios, the variance of the zeroth-order oracle scales with the magnitude of the gradient at the query point. Here, we show analogous results to Lemma 2 and Theorem 3 under the following relaxed assumption. Analogous results for other lemmas and theorems can be derived following similar steps. Assumption 2. The variance of the oracle gµ defined in (12) satisfies h i σ2 σ2 2 2 Eu ∥gµ (x) − ∇fµ (x)∥ ≤ 0 + 1 ∥∇f (x)∥ , (97) t t where σ0 , σ1 ≥ 0. Note that Assumption 1 is recovered by setting σ1 = 0 and σ0 = σ. Under this relaxed assumption, the proof of Lemma 2 proceeds as follows. Consider the recursion for EUk [∥ek+1 ∥2 ] derived in (52). With E[∥vk ∥2 ] ≤ σ 2 /t, equation (55) in the proof of Lemma 2 gives     2h2 σ 2 , EUk ∥ek+1 ∥2 ≤ ρ EUk ∥ek ∥2 + t where ρ = 1 − hβ + 2L21 h2 . Under (97), the additive term becomes state-dependent. Using the L1 -Lipschitz continuity of ∇f and ∇f (z ∗ ) = 0, we have 2

2

2

2

2

∥∇f (z̃k )∥ = ∥∇f (z̃k ) − ∇f (z ∗ )∥ ≤ L21 ∥z̃k − z ∗ ∥ ≤ 2L21 ∥ek ∥ + 2L21 ∥z̄k − z ∗ ∥ , where the last inequality uses ∥z̃k − z ∗ ∥2 ≤ 2∥ek ∥2 + 2∥z̄k − z ∗ ∥2 . Defining Σ1 := σ12 L21 , in place of (55) we obtain     2h2 σ02   4h2 Σ1 4h2 Σ1 2 2 E ∥ek ∥2 + + ∥z̄k − z ∗ ∥ . (98) EUk ∥ek+1 ∥ ≤ ρ + t t t Unlike the case of Assumption 1, bounding E[∥ek ∥2 ] alone is no longer sufficient, since (98) couples the error ek = z̃k − z̄k to the averaged trajectory z̄k . We therefore analyse the joint evolution of ak := E[∥ek ∥2 ] and bk := ∥z̄k − z ∗ ∥2 . Equation (98) gives ak+1 ≤ ρ′ ak +

2h2 σ02 4h2 Σ1 + bk , t t

(99)

where ρ′ := ρ + 4h2 Σ1 /t. For the bk recursion, recall the averaged dynamics (16): z̄k+1 = z̄k − h∇f (z̄k ) + qk . Expanding, we have 2

bk+1 = ∥z̄k − h∇f (z̄k ) + qk − z ∗ ∥ 2

= ∥z̄k − z ∗ − h∇f (z̄k )∥ + ∥qk ∥2 + 2 ⟨z̄k − z ∗ − h∇f (z̄k ), qk ⟩ p 2 ≤ (1 − c) bk + ∥qk ∥ , where c = 2hβL1 /(β + L1 ) ∈ (0, 1) is the GD contraction constant derived in [N+ 18, Thm 2.1.15]. Using Young’s inequality, for any ε > 0,
bk+1 ≤ (1 + ε)(1 − c) bk + 1 + 1ε ∥qk ∥2 . (100) From (46) and Jensen’s inequality, ∥qk ∥2 ≤ 2 bias2 + 2h2 L21 ak , where bias := hµL1 (n + 3)3/2 /2 denotes the smoothing bias term. Substituting into (100) yields

bk+1 ≤ (1 + ε)(1 − c) bk + 2 1 + 1ε h2 L21 ak + 2 1 + 1ε bias2 . (101) 23

Combining (99) and (101) and considering their worst-case setting, i.e., equality, the coupled system can be written in matrix form as       2h2 σ02 /t ak+1 ak =M + , (102) bk+1 bk 2(1 + 1/ε) bias2 where

 4h2 Σ1 ′ ρ   t .  M = (103)    1 2 2 2 1+ h L1 (1 + ε)(1 − c) ε If spectral radius of M is less than 1, the coupled system (102) converges to a neighbourhood of the origin determined by (I − M )−1 applied to the constant vector, and both supk ak and supk bk are uniformly bounded, which means that supk ∥qk ∥ is uniformly bounded. A sufficient condition for spectral radius of M being less than 1 is that both diagonal entries are strictly less than 1 and 

4h2 Σ1
(1 − ρ′ ) 1 − (1 + ε)(1 − c) > · 2 1 + 1ε h2 L21 , (104) t This follows from standard bounds on the spectral radius of a nonnegative matrix [HJ12, Ch 1, 5, and 8]. The condition ρ′ < 1 is satisfiable for any Σ1 > 0 by choosing t sufficiently large or 1 h sufficiently small. The condition (1 + ε)(1 − c) < 1 is satisfied for ε < 1−c − 1. The cross 4 2 condition (104) has a right-hand side of order O(h Σ1 L1 /t) versus a left-hand side of order O(h2 ), which is satisfiable for h small enough or t large enough. Once uniform boundedness of ak and bk is established, the perturbation bound takes the form √ sup ∥qk ∥ ≤ bias + hL1 sup ak , k

k

which is finite and controllable through the parameters h, µ, and t. The remainder of the ISS argument (Theorem 3) then applies without modification, yielding the same exponential decay rate and convergence to a neighbourhood whose radius depends on the perturbation bound. The key difference compared to Assumption 1 is that t must now be chosen large enough or h be chosen small enough to ensure that the spectral radius of M is less than 1, introducing a coupling between the variance parameter Σ1 = σ12 L21 and the GD contraction rate c. Remark 5. Setting h2 = 0 reduces the heavy ball method to gradient descent. In this case, ρHB reduces to ρGD = 1 − h1 β + 2h21 L21 , and the bound recovers the result of the GD case in Section 3.1. 3/2

Remark 6. The bound (26) consists of two terms: the first term, h1 µ L1 (n+3) , is the smoothing 2 bias and which depends on the smoothing parameter µ; the second term captures the variance of the ZO oracle amplified through the HB dynamics and can be manipulated through the number of samples t, the step size h1 , and the momentum parameter h2 (through ρHB ). Remark 7. Setting h2 = 0 reduces NAG to gradient descent. In this case, ρNAG reduces to ρGD = 1 − h1 β + 2h21 L21 , and the bound recovers the result of the GD case in Section 3.1. Remark 8. The bound (28) has the same two-term structure as the HB bound (26): the first term is the smoothing bias (depending on µ), and the second term captures the ZO oracle variance amplified through the dynamics (and depending on t, h1 , and h2 through ρNAG ). Remark 9. If, in the proof of Lemma 5, one bounds E[∥ek ∥] instead of E[∥ek ∥2 ], the resulting  1−c requirement becomes h ∈ 0, 2L1 , which implies λ > L1 and hence guarantees that the regularised function is L1 -strongly convex; in that case, the results of Section 3.1 apply directly. The bound obtained in Lemma 5 is considerably more general, as it covers regimes where the regularised objective is not strongly convex. Remark 10. Following a similar argument to the proof of Lemma 5, one can show that for h ∈ √ 2(1+c2 )−2c
0, , the same bound as in (36) holds but with 2L21 h2 in place of 8L21 h2 . This tighter 2L1 bound permits a nearly twice as large admissible range for h; however, the corresponding choice of λ satisfies λ > L21 , which is precisely the threshold at which the regularised objective becomes strongly convex. Consequently, λ cannot be made arbitrarily small, and the minimiser of the regularised problem cannot be brought arbitrarily close to that of the original problem. This reflects the wellknown trade-off between performance and robustness. 24

Remark 11. The same approach of Section 3.3 applies to the regularised problem (31) with HB or NAG: one can show that their ZO counterparts are perturbed versions of the corresponding FO algorithms with bounded and controllable perturbations, yielding results analogous to Theorem 6. The details are omitted for brevity.

E

Numerical examples complementary details

In this section, we present complementary experiments and results for the numerical examples given in Section 4. All the tests have been run on a Dell Latitude 7430 Laptop with a 12th Gen Intel Core i7 CPU. E.1

Additional experiments for the quadratic objective

This section complements the quadratic experiment in Section 4.1 with additional parameter studies. Effect of the smoothing parameter µ. We fix t = 1 and h = 10−5 and vary µ ∈ {10−1 , 4 × 10−3 , 10−6 , 10−8 }. The results are shown in Figure 4. For large µ, ZO-GD fails to converge; reducing µ restores convergence. However, once µ is sufficiently small, further reduction has negligible effect, since the dominant term in the perturbation bound no longer depends on µ. This is in contrast to h, which controls both terms in the perturbation bound. GD ZO-GD = 1e 01

101

x x

x x

GD ZO-GD = 1e 01

0

200

400

Iteration

600

800

1000

0.00

0.25

0.50

0.75

1.00

Wall-clock time (s)

1.25

2 × 100

1.50

1.75

GD ZO-GD = 4e 03 ZO-GD = 1e 06 ZO-GD = 1e 08

3 × 100

x x

3 × 100

x x

101

2 × 100

GD ZO-GD = 4e 03 ZO-GD = 1e 06 ZO-GD = 1e 08 0

200

400

Iteration

600

800

1000

0.00

0.25

0.50

0.75

1.00

Wall-clock time (s)

1.25

1.50

1.75

Figure 4: Parameter analysis: smoothing parameter. Effect of the number of sampled directions t. We fix µ = h = 10−5 and vary t ∈ {1, 5, 10, 20}. The results are shown in Figure 5. Increasing t reduces the perturbation, and ZO-GD converges to a neighbourhood of GD that shrinks with t. Parameter compensation. We revisit the cases in which GD converged, but ZO-GD diverged, demonstrating that the perturbation can always be controlled by adjusting the remaining parameters. First, consider the divergent case h = 10−4 , µ = 10−5 , t = 1 from Figure 1. As shown in the first row of Figure 6, increasing the number of sampled directions to t = 25 and reducing µ = 10−7 brings the perturbation within the admissible range. Second, consider the divergent case µ = 0.01, t = 1, h = 10−5 from Figure 4. Since µ and t do not interact directly in the perturbation bound, we instead reduce h to 10−7 . The result, shown in the second row of Figure 6, confirms that ZO-GD again tracks GD. 25

GD ZO-GD t = 1 ZO-GD t = 5 ZO-GD t = 10 ZO-GD t = 20

3 × 100

x x

x x

3 × 100

2 × 100

GD ZO-GD t = 1 ZO-GD t = 5 ZO-GD t = 10 ZO-GD t = 20 0

200

400

Iteration

600

800

2 × 100

1000

0

2

4

6

8

Wall-clock time (s)

10

Figure 5: Parameter analysis: number of sampled directions. GD h = 1e 04 Set 1: h = 10 4, = 10 7, t = 25

3 × 100

2 × 100

x x

x x

2 × 100

100

6 × 10 1

100

6 × 10 1

4 × 10 1

4 × 10 1 0

200

400

Iteration

3.2 × 100

600

800

1000

0

2.6 × 100

x x

2.8 × 100

2.6 × 100 2.4 × 100

2.4 × 100

2.2 × 100

2.2 × 100

10000

20000

30000

Iteration

40000

4

6

50000

8

Wall-clock time (s)

10

GD h = 1e 07 Set 2: h = 10 7, = 10 2, t = 1

3 × 100

2.8 × 100

0

2

3.2 × 100

GD h = 1e 07 Set 2: h = 10 7, = 10 2, t = 1

3 × 100

x x

GD h = 1e 04 Set 1: h = 10 4, = 10 7, t = 25

3 × 100

0

20

40

Wall-clock time (s)

60

80

Figure 6: Parameter analysis: compensation for divergent cases. Limitations of the number of sampled directions compensation. To demonstrate that increasing t alone cannot fully compensate for a large µ, we consider three configurations with µ = 5 × 10−1 fixed: (i) h = 10−5 , t = 1; (ii) h = 10−7 , t = 1; and (iii) h = 10−5 , t = 20. The results are shown in Figure 7. In case (i), ZO-GD diverges while GD converges. Reducing h in case (ii) shrinks the perturbation and restores convergence. Increasing t in case (iii) also reduces the perturbation, but since t does not affect the smoothing bias term in Lemma 2, it cannot fully compensate for a large µ. E.2

Additional experiments for neural network classification

This section complements the neural network experiment in Section 4.2 by providing a parameter sensitivity analysis. Parameter sensitivity. We conduct parameter sweeps for ZO-GD analogous to those in Appendix E.1. Figure 8 shows the effect of varying h, µ, and t individually while holding the other two fixed, for the network with 100, 609 parameters. The same phenomena observed for the quadratic objective are reproduced: when h is too large, ZO-GD diverges, and reducing h brings ZO-GD progressively closer to GD (Figure 8, left). Reducing µ restores convergence when it is initially too large, but further reduction has diminishing effect once the variance term dominates the perturbation bound (Figure 8, middle). Increasing the number of sampled directions t shrinks the neighbourhood, with ZO-GD converging closer to GD for larger t (Figure 8, right). These observations are consistent with the two-term structure of the perturbation bound in Lemma 5: h and µ control the smoothing bias term, while h and t control the variance term.

26

101

0

200

400

Iteration

600

800

x x

x x

GD h = 1e 05 Set 1: h = 10 5, = 3 × 10 2, t = 1

1000

0.00

GD h = 1e 07 Set 2: h = 10 7, = 3 × 10 2, t = 1

3.2 × 100 3 × 100

3 × 100

2.8 × 100

2.8 × 100

2.6 × 100

2.4 × 100

2.2 × 100

2.2 × 100 10000

20000

30000

40000

Iteration

0.50

0.75

1.00

1.25

Wall-clock time (s)

1.50

1.75

GD h = 1e 07 Set 2: h = 10 7, = 3 × 10 2, t = 1

2.6 × 100

2.4 × 100

0

0.25

3.2 × 100

x x

x x

GD h = 1e 05 Set 1: h = 10 5, = 3 × 10 2, t = 1

101

50000

0

40

60

80

Wall-clock time (s)

x x

3 × 100

x x

3 × 100

20

2 × 100

2 × 100

GD h = 1e 05 Set 3: h = 10 5, = 3 × 10 2, t = 20 0

200

400

Iteration

GD h = 1e 05 Set 3: h = 10 5, = 3 × 10 2, t = 20

600

800

1000

0

1

2

3

4

Wall-clock time (s)

5

6

7

Figure 7: Parameter analysis: limitations of the number of sampled directions compensation.

Effect of step size h ( =1e-07, t=5)

7.95 × 10 1

Effect of smoothing (h=0.001, t=5)

7.9 × 10 1

Training Loss

7.8 × 10 1 7.75 × 10 1 7.7 × 10 1

7.9 × 10 1

100 9 × 10 1

7.65 × 10 1

GD ( =0.002) ZO-GD ( =0.002) GD ( =0.001) ZO-GD ( =0.001)

7.6 × 10 1 0

20

60

Iteration k

7.875 × 10 1 7.85 × 10 1 7.825 × 10 1 GD (FO) ZO-GD (t=1) ZO-GD (t=3) ZO-GD (t=5) ZO-GD (t=10)

7.8 × 10 1

GD ( =0.0005) ZO-GD ( =0.0005) GD ( =0.0001) ZO-GD ( =0.0001)

40

Effect of mini-batch size t (h=0.001, =1e-07)

7.925 × 10 1

Training Loss

7.85 × 10 1

Training Loss

7.95 × 10 1

GD (FO) ZO-GD ( =5e-03) ZO-GD ( =1e-03) ZO-GD ( =1e-06) ZO-GD ( =1e-07)

8 × 10 1 80

100

7.775 × 10 1 0

20

40

60

Iteration k

80

100

0

20

40

60

Iteration k

80

100

Figure 8: Parameter sensitivity for ZO-GD on the neural network. Left: effect of h. Centre: effect of µ. Right: effect of t. Solid black: GD (FO); dashed coloured: ZO-GD.

27