Descending without gradients — how to minimize a function you can only ask, never differentiate.
Everything you know about optimization assumes you can compute — or at least approximate — a gradient. This lesson is about the awkward world where you can't: the objective is a simulator, a physical experiment, another neural network, or a loss you can only evaluate by running code someone else wrote. Zeroth-order (ZO) optimization — also called derivative-free or gradient-free optimization — is the art of minimizing a function when the only interface available is:
Why would anyone accept such a crippled interface? Because it is surprisingly common. Gradients require structure — a formula, a computation graph, a chain rule. Many objectives in engineering and machine learning have none:
Minimize drag of a wing shape: each evaluation runs a CFD solver that outputs one number. There is no \(\nabla\) of a CFD solver.
Adversarial attacks on a hosted classifier, or fine-tuning a 30-billion-parameter LLM whose gradients would not fit in GPU memory — you can only run the model forward and read the loss.
Tuning controller gains, chemical reactor temperatures, pricing: each "query" is a physical trial whose outcome is one noisy measurement.
Backpropagation needs differentiability. Objectives like "accuracy" or "F1 score" are step functions — you can evaluate them, but \(\nabla\) does not exist.
Optimization methods are usually classified by what they demand from the oracle. This ladder is worth memorizing, because the price of optimization scales directly with how high you can climb:
values only: \(f(x)\)
ZO-SGD, SPSA, CMA-ES, Bayesian optimization, Nelder–Mead
values + gradients: \(\nabla f(x)\)
SGD, Adam, L-BFGS — the workhorses of deep learning
+ curvature: \(\nabla^2 f(x)\)
Newton's method, interior-point, SQP
The guiding question of this lesson: the gradient is just one way to get a descent direction — can we manufacture one from function values alone, and what does it cost? The answer is yes, with a price that grows with dimension. We will build the estimator from scratch, prove its key property in three lines, assemble the algorithm, and see exactly where the price comes from.
Before replacing the gradient, it pays to ask precisely what role the gradient plays. Gradient descent repeatedly steps against the local slope:
where \(\eta > 0\) is the step size. Why does stepping along \(-\nabla f\) make sense? Because the gradient is the direction of steepest ascent: among all unit directions \(d\), the directional derivative \(\nabla f(x) \cdot d\) is minimized by \(d = -\nabla f(x)/\lVert \nabla f(x)\rVert\) — an immediate consequence of the Cauchy–Schwarz inequality, since \(\nabla f \cdot d \ge -\lVert \nabla f\rVert\) with equality exactly when \(d\) points against \(\nabla f\).
Now the reframe that unlocks everything: descent only needs a direction that is positively correlated with \(-\nabla f\) on average. It does not need the exact gradient at every step. Stochastic gradient descent (SGD) already lives in this world — it replaces \(\nabla f\) with a random, noisy \(\hat g\) satisfying \(\mathbb{E}[\hat g] \approx \nabla f\), and it works because the noise averages out over iterations. So our plan is bold but simple:
Start from the most obvious idea, because the sophisticated methods are refinements of it. The derivative is defined as a limit of slopes of secant lines, and a secant only needs two function values — which is exactly what the oracle gives us. Taylor's theorem tells us how good the approximation is: since \(f(x+\varepsilon) = f(x) + \varepsilon f'(x) + \tfrac{\varepsilon^2}{2} f''(\xi)\) for some \(\xi\) between \(x\) and \(x+\varepsilon\), the forward difference satisfies
and averaging the probes on both sides — the central difference — cancels the even-order terms of the two Taylor expansions (they enter with opposite signs while odd terms do not), leaving
The two probe points straddle \(x\) symmetrically, so their secant is nearly parallel to the tangent — the picture that makes this obvious:
In \(d\) dimensions the gradient is the vector of partial derivatives, and each partial derivative is a directional slope along its own coordinate axis \(e_i\). So one full gradient estimate probes every axis:
Why not just take \(\varepsilon \to 0\) and get exact derivatives? Because the oracle is not the Platonic function \(f\) — it is a measurement \(F(x) = f(x) + \xi\) with noise \(\xi\). Differentiating the measurement, the noise is amplified by \(1/\varepsilon\): subtracting two nearly equal numbers divides their random difference by a tiny \(\varepsilon\). Truncation error shrinks as \(\varepsilon\) shrinks; noise error grows. The total error of a forward difference is roughly \(\tfrac{\varepsilon}{2}\lvert f'' \rvert + \tfrac{\sigma}{\varepsilon}\) (with \(\sigma\) the noise standard deviation), which is minimized at \(\varepsilon \approx \sqrt{2\sigma/\lvert f''\rvert}\) — a finite, nonzero probe. Finite-difference optimization is therefore always a bias–variance compromise, never a limit process.
| \(\varepsilon\) | truncation bias (\(\approx\varepsilon\)) | noise error (\(\approx\sigma/\varepsilon\)) | typical total error |
|---|---|---|---|
| 1 | 1 | 0.01 | ≈ 1 |
| 0.1 | 0.1 | 0.1 | ≈ 0.14 ← optimal |
| 0.01 | 0.01 | 1 | ≈ 1 |
| 0.001 | 0.001 | 10 | ≈ 10 |
Forward difference of \(f(x)=x^2\) at \(x=1\), measured with noise \(\sigma=0.01\): true slope 2, estimate \(= 2 + \varepsilon + \Delta\xi/\varepsilon\).
Even with zero noise, coordinate differences have a cost problem that gets worse the more you need them: one gradient estimate costs \(2d\) queries, so a \(d = 10^6\) problem (a modest neural network) needs two million evaluations per iteration. The probes also explore blindly — each coordinate axis is examined separately, learning nothing about how the objective couples variables.
Here is the conceptual leap. The \(2d\)-probe scheme is wasteful because it insists on knowing every partial derivative separately. But we established in §2 that descent only needs some vector correlated with the gradient. So instead of measuring all \(d\) coordinate slopes, measure the slope along one random direction \(u\), then lay that slope back along \(u\):
Multiplying by \(u\) looks odd at first — why not just step along \(u\) with weight = slope? Why is the estimate of a gradient vector allowed to be a product like this? The answer is a beautiful identity, and it is the mathematical heart of this whole subject, so it deserves its own box.
Define the Gaussian smoothing of \(f\) (think: blurring a photo, kernel radius \(\varepsilon\)): \[ f_{\varepsilon}(x) \;=\; \mathbb{E}_{u \sim \mathcal N(0, I_d)}\bigl[\, f(x + \varepsilon u) \,\bigr]. \] Claim: \(\; \mathbb{E}\bigl[\hat g \mid x\bigr] = \nabla f_{\varepsilon}(x)\) — the estimator is an unbiased gradient of the blurred objective.
Why. We need one lemma (Stein's identity): for \(u \sim \mathcal N(0,I_d)\) and any smooth \(h\), \[ \mathbb{E}\bigl[u\, h(u)\bigr] = \mathbb{E}\bigl[\nabla_{\!u} h(u)\bigr], \] because component-wise \(\int u_i\, h\, e^{-\lVert u\rVert^2/2} \mathrm{d}u = \int \partial_i h\; e^{-\lVert u\rVert^2/2} \mathrm{d}u\) — integrate by parts; the boundary terms vanish since the Gaussian decays faster than \(h\) can grow. Now apply it to \(h(u) = \frac{f(x+\varepsilon u) - f(x-\varepsilon u)}{2\varepsilon}\), whose \(u\)-gradient is \(\nabla_{\!u} h = \frac{\varepsilon \nabla f(x+\varepsilon u) + \varepsilon \nabla f(x-\varepsilon u)}{2\varepsilon}\): \[ \mathbb{E}[\hat g] = \mathbb{E}\bigl[\nabla_{\!u} h(u)\bigr] = \mathbb{E}\Bigl[\tfrac{1}{2}\nabla f(x+\varepsilon u) + \tfrac{1}{2}\nabla f(x-\varepsilon u)\Bigr] = \nabla f_{\varepsilon}(x), \] where the last equality is differentiating \(f_\varepsilon(x) = \mathbb E[f(x+\varepsilon u)]\) under the integral (the Gaussian average of gradients is the gradient of the Gaussian average). So each \(\hat g\) is the true gradient of a nearby, slightly blurred problem — and \(\nabla f_\varepsilon \to \nabla f\) as \(\varepsilon \to 0\).
Read that figure again slowly, because it encodes the deepest intuition in this subject: ZO methods never descend the landscape you can see; they descend its blurred twin. The blur costs a little bias (below) and buys everything: one slope measurement per direction, in any dimension, with the same formula.
Two quantities control everything. (i) Blur bias: for \(L\)-smooth \(f\) (\(\lVert \nabla f(x)-\nabla f(y)\rVert \le L \lVert x - y\rVert\)), Nesterov's smoothing theorem gives \(\lvert f_\varepsilon(x) - f(x) \rvert \le \tfrac{L n \varepsilon^2}{2}\) and \(\lVert \nabla f_\varepsilon(x) - \nabla f(x) \rVert \le \tfrac{L (n+3)\,\varepsilon}{2}\) — small \(\varepsilon\) means "nearly optimizing \(f\)", large \(\varepsilon\) means optimizing an increasingly fictional landscape (and it caps how accurate you can ever get at fixed \(\varepsilon\)). (ii) Estimator variance: the same theorem bounds the second moment of the one-point variant \(\hat g_1 = \frac{f(x+\varepsilon u)-f(x)}{\varepsilon} u\): \[ \mathbb{E}\bigl[\lVert \hat g_1 \rVert^2\bigr] \;\le\; 2\,(d+4)\, \lVert \nabla f(x) \rVert^2 \;+\; \frac{\varepsilon^2 L^2 (d+6)^3}{2}. \] Unbiased, yes — but the scatter around the true gradient grows roughly linearly with \(d\) (the \((d+4)\lVert\nabla f\rVert^2\) term). One random direction in \(d\) dimensions is a very noisy summary of a \(d\)-dimensional vector: a single draw of \(u\) is nearly orthogonal to \(\nabla f\) when \(d\) is large, so most draws teach you little and only the average is right.
The estimator above uses two probes per step. A cheaper one-point version reuses the current value: \(\hat g_1 = \frac{f(x+\varepsilon u) - f(x)}{\varepsilon} u\) — one query per step, and the identity in the box holds for it verbatim (Nesterov's original form; that is exactly the \(h\) used there with the \(+\) branch only). When would you prefer paying for the second probe? Whenever measurements are noisy and comparable:
This is the common-random-numbers trick: in machine learning the "noise" \(\xi\) is the random minibatch — evaluate \(\mathcal{L}(\theta + \varepsilon u; B)\) and \(\mathcal{L}(\theta - \varepsilon u; B)\) on the same batch \(B\), and the batch's idiosyncrasies cancel in the difference. Spall's SPSA method (1992) is built entirely on this; MeZO (2023) uses the same pairing. If instead each query's noise is independent (different batch, different experimental day), pairing still helps but only partially.
Assemble the pieces and you get ZO-SGD (also called stochastic gradient descent without gradients, or the random-gradient method). It is refreshingly boring — the entire content of this lesson hides inside line 4:
Input: start x₀ · probe scale ε · step size η · budget K for k = 0, …, K−1: draw u ~ N(0, Id) # one random Gaussian direction F⁺ ← F(xk + εu); F⁻ ← F(xk − εu) # 2 function queries, paired noise ĝ ← ((F⁺ − F⁻) / (2ε)) · u # slope along u, laid back along u xk+1 ← xk − η·ĝ # plain SGD step return the best iterate seen
Three remarks before the theory. (1) The two knobs do different jobs: \(\eta\) is a standard SGD step size, while \(\varepsilon\) controls the bias–variance of the measurement — blur bias \(\sim Ld\varepsilon^2\) favors small \(\varepsilon\), noise amplification \(\sim \sigma/\varepsilon\) (and the \(\varepsilon^2 L^2 d^3\) variance term) favors larger \(\varepsilon\). (2) The step size must shrink with dimension: because \(\mathbb{E}\lVert \hat g\rVert^2\) grows \(\sim d\), the safe step is \(\eta \lesssim \tfrac{1}{(d+4)L}\) — about \(d\times\) smaller than gradient descent's \(\tfrac{1}{L}\). (3) Both knobs get annealed: keep \(\varepsilon\) fixed and you converge to a neighborhood of radius \(\sim Ld\varepsilon^2\) of the optimum, not to the optimum — so shrink \(\varepsilon\) and \(\eta\) together, exactly like a learning-rate schedule.
Now the honest accounting — what does all this cost compared to first-order methods? Convergence statements for nonconvex smooth problems target an ε-stationary point: \(\mathbb{E}\lVert \nabla f(\bar x)\rVert^2 \le \varepsilon_0\) (I switch to \(\varepsilon_0\) here for the target accuracy, to avoid clashing with the probe scale \(\varepsilon\)). The comparison, with everything verified against the primary sources:
| setting | first-order cost | zeroth-order cost | price |
|---|---|---|---|
| smooth convex, deterministic target \(f(\bar x)-f^{\star}\le\varepsilon_0\) |
\(\mathcal O(1/\varepsilon_0)\) | \(\mathcal O(d/\varepsilon_0)\) | factor \(d\) |
| smooth nonconvex, stochastic target \(\mathbb E\lVert\nabla f(\bar x)\rVert^2\le\varepsilon_0\) |
\(\mathcal O(1/\varepsilon_0^2)\) | \(\mathcal O(d/\varepsilon_0^2)\) | factor \(d\) |
| nonsmooth convex, deterministic target \(f(\bar x)-f^{\star}\le\varepsilon_0\) |
\(\mathcal O(1/\varepsilon_0^2)\) | \(\mathcal O(d/\varepsilon_0^2)\) | factor \(d\) |
One sentence summarizes the table: ZO methods need at most a factor \(d\) more iterations than their first-order counterparts — Nesterov & Spokoiny's headline result — and the same factor shows up whether the problem is smooth or not, convex or not. The \(d\) is not an artifact of sloppy analysis, either: lower bounds (Jamieson et al. 2012; Duchi et al. 2015) show that with one-point feedback you genuinely cannot beat \(\Omega(d)\), because each query returns one scalar while the gradient has \(d\) unknowns — you are estimating a \(d\)-dimensional vector from single random projections, and information theory charges \(\Omega(d)\) samples for that. (Two-point feedback is the interesting exception: with paired queries, some settings — e.g. stochastic nonsmooth convex — recover first-order dimension dependence.)
Where does the factor \(d\) come from, mechanically? From the variance term we computed in §4. Each step is an unbiased SGD step, but with \(\mathrm{Var}[\hat g] \sim d \lVert \nabla f\rVert^2\); the trajectory only trusts the average of the estimates, and averaging \(T\) independent estimates shrinks noise like \(1/T\). To wash out a \(d\times\) noisier signal you need \(\sim d\times\) more steps — the same arithmetic that makes Monte Carlo converge as \(1/\sqrt{T}\). Watch it happen on an actual run:
Let's run one full ZO-SGD step with numbers small enough to check on paper. Minimize the (intentionally ill-conditioned) bowl
starting from \(x_0 = (2,\,1)\), with probe scale \(\varepsilon = 0.1\), step size \(\eta = 0.05\), and one draw of the random direction (I use a unit vector to keep the arithmetic tidy — a genuine Gaussian draw has random length, and the averaging works out identically): \(u = (0.6,\, 0.8)\).
| quantity | computation | value |
|---|---|---|
| current point & value | \(f(x_0) = \tfrac12(2^2 + 10\cdot 1^2)\) | \(x_0=(2,1),\ f = 7.000\) |
| probe plus | \(x_0 + \varepsilon u = (2.06,\ 1.08)\), \(f = \tfrac12(2.06^2 + 10\cdot 1.08^2)\) | 7.954 |
| probe minus | \(x_0 - \varepsilon u = (1.94,\ 0.92)\), \(f = \tfrac12(1.94^2 + 10\cdot 0.92^2)\) | 6.114 |
| slope along \(u\) | \((7.954 - 6.114)\,/\,0.2\) — check: \(\nabla f(x_0)\cdot u = 1.2 + 8 = 9.2\) ✓ | 9.20 |
| estimate | \(\hat g = 9.2\cdot(0.6,\ 0.8)\) | \((5.52,\ 7.36)\) |
| true gradient | \(\nabla f(x_0) = (2,\ 10)\); alignment \(\cos\angle(\hat g, \nabla f) = 0.90\) | \((2,\ 10)\) |
| step | \(x_1 = x_0 - 0.05\,\hat g = (2 - 0.276,\ 1 - 0.368)\) | \((1.724,\ 0.632)\) |
| new value | \(f(x_1) = \tfrac12(1.724^2 + 10 \cdot 0.632^2)\) | 3.483 |
One step, two queries, \(7.00 \to 3.48\). Notice the estimate \((5.52, 7.36)\) is not the gradient \((2,10)\) — it overshoots in \(x\) and undershoots in \(y\) — yet it still descends, because what matters is the average over many draws. To see how much a single bad draw can hurt, redo the computation with the mirror direction \(u' = (0.6, -0.8)\): the slope is \(\nabla f\cdot u' = -6.8\), giving \(\hat g' = (-4.08,\ 5.44)\) — the \(x\)-component even points the wrong way. But alignment with the true gradient is still positive (\(\cos = 0.67\)), and the step lands at \((2.204,\ 0.728)\) with \(f = 5.09 < 7\): descent anyway. That is the stochastic-approximation bargain in miniature — every step is locally unreliable, the average is not:
Reading about noise is not the same as feeling it. Below is a live ZO-SGD run on a rippled, nonconvex objective you cannot differentiate into submission:
Every step draws a fresh random direction and makes two measured queries \(F = f + \sigma\xi\). Try these experiments, in order:
amber = ZO-SGD · teal dashed = gradient descent (for reference, it “cheats” by using ∇f)
ZO-SGD is the trunk; decades of work grew branches that attack its weaknesses (variance, conditioning, query budget) or exploit special structure. The map:
ZO-SVRG, SpiderZO: occasionally compute a (noisy) full gradient estimate, then reuse it to correct the random directions — same trick as SVRG in first-order world. This is how the factor \(d\) gets clawed back to \(d^{1/2}\) or better on structured problems.
ZO-Adam (Chen et al. 2019), momentum, averaging: since ZO-SGD is literally SGD, every SGD trick applies — but they now also fight estimator noise, so gains are larger.
The engineering workhorse: directions \(u = \pm 1\) per coordinate (Rademacher instead of Gaussian), 2 queries per gradient regardless of \(d\), paired measurements. Decades of use tuning controllers, queues, and industrial processes.
A population of perturbations = several directional probes per step; "fitness" ranks them (rank-shaped weights), antithetic pairs \(\pm u\) cut variance. Recognize it: randomized smoothing plus resampling, dressed in biology vocabulary (Salimans et al. 2017 made the connection explicit).
Nelder–Mead simplex, pattern/Hooke–Jeeves search: deterministic, derivative-free, no variance to manage — robust for small \(d\), but no complexity theory and poor scaling. The classical pre-randomness era.
A different regime: when each query costs hours (drug assay, hyperparameter search) and the budget is tens-to-hundreds of queries with \(d \lesssim 20\), model \(f\) globally (Gaussian process) and choose each query to be maximally informative — don't descend locally at all.
The most striking modern demonstration that this 1990s machinery matters: MeZO (Malladi et al., NeurIPS 2023) fine-tunes language models with two forward passes per step and no backpropagation. The insight is that for huge models the binding constraint is not time but memory: backprop must cache activations and gradients (\(\approx 12\times\) inference memory), while MeZO's two-point estimator only needs the loss difference \(\mathcal L(\theta + \varepsilon u) - \mathcal L(\theta - \varepsilon u)\) and the direction \(u\) — which it never even stores: it re-derives \(u\) from a saved random seed. Result: a 30-billion-parameter model tuned on a single 80 GB GPU (backprop manages 2.7B), accuracy within ~1% of backprop fine-tuning on most tested tasks, checkpoints of < 0.1 MB, and — since nothing needs derivatives — it directly optimizes non-differentiable objectives like accuracy and F1.
But wait — the theory of §5 says ZO needs \(d\times\) more steps, and \(d = 30\) billion! Why does this work at all? Because that pessimism assumes the worst case \(\mathrm{Var}[\hat g] \sim d \lVert\nabla f\rVert^2\). MeZO's theory shows that near a well-pre-trained checkpoint the loss landscape has low effective rank: the gradients that matter live in a low-dimensional subspace, so the estimator's variance is governed by that intrinsic dimension, not by the parameter count. Lesson: complexity bounds are statements about problem classes; structure can quietly move your problem into a much friendlier class.
| your situation | reach for |
|---|---|
| Autodiff works and memory is fine | First-order methods. ZO is strictly worse — it injects variance for information you could have had for free. |
| Gradients exist but their memory is the wall (huge models) | MeZO / ZO-Adam — trade \(d\times\) steps for inference-level memory. |
| Black-box \(f\), moderate \(d\), queries are cheap-ish, possibly non-differentiable | ZO-SGD / ZO-Adam, evolution strategies (CMA-ES) |
| Noisy system to be tuned robustly (controllers, queues) | SPSA — built for exactly this |
| Each query very expensive, \(d \lesssim 20\), budget ≤ a few hundred | Bayesian optimization |
| Small \(d\), want determinism and robustness | Nelder–Mead, pattern search |
Figures' trajectories and curves are illustrative simulations computed for presentation; the complexity statements are quoted from the references above. This page is fully self-contained: MathJax and its fonts are embedded, and it renders identically offline.