Trust the Direction, Search the Step: Zero-and-First-Order Methods for LLM Fine-Tuning

Trust the Direction, Search the Step: Zero-and-First-Order Methods for LLM Fine-Tuning

arXiv:2610.02190v1 [cs.LG] 1 Oct 2026

Cristian McGee University of Central Florida Orlando, FL, USA

El Houcine Bergou Mohammed VI Polytechnic University Ben Guerir, Morocco

Aritra Dutta University of Central Florida Orlando, FL, USA

Abstract Step-size selection remains a central challenge in large-scale neural network optimization; conservative steps slow convergence, while aggressive steps can destabilize it. We combine Zero-and-First-Order optimization (ZFO) and propose a framework that decouples direction selection from step-size. ZFO uses a trusted first-order optimizer to determine the direction and performs zeroth-order evaluations only along this one-dimensional subspace to choose how far to move. Using the current gradient information and two additional objective function evaluations, ZFO instances construct a local model of the objective function along the proposed direction and select a curvature-aware step within a bounded search interval. This yields an adaptive step-selection mechanism that costs less than a full line search. We provide theoretical guarantees to show that shared-sample evaluations produce reliable finite-difference curvature estimates, that the induced local model selects a near-optimal step along the search interval, and that ZFO converges to a neighborhood of a stationary point. Across the evaluated settings, language models and datasets, ZFO frequently improves optimization and final performance relative to fixed-step first-order baselines, with the magnitude and preferred local model depending on the objective. Our code is publicly available at: https://github.com/nizswan/Zeroth-First-Order-Framework.

1

Introduction

Gradient-based first-order (FO) optimization methods are the main workhorse of modern deep neural network training [14, 25, 31]. During training, gradient information provides the direction to move along, which, coupled with a step size (or learning rate), updates the model parameters. While the gradient direction is generally informative, the step size used along that direction is typically chosen through schedules, heuristics, or hyperparameter tuning [18, 30, 48, 49]. For larger networks, whose local geometry can change rapidly during adaptation, a poorly chosen step size may lead to unstable training [13]. Despite these step-size instabilities, FO methods remain the standard choice for neural network adaptation because the FO gradient provides strong directional guidance. In contrast, zeroth-order (ZO) methods [5–7, 16, 17, 26, 34, 38] do not compute gradients. They construct approximate gradient updates using only objective function evaluations, often via finite differences along random directions. Although ZO methods reduce computational overhead, gradient approximations are often noisy and high-variance in large-dimensional settings [2, 23]. However, when FO provides the update direction, the high-dimensional problem becomes a univariate step40th Conference on Neural Information Processing Systems (NeurIPS 2026).

(a) First-Order Direction

(b) Zeroth-Order Step Selection

zt = gt =kgt k

Át (s) = Ft (µt + szt )

Ft (µ) Objective

Át (§²)

µt µt + ®gt

I t = [0; Rt ]

ZFO

FO

Surrogate Át (s) Local model M(s)

µt + 1

zt ¡²

0

+² sFO

sM

Rt

Step length s

Given: Át (0); Át0 (0) Extra Evaluations: Át (¡²); Át (+²) d 00 000 d Át (0); Át (0) ⟶ sM 2 arg max M(s) ⟶ µt + 1 = µt + sM zt s2I t

Figure 1: The Zeroth First -Order (ZFO) Framework. Left: A first-order optimizer proposes the direction zt at θt , which restricts the objective to the line θt + szt . Right: Along this line, θt + szt , the forward pass gives ϕt (0), the backward pass gives ϕ′t (0), and two symmetric probes at ±ϵ supply curvature estimates. ZFO fits a local model M (s) and takes the optimized step sM over the interval [0, Rt ]; see Algorithm 1.

selection problem; see Figure 1. Consequently, computationally inexpensive ZO evaluations become substantially informative for modeling local function behavior. This raises a natural question: If FO methods are strong at choosing an update direction, can ZO evaluations be used to decide how far to move along that direction? We investigate this question by using ZO evaluations not to replace gradients, but to select reliable step sizes. This idea is related to recent work on generalized Newton step selection GeN [9], which uses function evaluations along a first-order direction to construct a curvature-informed step. In contrast, ZFO takes a broader view: rather than committing to a particular step construction, we formulate fixed-budget directional step selection as a local modeling problem. This perspective is not restricted to a specific optimization setup; in principle, the framework applies whenever an approximate gradient direction is available, but its appropriate magnitude is uncertain, a pattern that recurs in existing step-size rules [33]. This uncertainty becomes especially relevant in noisy or stochastic optimization settings. Large language model (LLM) fine-tuning is a natural setting for this hybridization. Modern fine-tuning pipelines construct minibatches of prompt-response examples and optimize a surrogate objective on these fixed examples. Once a minibatch is formed, nearby parameter evaluations can reuse the same data and require only additional forward evaluations of the perturbed model. This makes ZO probes especially useful—not as a replacement for gradients, but as a cheap mechanism for modeling objective variation along a trusted FO direction. In contrast, fully ZO LLM fine-tuning methods use forward evaluations to estimate high-dimensional gradients [26, 34]. ZFO exploits this structure by assigning direction selection to a base FO optimizer, such as AdamW [31] or Muon [25], while using a small number of ZO function evaluations to build a univariate local model of the objective function and select the step size along that direction. Taken together, the ZFO framework acts as a lightweight step-selection wrapper around a standard FO optimizer. This is particularly valuable in LLM fine-tuning, where local geometry can shift rapidly and poorly tuned learning rates can cause unstable or slow training [27, 50]. Problem Setup. Let fine-tuning an LLM be represented by a conditional distribution πθ (y | x), where x denotes an input prompt, y denotes a model response, and θ ∈ Rd are the model parameters. We initialize the model from a pretrained reference model and further optimize it using task-specific feedback. We summarize the main notation used throughout the paper in §D; see Table 14. At a given iteration t, a FO fine-tuning algorithm constructs a minibatch of information, denoted by St . This minibatch contains all quantities needed to define the current training objective. Depending on the method, St may include prompts, sampled responses, supervised labels, etc. Conditioned on St , the fine-tuning method defines a sampled surrogate objective optimized during the current update 2

Ft (θ) = F (θ; St ). Different LLM fine-tuning algorithms instantiate this surrogate differently. Our framework does not require a particular choice of Ft , and only assumes that the surrogate can be evaluated on the fixed minibatch St . This fixed-batch viewpoint is standard in LLM fine-tuning [22, 26, 34, 46]. Once St is constructed, it is held fixed within the update. Evaluating Ft at a nearby parameter value therefore does not require constructing a new batch or collecting new feedback; it only requires re-scoring the same prompt-response pairs under the perturbed model. A standard fine-tuning step computes a stochastic gradient of the sampled surrogate at the current parameters, and passes it through an optimizer such as AdamW. If we denote the resulting update direction by gt , the usual first-order gradient ascent update is: θt+1 = θt + αgt , where α > 0 is the learning rate. This update contains two decisions: the direction gt and the step length controlled by α. The chosen fine-tuning surrogate and the optimizer determine the direction. The step length, however, is usually specified by a learning-rate schedule or hyperparameter tuning [32]. In LLM fine-tuning, this scalar choice can be sensitive; overly tiny steps may not entirely use the information in the current batch, while aggressive steps can destabilize training [51]. We improve the step-length decision while leaving the optimizer direction unchanged and formulate all objectives as maximization problems. Let gt ̸= 0, and define the normalized direction gt zt = . (1) ∥gt ∥ We restrict the current surrogate to the one-dimensional path along this direction: ϕt (s) = Ft (θt + szt ), for s ∈ R. The standard optimizer step corresponds to sFO = α∥gt ∥. Thus, choosing the learning rate is equivalent to choosing a scalar displacement along the fixed-batch improvement of ϕt . Rather than relying only on a global learning-rate schedule, we ask whether the current surrogate Ft supports a better step along the direction already selected by the optimizer. We search over a bounded interval It = [0, Rt ], Rt = βα∥gt ∥, where β > 0 controls the search radius relative to the baseline optimizer step. When β = 1, the interval ends at the usual first-order update; larger values allow the method to consider longer steps while remaining constrained to the same direction: s⋆t ≈ arg max ϕt (s). (2) s∈It

A full line search would evaluate ϕt at many candidate values of s, which is costly for LLMs. Instead, we use a fixed-budget step-selection procedure: the current value ϕt (0) is obtained from the usual forward evaluation of the surrogate on St , and the backward pass gives the directional derivative ϕ′t (0) = ⟨∇Ft (θt ), zt ⟩ . We then evaluate the same fixed-batch surrogate at two symmetric perturbations along the update direction: ϕt (ϵ) = Ft (θt + ϵzt ), ϕt (−ϵ) = Ft (θt − ϵzt ), for a small ϵ > 0. Together, these quantities {ϕt (0), ϕ′t (0), ϕt (ϵ), ϕt (−ϵ)} , define a local model of the surrogate improvement along the optimizer direction. The selected step st is obtained by maximizing this model over It , and the parameters are updated as θt+1 = θt + st zt . In this formulation, the base LLM fine-tuning algorithm determines the surrogate Ft and the update direction gt , while our method determines how far to move along that direction; see Figure 1. The procedure is therefore compatible with a broad class of LLM fine-tuning objectives. Due to space limitations, we defer the related work to §A. Our contributions. We develop ZFO as a general framework that separates first-order direction selection from low-query zeroth-order step selection. This formulation allows different local stepselection models to be developed and analyzed within the same optimization framework. Below, we summarize our main contributions: (i) A ZFO optimization framework. We propose the general ZFO framework in which a standard first-order optimizer proposes an update direction, and a small number of zeroth-order evaluations are used to choose an improved step size along this direction. Within this framework, we develop four practical ZFO instances based on Taylor and Padé local models. While existing methods such as GeN is a one step-selection rule, our ZFO is a fixed-budget directional modeling framework. GeN can be interpreted within the ZFO framework as a particular fixed-budget directional step-selection.

3

(ii) Theoretical justification of adaptive step selection. Our theoretical analysis of the onedimensional objective follows from restricting the training objective to the first-order search direction. We show that the finite-difference curvature estimates with common random numbers (CRNs) concentrate around their population counterparts with high probability. We further prove that maximizing the resulting local model yields a step whose objective value is near-optimal within a bounded search region. Finally, we establish a global stationarity guarantee showing that ZFO converges to a stationary neighborhood whose size is controlled by the accumulated step-selection error. (iii) Empirical validation across scales. We evaluate ZFO on both controlled low-dimensional problems and LLM fine-tuning tasks. The experiments show that ZFO is robust to the choice of initial learning rate and can match or outperform standard first-order and zeroth-order baselines on most benchmarks. These results demonstrate that lightweight zeroth-order step adaptation can improve the stability and reliability of first-order fine-tuning without incurring the cost of full line search or extensive hyperparameter tuning.

2

Methodology

We fix an iteration t and suppress the subscript t for readability. We use this convention throughout the section. Local derivative estimates. We recall that the forward-backward pass provides ϕ(0) and the directional derivative, ϕ′ (0) = ⟨∇F (θ), z⟩. ZFO adds two symmetric fixed-batch evaluations, ϕ(ϵ) = F (θ + ϵz), and ϕ(−ϵ) = F (θ − ϵz). We estimate the second directional derivative by the central finite difference c′′ (0) = ϕ(ϵ) − 2ϕ(0) + ϕ(−ϵ) . ϕ ϵ2 Using the available first-order derivative, the same two probes also yield 3 ′′′ (0) = ϕc (ϕ(ϵ) − ϕ(−ϵ) − 2ϵϕ′ (0)) . ϵ3 ′′′ (0). ZFO uses these quantities to c′′ (0), and c = ϕc For compactness, define: a = ϕ′ (0), b = ϕ construct a local approximation of ϕ (see Figure 1), considering both polynomial models derived from Taylor expansions and rational models derived from Padé approximations. (a) Taylor models. The second-order Taylor model is T2 (s) = ϕ(0) + as + 12 bs2 . Since T2 is quadratic, the maximizer is obtained by checking the endpoints of I and, when feasible, the stationary point s = −a/b. Similarly, the third-order Taylor model is T3 (s) = ϕ(0) + as + 21 bs2 + 16 cs3 . The candidate maximizers are the endpoints of I and the real roots inside I of a + bs + 12 cs2 = 0. (b) Padé models. The Padé variants approximate ϕ using rational functions. The second-order Padé model is as b P2 (s) = ϕ(0) + , q2 = − , 1 + q2 s 2a Similarly, the third-order Padé model is as + ρs2 c b P3 (s) = ϕ(0) + , where q3 = − , ρ = + aq3 . 1 + q3 s 3b 2 If a Padé denominator has a pole inside I, or if its coefficients are ill-defined, ZFO rejects the Padé step and falls back to the quadratic Taylor step. This safeguard avoids unstable rational extrapolation near singularities. Model selection and update. Each ZFO variant corresponds to one model choice from {T2 , T3 , P2 , P3 }. For a chosen model M in this set, ZFO computes the step by solving sM ∈ arg max M (s). s∈I

(3)

All ZFO variants use the same direction, surrogate, and feasible interval. They differ only in how they convert two additional evaluations into a scalar step length. GeN [9] also fits within this abstraction as a specific quadratic ZFO instance. GeN uses two function evaluations along the first-order direction to construct a generalized Newton step. Our instances reuse the directional derivative already available from the backward pass so that the same two probes recover both second- and third-order directional 4

Algorithm 1 Zeroth First -Order (ZFO) Framework Require: θ0 , base optimizer, learning rate α, ϵ > 0, β > 0, model M ∈ {T2 , T3 , P2 , P3 } 1: for t = 0, 1, 2, . . . do 2: Compute the direction gt using the base optimizer ▷ First-order direction 3: Set zt = gt /∥gt ∥ and It = [0, βα∥gt ∥] t (−ϵ) , ct = ϵ33 (ϕt (ϵ) − ϕt (−ϵ) − 2ϵat ) 4: Compute: at = ϕ′t (0), bt = ϕt (ϵ)−2ϕtϵ(0)+ϕ 2 5: Construct the chosen local model M using at , bt , ct 6: Select st = sM ∈ arg max M (s), s.t. s ∈ It ▷ Zeroth-order step selection 7: if M is a Padé model and has a pole in It then 8: Replace M by T2 , recompute sM , and set st = sM 9: end if 10: Update θt+1 = θt + st zt ▷ Parameter update 11: end for

information and support Taylor and Padé models. GeN uses curvature-based rejection as its safeguard, while the instances studied here optimize over a bounded interval and fall back when a Padé model is ill-defined. ZFO’s novelty is the fixed-budget local-model abstraction, its family of Taylor/Padé instantiations, and the associated unified analysis and empirical validation. Algorithm 1 summarizes our framework.

3

Theoretical Analysis

This section provides a theoretical justification for the ZFO step-selection framework. The analysis has three components. First, we show that shared-batch evaluations yield reliable finite-difference estimates of directional curvature; see §3.1 (i). Second, conditional on the sampled examples, we show that maximizing the local ZFO model returns a step whose objective value is close to the best value attainable over the search interval; see §3.1 (ii). Third, we connect this local step-selection property to a global stationarity guarantee for the ZFO update sequence; see §3.2. We omit the iteration index t when no ambiguity arises. 3.1

Local Convergence Analysis

