Optimal Stochastic Bilevel Optimization with First-Order Oracles

1 Introduction

• With second-order information, O(ϵ^-4) oracle complexity is achievable; standard-smoothness first-order methods achieve O(ϵ^-6).Higher-order smoothness in y improves a prior first-order upper bound to eO(ϵ^-4-2/p), interpolating between those rates.Source: §1 · pp. 1–2
• Existing lower bounds do not characterize how first-order stochastic bilevel complexity depends on p, leaving the 4+2/p exponent’s optimality unresolved.The motivating question is the optimal complexity under a first-order oracle and how higher-order smoothness in the lower variable changes it.Source: §1 · p. 2
• MRT-FD is a single-loop first-order method requiring order-p smoothness only for g in the lower variable, for any p ≥ 1.It tracks the lower-level solution and auxiliary response while updating the upper-level variable, using order-p finite differences for second-order derivative actions.Source: §1 · p. 2
• The lower bound is Ω(ϵ^-4-2/p), matching MRT-FD’s fixed-p upper dependence and establishing the ε exponent as an oracle barrier.The construction embeds a hard nonconvex chain into a bilevel problem using a sinusoidal coupling and scalar lower-level variable.Source: §1 · p. 2

1.1 Related Work

• Prior first-order approaches include value-function and penalty reformulations; SGHA achieves O(ϵ^-6) but requires additional boundedness assumptions.Classical bilevel methods such as stocBiO and SUSTAIN use Hessian- and Jacobian-vector products, unlike MRT-FD’s first-order oracle.Source: §1.1 · pp. 2–3
• FdeHBO is structurally closest to MRT-FD, but relies on stronger stochastic assumptions and tracks a different auxiliary error.FdeHBO reuses samples and uses variance reduction; MRT-FD uses fresh independent samples and controls distance to z⋆(x).Source: §1.1 · p. 3
• Chen et al. (2026b) obtain eO(ϵ^-4-2/p) with smoothness on both f and g in y; this paper imposes the additional order-p condition only on g.The paper’s lower bound complements prior bounds in other bilevel settings by characterizing the effect of higher-order smoothness in the lower variable.Source: §1.1 · p. 3

2 Problem Setting and Oracle Model

• The goal is an approximately stationary point of the hyperobjective Φ, with ε-stationarity defined by E∥∇Φ(xR)∥² ≤ ε².The lower-level objective is strongly convex, while the hyperobjective may be nonconvex.Source: §2 · p. 3

2.1 Assumptions

• Base assumptions include a jointly Lipschitz gradient for f, strong convexity of g(x,·), Lipschitz derivative conditions, bounded ∥∇yf∥, and a lower-bounded Φ.The setup is smooth nonconvex–strongly-convex bilevel optimization.Source: §2.1 · p. 4
• The auxiliary updates are projected onto a ball whose radius ensures both the current-point and target linear-system solutions lie within it.For analysis, minibatch estimates average m fresh oracle calls at a fixed query point.Source: §2.1 · p. 4, §2.2 · p. 5
• For fixed p ≥ 1, order-p smoothness is imposed on gradients of g with respect to y, with no separate Lipschitz conditions required at intermediate derivative orders.The additional higher-order condition concerns g only; the p=1 case follows from the base regularity condition.Source: §2.1 · p. 4

2.2 Stochastic gradient oracle

• Each oracle call estimates the x- and y-gradients of both f and g; the estimates are unbiased with bounded variance.Calls use fresh randomness and are conditionally independent across query points fixed by a common history, but estimates within one call need not be independent.Source: §2.2 · pp. 4–5

3.1 Estimating matrix–vector products from first-order queries

• Order-p interpolation stencils combine gradients along direction z to estimate Hz and Jz from first-order queries, using at most p+1 query locations.The stencil moment conditions preserve the first-order Taylor term while canceling polynomial terms of orders 0 and 2 through p.Source: §3.1 · p. 5