Shared-batch common-random-numbers (CRN). ZFO evaluates all local objective values on the same minibatch. At a fixed iteration, the examples are sampled once and then reused across the evaluations {ϕ(−ϵ), ϕ(0), ϕ(ϵ)} . After conditioning on these shared examples, ϕ becomes a deterministic function of the scalar step length, s. This conditional determinism is the basis of the local approximation analysis. Let τ = {τ1 , . . . , τn } denote the shared sampled examples used under CRN. For each example τi , define the per-sample directional objective ϕi (s) = F (θ + sz; τi ). The Pn empirical shared-batch CRN objective is ϕ(s) = n1 i=1 ϕi (s). Definition 1 (Population directional objective). The population counterpart of the empirical CRN objective is ψ(s) = Eτ [ϕi (s)], where the expectation is taken with respect to the distribution used to generate the empirical objective. (i) Local smoothness assumptions. The first group of assumptions ensures that the directional objective admits Taylor expansions with controlled remainders. Assumption 1 (Pathwise smoothness). For each example τi , the per-sample directional objective ϕi is five times continuously differentiable on the interval [−ϵ, R]: ϕi ∈ C 5 ([−ϵ, R]). Assumption 2 (Higher-order derivative bounds). There exist constants W, M, B, L > 0 such that, (4) (5) for every i and every s ∈ [−ϵ, R], |ϕ′′i (s)| ≤ W, |ϕ′′′ i (s)| ≤ M, |ϕi (s)| ≤ B, and |ϕi (s)| ≤ L. (ii) Sampling assumptions. Population concentration is analyzed at a fixed iterate and direction under the independent-example assumptions below. Algorithm 1 reuses its batch across probes, but also computes the direction from that batch; CRN alone does not establish the independence needed by this concentration argument. Thus Theorem 1 is a fixed-direction result, whereas the approximation guarantee in Theorem 2 is pathwise once the batch and direction are fixed. Bounded directional derivatives imply the sub-Gaussian condition by Hoeffding’s lemma; see (4). 5

Assumption 3 (CRN). For each index i, the same example τi is used to evaluate {ϕi (−ϵ), ϕi (0), ϕi (ϵ)}. Moreover, for i ̸= j, the examples τi and τj are sampled independently. Assumption 4. (Sub-Gaussian Deviation) For each ϕi , the directional deviations at s = 0 are sub-gaussian up to third order, i.e., there exists constant σ such that, for every k ∈ {0, 1, 2, 3} and 2 2 (k) all λ ∈ R, E[exp(λ[ϕi (0) − ψ (k) (0)])] ≤ exp( λ 2σ ). Theorem 1 (CRN curvature concentration). at q Suppose Assumptions 1–4 hold. Then, with probability q 2 log(4/δ) 2 log(4/δ) Bϵ2 Lϵ2 ′′ ′′′ ′′ ′′′ c c least 1 − δ, ϕ (0) − ψ (0) ≤ +σ , and ϕ (0) − ψ (0) ≤ +σ . 12

n

20

n

Theorem 1 shows that the finite-difference curvature estimates used by ZFO concentrate around their population counterparts. Each bound contains two terms. The first term is a deterministic finitedifference bias controlled by the perturbation scale ϵ and the relevant higher-order √ derivative bound. The second term is a stochastic sampling term that decreases at the classic rate 1/ n. Thus, under CRN, the two additional objective evaluations provide statistically meaningful curvature information. After conditioning on the shared examples, the empirical objective ϕ is deterministic. Let sM ∈ arg maxs∈[0,R] T2 (s) denote the model-selected step by T2 , and let sϕ ∈ arg maxs∈[0,R] ϕ(s) denote the best step for the empirical CRN objective over the same interval. Theorem 2 (Second-order Taylor approximation bound). Suppose Assumptions 1 and 2 hold. Then the 3 2 2 step selected by maximizing the second-order Taylor model satisfies ϕ(sM ) ≥ ϕ(sϕ )− M3R − Bϵ12R . Theorem 2 shows that, once the examples are fixed, ZFO selects a step that is close to the best step 3 available in the search interval. The first error term, M3R , is the deterministic Taylor approximation 2 2 error over the interval. The second error term, Bϵ12R , is the error caused by estimating curvature through finite differences. The guarantee is strongest when the search radius R is moderate, the perturbation scale ϵ is small, and the directional objective has limited higher-order variation. The same proof strategy extends to the third-order Taylor model. For the Padé variants, the local model is rational rather than polynomial. These models can better capture nonlinear curvature changes over a finite interval. Their approximation guarantees follow from the same maximizer-comparison argument, replacing the Taylor remainder by a rational residual term. We defer the analysis of these variants to §B.1. 3.2

Global convergence Analysis

We now connect the local analysis to the true optimization objective J. Recall that ZFO optimizes on the surrogate Ft , whose expectation under the current sampling distribution is F̄t (θ + szt ) = ψ(s). For instance, standard policy optimization [53] and its modern extensions [45, 46] construct surrogates such that ∇J(θt ) = ∇F̄ (θt ) at the current-iterate. Define the expected directional alignment quantity qt = E[⟨∇J(θt ), zt ⟩]. We start by stating two main classical assumptions that we use in the sequel. Assumption 5 (Smooth and bounded objective). The true objective J : Rd → R is LJ -smooth and bounded above by J ⋆ . That is, for all x, y, J(y) ≥ J(x) + ⟨∇J(x), y − x⟩ − L2J ∥y − x∥2 , and J(θt ) ≤ J ⋆ for all iterations t. Additionally, the expected surrogate F̄ is also LF̄ -smooth. Assumption 6 (Aligned update direction). There exists µ > 0 such that, for all t, qt = Et [⟨∇J(θt ), zt ⟩] ≥ µ∥∇J(θt )∥. Assumption 6 is natural since gradient methods estimate an ascent direction for the objective [46]. Although optimizers such as AdamW rescale this direction, the resulting update remains positively oriented with ∇J(θt ) in expectation [3], supporting the assumption of µ > 0. We start with a Lemma that bounds the step-selection error at iteration t of the objective: maxs∈It J(θt + szt ) − J(θt + st zt ). This quantity measures the loss in objective value by using the ZFO-selected step st instead of the best step available in the interval It for the true objective J. Lemma 1 (Surrogate vs. Empirical Difference). Let Assumptions 1-5 hold, then maxs∈It J(θt + szt ) − J(θt + st zt ) ≤ ∆t , holds with probability r at least 1 − δ, where 2 log(4/δ) ∆t := (LJ + LF̄ )Rt2 + 2σ(1 + Rt ) + 2W Rt2 + ϕ(sϕ ) − ϕ(sM ). n 6

Table 1: Evaluation performance for baseline first-order (FO) fine-tuning (using AdamW) and zeroth-order (ZO) fine-tuning (using MeZO) against ZFO methods across models and reasoning benchmarks. All runs report the final evaluation and mean results over 3 seeds. Green is the best and red is the second best result. Model

Dataset

Baseline

ZFO Methods (This Paper)

AdamW (FO) MeZO (ZO)

Qwen-2.5-Math-1.5B

Taylor2

Taylor3

Padé2

Padé3 80.97 ± 0.73

GSM8K

78.24 ± 0.59

0.00

81.00 ± 0.11

75.56 ± 0.24

80.95 ± 0.49

MATH

31.77 ± 1.63

22.01 ± 3.93

40.10 ± 1.19

20.18 ± 3.03

37.50 ± 1.79

39.58 ± 0.45

SVAMP

88.11 ± 0.69

68.67 ± 10.33

90.78 ± 0.96

89.11 ± 0.51

90.89 ± 1.02

90.33 ± 2.19

AsDiv

91.45 ± 0.36

64.72 ± 3.51

94.09 ± 0.29

88.27 ± 0.35

93.04 ± 0.22

93.46 ± 0.18

OpenBookQA

26.07 ± 2.53

8.00 ± 10.05

37.80 ± 17.85

39.73 ± 20.99

60.67 ± 0.95

64.07 ± 1.22

SVAMP

55.44 ± 4.91

22.44 ± 3.89

55.67 ± 7.86

39.33 ± 0.88

64.22 ± 0.19

65.11 ± 5.39

AsDiv

56.51 ± 18.76

10.15 ± 4.88

65.42 ± 2.59

55.85 ± 2.55

65.58 ± 3.74

69.27 ± 1.92

OpenBookQA

51.20 ± 1.11

20.33 ± 10.90

58.27 ± 13.62

51.27 ± 0.95

49.40 ± 0.20

52.07 ± 0.23

Gemma-2-2B

SVAMP

52.00 ± 1.76

14.44 ± 5.93

59.56 ± 0.51

47.22 ± 2.01

59.00 ± 1.76

58.22 ± 3.69

GSM8K

4.85 ± 5.47

1.30 ± 0.98

4.02 ± 3.98

9.00 ± 6.68

5.05 ± 6.94

1.72 ± 0.88

Llama-3.2-1B

AsDiv

16.60 ± 28.12

3.33 ± 0.48

13.72 ± 23.17

44.77 ± 1.42

25.23 ± 21.97

52.16 ± 4.75

OpenBookQA

25.53 ± 0.31

0.00

34.40 ± 1.56

31.53 ± 0.81

34.60 ± 1.04

25.64 ± 0.23

Phi-2

This lemma quantifies the loss in objective value when the ZFO step st , selected through the local model, is used instead of the best step in the interval It . The error term ∆t combines smoothness mismatch, finite-sample concentration error, approximation error, and the model-selection gap ϕ(sϕ ) − ϕ(sM ). E.g., when the local model is the second-order Taylor model T2 , using Theorem 2 we get r M R3 Bϵ2 R2 2 log(4/δ) 2 + 2W Rt2 + + . ∆t ≤ (LJ + LF̄ )Rt + 2σ(1 + Rt ) n 3 12 Finally, Theorem 3 converts the local quality of the ZFO step-selection rule into a global stationarity guarantee. Theorem 3 (Stationarity of the ZFO update). Let Assumptions 5–6 hold and qt /LJ ∈ It . Then, for every T > 0, with probability at least 1 − δ, we have T −1 X ˜T 2LJ (J ⋆ − J(θ0 )) 2LJ ∆ ˜T = 1 E[∆t ]. min E∥∇J(θt )∥2 ≤ + , where ∆ 0≤t≤T µ2 T µ2 T t=0 ⋆