3.2 Algorithm

• MRT-FD updates x, y, and z simultaneously, tracking the lower solution and auxiliary targets instead of solving inner problems to high accuracy at each outer step.The y step tracks y⋆(x), while the auxiliary update tracks the linear-system solution that approaches z⋆(x) as y approaches y⋆(x).Source: §3.2 · pp. 6–7

3.3 Convergence guarantee

• A higher smoothness order permits a larger finite-difference radius at the same bias, reducing variance amplification and the required minibatch size.The bias scales as τ^(2p), while the variance amplification scales as 1/τ².Source: §3.1 · p. 6, §3.3 · p. 7
• For p=1, the rate is O(eps^-6); for p≥2, it matches Chen et al. (2026b)'s O(eps^-4-2/p) dependence without the logarithmic factor.Under the theorem's parameter choices, the output also satisfies E||∇Φ(x_R)||² ≤ eps² for any fixed finite p.Source: §3.3 · pp. 7–8

4.1 A stochastic nonconvex minimization hard instance

• The hard chain exposes at most one new coordinate per gradient along prefix-supported queries, so reaching its last coordinate is necessary for a small gradient.The random rotation hides chain directions from adaptive algorithms.Source: §4.1 · pp. 8–9
• The comparison oracle strengthens the stochastic-gradient oracle, and its extra scalar response supports simulation of the lower-gradient oracle in the bilevel construction.Subtracting its third component from its first recovers the stochastic-gradient estimate in one call.Source: §4.1 · p. 9, §4.2 · p. 10

4.2 A bilevel hard instance

• A sinusoidal coupling with a scalar lower-level variable lifts the hard minimization instance to a bilevel problem while keeping the hyperobjective close to the original objective.One comparison-oracle call simulates one primitive bilevel-oracle call, including joint randomness.Source: §4.2 · p. 10
• Choosing the reveal probability ϖ = Θ(δ²ε²) keeps primitive-oracle variance bounded and transfers the minimization lower bound to the bilevel objective.The construction retains the relevant hypergradient information while attenuating the primitive signal.Source: §4.2 · p. 10

4.3 Lower bound and minimax complexity

• For every fixed finite p, the lower bound makes the upper bound tight in its ε dependence, with complexity scaling as ε^-4-2/p.The hard instances have a scalar lower-level variable; their upper-level dimension can be O(ε^-6-2/p log(1/ε)).Source: §4.3 · p. 11

AI use statement

• The authors report using GPT-5.5 for literature and reference checking, drafting and revision, selected computations, and proof checking.They state that the algorithmic design and main convergence proof were developed by the authors.Source: AI use statement · p. 11

A Proof of Theorem 3.1

• The proof fixes z0 = 0 and organizes the guarantee around constants, regularity, update control, and a summed Lyapunov inequality.Its assumptions are 2.1–2.3, with fixed initial x0 and y0.Source: A Proof of Theorem 3.1 · p. 14

A.1 Constants and parameter choices

• The appendix chooses proof constants and step-size parameters from problem parameters and p, independently of ε, batch size, and iteration count.The auxiliary projection radius is Rz = 4q0.Source: §A.1 · p. 14
• The proof combines lower-solution and auxiliary tracking errors into Tk, then adds them to the objective gap to form Vk.An upper bound suffices for the parameter rule; the algorithm need not query the exact lower solution or objective value.Source: §A.1 · pp. 14–15

A.2 Geometry and errors from the lower iterate

• Strong convexity gives a unique lower solution, and the appendix establishes the hypergradient identity and Lipschitz continuity of the hypergradient.The lower solution and auxiliary target have regularity bounds under Assumption 2.1.Source: §A.2 · p. 15

A.3 Finite differences and conditional sampling