0 )) The upper bound in Theorem 3 has two terms. The first term, 2LJ (Jµ2−J(θ , is the usual optimization T

˜

term and decays at rate 1/T . The second term, 2LµJ 2∆T , is the price paid for imperfect step selection. Thus, ZFO does not necessarily converge exactly to a stationary point unless the average step-selection error vanishes. Instead, it converges to a neighborhood of stationarity whose size is controlled by the ¯T. error ∆ The three results establish the theoretical foundation of ZFO. The CRN concentration theorem establishes the statistical reliability of the finite-difference curvature estimates; the deterministic approximation theorem shows that maximizing the ZFO model yields a near-optimal local step; and the global stationarity theorem shows that controlled step-selection error preserves convergence to stationarity. Together, these results justify the central mechanism of ZFO: a small number of fixed-batch objective evaluations can recover useful local curvature information and improve step selection over a fixed first-order update while retaining a principled convergence interpretation.

4

Experiments

Organization of the Results. We organize our experimental results as follows: In the main paper, (i) we compare ZFO against first-order fine-tuning with AdamW and zeroth-order fine-tuning with MeZO across several LLMs and datasets in Table 1; (ii) we conduct a focused Qwen-2.5-Math-1.5B sub-study over additional reasoning benchmarks in Table 2; (iii) we evaluate ZFO on classical optimization problems to assess its benefits beyond LLM fine-tuning in Figure 2; (iv) we summarize the ablation studies in the appendix; and (v) we report wall-clock and memory overheads in Fig7

Table 2: Qwen-2.5-Math-1.5B study on ZFO methods in our settings via additional reasoning benchmarks. Results report evaluation performance averaged over 3 seeds for FO and ZFO methods. Method

ARC-Challenge

StrategyQA

FOLIO

CODAH

AdamW (FO)

49.29 ± 16.08

68.56 ± 0.91

34.48 ± 1.13

57.91 ± 3.24

Taylor2 Taylor3 Padé2 Padé3

58.10 ± 1.33 54.24 ± 0.94 53.58 ± 2.77 57.05 ± 1.24

69.63 ± 1.82 61.67 ± 1.46 65.55 ± 4.43 70.06 ± 0.61

36.6 ± 3.26 34.48 ± 1.13 36.27 ± 2.73 36.44 ± 4.53

62.83 ± 3.14 55.28 ± 1.66 53.96 ± 17.72 41.85 ± 20.58

ures 3 and 4 to quantify the computational cost of ZFO. We defer the remaining experiments and expanded diagnostics to the Appendix C.2, where (vi) we provide the complete six-problem classical optimization study in Table 8 and Figure 5, (vii) we validate the local ZFO approximants on analytical functions in Figure 6 and the step-selection mechanism on controlled toy RLVR environments in Figure 7, (viii) we provide learning-rate sensitivity curves in Figure 8 and study the interaction between α and β in Figure 9, including the complementary β = 1 control in Table 12, (ix) we analyze the frequency of the steps selected by each ZFO method in Table 13, (x) we summarize cross-setting improvements in Table 9 and compare against additional first-order baselines (Prodigy [37] and AdamW [31]) in Table 10, and (xi) we evaluate the sequential ZFO-Seq against line-search baselines (Armijo [39], BB [4], Polyak [42], and PLS [33]) in Table 11. Unless otherwise stated, all reported results are averaged over three seeds. The LLM experiments use the implementation described in §C.1.1. Additional experimental configuration details are provided in §C.1.2. RLVR, LLMs, and Datasets used for Fine-tuning. We report the performance of various highdimensional models configured for an RLVR setup using the GRPO policy gradient [47] and the AdamW optimizer [31], including Qwen-2.5-MATH-1.5B [54], Phi-2 [24], Gemma-2-2B [52], and Llama-3.2-1B [19] on the datasets GSM8K [12], SVAMP [41], AsDiv [35], OpenBookQA [36], and a (512/256 train/evaluation) subset of MATH [21]; Table 3 summarizes the model and dataset details. Learning rate scheduler. For a fair comparison, we fix the learning rate α to the value that performs the best for the FO baseline, chosen from the set {10−6 , 5 × 10−6 , 10−5 }, and consequently select the best bound β from {3, 5, 10}; the remainder of the configuration details are in §C.1.2. 4.1

ZFO for LLM fine-tuning and classic optimization

(i) Does directional ZO step selection improve FO optimization? We show the performance of the ZFO framework on LLM fine-tuning in Table 1, where at least one ZFO variant improves on the first-order fine-tuning (FO FT) baseline for every model-dataset pair, indicating that ZFO’s benefits are not confined to a single model family or benchmark. We observe the largest gains where the FO baseline is relatively weak or highly variable, including SVAMP for Phi-2 and Gemma-2-2B, as well as AsDiv for Llama and MATH and OpenBookQA for Qwen. When the FO baseline is already strong, such as Qwen on SVAMP, ZFO yields smaller but still positive improvements in mean performance. Although the best-performing ZFO variant changes across tasks and models, Padé3 provides a strong representative instance of the framework, improving over FO in 14 of the 16 LLM settings summarized in Table 9. Padé3 is a strong representative default in our experiments, although the preferred local model remains task-dependent. (ii) Qwen Sub-Study. Table 2 extends the main results by isolating Qwen-2.5-Math-1.5B on four reasoning benchmarks: ARC-Challenge [11], StrategyQA [15], FOLIO [20], and CODAH [10]. At least one ZFO method outperforms the FT baseline on every dataset, with the largest improvements on ARC-Challenge (↑8.81 over FO FT) and CODAH (↑4.92 over FO FT). Even under the best learning rate configurations, the main results and the Qwen sub-study indicate that FO FT exhibits seed variance under our settings. (iii) Classic Optimization. To show that ZFO is not intrinsically tied to RLVR and to examine ZFO’s performance beyond LLM fine-tuning, Figure 2 presents three classic optimization problems; see Table 8 for the comprehensive summary of the six classic optimization problems. At least one ZFO variant achieves a substantially lower median objective than its corresponding FO baseline on every displayed problem after 100 updates. Taylor3 –AdamW yields the strongest results on nonlinear least squares and Beale, while Padé3 –AdamW performs best on Rosenbrock. Consistent with our 8

Median log10 objective

(a) Nonlinear least squares

(b) Rosenbrock

(c) Beale 2

0

0

0 −2

−2 −5

−4

−4 −6

0

25

50

75

100

0

25

Gradient updates

50

75

100

0

25

Gradient updates

FO

Taylor2

Taylor3

Pade2

AdamW directions

50

75

100

Gradient updates Pade3

GeN

SGD directions

Figure 2: Results on classic optimization problems. We plot the median log10 objective over three perturbed

80000

2

60000

1 0.20

0.17

0

GSM8K

MATH

Mean Relative Percentages Across Datasets

0.17

0.17

AsDiv

SVAMP

94.8%

Iteration Wall Time FO-only Time ZFO Extra Time

5.2% 2.8% 0

0.17

OpenBookQA

20

40

60

Relative Percentage

FO ZFO overhead Probe-only overhead Trajectory rollout

+11.0% 53,766

59,655

40000 20000 0

Timing bars

Percentage bars 88.1% 80 100

GPU Memory Required (MB)

3

GPU Memory Required (MB)

Time (s)

initializations. GeN has a similar computational budget as ZFO, as both methods perform two additional function evaluations per iteration. We study the extended classic optimization problems in Table 8 and Figure 5. Per-Dataset Timing Breakdown Qwen Memory Allocation Qwen GSM8K

50000 40000 30000 20000 10000 0

FO

ZFO

+13.8% 34,177

FO

38,891

ZFO

Llama AsDiv

42.3% 95.9% 100.0% 88.1% 66.3% 82.8% 85.3%

100% 50%

+19.1% 55,451

FO

66,054

ZFO

Trajectory rollout Policy forward Reference forward Backward Optimizer step Stored direction z Final blend 0%

80000 60000 40000 20000 0

Phi-2 OpenBookQA

Figure 3: Wall-clock time comparison. ZFO step selection Figure 4: GPU memory requirement comadds little wall-clock overhead. Across five benchmarks, trajec- parisons. by FO FT vs ZFO FT across three tory rollout dominates iteration time, while the additional ZFO configurations. The top right panel depicts computation contributes only a small fraction of the time that the relative memory allocation for FO and FO-only updates cost. ZFO FT.

LLM experiments, the best approximant varies across problems, indicating that the gains from step selection depend on the local objective. (iv) Ablation Studies. We provide additional ablation studies in Appendix C.2 to better understand when and how ZFO improves step selection. Figure 8 compares ZFO against AdamW across fixed learning rates, showing that gains in Table 1 cannot generally be recovered by simply changing the FO learning rate. Figure 9 further studies the interaction between the learning rate α and search bound β, where increasing β is generally beneficial at smaller learning rates but can become unstable once α is large, showing that the two parameters are complementary but not interchangeable. Figure 7 provides a clear β = 1 control. Table 13 shows that the four ZFO instances use the search interval differently, ranging from primarily interior steps for Taylor3 to mostly endpoint decisions for Padé2 . 4.2

Computational Cost

The relative computational overhead of ZFO is small. A standard iteration involves 4 key components: a rollout call for the trajectory τt , a reward evaluation R(τt ), a forward pass, and a backward pass. In contrast, our ZFO methods leave all 4 operations untouched and only call two additional forward-pass probes. Since these probes are only forward passes, the overhead is small and comparable to the forward-only cost structure of ZO methods such as MeZO [34]. From this perspective, we present both wall-clock time costs in Figure 3 and memory overhead costs in Figure 4. We train the models on an NVIDIA H100 80GB GPU. Figure 3 decomposes per-iteration runtime for Qwen-2.5-MATH-1.5B on the benchmarks reported in Table 1. ZFO adds ∼ 5.2% relative wall-clock overhead in this configuration; ZFO incurs real memory overhead. 9

Figure 4 reports peak GPU memory across three model–dataset pairs, alongside a per-phase breakdown for Qwen-2.5-MATH-1.5B. The overhead ranges from 11.0% (Qwen/GSM8K) to 19.1% (Phi-2/OpenBookQA). We attribute this range to configuration differences; ZFO probes are absorbed with a smaller relative increase, whereas configurations with a small fixed group size and a larger batch see the overhead grow more. The memory allocation panel shows that the dominant cost remains the reference forward pass, while storing the update direction and applying the final blend sit below the existing bottleneck; hence, ZFO does not introduce a new peak; it widens the existing one. We note that the memory-saving trick MeZO relies on does not transfer here, since our update direction is the FO gradient rather than a random direction.

5

Limitations, Discussion, and Conclusion

What can two probes recover? The backward pass already supplies the local slope; two symmetric probes add estimates of curvature and its variation. Theorem 2 translates local-model error into a bound on the value lost relative to the best step in I. A wider interval permits larger updates but increases approximation error. Bounds and Padé fallbacks limit extrapolation and avoid poles, but do not guarantee that every step improves the surrogate. Theorem 3 then connects step-selection accuracy to convergence. Under the stated assumptions, ZFO therefore retains a stationarity guarantee even with approximate steps. Padé2 Behavior. Padé2 provides a useful diagnostic because, as shown in the toy validation in §C, it behaves the most like FO since its updates select the endpoint of the search region unless a pole is detected, in which case it falls back to Taylor2 . This makes Padé2 a conservative reference for distinguishing whether ZFO gains come from expanding the search region alone or intelligent step selection. In Tables 1 and 2, the richer surrogates frequently outperform Padé2 (e.g., Qwen on MATH, ARC-Challenge, and CODAH). Conversely, when Padé2 improves over FO, then the increased search region contributes in beneficial ways that are not obtainable simply by scaling α. Beyond the Four Instances. We view Taylor2 , Taylor3 , Padé2 , and Padé3 as initial ZFO instances rather than the strongest possible realizations of the framework. Each method is a simple analytic model built from lightweight probes, yet these lightweight probes already yield measurable gains in evaluation performance and training stability across several benchmarks. The main opportunity ZFO provides is not limited to the specific empirical gains reported in this paper, but instead lies in the observation that a meaningful curvature signal can be extracted along the FO direction at low additional cost. Conclusion. Step-size selection is a classic problem in large-scale optimization. The direction selected by a first-order optimizer is valuable, but the step length remains uncertain. We show that two additional function evaluations along that direction are sufficient to construct useful local models of the one-dimensional objective. This creates a general fixed-query step-selection framework, ZFO. In contrast to ZO methods, ZFO does not estimate a gradient with zeroth-order queries; it estimates how the objective changes along a gradient direction whose orientation has already been supplied by a first-order optimizer. Existing methods have already explored function-based step selection along a first-order direction, most notably GeN. Our contribution elevates this idea from a single-step rule to a fixed-query directional modeling framework—the same two probes can support multiple local models, including higher-order Taylor and rational Padé approximations, while the analysis separates finite-difference estimation error, local-model approximation error, and global optimization error. Therefore, ZFO is not primarily a new Taylor or Padé step-size formula. It is a fixed-query framework for exploiting zeroth-order information along an already trusted first-order direction, together with a theory that explains the statistical and optimization error introduced by this separation. In that sense, GeN is a special case of the broader ZFO abstraction, rather than an unrelated baseline. In the same spirit, we encourage the practitioners to add their favorite functional approximators to the ZFO framework. Across our experiments, ZFO often improves step selection over the fixed FO update. The quality of the selected step remains dependent on the local model, search radius, and alignment of the underlying first-order direction. Acknowledgment. During this work, Aritra Dutta was partially supported by the Florida Department of Health Grant AWD00007072 and the National Science Foundation Grant 2321986. 10

References [1] Shihab Ahmed, El Houcine Bergou, Yue Wang, and Aritra Dutta. Stabilizing policy gradient methods via reward profiling. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 40, pages 19560–19568, 2026. [2] Krishnakumar Balasubramanian and Saeed Ghadimi. Zeroth-order (non)-convex stochastic optimization via conditional gradient and gradient updates. Advances in Neural Information Processing Systems, 31, 2018. [3] Lukas Balles and Philipp Hennig. Dissecting Adam: The Sign, Magnitude and Variance of Stochastic Gradients. In Jennifer Dy and Andreas Krause, editors, Proceedings of the 35th International Conference on Machine Learning, volume 80 of Proceedings of Machine Learning Research, pages 404–413. PMLR, 2018. [4] Jonathan Barzilai and Jonathan M Borwein. Two-point step size gradient methods. IMA journal of numerical analysis, 8(1):141–148, 1988. [5] El Houcine Bergou, Eduard Gorbunov, and Peter Richtárik. Stochastic three points method for unconstrained smooth minimization. SIAM Journal on Optimization, 30(4):2726–2749, 2020. [6] Adel Bibi, El Houcine Bergou, Ozan Sener, Bernard Ghanem, and Peter Richtárik. A stochastic derivative-free optimization method with importance sampling: Theory and learning to control. Proceedings of the 34th AAAI Conference on Artificial Intelligence, 2020. [7] Soumia Boucherouite, Grigory Malinovsky, Peter Richtárik, and El Houcine Bergou. Minibatch stochastic three points method for unconstrained smooth minimization. Proceedings of the AAAI Conference on Artificial Intelligence, 38(18), 2024. [8] Richard P Brent. Algorithmsfor minimization without derivatives. Prentiee-Hall, Englewood Cliffs, New Jersey, 1973. [9] Zhiqi Bu and Shiyun Xu. Gradient descent with generalized newton’s method. In The Thirteenth International Conference on Learning Representations, 2025. [10] Michael Chen, Mike D’Arcy, Alisa Liu, Jared Fernandez, and Doug Downey. CODAH: An adversarially-authored question answering dataset for common sense. In Proceedings of the 3rd Workshop on Evaluating Vector Space Representations for NLP, pages 63–69, 2019. [11] Peter Clark, Isaac Cowhey, Oren Etzioni, Tushar Khot, Ashish Sabharwal, Carissa Schoenick, and Oyvind Tafjord. Think you have solved question answering? try ARC, the AI2 reasoning challenge. arXiv preprint arXiv:1803.05457, 2018. [12] Karl Cobbe, Vineet Kosaraju, Mohammad Bavarian, Mark Chen, Heewoo Jun, Lukasz Kaiser, Matthias Plappert, Jerry Tworek, Jacob Hilton, Reiichiro Nakano, et al. Training verifiers to solve math word problems. arXiv preprint arXiv:2110.14168, 2021. [13] Jeremy Cohen, Simran Kaur, Yuanzhi Li, J Zico Kolter, and Ameet Talwalkar. Gradient descent on neural networks typically occurs at the edge of stability. In International Conference on Learning Representations, 2021. [14] Aritra Dutta, El Houcine Bergou, Soumia Boucherouite, Nicklas Werge, Melih Kandemir, and Xin Li. Demystifying the Myths and Legends of Nonconvex Convergence of SGD. arXiv preprint arXiv:2310.12969, 2023. [15] Mor Geva, Daniel Khashabi, Elad Segal, Tushar Khot, Dan Roth, and Jonathan Berant. Did aristotle use a laptop? a question answering benchmark with implicit reasoning strategies. Transactions of the Association for Computational Linguistics, 9:346–361, 2021. [16] Saeed Ghadimi and Guanghui Lan. Stochastic first-and zeroth-order methods for nonconvex stochastic programming. SIAM journal on optimization, pages 2341–2368, 2013. [17] Eduard Gorbunov, Adel Bibi, Ozan Sener, El Houcine Bergou, and Peter Richtarik. A stochastic derivative free optimization method with momentum. In International Conference on Learning Representations, 2020. 11

[18] Priya Goyal, Piotr Dollár, Ross Girshick, Pieter Noordhuis, Lukasz Wesolowski, Aapo Kyrola, Andrew Tulloch, Yangqing Jia, and Kaiming He. Accurate, large minibatch SGD: Training ImageNet in 1 hour. arXiv preprint arXiv:1706.02677, 2017. [19] Aaron Grattafiori, Abhimanyu Dubey, Abhinav Jauhri, Abhinav Pandey, Abhishek Kadian, Ahmad Al-Dahle, Aiesha Letman, Akhil Mathur, Alan Schelten, Alex Vaughan, et al. The llama 3 herd of models. arXiv preprint arXiv:2407.21783, 2024. [20] Simeng Han, Hailey Schoelkopf, Yilun Zhao, Zhenting Qi, Martin Riddell, Wenfei Zhou, James Coady, David Peng, Yujie Qiao, Luke Benson, et al. Folio: Natural language reasoning with first-order logic. In Proceedings of the 2024 Conference on Empirical Methods in Natural Language Processing, pages 22017–22031, 2024. [21] Dan Hendrycks, Collin Burns, Saurav Kadavath, Akul Arora, Steven Basart, Eric Tang, Dawn Song, and Jacob Steinhardt. Measuring mathematical problem solving with the MATH dataset. In Thirty-fifth Conference on Neural Information Processing Systems Datasets and Benchmarks Track (Round 2), 2021. [22] Edward J Hu, Yelong Shen, Phillip Wallis, Zeyuan Allen-Zhu, Yuanzhi Li, Shean Wang, Lu Wang, Weizhu Chen, et al. LoRA: Low-rank adaptation of large language models. ICLR, 1 (2):3, 2022. [23] Kevin G Jamieson, Robert Nowak, and Ben Recht. Query complexity of derivative-free optimization. Advances in neural information processing systems, 25, 2012. [24] Mojan Javaheripi, Sébastien Bubeck, Marah Abdin, Jyoti Aneja, Caio Cêsar Teodoro Mendes, Weizhu Chen, Allie Del Giorno, Ronen Eldan, Sivakanth Gopi, et al. Phi-2: The surprising power of small language models. Microsoft Research Blog, 2023. [25] Keller Jordan, Yuchen Jin, Vlado Boza, Jiacheng You, Franz Cesista, Laker Newhouse, and Jeremy Bernstein. Muon: An optimizer for hidden layers in neural networks, 2024. [26] Bumsu Kim, Daniel McKenzie, Hanqin Cai, and Wotao Yin. Curvature-aware derivative-free optimization. Journal of Scientific Computing, 103(2):43, 2025. [27] Yu-Ang Lee, Ching-Yun Ko, Pin-Yu Chen, and Mi-Yen Yeh. Learning rate matters: Vanilla loRA may suffice for LLM fine-tuning. In ICML’26 workshop on CoLoRAI - The 2nd Workshop on Connecting Low-rank Representations in AI, 2026. [28] Zhe Li, Bicheng Ying, Zidong Liu, Chaosheng Dong, and Haibo Yang. Reconciling HessianInformed Acceleration and Scalar-Only Communication for Efficient Federated Zeroth-Order Fine-Tuning. arXiv preprint arXiv:2506.02370, 2025. [29] 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. In The Thirty-ninth Annual Conference on Neural Information Processing Systems, volume 38, pages 21139–21169, 2026. [30] Ilya Loshchilov and Frank Hutter. SGDR: Stochastic gradient descent with warm restarts. In International Conference on Learning Representations, 2017. [31] Ilya Loshchilov and Frank Hutter. Decoupled weight decay regularization. In International Conference on Learning Representations, 2019. [32] Dougal Maclaurin, David Duvenaud, and Ryan Adams. Gradient-based hyperparameter optimization through reversible learning. In International conference on machine learning, pages 2113–2122. PMLR, 2015. [33] Maren Mahsereci and Philipp Hennig. Probabilistic line searches for stochastic optimization. Journal of Machine Learning Research, 18(119):1–59, 2017. [34] 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. 12

[35] Shen-Yun Miao, Chao-Chun Liang, and Keh-Yih Su. A diverse corpus for evaluating and developing English math word problem solvers. In Proceedings of the 58th annual meeting of the Association for Computational Linguistics, pages 975–984, 2020. [36] Todor Mihaylov, Peter Clark, Tushar Khot, and Ashish Sabharwal. Can a suit of armor conduct electricity? a new dataset for open book question answering. In Proceedings of the 2018 conference on empirical methods in natural language processing, pages 2381–2391, 2018. [37] Konstantin Mishchenko and Aaron Defazio. Prodigy: an expeditiously adaptive parameter-free learner. In Proceedings of the 41st International Conference on Machine Learning, 2024. [38] Yurii Nesterov and Vladimir Spokoiny. Random gradient-free minimization of convex functions. Foundations of Computational Mathematics, 17(2):527–566, 2017. [39] Jorge Nocedal and Stephen J Wright. Numerical optimization. Springer, 2nd edition, 2006. [40] Courtney Paquette and Katya Scheinberg. A stochastic line search method with expected complexity analysis. SIAM Journal on Optimization, 30(1):349–376, 2020. [41] Arkil Patel, Satwik Bhattamishra, and Navin Goyal. Are NLP models really able to solve simple math word problems? In Proceedings of the 2021 conference of the North American chapter of the association for computational linguistics: human language technologies, pages 2080–2094, 2021. [42] Boris T. Polyak. Minimization of nonsmooth functionals. USSR Computational Mathematics and Mathematical Physics, 9(3):14–29, 1969. [43] Junbin Qiu, Zhengpeng Xie, Xiangda Yan, Yongjie Yang, and Yao Shu. Zeroth-Order Optimization is Secretly Single-Step Policy Optimization. arXiv preprint arXiv:2506.14460, 2025. [44] Ruizhong Qiu and Hanghang Tong. Gradient compressed sensing: A query-efficient gradient estimator for high-dimensional zeroth-order optimization. arXiv preprint arXiv:2405.16805, 2024. [45] John Schulman, Sergey Levine, Pieter Abbeel, Michael Jordan, and Philipp Moritz. Trust region policy optimization. In International conference on machine learning, pages 1889–1897. PMLR, 2015. [46] John Schulman, Filip Wolski, Prafulla Dhariwal, Alec Radford, and Oleg Klimov. Proximal policy optimization algorithms. arXiv preprint arXiv:1707.06347, 2017. [47] Zhihong Shao, Peiyi Wang, Qihao Zhu, Runxin Xu, Junxiao Song, Xiao Bi, Haowei Zhang, Mingchuan Zhang, YK Li, Yang Wu, et al. Deepseekmath: Pushing the limits of mathematical reasoning in open language models. arXiv preprint arXiv:2402.03300, 2024. [48] Leslie N Smith. Cyclical learning rates for training neural networks. In 2017 IEEE winter conference on applications of computer vision (WACV), pages 464–472. IEEE, 2017. [49] Leslie N Smith. A disciplined approach to neural network hyper-parameters: Part 1–learning rate, batch size, momentum, and weight decay. arXiv preprint arXiv:1803.09820, 2018. [50] Jacob Mitchell Springer, Sachin Goyal, Kaiyue Wen, Tanishq Kumar, Xiang Yue, Sadhika Malladi, Graham Neubig, and Aditi Raghunathan. Overtrained language models are harder to fine-tune. arXiv preprint arXiv:2503.19206, 2025. [51] Joe Suk and Yaqi Duan. On the optimization dynamics of RLVR: Gradient gap and step size thresholds. arXiv preprint arXiv:2510.08539, 2025. [52] Gemma Team, Morgane Riviere, Shreya Pathak, Pier Giuseppe Sessa, Cassidy Hardin, Surya Bhupatiraju, Léonard Hussenot, Thomas Mesnard, Bobak Shahriari, Alexandre Ramé, et al. Gemma 2: Improving open language models at a practical size. arXiv preprint arXiv:2408.00118, 2024. 13

[53] Ronald J Williams. Simple statistical gradient-following algorithms for connectionist reinforcement learning. Machine learning, 8(3):229–256, 1992. [54] An Yang, Beichen Zhang, Binyuan Hui, Bofei Gao, Bowen Yu, Chengpeng Li, Dayiheng Liu, Jianhong Tu, Jingren Zhou, Junyang Lin, Keming Lu, Mingfeng Xue, Runji Lin, Tianyu Liu, Xingzhang Ren, and Zhenru Zhang. Qwen2.5-Math Technical Report: Toward Mathematical Expert Model via Self-Improvement, 2024. [55] Yihua Zhang, Pingzhi Li, Junyuan Hong, Jiaxiang Li, Yimeng Zhang, Wenqing Zheng, Pin-Yu Chen, Jason Lee, Wotao Yin, Mingyi Hong, et al. Revisiting Zeroth-Order Optimization for Memory-Efficient LLM Fine-Tuning: A Benchmark. In International Conference on Machine Learning, volume 235, pages 59173–59190, 2024. [56] Yanjun Zhao, Sizhe Dang, Haishan Ye, Guang Dai, Yi Qian, and Ivor Tsang. Second-Order Fine-Tuning without Pain for LLMs: A Hessian-Informed Zeroth-Order Optimizer. In The Thirteenth International Conference on Learning Representations, 2025.

14

Contents 1

Introduction

1

2

Methodology

4

3 Theoretical Analysis

4

5

5

3.1

Local Convergence Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5

3.2

Global convergence Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

Experiments

7

4.1

ZFO for LLM fine-tuning and classic optimization . . . . . . . . . . . . . . . . .

8

4.2

Computational Cost . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

9

Limitations, Discussion, and Conclusion

10

A Related work

15

B Addendum to Theoretical Analysis

16

B.1 ZFO Maximum Error . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

18

B.2 ZFO Global Stationarity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

23

B.3 Finite Difference Lemmas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

24

C Addendum to Experimental Results

26

C.1 Configuration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

26

C.1.1

Surrogate Objective Formulation . . . . . . . . . . . . . . . . . . . . . . .

26

C.1.2

Empirical Configurations . . . . . . . . . . . . . . . . . . . . . . . . . . .

28

C.2 Additional Insights . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

29

D Notation

35

Organization. We organize the Appendix as follows: In Section A, we detail the related work. In Section B, we quote the general assumptions and provide detailed proofs of our theoretical result in the main paper. Finally, in Section C, we discuss the configuration details for our experiments in §C.1 and provide additional performance insights in §C.2.

A

Related work

We position ZFO within three lines of work that motivate or contrast with our approach. Together, these related works inspire ZFO’s central premise that cheap query calls can help model local objective behavior along the gradient direction. Policy Optimization Stability & Update Regulation. Stability in policy optimization has motivated a range of update-regulation strategies. In policy optimization, methods such as TRPO [45] enforce a KL-constrained trust region, while PPO [46] replaces explicit constraints with a clipped surrogate objective to stabilize updates. More recent approaches, such as GRPO [47], further emphasize that training stability is highly sensitive to the way policy updates are regulated. Related works also consider wrapper-based stabilization strategies for policy-gradient methods that regulate whether an update should be accepted [1]. Our work is motivated by the same stability concern; however, 15

we approach it from a different perspective by explicitly modeling the objective along the update direction path. Zeroth-Order Optimization. Our work is also related to ZO optimization for fine-tuning large language models. Methods such as MeZO and CARS use central finite differences along random directions to approximate gradient or curvature behavior [26, 34]; particularly, MeZO uses finite differences along random directions to approximate the gradient, and CARS uses a second-order finite difference to approximate the Hessian. Extensions of MeZO include curvature, sparsity, and queryefficient estimators. Hi-ZOO incorporates a diagonal Hessian estimate at the cost of an extra forward pass [56]. Sparse-MeZO restricts updates to noise-resistant weights to enable larger step sizes [29]. GraCe leverages gradient sparsity to achieve a dependence on dimension in query complexity [44]. HiSo carries Hessian-informed ZO into the federated setting [28]. Complementary benchmarking has mapped the landscape of ZO optimizers across fine-tuning paradigms [55]. Recent theoretical work further establishes that ZO optimization is equivalent to single-step policy optimization [43], suggesting a natural connection between zeroth-order methods and RLVR that ZFO exploits. Line-Search & Adaptive Step Size. Broadly, our method is similar to line-search and adaptive learning-rate methods, including Polyak step size [42], Barzilai-Borwein [4], and stochastic variants such as Probabilistic line-search [33] and Stochastic Armijo line-search [40]. These approaches adapt step sizes by repeatedly evaluating candidate points along a univariate objective, often requiring multiple queries. In contrast, our method uses only two symmetric probes to approximately recover higher-order directional information to construct a local model of the objective. This results in a fixed-budget alternative to classical line-search methods. Generalized Newton Step Selection. Among closely related recent methods, generalized Newton step selection or GeN [9] has already explored function-based step selection along a first-order direction, providing an early example of combining FO directions with low-query step-size adaptation. Our contribution is to elevate this idea from a single-step rule to a fixed-query directional modeling framework: the same two probes can support multiple local models, including higher-order Taylor and rational Padé approximations, while the analysis separates finite-difference estimation error, local-model approximation error, and global optimization error. GeN is a special case of the broader ZFO abstraction, rather than an unrelated baseline.

B

Addendum to Theoretical Analysis

This section provides the theoretical foundations needed to characterize when, why, and how ZFO can reliably improve step selection. We first establish concentration guarantees for the finite-difference curvature estimates, then bound how closely the resulting models recover the best step within their local interval, and we finally connect these local guarantees to global stationarity of the ZFO update. Together, these results separate the sources of error introduced by stochastic estimation, local function approximation, and step selection, providing a theoretical justification for the mechanisms used by ZFO. To justify using finite-difference curvature information for ZFO step selection, we first ask whether the second- and third-order directional derivatives estimated from a shared batch can reliably approximate their population counterparts. Hoeffding’s Inequality. Let X be a random variable such that a ≤ X ≤ b almost surely. Then Hoeffding’s lemma gives  2  λ (b − a)2 E [exp (λ(X − E[X]))] ≤ exp , for all λ ∈ R. (4) 8 Thus, every bounded random variable is sub-Gaussian with parameter σ = (b − a)/2. Let X1 , . . . , Xn be independent sub-Gaussian random variables with parameter σ > 0, i.e.,   2 2 E eλXi ≤ eλ σ /2 , ∀λ ∈ R. Then the sub-Gaussian concentration inequality ! r gives: n X 1 2 log(2/ζ) P Xi ≤ σ ≥ 1 − ζ. n i=1 n 16

(5)

Theorem 1 (CRN curvature concentration). at q Suppose Assumptions 1–4 hold. Then, with probability q 2 log(4/δ) 2 log(4/δ) Bϵ2 Lϵ2 ′′ ′′′ ′′ ′′′ c c least 1 − δ, ϕ (0) − ψ (0) ≤ +σ , and ϕ (0) − ψ (0) ≤ +σ . 12

n

20

n

Proof. Fix an index i ∈ {1, . . . , n}. By Assumption 1 and the second-order central finite-difference expansion, there exists ξˆi ∈ (−ϵ, ϵ) such that ϕi (ϵ) − 2ϕi (0) + ϕi (−ϵ) ϵ2 (4) ˆ ′′ = ϕ (0) + ϕ (ξi ). i ϵ2 12 i Using Assumption 2, the remainder satisfies ϵ2 (4) ˆ Bϵ2 ϕi (ξi ) ≤ . 12 12 Averaging the above identity over the shared batch gives n n X 1 X ϵ2 (4) ˆ c′′ (0) = 1 ϕ′′i (0) + ϕ (ξi ). ϕ n i=1 n i=1 12 i Therefore, by the triangle inequality, n

X Bϵ2 c′′ (0) − ψ ′′ (0) ≤ 1 (ϕ′′i (0) − ψ ′′ (0)) + . ϕ n i=1 12 By Assumptions 3 and 4, the random variables ϕ′′i (0) − ψ ′′ (0), i = 1, . . . , n, are independent, centered, and sub-Gaussian with parameter σ. Applying Inequality 5 with ζ = δ/2 yields, with probability at least 1 − δ/2, r n 1 X ′′ 2 log(4/δ) ′′ (ϕ (0) − ψ (0)) ≤ σ . n i=1 i n Consequently, r 2 Bϵ 2 log(4/δ) c′′ (0) − ψ (0) ≤ ϕ +σ . 12 n We now repeat the same argument for the third derivative estimate. By the third-order finite-difference expansion, for each i there exists ξ˜i ∈ (−ϵ, ϵ) such that ϵ2 (5) ˜ 3 (ϕi (ϵ) − ϕi (−ϵ) − 2ϵϕ′i (0)) ′′′ = ϕ (0) + ϕ (ξi ). i ϵ3 20 i Assumption 2 gives Lϵ2 ϵ2 (5) ˜ ϕi (ξi ) ≤ . 20 20 Averaging over the shared batch and applying the triangle inequality gives n 1 X ′′′ Lϵ2 ′′′ (0) − ψ ′′′ (0) ≤ ϕc (ϕi (0) − ψ ′′′ (0)) + . n i=1 20 ′′

′′′ Again, by Assumptions 3 and 4, the centered variables ϕ′′′ i = 1, . . . , n, are i (0) − ψ (0), independent and sub-Gaussian with parameter σ. Applying Inequality 5 with ζ = δ/2 yields, with probability at least 1 − δ/2, r n 1 X ′′′ 2 log(4/δ) ′′′ (ϕi (0) − ψ (0)) ≤ σ . n i=1 n

Thus,

r Lϵ2 2 log(4/δ) ′′′ ′′′ c ϕ (0) − ψ (0) ≤ +σ . 20 n Finally, taking a union bound over the two concentration events shows that both inequalities hold simultaneously with probability at least 1 − δ. Therefore,r 2 c′′ (0) − ψ ′′ (0) ≤ Bϵ + σ 2 log(4/δ) , ϕ 12 n 17

and ′′′ (0) − ψ ′′′ (0) ϕc

Lϵ2 ≤ +σ 20

r

2 log(4/δ) . n

This completes the proof. B.1

ZFO Maximum Error

Theorem 2 (Second-order Taylor approximation bound). Suppose Assumptions 1 and 2 hold. Then the 2 2 3 step selected by maximizing the second-order Taylor model satisfies ϕ(sM ) ≥ ϕ(sϕ )− M3R − Bϵ12R . Proof. By Assumption 1 and Taylor’s theorem, there exists ξs such that ϕ′′ (0) ϕ(3) (ξs ) ϕ(s) = ϕ(0) + sϕ′ (0) + s2 + s3 . 2 6 Similarly, for some ξˆs , the approximation T2 admits the expansion, ϕ′′ (0) ϕ(4) (ξˆs ) T2 (s) = ϕ(0) + sϕ′ (0) + s2 + ϵ2 s2 . 2 24 Subtracting these identities gives ϕ(4) (ξˆs ) ϕ(3) (ξs ) − ϵ2 s2 . ϕ(s) − T2 (s) = s3 6 24 Therefore, by the triangle inequality and Assumption 2 we have, M B M R3 ϵ2 BR2 ϕ(3) (ξs ) ϕ(4) (ξˆs ) + ϵ2 s2 ≤ s3 + ϵ2 s2 ≤ + . |ϕ(s) − T2 (s)| ≤ s3 6 24 6 24 6 24 Consequently, |ϕ(s) − T2 (s)| ≤

M R3 ϵ2 BR2 + =: φ. 6 24

We have sϕ := arg max ϕ(s). s∈[0,R]

By the definitions of sϕ and sM , ϕ(sϕ ) ≥ ϕ(sM ) and T2 (sM ) ≥ T2 (sϕ ). Chaining these inequalities with the uniform bound φ yields ϕ(sM ) ≥ T2 (sM ) − φ ≥ T2 (sϕ ) − φ ≥ ϕ(sϕ ) − 2φ. Therefore, M R3 ϵ2 BR2 ϕ(sM ) ≥ ϕ(sϕ ) − − . 3 12

Additional Model Error Bounds. To understand how accurately ZFO’s local models can choose a step, we next ask how much objective value may be lost when maximizing a finite-difference Taylor or Padé approximation. The proof technique used in Theorem 2 generalizes naturally to higher-order ZFO models. This allows us to characterize the approximation gap of Taylor3 and Padé over the search region [0, R]. Corollary 1 establishes a similar bound for Taylor3 as Theorem 2 did for Taylor2 . Corollaries 2 and 3 adapt a similar structure through the residual r(s) := D(s)ϕ(s)−N (s). All three results yield finite-difference and model truncation errors. Padé3 alternative forms. Third-order matching does not uniquely determine the Padé form, since the numerator and denominator degrees can be allocated differently. The construction above uses a quadratic numerator and linear denominator, corresponding to a 2/1 approximant. Below, we consider the complementary 1/2 form, with a linear numerator and quadratic denominator, which is also used in our LLM experiments. Both constructions match the local expansion through third order under their respective nondegeneracy conditions, but generally produce different performance curves. These choices illustrate the flexibility of ZFO’s rational surrogate construction. Definition 2. Let Pk (s) be the k-th order Padé approximation of ϕ(s) over [0, R], where Pk (s) = N (s) 2 D(s) , N (s) = a0 + a1 s, with D(s) = 1 + d1 s for k = 2 and D(s) = 1 + d1 s + d2 s for k = 3, and R = sFO β, such that β > 0. The coefficients of N (s) and D(s) are chosen so that Pk (s) 18

agrees with the k-th order Taylor model Tk (s) when Taylor expanded. Define the Padé residual r(s) := D(s)ϕ(s) − N (s). When we use the Padé2 and Padé3 models we assume the following: Assumption 7. The Padé denominator is bounded away from zero on [0, R], i.e., there exists η > 0 such that |D(s)| ≥ η, ∀s ∈ [0, R], for both Padé2 and Padé3 . Assumption 8. The Padé residual has bounded higher-order derivatives, i.e., there exist Q, P > 0 such that |r(3) (s)| ≤ Q (Padé2 ), |r(4) (s)| ≤ P (Padé3 ), ∀s ∈ [0, R]. Assumption 7 is consistent with our safeguard of Padé3 falling back to Taylor2 when encountering a pole. Corollary 1 (Third-order Taylor approximation bound). Suppose Assumptions 1-2 hold. Then the step 4 Bϵ2 R2 Lϵ2 R3 selected by maximizing the third-order Taylor model satisfies ϕ(sM ) ≥ ϕ(sϕ )− BR 12 − 12 − 60 . Proof. By Assumption 1 and Taylor’s theorem, there exists ξs such that ϕ′′ (0) ϕ′′′ (0) ϕ(4) (ξs ) + s3 + s4 . ϕ(s) = ϕ(0) + sϕ′ (0) + s2 2 6 24 Similarly, the approximation T3 admits the expansion, for some ξˆs and ξ˜s , ϕ′′ (0) ϕ′′′ (0) ϕ(4) (ξˆs ) ϕ(5) (ξ˜s ) T3 (s) = ϕ(0) + sϕ′ (0) + s2 + s3 + ϵ2 s2 + ϵ2 s3 . 2 6 24 120 Subtracting these identities gives ϕ(4) (ξs ) ϕ(4) (ξˆs ) ϕ(5) (ξ˜s ) ϕ(s) − T3 (s) = s4 − ϵ2 s2 − ϵ2 s3 . 24 24 120 Therefore, by the triangle inequality and Assumption 2, ϕ(4) (ξˆs ) ϕ(5) (ξ˜s ) B B L ϕ(4) (ξs ) + ϵ2 s2 + ϵ2 s3 ≤ s4 + ϵ2 s2 + ϵ2 s3 |ϕ(s)−T3 (s)| ≤ s4 24 24 120 24 24 120 ≤

BR4 ϵ2 BR2 ϵ2 LR3 + + . 24 24 120

Consequently, |ϕ(s) − T3 (s)| ≤

ϵ2 BR2 ϵ2 LR3 BR4 + + =: φ. 24 24 120

Define sϕ := arg max ϕ(s). s∈[0,R]

By the definitions of sϕ and sM , ϕ(sϕ ) ≥ ϕ(sM ) and T3 (sM ) ≥ T3 (sϕ ). Chaining these inequalities with the uniform bound φ yields ϕ(sM ) ≥ T3 (sM ) − φ ≥ T3 (sϕ ) − φ ≥ ϕ(sϕ ) − 2φ. Therefore, ϵ2 BR2 ϵ2 LR3 BR4 − − . ϕ(sM ) ≥ ϕ(sϕ ) − 12 12 60 Corollary 2 (Padé2 approximation bound). Let Assumptions 1-2 and 7-8 hold. Then the step selected 3 2 2 by maximizing the Padé2 model satisfies ϕ(sM ) ≥ ϕ(sϕ ) − Bϵ12ηR − QR 3η . Proof. Since P2 (s) = N (s)/D(s), we have D(s)ϕ(s) − N (s) r(s) ϕ(s) − P2 (s) = = . D(s) D(s) By construction, P2 (s) agrees with the second-order Taylor model T2 (s) when Taylor expanded. Since ϕ(0) and ϕ′ (0) are exact, this gives r(0) = 0, r′ (0) = 0. 19

Notice that r(s) = D(s)(ϕ(s) − P2 (s)) = (1 + d1 s)

! c′′ (0) ϕ′′ (0) − ϕ s2 + O(s3 ) . 2

Therefore, the second-order coefficient satisfies c′′ (0) r′′ (0) ϕ′′ (0) − ϕ = , 2 2 and using the finite-difference error from Lemma 2 yields Bϵ2 r′′ (0) ≤ . 2 24 Taylor expanding r(s) gives

s2 ′′ s3 r (0) + r(3) (ξs ) 2 6 for some ξs ∈ [0, s]. Therefore, by Assumption 8, Bϵ2 s2 Qs3 |r(s)| ≤ + . 24 6 Since |D(s)| ≥ η by Assumption 7,   QR3 1 Bϵ2 R2 =: φ. + |ϕ(s) − P2 (s)| ≤ η 24 6 r(s) =

By the definition of sM = arg maxs∈[0,R] P2 (s), P2 (sM ) ≥ P2 (sϕ ). Chaining this inequality with the uniform bound φ yields ϕ(sM ) ≥ P2 (sM ) − φ ≥ P2 (sϕ ) − φ ≥ ϕ(sϕ ) − 2φ. Substituting the value of φ gives QR3 Bϵ2 R2 − . ϕ(sM ) ≥ ϕ(sϕ ) − 12η 3η Hence, the result. Corollary 3 (Padé3 approximation bound). Let Assumptions 1-2 and 7-8 hold. Then the step selected 2 3 2 2 2 3 R R4 by maximizing the Padé3 model satisfies ϕ(sM ) ≥ ϕ(sϕ ) − Bϵ12ηR − Lϵ60ηR − |d1 |Bϵ − P12η . 12η Proof. Since P3 (s) = N (s)/D(s), we have r(s) D(s)ϕ(s) − N (s) = . ϕ(s) − P3 (s) = D(s) D(s) By construction, P3 (s) agrees with the third-order Taylor model T3 (s) when Taylor expanded. Since ϕ(0) and ϕ′ (0) are exact, this gives r(0) = 0, r′ (0) = 0. Notice that ! ϕ′′ (0) − ϕ̂′′ (0) 2 ϕ′′′ (0) − ϕ̂′′′ (0) 3 4 2 r(s) = D(s)(ϕ(s)−p3 (s)) = (1+d1 s+d2 s ) s + s + O(s ) . 2 6 Therefore, the second-order coefficient satisfies r′′ (0) ϕ′′ (0) − ϕ̂′′ (0) = , 2 2 and using the finite-difference error yields r′′ (0) Bϵ2 ≤ . 2 24 Similarly, the third-order coefficient satisfies r′′′ (0) ϕ′′′ (0) − ϕ̂′′′ (0) ϕ′′ (0) − ϕ̂′′ (0) = + d1 . 6 6 2 Consequently, r′′′ (0) Lϵ2 Bϵ2 ≤ + |d1 | . 6 120 24 20

Taylor expanding r(s) gives s2 ′′ s3 s4 r (0) + r′′′ (0) + r(4) (ξs ) 2 6 24 for some ξs ∈ [0, s]. Therefore, by Assumption 8,  2  Lϵ Bϵ2 P s4 Bϵ2 s2 s3 + + + |d1 | . |r(s)| ≤ 24 120 24 24 r(s) =

Since |D(s)| ≥ η by Assumption 7,   2   1 Bϵ2 R2 Lϵ Bϵ2 P R4 |ϕ(s) − P3 (s)| ≤ + + |d1 | R3 + =: φ. η 24 120 24 24 By the definition of sM = arg maxs∈[0,R] P3 (s), P3 (sM ) ≥ P3 (sϕ ). Chaining this inequality with the uniform bound φ yields ϕ(sM ) ≥ P3 (sM ) − φ ≥ P3 (sϕ ) − φ ≥ ϕ(sϕ ) − 2φ. Substituting the value of φ gives Bϵ2 R2 Lϵ2 R3 |d1 |Bϵ2 R3 P R4 ϕ(sM ) ≥ ϕ(sϕ ) − − − − . 12η 60η 12η 12η

Together, these three corollaries and Theorem 2 quantify how the additional curvature information used by ZFO models affects the worst-case approximation gap. The Taylor3 bound replaces the O(R3 ) term in Theorem 2 with an O(R4 ) term. Hence the relationship between the errors of Taylor2 and Taylor3 can be most significantly characterized by the power of their search regions R. The Padé2 bound remains closer in order to the Taylor2 bound, but introduces the residual-controlled QR3 4 3η term. The Padé3 bound exhibits a similar O(R ) term as Taylor3 and also uses a residual error term. Overall, the two corollaries reinforce the same conclusion as Theorem 2; when the search region is moderate and the perturbation scale is small, then ZFO maximizers achieve values close to the best step within a local trust region. ZFO Dominance. The following theorem is referenced in the main text as an appendix result. As a supplementary result beyond the main-text analysis, we next ask whether the ZFO-selected Taylor2 step can be guaranteed to perform at least as well as the original first-order step under suitable local curvature conditions. Under Assumption 1 and Assumptions 9-11, we compare the Taylor2 step with the typical FO step by bounding the curvature approximation error, imposing strong concavity on [0, R], and requiring the Taylor-model maximizer to be an interior point. The theorem gives conditions to guarantee ϕ(sM ) ≥ ϕ(sFO ). Assumption 9. Let the difference between the second derivative of ϕ and the second order approximation of ϕ at s = 0 be bounded, i.e., there exists D > 0 such that |ϕ′′ (s) − ϕ(ϵ)−2ϕ(0)+ϕ(−ϵ) |≤ ϵ2 D; ∀s ∈ [0, R]. Assumption 10. Let ϕ be ν-strongly concave on [0, R], i.e., ϕ(u) ≤ ϕ(s) + ϕ′ (s)(u − s) − ν2 (u − s)2 ; for all s, u ∈ [0, R] where ν > 0. 2

′

ϵ ϕ (0) Assumption 11. Let the maximum of T2 (s) be an interior point, i.e., − ϕ(ϵ)−2ϕ(0)+ϕ(−ϵ) ∈ [0, R] and ϕ(ϵ) − 2ϕ(0) + ϕ(−ϵ) < 0.

Theorem 4 (Second-order Taylor win guarantee). If Assumption 1 and Assumptions 9-11 hold, and if sFO ∈ [0, R], then the following relation holds between ϕ(sM ) and ϕ(sFO ): 2DsM if |sFO − sM | = 0 or |sFO − sM | ≥ , then ϕ(sM ) ≥ ϕ(sFO ). ν Proof. We have

1 c′′ T2 (s) = ϕ(0) + ϕ′ (0)s + ϕ (0)s2 , 2

which gives c′′ (0)s. T2′ (s) = ϕ′ (0) + ϕ 21

By Assumption 11, there exists sM ∈ [0, R] such that T2′ (sM ) = 0, which gives c′′ (0)sM = 0, ϕ′ (0) + ϕ and therefore ϕ′ (0) sM = − . c′′ (0) ϕ Notice that Z s  ′ ′ ′ ′ ′′ c c′′ (0) du. ϕ (s) − T2 (s) = ϕ (s) − ϕ (0) − ϕ (0)s = ϕ′′ (u) − ϕ 0

Therefore, by Assumption Z 9,

sM 

|ϕ′ (sM )−T2′ (sM )| =

Z sM Z sM  c′′ (0) du ≤ c′′ (0) du ≤ Ddu = DsM . ϕ′′ (u) − ϕ ϕ′′ (u) − ϕ 0

0

0

Since T2′ (sM ) = 0, this gives

|ϕ′ (sM )| ≤ DsM .

(⋆)

Let d = sFO − sM . If |d| = 0, then sFO = sM , and hence ϕ(sM ) = ϕ(sFO ). Now suppose |d| > 0. By Assumption 10, ν ϕ(sFO ) ≤ ϕ(sM ) + ϕ′ (sM )d − d2 , 2 which gives ν ϕ(sFO ) − ϕ(sM ) ≤ ϕ′ (sM )d − d2 . 2 Furthermore, by (⋆), ϕ′ (sM )d ≤ |ϕ′ (sM )||d| ≤ DsM |d|. Therefore,  ν ν  ϕ(sFO ) − ϕ(sM ) ≤ DsM |d| − d2 = |d| DsM − |d| . 2 2 Thus, if 2DsM , |d| ≥ ν then ν DsM − |d| ≤ 0, 2 and therefore ϕ(sFO ) − ϕ(sM ) ≤ 0. Hence ϕ(sM ) ≥ ϕ(sFO ). Finally, since d = sFO − sM , the sufficient condition is 2DsM . |sFO − sM | ≥ ν If sFO = αϕ′ (0), then ! 1 ϕ′ (0) ′ |sFO − sM | = αϕ (0) − − = |ϕ′ (0)| α + . ′′ ′′ c c ϕ (0) ϕ (0) ′

ϕ (0) , this gives the equivalent condition Using sM = − c ′′ ϕ (0)

α+

1 2D ≥− . ′′ c c ϕ (0) ν ϕ′′ (0)

M Therefore, if |sFO − sM | = 0 or |sFO − sM | ≥ 2Ds ν , then ϕ(sM ) ≥ ϕ(sFO ).

22

B.2

ZFO Global Stationarity

To connect the local ZFO step-selection guarantees to overall optimization behavior, we finally ask whether repeatedly applying these approximate directional steps still drives the training process toward a stationary point of the population objective. ZFO is computed on the sampled surrogate Ft (θ) = F (θ; St ), where St is generated at the current iterate θt and then fixed for nearby perturbations. We define the expected surrogate obtained under the current sampling distribution F̄t (θ) := ES [F (θ; S)]. In our theoretical analysis, the true objective J and the expected surrogate F̄t have the same gradient ∇F̄t (θ) = ∇J(θt ). In this formulation Ft is used for parameter updates, while its expectation F̄t develops a reliable objective signal. Lemma 1 (Surrogate vs Empirical Difference). Let Assumptions 1-5 hold, then maxs∈It J(θt + szt ) − J(θt + st zt ) ≤ ∆t , holds r with probability at least 1 − δ, where 2 log(4/δ) ∆t := (LJ + LF̄ )Rt2 + 2σ(1 + Rt ) + 2W Rt2 + ϕ(sϕ ) − ϕ(sM ). n Proof. Let jt (s) := J(θt + szt ), so that ∆t ≥ jt (s∗ ) − jt (sM ), where s∗ maximizes jt (s) over s ∈ It . Since ∇J(θt ) = ∇F̄ (θt ), we have jt′ (0) = ψt′ (0). Therefore, for any s ∈ It , Z s

jt (s) − jt (0) = sjt′ (0) +

(s − u)jt′′ (u)du.

0

Similarly, ψt (s) − ψt (0) = sψt′ (0) +

Z s

(s − u)ψt′′ (u)du.

0

Subtracting these identities gives

Z s

[jt (s) − jt (0)] − [ψt (s) − ψt (0)] =

(s − u)[jt′′ (u) − ψt′′ (u)]du.

0

By Assumption 5, jt′′ (u) ≤ LJ and ψt′′ (u) ≤ LF̄ , and therefore, Z s by the triangle inequality, LJ + LF̄ 2 (s − u)du = s . [jt (s) − jt (0)] − [ψt (s) − ψt (0)] ≤ (LJ + LF̄ ) 2 0 Consequently, LJ + LF̄ ∗ 2 jt (s∗ ) − jt (0) ≤ (s ) + ψt (s∗ ) − ψt (0), 2 and similarly, LJ + LF̄ jt (sM ) − jt (0) ≥ ψt (sM ) − ψt (0) − (sM )2 . 2 Combining these inequalities yields max J(θt + szt ) − J(θt + sM zt ) ≤ ψt (s∗ ) − ψt (sM ) + (LJ + LF̄ )Rt2 . s∈It

We now bound the difference ψt (s∗ ) − ψt (sM ). By the triangle inequality, ψt (s∗ ) − ψt (sM ) ≤ |ψt (s∗ ) − ϕt (s∗ )| + ϕt (s∗ ) − ϕt (sM ) + |ϕt (sM ) − ψt (sM )|. By the definition of sϕ , ϕt (s∗ ) − ϕt (sM ) ≤ ϕt (sϕ ) − ϕt (sM ). Furthermore, |ψt (s∗ ) − ϕt (s∗ )| + |ϕt (sM ) − ψt (sM )| ≤ 2 max |ψt (s) − ϕt (s)|. s∈It

To bound maxs∈It |ψt (s) − ϕt (s)|, Taylor expanding gives |ψt (s) − ϕt (s)| = (ψt − ϕt )(0) + s(ψt′ − ϕ′t )(0) +

(ψt′′ − ϕ′′t )(ξ) 2 s 2

≤ |(ψt − ϕt )(0)| + |(ψt′ − ϕ′t )(0)|Rt + W Rt2 . Therefore, by Assumption 4, r 2 log(4/δ) max |ψt (s) − ϕt (s)| ≤ σ(1 + Rt ) + W Rt2 . s∈It n 23

Consequently, r ∆t := (LJ + LF̄ )Rt2 + 2σ(1 + Rt )

2 log(4/δ) + 2W Rt2 + (ϕt (sϕ ) − ϕt (sM )). n

Theorem 3 (Stationarity of the ZFO update). Suppose Assumptions 5–6 hold. Assume also that qt /LJ ∈ It . Then, for every T > 0, with probability at least 1 − δ, T −1 X ˜T 2LJ (J ⋆ − J(θ0 )) 2LJ ∆ ˜T = 1 min E∥∇J(θt )∥2 ≤ , where ∆ + E[∆t ]. 0≤t≤T µ2 T µ2 T t=0 Proof. Since LqtJ ∈ It , the comparison step exists within the relevant search interval. Therefore,   qt max J(θt + szt ) ≥ J θt + zt . s∈It LJ By Assumption 5,   qt qt LJ qt2 J θt + zt ≥ J(θt ) + ⟨∇J(θt ), zt ⟩ − . LJ LJ 2 L2J Taking expectation over the randomness used to construct zt , with the history before iteration t held fixed, and using qt = t [⟨∇J(θt ), z t ⟩], gives  E qt LJ qt2 q2 q2 Et J θt + zt = J(θt ) + t . ≥ J(θt ) + t − 2 LJ LJ 2 LJ 2LJ Consequently,   q2 Et max J(θt + szt ) ≥ J(θt ) + t . s∈It 2LJ By construction, J(θt+1 ) = J(θt + st zt ) ≥ max J(θt + szt ) − ∆t . s∈It

Taking Et on both sides therefore gives

qt2 − Et [∆t ]. 2LJ Assumption 6 gives qt2 ≥ µ2 ∥∇J(θt )∥2 , and therefore µ2 ∥∇J(θt )∥2 ≤ Et [J(θt+1 )] − J(θt ) + Et [∆t ]. 2LJ Taking full expectation yields µ2 b t, E∥∇J(θt )∥2 ≤ E[J(θt+1 )] − E[J(θt )] + ∆ 2LJ b t = E[∆t ]. Summing from t = 0 to T and using Assumption 5 yields where ∆ T T X µ2 X b t. E∥∇J(θt )∥2 ≤ J ⋆ − J(θ0 ) + ∆ 2LJ t=0 t=0 Et [J(θt+1 )] ≥ J(θt ) +

Finally, T

min E∥∇J(θt )∥2 ≤

0≤t≤T

˜T = 1 and therefore, using ∆ T

PT

t=0 ∆t ,

b

min E∥∇J(θt )∥2 ≤

0≤t≤T

B.3

1X E∥∇J(θt )∥2 , T t=0

˜T 2LJ (J ⋆ − J(θ0 )) 2LJ ∆ + . µ2 T µ2

Finite Difference Lemmas

Finite-Difference Theorems and Proofs. We provide the proofs for two finite-difference approximations used by ZFO in the form of Lemmas. 24

Lemma 2. Let ϕ(s) be four times differentiable in an ϵ-neighborhood, then: c′′ (0) := ϕ(ϵ) − 2ϕ(0) + ϕ(−ϵ) = ϕ′′ (0) + O(ϵ2 ) ϕ ϵ2 Proof. Expanding ϕ(s) by its Taylor Series expansion gives ϕ(s) = ϕ(0) + sϕ′ (0) + s2

ϕ′′ (0) ϕ′′′ (0) ϕ(4) (ξs ) + s3 + s4 , 2 6 24

ϕ(ϵ) = ϕ(0) + ϵϕ′ (0) + ϵ2

ϕ′′′ (0) ϕ(4) (ξϵ ) ϕ′′ (0) + ϵ3 + ϵ4 2 6 24

ϕ(−ϵ) = ϕ(0) − ϵϕ′ (0) + ϵ2

ϕ′′ (0) ϕ′′′ (0) ϕ(4) (ξ−ϵ ) − ϵ3 + ϵ4 . 2 6 24

hence

and

Adding the two gives ϕ(ϵ) + ϕ(−ϵ) = 2ϕ(0) + ϵ2 ϕ′′ (0) + ϵ4

ˆ ϕ(4) (ξ) , 12

ϕ(ϵ) − 2ϕ(0) + ϕ(−ϵ) = ϵ2 ϕ′′ (0) + ϵ4

ˆ ϕ(4) (ξ) , 12

subtracting 2ϕ(0) gives

and dividing by ϵ2 gives (4) ˆ ϕ(ϵ) − 2ϕ(0) + ϕ(−ϵ) (ξ) ′′ 2ϕ = ϕ′′ (0) + O(ϵ2 ). = ϕ (0) + ϵ ϵ2 12

The equality above represents a natural second-order central finite difference; we prove the asymmetric finite difference for Lemma 3 below, Lemma 3. Let ϕ(s) be five times differentiable in an ϵ-neighborhood, then: ϕ(ϵ) − ϕ(−ϵ) − 2ϵϕ′ (0) ′′′ (0) := 3 ϕc = ϕ′′′ (0) + O(ϵ2 ). ϵ3 Proof. Consider the new expansion of ϕ(s) ϕ(s) = ϕ(0) + sϕ′ (0) + s2

ϕ′′ (0) ϕ′′′ (0) ϕ(4) (0) ϕ(5) (ξs ) + s3 + s4 + s5 . 2 6 24 120

Then taking the difference between ϕ(ϵ) and ϕ(−ϵ) returns ϕ(ϵ) − ϕ(−ϵ) = 2ϵϕ′ (0) + ϵ3

ˆ ϕ′′′ (0) ϕ(5) (ξ) + ϵ5 . 3 60

ϕ(ϵ) − ϕ(−ϵ) − 2ϵϕ′ (0) = ϵ3

ˆ ϕ′′′ (0) ϕ(5) (ξ) + ϵ5 , 3 60

Using ϕ′ (0) gives

3

and dividing by ϵ3 returns 25

Table 3: Summary of the models and datasets used in the large-model experiments. Parameter counts denote the total number of model parameters. Dataset sizes are reported as the number of training/evaluation examples used in our experiments. MATH uses the 512/256 train/evaluation subset described in the main text. Study

Model

Model Parameters

Dataset Size (Train / Eval)

Qwen-2.5-Math-1.5B

1.54B

Phi-2

2.70B

SVAMP AsDiv OpenBookQA

Gemma-2-2B

2.61B

SVAMP

Llama-3.2-1B

1.23B

GSM8K AsDiv OpenBookQA

7,473 / 1,319 1,844 / 461 4,957 / 500

Qwen-2.5-Math-1.5B

1.54B

ARC-Challenge StrategyQA FOLIO CODAH

1,119 / 1,172 2,061 / 229 1,001 / 204 2,498 / 278

Main Study

Qwen Sub-Study

Dataset GSM8K MATH SVAMP AsDiv OpenBookQA

7,473 / 1,319 512 / 256 700 / 300 1,844 / 461 4,957 / 500 700 / 300 1,844 / 461 4,957 / 500 700 / 300

ˆ 3(ϕ(ϵ) − ϕ(−ϵ) − 2ϵϕ′ (0)) ϕ(5) (ξ) = ϕ′′′ (0) + ϵ2 = ϕ′′′ (0) + O(ϵ2 ). 3 ϵ 20

C

Addendum to Experimental Results

This addendum collects supplementary diagnostics that complement the main results. We first validate the ZFO step-selection mechanism on two low-dimensional toy RLVR environments, then discuss the local-approximation behavior of each ZFO method on a panel of analytic test functions, and finally report cross-setting win ratios that summarize how consistently each variant improves over the FO baseline. C.1

Configuration

This section specifies the objective and implementation choices used throughout our experiments. First, we define the GRPO surrogate, optimized by each method, then report the hyperparameter configurations used for the main large-model runs, Qwen sub-study, and toy validation. C.1.1

Surrogate Objective Formulation

GRPO Surrogate. For this paper, the surrogate Ft instantiates the GRPO objective. Let Pt denote the set of prompts in the current batch, with |Pt | the batch size and G the number of examples Parameter

Setting

Optimizer Policy Temperature Top-p KL coefficient β ⋆ Perturbation scale ϵ Seed

AdamW grpo 0.8 0.95 0.02 1e-3 {100, 101, 102}

Table 4: Fixed experimental settings used across all large language RLVR-ZFO runs.

Config

Batch

Group

Epochs

Max Len.

Output Len.

a

4

4

10

512

192

b

2

2

8

384

192

c

8

2

8

384

192

d

4

2

5

256

192

Table 5: Training configuration types used across runs. Batch represents batch size, Group represents group size, Max Len. and Output Len. represent Max Length and Maximum Output Length. 26

Table 6: Hyperparameter Configuration. Learning rates, bounds, and configuration types used for the main empirical study. FO LR and ZO LR denote the learning rates used for first-order and zeroth-order fine-tuning, respectively. BT2, BT3, BP2, and BP3 represent the selected bounds for Taylor2, Taylor3, Padé2, and Padé3, respectively. Learning rate search is over {1e − 6, 5e−6, 1e−5} and bound search is over {3, 5, 10}. ALL indicates that every lr for MeZO produced the same result. Configuration

Model

Dataset

Qwen-2.5-Math-1.5B

GSM8K MATH SVAMP AsDiv OpenBookQA

5e-6 5e-6 1e-5 5e-6 1e-5

ALL 1e-6 5e-7 5e-7 1e-7

10 5 3 5 10

5 5 5 5 5

10 5 3 5 5

10 5 3 5 5

a c a c a

Phi-2

SVAMP AsDiv OpenBookQA

5e-6 5e-6 5e-6

5e-7 1e-7 1e-7

3 3 3

5 3 5

5 5 3

3 3 3

a d d

Gemma-2-2B

SVAMP

1e-6

5e-9

3

3

3

3

b

Llama-3.2-1B

GSM8K AsDiv OpenBookQA

5e-6 5e-6 1e-6

1e-7 5e-8 ALL

3 5 5

3 5 5

3 3 5

3 3 3

a c b

FO LR ZO LR BT2 BT3 BP2 BP3 Config Type

Table 7: Learning rate, bound, and configuration type used for the Qwen sub-study on additional reasoning benchmarks. Learning rate search is over {5e-6, 1e-5} and bound search is over {3, 5}. Model

Dataset

Qwen-2.5-Math-1.5B

ARC-Challenge StrategyQA FOLIO CODAH

Configuration Learning Rate BT2 BT3 BP2 BP3 Config Type 1e-5 1e-5 1e-5 1e-5

3 5 3 5

5 5 3 5

3 3 3 5

3 5 3 5

c c c c

sampled per prompt. For each prompt i ∈ {1, . . . , |Pt |} and each sampled example j ∈ {1, . . . , G}, let τi,j,t denote the sampled example, |τi,j,t | its length, and (xi,j,v , yi,j,v ) the state-action pair at token position v. The surrogate is |Pt | G 1 XX e Ai,j,t log πθ (τi,j,t ) − κ · KL(πθ ∥πref ), Ft (θ; St ) = |Pt |G i=1 j=1 ei,j,t is the group-relative advantage given as where κ > 0 is the KL coefficient, A G X ei,j,t = ρ(τi,j,t ) − 1 A ρ(τi,k,t ), G k=1

and ρ is the verifiable reward given by ρ(τi,j,t ) = 1τi,j,t ∈Ct , with Ct denoting the set of sampled examples deemed correct by the verifier at iteration t. The KL penalty is computed at the token level and averaged across the batch: |τi,j,t | |Pt | G X  1 XX 1 KL(πθ ∥πref ) = log πθ (yi,j,v | xi,j,v ) − log πref (yi,j,v | xi,j,v ) . |Pt |G i=1 j=1 |τi,j,t | v=1 Implementation-Level Interpretation of α and β. In the theory, we describe the ZFO search interval using the effective radius Rt = βsFO , where sFO represents the scale of the base first-order step. Conceptually, this makes the learning rate α and the bound β nearly interchangeable as step-scale controls. Under an idealized linear scaling view, replacing (α, β) with another pair having the same product αβ would therefore produce a comparable search radius. 27

In implementation, however, the two quantities enter at different stages. The base optimizer first forms its tentative AdamW update using learning rate α. If θt is the parameter vector before this tentative update and θbt is the parameter vector after it, then the realized optimizer displacement is dt = θbt − θt . ZFO then searches along this displacement through parameters θt + adt with a ∈ [0, β]. Thus the implemented effective search radius is Rt = β∥dt ∥. Therefore, matched-product pairs such as (α, β) = (10−6 , 5) and (5 × 10−6 , 1) should be understood as equivalent in nominal search scale, but not necessarily identical in implementation. They can produce different results because α affects the optimizer displacement dt before ZFO constructs its search line. LLM implementation. GRPO supplies the optimizer direction and first-order coefficient. Symmetric probes evaluate a separate reward-minus-KL score on fixed sampled sequences, with rewards held constant. The probe values and GRPO slope jointly define the local model used for step selection. The theory analyzes an idealized ZFO setting in which the directional derivative and function probes are evaluations of the same sampled surrogate. The RLVR implementation approximates this setting by holding sampled responses and rewards fixed during the probes.

C.1.2

Empirical Configurations

All large-model runs share the fixed hyperparameters listed in Table 4. We hold these values constant so that observed differences in performance can be attributed to the update method as opposed to incidental variance in sampling, regularization, or optimization settings. Each reported result is averaged over the three seeds {100, 101, 102}, and the standard deviations reported in Tables 1 and 2 are based on these various runs. We evaluate the models after each epoch, and all runs use the same early-stopping rule; we terminate the training when evaluation accuracy is unchanged for three consecutive epochs. Otherwise, we continue training until the maximum epoch count specified for the corresponding configuration. Per-run hyperparameters differ across model-dataset pairs and include the learning rate α, the bound hyperparameter β, and the selected configuration from Table 5. We summarize four configuration types in Table 5 and represent a performance-versus-cost trade-off. Config a is the default, with the largest group size and longest context length. Configs b–d are progressively cheaper alternatives selected when a is infeasible based on either memory or time constraints. Config b reduces both batch and group sizes for memory-constrained runs, c trades group size for a larger batch when throughput is more valuable than rollout diversity, and d additionally shortens the context and reduces epochs. These per-run details are reported in Table 6 for the main results, including the FO learning rate, ZO learning rate, and each ZFO bound β. The FO learning rate is selected from the sweep {1e-6, 5e-6, 1e-5 }, and the corresponding bound for each ZFO method is then chosen from a validated sweep over {3, 5, 10}. The ZO learning rate is tuned independently from the sweep {1e−8, 5e−8, 1e−7, 5e−7}, since ZO methods generally require a smaller effective step scale than FO updates. The performance reported is the mean verifiable reward on the evaluation set. The Qwen sub-study from Table 2 follows a similar protocol but fixes config c across runs (a memory-conservative choice given the sub-study sweeps over four datasets), restricts the learning rate sweep to {5e-6, 1e-5 }, and limits the bound sweep to {3, 5}. The Qwen sub-study configurations are reported in Table 7. The toy models in Figure 7 use a lightweight configuration designed to validate the ZFO step-selection mechanism in a controlled low-dimensional setting. Both environments share a context dimension of 6, a sequence length of 4, and a two-layer policy network of hidden width 48. Each run uses a dataset of 200 contexts trained for 50 steps with batch size 16, group size G = 4, perturbation scale ϵ = 10−2 , and KL coefficient of 0.1. The bound hyperparameter β is fixed to 1, so the ZFO search interval collapses to the FO typical step. Table 3 depicts the models and datasets used throughout our large-model experiments, along with their corresponding parameter counts and dataset sizes. We consider models ranging from 1.23B to 2.7B parameters across a diverse set of mathematical reasoning, commonsense reasoning, and knowledge-based benchmarks. We additionally report the datasets used in our Qwen sub-study to provide a complete summary of the main experimental settings considered. 28

Table 8: Compact summary of the Classical Optimization Problems used in Figure 5. Before optimization, each listed initial position θ0 is perturbed independently per seed as θ0,s = θ0 + ξs , where ξs ∼ N (0, 0.04I) applies independent Gaussian noise (standard deviation 0.2) to every coordinate, drawn from a dedicated, task-specific random stream and shared across all FO, ZFO, and GeN methods evaluated on that task–seed pair. Problem Linear least squares Nonlinear least squares Logistic regression Low-rank matrix factorization Rosenbrock Beale

Objective 1 ∥Aθ − y∥2 2n
1 P bxi + c − y 2 mina,b,c i i ae 2N  λ 1P  ⊤ ∥θ∥2 log(1 + exi θ ) − yi x⊤ minθ i θ + n i 2 1 λ minU,V ∥U V ⊤ − M ∥2F + ∥θ∥2 2mn 2|θ| minx1 ,x2 100(x2 − x21 )2 + (1 − x1 )2
2 P minx1 ,x2 3k=1 ck − x1 + x1 xk2 minθ

Initial position θ0 = 020 (a0 , b0 , c0 ) = (0.4, 0.8, 1.0) θ0 = 016 (U0 , V0 ) ∼ N (0, 0.152 ) (x1 , x2 )0 = (−1.2, 1.0) (x1 , x2 )0 = (−1.0, 1.0)

Our MeZO baseline uses a forward-only, two-sided finite-difference update along an unnormalized Gaussian perturbation direction. At each iteration, sampled responses, group-centered reward advantages, and reference-policy log probabilities are held fixed across the positive and negative perturbations. The evaluated objective combines advantage-weighted completion log probabilities with a nonnegative k3 KL penalty, averaged over completion tokens within each sequence. Both terms use completion-only masking, and responses are decoded from the end of the padded prompt. Objective evaluations run in evaluation mode, with log-ratio clipping for numerical stability. The perturbation direction is regenerated from the same random seed for both probes and the parameter update, avoiding backpropagation.

C.2

Additional Insights

Optimization Beyond Neural Network Training. To evaluate whether ZFO extends beyond neural network optimization, we consider six classical optimization problems: linear least squares, nonlinear least squares, logistic regression, low-rank matrix factorization, Rosenbrock, and Beale. A summary of the problems can be found in Table 8. We evaluate Taylor2, Taylor3, Padé2, and Padé3 with both SGD and AdamW-style directions, alongside their first-order baselines and GeN [9] and present the results in Figure 5. Every method receives an independent search over the 18 learning rates, {1, 2, 5} × 10−5,...,0 , while ZFO methods use a fixed bound β = 3 and perturbation radius ϵ = 10−3 . GeN uses γ = 0.9 and Φ = 1 for fair comparison against the ZFO methods. Each experiment runs for 100 gradient updates over three seeds. The results show that ZFO’s benefits are substantial but depend on both the local model and direction. Taylor3 –AdamW attains the lowest final median objectives on nonlinear least squares and Beale, whereas Padé3 –AdamW performs best on Rosenbrock. Several ZFO and GeN variants also accelerate matrix factorization relative to ordinary SGD. On linear least squares, most adaptive methods improve moderately over their corresponding first-order baselines, whereas logistic regression converges to essentially the same objective across methods. Some combinations also stall or underperform, particularly on the nonconvex objectives. Thus, the experiment supports ZFO’s applicability beyond LLMs and its ability to yield order-of-magnitude improvements in favorable geometries. All Roads Lead to ϕ. Figure 6 examines the four ZFO surrogates, Taylor2 , Taylor3 , Padé2 , and Padé3 , on a panel of locally analytic functions. In each panel, f (x0 ) and f ′ (x0 ) are computed exactly, while the second and third derivatives are approximated by the same finite differences discussed in this paper. The top row consists of functions that are either sub-linear or polynomial on a tight interval; hence the approximations are nearly indistinguishable from the actual function. The bottom row deliberately widens the window by introducing functions that do not have sub-linear properties. Across both rows, the consistent takeaway is that all four interpolants generally agree with a bit of variance on super-polynomial functions. 29

Optimization beyond LLM Fine-tuning Linear LSE

−2 −3

Matrix Factorization

0 −2 −4

20

40 60 Update Iteration

SGD ZFO-T2-SGD

−4

80

100

ZFO-T3-SGD ZFO-P2-SGD

0.4 0.2 0.0 −0.2 −0.4

Update Iteration

Update Iteration

Rosenbrock

Beale 4

2.5 0.0 −2.5 −5.0 −7.5 −10.0

0.6

Median log10 Objective

Median log10 Objective

Median log10 Objective

2

0

−2

5.0

4

−6

0

−6

Update Iteration

Logistic

0.8 Median log10 Objective

0 −1

−4

Nonlinear LSE

2

1

Median log10 Objective

Median log10 Objective

2

0

20

40 60 Update Iteration

ZFO-P3-SGD GeN-SGD

80

AdamW ZFO-T2-AdamW

100

2 0 −2 −4 −6 −8

0

20

ZFO-T3-AdamW ZFO-P2-AdamW

40 60 Update Iteration

80

100

ZFO-P3-AdamW GeN-AdamW

Figure 5: Median log10 objective over three seeds for six optimization problems. Solid curves use AdamW directions, and dashed curves use SGD directions. Each method uses its best learning rate from the same 18-point search. Lower values indicate better optimization.

Figure 6: ZFO approximations on a variety of locally analytic functions. Each panel plots Taylor2 , Taylor3 , Padé2 , and Padé3 around an x0 . Top row consists of panels where all four interpolants are visually indistinguishable from the actual function. The bottom row consists of panels that expose the qualitative differences between the Taylor and Padé approximants.

Toy Model Validation. We first validate the ZFO framework on two low-dimensional toy RLVR environments shown in Figure 7. The first toy model, Grouped Sequence Bandit, is a sequencedecision task that requires a learning algorithm to choose a short sequence of tokens to match a hidden target sequence. Specifically, the model samples a complete sequence based on a small input and updates the parameters based on the reward feedback. The second model, Grouped Path-Matching, uses a similar idea, but makes the reward structure more dependent on the entire sequence path. Instead of rewarding independent tokens, the policy receives a higher reward when its sampled actions remain aligned with the hidden target path. Both environments use a stochastic policy trained with grouped advantages and KL regularization against a frozen reference, replicating the surrogate structure of Ft (θ; St ) used in our larger experiments. Each method is averaged over 5 seeds with the bound hyperparameter β = 1; the remainder of the configuration can be found in §C.1.2. 30

Figure 7: Toy RLVR validation on the Grouped Sequence Bandit and Grouped Path-Matching environments. ZFO methods improve over FO and ZO baselines under the same search scale β = 1, showing that local curvature models can select more effective steps along update directions.

Figure 8: Stability of ZFO. Empirical evaluation curves demonstrating ZFO’s capacity to stabilize and improve over fine-tuning with AdamW (FO FT). Demonstrations involve AdamW FT (FO FT) curves of multiple learning rates.

Figure 7 reports the training curves and final performance across both environments. Across both toy environments, the ZFO methods improve over the FO and ZO baselines, supporting the central claim that additional function evaluations can improve step-size selection along a fixed direction. Importantly, these improvements occur with the same search scale, since all methods use β = 1. Thus, the ZFO methods are not benefiting from a larger trust region; rather, they are using local curvature behavior to select more effective intermediate steps within the same update interval. This is most visible in the separation between Taylor2 , Taylor3 , and Padé3 relative to FO, showing that the additional curvature information can meaningfully alter the chosen step even when the learning rate is fixed. The behavior of Padé2 further supports this interpretation; among the ZFO methods, Padé2 is the most restricted in the sense it contains only a linear numerator and denominator. Consequently, when β = 1, Padé2 often selects a step close to the standard FO update, making it the ZFO variant most similar to FO. This trend is reflected in the final performance summaries with Taylor3 performing the best in both environments – Taylor2 and Padé3 trailing close behind. Overall, these toy results provide validation of the ZFO step-selection mechanism before moving to higher-dimensional models. Learning Rate Sensitivity Curves. Figure 8 establishes that ZFO’s gains rise from intelligent step selection rather than a hidden learning-rate effect. In each panel, we plot ZFO trajectories alongside FO FT runs surrounding the base learning rate to compare against ZFO runs. The three panels highlight complementary behaviors; e.g., Qwen on OpenBookQA shows ZFO trajectories outperforming FO runs with the same nominal search regions, reinforcing the extremely large margin in Table 1. Phi-2 on Svamp shows ZFO trajectories steadily improving over the FO baselines. Llama on AsDiv shows that the most conservative FO baseline eventually destabilizes during training, whereas ZFO methods remain stable. Together, these curves indicate that ZFO occupies a regime that FO cannot reach by fixed learning rates. Bound vs. Learning Rate Ablation. Figure 9 reports the full ablation over the learning rate α and search bound β that was summarized in the main paper. We sweep both hyperparameters across all four ZFO methods on Qwen-2.5-Math-1.5B, evaluating after 10 epochs on subsets (512/256 train/eval) of GSM8K and SVAMP. As discussed in the main paper, at smaller learning rates, increasing β produces near-monotonic improvements, while at larger learning rates the gains saturate as the search region R grows too large. Comparing (α = 10−6 , β = 5) against (α = 5 × 10−6 , β = 1) also shows that configurations with the same nominal search region yield different final performance, 31

Figure 9: Ablation study inspecting the relationship between the learning rate α and search hyperparameter β on the reward rate. We report the performance over a training set of 512 elements and an evaluation set of 256 elements on the Qwen-2.5-MATH-1.5B model and datasets, GSM8K and SVAMP, on the top and bottom row, respectively. From left to right, the panels represent all four ZFO methods (Taylor2 , Taylor3 , Padé2 , Padé3 ) presented in this paper. Method

Main Results Wins

Qwen Sub-Study Wins

Win Ratio

Taylor2 Taylor3 Padé2 Padé3

10/12 6/12 11/12 11/12

4/4 1/4 2/4 3/4

14/16 7/16 13/16 14/16

Table 9: Win ratios for each ZFO method across the main results and Qwen sub-study. A win denotes a setting in which the method achieves a higher mean accuracy than the corresponding FO baseline.

indicating that α and β are complementary rather than interchangeable. See §C.1.1 for the formal implementation of R. Examining Figure 9 closely, we observe the high-reward region forms a diagonal band that shifts toward smaller β as α grows. At α = 10−6 and α = 5 × 10−6 , performance improves steadily as β increases from 1 to 100 across all four methods. At α = 10−5 , the best results pull back toward small β. Combining the largest bound with α = 10−5 collapses the reward to near-zero for almost all methods. Cross-Setting Improvement Frequency. As an additional appendix diagnostic, we summarize how often each ZFO variant improves over the corresponding FO baseline across the reported evaluation settings in Table 9. Across the 12 main result settings and the 4 Qwen sub-study settings, Taylor2 , Taylor3 , Padé2 , and Padé3 improve over FO in 14/16, 7/16, 13/16, and 14/16 settings, respectively. This supports the view that ZFO is driven by a recurring advantage in adaptive step selection across settings. Additional first-order baselines. To compare ZFO with different learning-rate adaptation strategies, we additionally evaluate AdamW with a warmup-cosine schedule [30] and the parameter-free Prodigy optimizer [37] in Table 10. We preserve the Qwen-2.5-Math-1.5B model, dataset-specific training configurations, GRPO objective, evaluation protocol, and seeds used in the corresponding main and sub-study experiments. For AdamW with cosine decay, the learning rate selected for the constant-rate FO baseline is used as the peak rate; it gets reached by linear warmup over the first 10% of optimizer steps and then decayed to zero according to a cosine schedule. Prodigy instead uses its recommended configuration with adjustment factor 1.0, dcoef = 1.0, decoupled weight decay 0.01, and no external learning-rate schedule. To contextualize these baselines within the range of behavior exhibited by ZFO, the table also reports the lowest- and highest-mean ZFO variants for each dataset. 32

Table 10: Comparison with additional first-order baselines. AdamW with linear warmup and cosine decay and Prodigy are compared with the worst- and best-performing ZFO variants for each Qwen-2.5-Math-1.5B dataset. The ZFO extrema are selected by mean final evaluation performance across Taylor2 , Taylor3 , Padé2 , and Padé3 ; the selected variant is shown in parentheses. The best result in each row is bold. Model

Dataset

Qwen-2.5-Math-1.5B

GSM8K SVAMP AsDiv OpenBookQA ARC-Challenge StrategyQA FOLIO CODAH

AdamW + Cosine

Prodigy

Worst ZFO

Best ZFO

72.88 ± 0.54 86.00 ± 0.33 81.98 ± 0.52 59.60 ± 0.87 55.15 ± 1.11 66.28 ± 0.51 34.48 ± 1.13 48.32 ± 14.52

64.32 ± 0.90 83.22 ± 0.19 80.82 ± 0.30 38.73 ± 0.31 30.80 ± 0.15 60.80 ± 0.95 34.31 ± 0.49 30.58 ± 1.57

75.56 ± 0.24 (Taylor3 ) 89.11 ± 0.51 (Taylor3 ) 88.27 ± 0.35 (Taylor3 ) 37.80 ± 17.85 (Taylor2 ) 53.58 ± 2.77 (Padé2 ) 61.67 ± 1.46 (Taylor3 ) 34.48 ± 1.13 (Taylor3 ) 41.85 ± 20.58 (Padé3 )

81.00 ± 0.11 (Taylor2 ) 90.89 ± 1.02 (Padé2 ) 94.09 ± 0.29 (Taylor2 ) 64.07 ± 1.22 (Padé3 ) 58.10 ± 1.33 (Taylor2 ) 70.06 ± 0.61 (Padé3 ) 36.60 ± 3.26 (Taylor2 ) 62.83 ± 3.14 (Taylor2 )

Table 10 shows that the strongest ZFO variant outperforms both AdamW with cosine decay and Prodigy across all eight datasets on Qwen. AdamW with cosine decay remains the stronger of the two additional FO baselines, while Prodigy consistently trails both the scheduled AdamW baseline and the best ZFO method. However, the worst-performing ZFO variant falls below AdamW with cosine decay on several datasets, indicating that the advantage does not arise uniformly across every local model. Overall, these results suggest that ZFO can improve over many step-adaptive baselines, while also signaling that the choice of local approximation remains important. ZFO Sequential Line Search. The four ZFO methods in the main paper select a step using a local model constructed from a fixed budget of two probes. ZFO-Seq extends this design to a sequential budget of k > 2 fixed-batch evaluations. The first two probes construct the initial model M ∈ {T2 , T3 , P2 , P3 } as in Algorithm 1, and the remaining k − 2 probes are spent refining the model’s proposal. At each refinement step, the current candidate maximizer is evaluated on the shared minibatch St , and a new candidate is proposed by fitting a parabola through the best evaluated point in It and its two nearest evaluated neighbors, aligning with safeguarded successive parabolic interpolation [8]. The returned step is the best evaluated point rather than a model maximizer. Every candidate must survive an actual evaluation of ϕt before it can be selected. Since (0, ϕt (0)) ∈ Et by construction, ZFO-Seq guarantees ϕt (st ) ≥ ϕt (0) at every iteration, a monotone non-decrease of the empirical surrogate. This replaces the model-selection gap ϕ(sϕ ) − ϕ(sM ) with the gap attained by the best of k evaluated points, which shrinks under refinement. Algorithm 2 summarizes the procedure. Table 11 reports the evaluation accuracy for a variety of line-search methods on the Qwen-2.5-Math1.5B model (seed 100, configuration C, learning rate 1e-5, search bound β = 3). Each ZFO-Seq model used P3 as its default initialization. Algorithm 2 : ZFO-Seq Require: θ0 , base optimizer, learning rate α, ϵ > 0, β > 0, budget k ≥ 3, initial model M ∈ {T2 , T3 , P2 , P3 } 1: for t = 0, 1, 2, . . . do 2: Algorithm 1 Inner Loop 3: for j = 1, . . . , k − 2 do 4: Evaluate ϕt (s(j) ) on St and add (s(j) , ϕt (s(j) )) to Et 5: if j < k − 2 then 6: Fit a parabola M through the best point in Et and its two nearest neighbors 7: Propose s(j+1) ∈ arg maxs∈It M(s) 8: end if 9: end for 10: Select st ∈ arg max(s,ϕt (s))∈Et , ϕt (s) 11: Update θt+1 = θt + st zt 12: end for Table 11 compares ZFO-Seq against several established line-search methods. ZFO-Seq achieves the strongest result on CODAH, while remaining within approximately two points of the bestperforming method on both ARC-Challenge and SVAMP. Armijo and Barzilai-Borwein remain 33

Table 11: Comparison with line-search algorithms. Extended line-search results comparing the basic ZFO-Seq method against Armijo, Barzilai-Borwein, Polyak, and a simplified probabilistic line search (PLS) on the Qwen Model. The best result is in bold, and the second best result is underlined. Dataset ZFO-Seq Armijo BB Polyak PLS CODAH ARC-Challenge SVAMP

55.47 55.86 83.59

32.81 57.81 85.55

52.73 54.69 85.55

30.47 35.16 82.81

31.25 35.16 76.17

Table 12: Qwen-2.5-Math-1.5B sub-study ablation with β = 1. Results are mean ± standard deviation, over three seeds. All datasets are trained with 3 epochs (ARC-Challenge, StrategyQA; config a’), (ASDiv; config c’). Under this 3-epoch configuration, reducing the ARC FO learning rate from 10−5 to 5 × 10−6 improves accuracy from 31.71 ± 8.11 to 55.12 ± 0.60. The ASDiv FO result uses the third-epoch checkpoint of the matched FO runs. Dataset ARC-Challenge StrategyQA ASDiv

Base LR −5

1 × 10 1 × 10−5 5 × 10−6

AdamW (FO)

Taylor2

Taylor3

Padé2

Padé3

31.71 ± 8.11

56.77 ± 0.61

56.51 ± 0.36

55.49 ± 0.98

49.86 ± 0.27

59.24 ± 10.14

70.69 ± 5.10

74.87 ± 0.37

67.25 ± 13.64

48.18 ± 0.44

82.08 ± 0.08

81.82 ± 0.26

78.05 ± 1.41

81.87 ± 0.26

81.92 ± 0.46

slightly stronger on these latter settings, whereas Polyak and PLS generally trail ZFO-Seq. Thus, ZFO-Seq is competitive with classical line-search and adaptive step-size methods across the three evaluated settings. Unit-Bound Ablation. To determine whether ZFO benefits require a search interval larger than the standard first-order update, we additionally evaluate the four ZFO models with the bound fixed to β = 1 on Qwen-2.5-Math-1.5B. Under this setting, ZFO may shorten or reshape the AdamW update but cannot select a step beyond the endpoint corresponding to the nominal FO displacement. ARC-Challenge and StrategyQA use three-epoch versions of configuration a, while ASDiv uses a three-epoch version of configuration c; we denote these shortened configurations by a’ and c’, respectively. Their batch size, group size, context length, and maximum output length otherwise remain unchanged. All methods use the fixed RLVR settings in Table 4, including perturbation scale ϵ = 10−3 , and results are averaged over seeds {100, 101, 102}. Table 12 reports the resulting final evaluation accuracy. Table 12 shows that ZFO can modify training behavior even when it is not permitted to move beyond the nominal FO update. On ARC-Challenge, all four ZFO variants outperform AdamW at the shared learning rate of 10−5 ; however, reducing the separately tuned FO learning rate to 5 × 10−6 raises its accuracy to 55.12, so the ARC result is best interpreted as ZFO recovering a competitive step scale. On StrategyQA, Taylor3 performs best at 74.87, followed by Taylor2 and Padé2 , whereas Padé3 falls below FO. On ASDiv, FO remains narrowly best, although Taylor2 , Padé2 , and Padé3 remain within 0.26 points. Thus, β = 1 is sufficient for gains in some settings, but we argue that the ZFO methods’ largest gains arise from the ability to safely navigate larger learning rates. Location of ZFO Steps. We next examine where the ZFO methods select steps within the search interval It = [0, R]. This diagnostic distinguishes between selecting no update (st = 0), the largest permitted update (st = R), and an interior step. Because numerical solutions may lie arbitrarily close to an endpoint, we report both the exact interior frequency, 0 < st < R, and a 1%-trimmed interior frequency, 0.01 < st /R < 0.99. We also report the median normalized step st /R. Statistics are computed within each run. Table 13 summarizes the β = 1 Family-B cohort, consisting of 48 runs across four configurations, and the fixed-β = 3 cohort, consisting of 108 runs across nine model-dataset pairs. Table 13 reveals distinct step-selection behaviors across the four ZFO models. Taylor3 makes the most frequent use of the interval interior under both bounds. Its exact interior frequency increases from 65.83% at β = 1 to 80.77% at β = 3, while its 1%-trimmed interior frequency rises more substantially from 23.88% to 70.28%. The difference between the two interior measures at β = 1, together with a median normalized step of 0.00156, indicates that many of its nominally interior selections lie close to the zero endpoint. At β = 3, however, most selections remain interior even 34

Table 13: Distribution of ZFO-selected steps over the search interval [0, R]. Frequencies are macro-averaged across runs. Exact and 1%-trimmed interior selections are reported separately, and median steps are normalized by R. For Padé3 , fallback frequency denotes activation of the quadratic safeguard. The largest interior frequency for each value of β is bold. Median st /R

Fallback frequency (%)

57.66 33.31 56.59 20.27

1.00000 0.00156 1.00000 0.00000

0.61 0.67 0.76 12.49

64.59 15.13 72.35 64.36

1.00000 0.02287 1.00000 1.00000

< 0.01

Selected-step frequency (%)

β

Method st = 0

0 < st < R

0.01 < st /R < 0.99

st = R

1

Taylor2 Taylor3 Padé2 Padé3

1.37 0.86 43.41 76.27

40.97 65.83 0.00 3.45

3.76 23.88 0.00 2.55

3

Taylor2 Taylor3 Padé2 Padé3

5.04 4.10 27.65 23.41

30.37 80.77 0.00 12.23

15.87 70.28 0.00 6.93

0.01

< 0.01 54.95

after trimming, showing that Taylor3 frequently chooses steps that are meaningfully separated from both boundaries. The other models rely more heavily on endpoint selections. Taylor2 most often selects the full step st = R, doing so in 57.66% of updates at β = 1 and 64.59% at β = 3. Padé2 never selects an interior point in either cohort and therefore acts empirically as a binary rule between st = 0 and st = R. Padé3 also remains largely endpoint-driven. The median of textttPadé3 shifts from zero at β = 1 to the upper endpoint at β = 3. Its fallback frequency simultaneously increases from 12.49% to 54.95%, indicating that the quadratic safeguard becomes substantially more important over the larger search interval. Overall, the results show that the cubic Taylor model uses higher-order information primarily to select intermediate step lengths, whereas the remaining models more often use their local approximations to choose between suppressing and taking the full permitted update.

D

Notation Table 14: Table of notations. Notation

Definition

∥·∥

Euclidean norm.

πθ (y | x)

Conditional language model with parameters θ, prompt x, and response y.

θt

Model parameters at iteration t.

St

Fixed sampled information used during iteration t, such as prompts, responses, labels, rewards, or other quantities needed to define the current surrogate.

Ft (θ) = F (θ; St )

Fixed-batch surrogate objective optimized during iteration t.

J(θ)

Population or global objective used in the stationarity analysis.

∇Ft (θt )

Gradient of the fixed-batch surrogate at the current parameters.

gt

Update direction produced by the base first-order optimizer.

zt

Normalized search direction, zt = gt /∥gt ∥.

α

Base learning rate used by the first-order optimizer.

β

Bound hyperparameter controlling the ZFO search radius relative to the base optimizer step.

sFO

Scalar length of the standard first-order step, typically sFO = α∥gt ∥.

Rt

Radius of the scalar ZFO search interval. In the simplified theory, Rt = βα∥gt ∥.

It

ZFO search interval, It = [0, Rt ].

ϕt (s)

One-dimensional restriction of the surrogate along the update direction, ϕt (s) = Ft (θt + szt ).

ϕ′t (0)

Directional derivative of Ft along zt , given by ⟨∇Ft (θt ), zt ⟩. Continued on next page

35

Table 14: Table of notations continued. Notation

Definition

ϵ

Probe radius used for the symmetric ZO evaluations ϕt (ϵ) and ϕt (−ϵ).

′′ (0) c ϕ t

Estimate of the second directional derivative.

′′′ (0) ϕc t

Estimate of the third directional derivative.

at , bt , ct

′′ (0), and c = ϕ ′′′ (0). c c at = ϕ′t (0), bt = ϕ t t t

T2 (s)

Quadratic Taylor model of ϕt built from ϕt (0), at , and bt .

T3 (s)

Cubic Taylor model of ϕt built from ϕt (0), at , bt , and ct .

P2 (s)

Second-order Padé model used by Padé2 .

P3 (s)

Third-order Padé model used by Padé3 .

M (s)

Chosen local ZFO model, with M ∈ {T2 , T3 , P2 , P3 }.

sM

Scalar step selected by maximizing the chosen local model over the feasible interval, sM ∈ arg maxs∈It M (s).

sϕ

Best scalar step for the true directional surrogate over the feasible interval, sϕ ∈ arg maxs∈It ϕt (s).

θt+1

ZFO-updated parameters, θt+1 = θt + sM zt .

N (s), D(s)

Numerator and denominator of a Padé approximation.

r(s)

Padé residual, defined as r(s) = D(s)ϕ(s) − N (s).

η

Lower bound on the Padé denominator over the search interval, |D(s)| ≥ η.

M, B, L

Derivative-bound constants used in the Taylor and finite-difference approximation bounds.

Q, P

Higher-order residual-bound constants used in the Padé approximation bounds.

ϕi,t (s)

Per-sample or per-iterate directional surrogate used in the CRN curvature analysis.

ψt (s)

Population counterpart of the directional surrogate.

n

Number of shared samples/trajectories used in CRN training.

σ

Sub-Gaussian parameter.

δ

Failure probability in high-probability concentration bounds.

∆t ¯T ∆

Step-selection error term controlling the gap between the selected ZFO step and the best step over It . ¯ T = 1 P T ∆t . Average step-selection error over T iterations, ∆ t=0 T

LJ , LF̄

Lipschitz constants.

µ

Alignment constant between the optimizer direction and the population objective gradient.

36

Related documents

• Trust the Direction, Search the Step: Zero-and-First-Order Methods for LLM Fine-Tuning#1157538
• AIM-ZO: Activation-Informed Subspace Maintenance for Zeroth-Order LLM Fine-Tuning#1129851
• SubZero+: Efficient Zeroth-Order LLM Fine-Tuning via Large Learning Rates#616813
• SubZero+: Efficient Zeroth-Order LLM Fine-Tuning via Large Learning Rates#616313
• Runtime-Orchestrated Second-Order Optimization for Scalable LLM Training#194280