• Normalized order-p finite differences have controlled approximation error without requiring a (p + 1)st derivative.Fresh endpoint batches yield centered noise with bounded variance; distinct nodes use distinct calls.Source: §A.3 · pp. 16–17
• Algorithm 1 uses separate fresh batches for center and node averages, with the outer finite difference using zk+1 and the auxiliary finite difference using zk.A final fresh batch supplies the lower-gradient update; the output index is drawn uniformly and independently of oracle calls.Source: §A.3 · p. 17

Here dz

• The appendix verifies finite second moments for the iterates, tracking errors, outer gradient, objective, and potential at every finite iteration.This follows by induction from linearly growing gradients, bounded auxiliary iterates, and uniformly bounded noise moments.Source: Here dz · p. 18

A.4 One-step descent and tracking estimates

• The projected auxiliary update contracts its error at the current outer point under the specified step size.The subsequent estimates control outer-direction error and objective descent.Source: §A.4 · pp. 18–19
• After the outer iterate moves, both lower-solution and auxiliary-target tracking errors are bounded, and their weighted sum yields a tracking inequality.The proof accounts for target motion as well as update noise.Source: §A.4 · pp. 19–20

A.5 Lyapunov descent

• Adding objective descent to the weighted tracking inequality gives a uniform one-step Lyapunov bound under Assumptions 2.1–2.3.The step sizes do not depend on the target accuracy.Source: §A.5 · p. 20

• E∥∇Φ(xR)∥² ≤ ϵ² for the algorithm's output.The finite-time output bound holds for 0 < τ ≤ 1 and positive integers m and K.Source: §A.6 · p. 21
• O(ϵ^-4-2/p) stochastic gradient queries suffice for fixed p, problem parameters, and initial gap.At p = 1 and p = 2, the exponents are 6 and 5, respectively.Source: §A.6 · p. 22
• The constants are not optimized in the condition number or uniformly over p.Source: §A.6 · p. 22

B.1.1 The Carmon–Arjevani chain

• The hard-chain lemma provides numerical constants Δc, L0, and G∞ that exist for every chain length T.The passage describes these constants as numerical and says their precise values are immaterial for scaling.Source: §B.1.1 · p. 22

B.1.2 A smooth decomposition for the gradient oracle

• The oracle decomposition expresses CT as a deterministic base minus a correction, whose Bernoulli-retained form has expectation CT and whose second-order jet comes from at most one chain summand.The correction has uniform bounds independent of chain length.Source: §B.1.2 · pp. 23–24

B.2 Construction and scaling of the smooth pair

• The compressed hard-chain and correction bounds are used to establish properties of the scaled pair uniformly over U, T, and ε.The construction defines a base, correction, and nonconvex objective before stating the scaled-pair proposition.Source: §B.2 · pp. 26–28

B.3 Finite-channel random rotation

• The finite-channel random-rotation argument covers two ambient-vector responses and scalar information, arbitrary adaptive queries and outputs, and uses projection-accumulation and posterior-invariance lemmas.The analysis also includes sequential posterior-stabilizer invariance.Source: §B.3 · pp. 29–31
• Conditioned on the past transcript and current query, the residual component of u_j is rotation-invariant on a sphere of determined radius, or is zero when that radius vanishes.Sequential posterior kernels establish invariance under orthogonal maps fixing the conditioned span pointwise; normalization is used only when the residual radius is positive.Source: §B.3 · pp. 32–33
• For d > n + T, Lemma B.12 bounds the chance that one query guesses a residual direction, using its conditionally uniform spherical law.The current query is measurable under the conditioning, and its projection onto the residual subspace has rank at most one.Source: §B.3 · pp. 33–34
• Under d > N + T, Lemma B.13 bounds the finite-channel zero-chain’s progress without dependence on the number of vector channels.The proof combines a rare-progress estimate with geometric guessing; progress on the relevant event requires an open gate.Source: §B.3 · pp. 34–35
• Corollary B.14 extends the geometric argument to a bounded terminal output measurable from the transcript, treating it as an additional point without another oracle call.The proof leaves the rare-progress estimate unchanged and applies the geometric argument with N + 1 points.Source: §B.3 · p. 35
• The finite-channel bounds also hold when queries use a bounded Borel soft projection and responses undergo known postprocessing independent of U.An algorithm using the actual oracle can be simulated with the canonical oracle, applying the postprocessing to its responses.Source: §B.3 · pp. 35–36

• The construction draws U uniformly from Ortho(d, T), independently of the algorithm seed and oracle draws.The resulting functions are stated to be smooth, with Lipschitz, value, and gap bounds following from (115)–(117).Source: §B.4 · p. 36
• One comparison call at an arbitrary x is simulated by querying (148) at the transformed point, and the compressed final output costs no oracle call.The transformations are jointly continuous and independent of U.Source: §B.4 · p. 36
• The section applies the comparison-oracle results and dimension choice to algorithms using at most Ns calls; c = 1/8 is admissible.Source: §B.4 · p. 36

B.4.1 Variance and the ordinary first-order lower bound

• Choosing ϖ = 2(1 + C²)ϵ² ≤ 1/2 gives a family with variance at most one; on a fixed rotation, E||∇P_U(x_out)||² ≥ 3ϵ² > ϵ², and T/ϖ = Θ(ϵ⁻⁴).The stochastic-gradient law is realized by r_A − r_∇u, with one comparison call per gradient call.Source: §B.4.1 · pp. 36–37

B.5 The bilevel construction

• The bilevel construction takes x ∈ R^d and a scalar lower-level variable y ∈ R, with 0 < ϵ ≤ δ ≤ 1.The appendix states that the dependence on δ is suppressed.Source: §B.5 · p. 37

B.5.1 Primitive oracle: simulation and variance

• Each primitive call returns four conditionally unbiased estimates with joint conditional variance at most Cvarδ²ϵ²/ϖ.One primitive call is a known measurable transformation of one comparison call, reproducing the full joint law including the common Bernoulli draw.Source: §B.5.1 · p. 37

B.5.2 Lower solution and outer gradient

• For |s|/δ≤1/4, Yδ(s) is continuously differentiable; the proof also bounds |Yδ(s)| by |s|.Source: §B.5.2 · p. 38
• The lower-level objective is strongly convex globally, with curvature at least 3/4.Source: §B.5.2 · p. 38, §B.5.3 · p. 39
• If ∥∇PU(x)∥≥2ϵ, then ∥∇ΦU(x)∥≥2ϵ−C∇ϵ^3≥3ϵ/2 for sufficiently small ϵ.The bound holds uniformly over ϵ≤δ≤1.Source: §B.5.2 · p. 39
• The mixed-derivative term in the hypergradient retains the correction ∇uU even though the primitive signal is smaller.Source: §B.5.2 · p. 39

B.5.3 Base and higher-order regularity

• The constructed family satisfies the base and order-p regularity assumptions with numerical bounds independent of T, d, U, ϵ, and p.Source: §B.5.3 · p. 40
• The admissible primitive oracle satisfies Assumption 2.3 with σ2≤1/2.Source: §B.5.3 · p. 40

• The construction sets δ=ϵ^(1/p) and ϖ=crrδ^2ϵ^2, with T taken from (110).Source: §B.6 · p. 40
• For sufficiently small ϵ0, T≥(cT/2)ϵ^-2 and T/(8ϖ)≥2, uniformly in p.Source: §B.6 · p. 40
• The primitive oracle simulates each algorithm call to Rϖ and preserves the output law, allowing Proposition 4.1's lower bound to be applied.Source: §B.6 · p. 40
• With probability at least 3/4 on a fixed instance, the algorithm's output satisfies ∥∇ΦU(xout)∥≥3ϵ/2.The probability is over the algorithm and oracle randomness.Source: §B.6 · p. 41
• The fixed instance lies in Pp(M0), with a common sufficiently small ϵ0 that is independent of p.Source: §B.6 · p. 41