Stochastic Bilevel Optimization with Heavy-Tailed Noise

In plain words

Some optimization problems are nested: an outer choice, such as a regularization strength or a weight per data source, is judged by how well an inner model does once it has been fitted to that choice. Refitting the inner model for every outer choice is what makes these problems expensive, and the gradient estimates used to do it are noisy. Here the noise is heavy-tailed: now and then a sample produces a gradient that is far larger than usual.

The paper's answer is a two-level loop. The inner fit takes capped steps, so one very large sample cannot throw it off. The outer update uses only the direction of its gradient estimate, not its size, so every outer step has the same length. The same idea is applied to a related two-player problem, where one variable is maximized and the other minimized.

The work is theoretical. It bounds how many noisy gradient calls the method needs to reach an approximate stationary point, and it reports three experiments. The bounds match the best known results in the bounded-variance case, and they are stated for heavy-tailed noise too.

TL;DR

• Problem: smooth bilevel optimization with a strongly convex lower level, where the stochastic gradient noise has only a finite p-th central moment, p in (1, 2].
• Method: N2SBA, a nested loop with clipped SGD on the inner subproblems and a normalized outer gradient step on a penalized surrogate §3.1.
• Result: an approximate stationary point with SFO complexity O~(kappa^((7p-3)/(p-1)) sigma^(p/(p-1)) eps^(-(4p-2)/(p-1))) in expectation and with high probability; at p = 2 this is O~(kappa^11 sigma^2 eps^-6) §3.2 tab. 1.
• Minimax variant N2SGDA reaches O~(kappa^((2p-1)/(p-1)) sigma^(p/(p-1)) eps^(-(3p-2)/(p-1))) §4.
• Status: theory-first; experiments are small and shown only as figures §5.

Problem

Bilevel problem: minimize phi(x) = f(x, y*(x)) where y*(x) is the unique minimizer of a strongly convex lower-level function g(x, .) §1 eq. 1. Both f and g are available only through stochastic gradients eq. 2. The paper targets the nonconvex-strongly-convex case, with noise assumed to have a finite p-th central moment (the p-BCM assumption §2) instead of finite variance. The authors state that existing stochastic bilevel methods analyzed under bounded variance do not carry over to this noise model §1 §3.1.

Methods that avoid Hessian-vector products build a penalized surrogate L*_lambda(x) and run stochastic gradient descent on it eq. 3 (Kwon et al. 2023). Standard SGD-based analyses need bounded variance; the authors argue that heavy-tailed noise breaks this §3.1.

Method

Penalized reformulation. For a penalty lambda > 0 the surrogate is the minimum over y of f(x, y) + lambda (g(x, y) - min_z g(x, z)) eq. 3. Its gradient needs only first-order information at two points: y*_lambda(x) and y*(x) §3.1. The paper notes that lambda = Theta(kappa^3 ell / eps) makes an eps-stationary point of the surrogate an O(eps)-stationary point of phi §3.1.

Outer and inner loops alg. 1:

1. Inner loops, one for y and one for z, run clipped SGD on the strongly convex subproblems f + lambda g (for y) and lambda g (for z) §3.1 eq. 4. The inner iterates are warm-started from the previous outer step alg. 1.
2. The outer gradient estimate g_x averages M mini-batch samples of the penalized gradient evaluated at (x_t, y_t) and (x_t, z_t) §3.1.
3. The outer step is normalized: x_{t+1} = x_t - eta_x g_x / ||g_x|| alg. 1. Only the direction of g_x is used.

Clipping operator §3: clip(g, tau) rescales g to norm at most tau, and returns 0 when g = 0.

# clip(g, tau): rescale g so its norm is at most tau  [§3]
def clip(vector, level):
    size = norm(vector)
    if size == 0:
        return zeros_like(vector)
    return min(1.0, level / size) * vector

Algorithm 1, N2SBA alg. 1:

# N2SBA: nested-loop normalized stochastic bilevel approximation  [alg. 1]
# Inputs: x0, y0 (also used as z0), lam, eta_x, eta_y[t], eta_z[t],
#         clip_y[t][k], clip_z[t][k], T (outer), K (inner), M (batch)
# Output: x_hat, drawn uniformly from x_1 .. x_T  [alg. 1]

def n2sba(x0, y0, lam, eta_x, eta_y, eta_z, clip_y, clip_z, T, K, M):
    # Inner iterates carry over between outer steps  [alg. 1, l. 330-331]
    y_hat = y0
    z_hat = y0
    x = x0
    iterates = []
    for t in range(T):
        for k in range(K):
            # Draw xi'_t, zeta'_t and zeta_t  [alg. 1, l. 324]
            f_sample = sample_f()
            g_sample_prime = sample_g()
            g_sample = sample_g()
            # Clipped step toward y*_lambda(x_t)  [alg. 1, l. 325-326]
            grad_y = (
                grad_y_F(x, y_hat, f_sample)
                + lam * grad_y_G(x, y_hat, g_sample_prime)
            )
            y_hat = y_hat - eta_y[t] * clip(grad_y, clip_y[t][k])
            # Clipped step toward y*(x_t)  [alg. 1, l. 327-328]
            grad_z = lam * grad_y_G(x, z_hat, g_sample)
            z_hat = z_hat - eta_z[t] * clip(grad_z, clip_z[t][k])
        # Warm start for the next outer step  [alg. 1, l. 330-331]
        y = y_hat
        z = z_hat
        # Mini-batch estimate of the surrogate gradient  [alg. 1, l. 332-333]
        grad_x = 0
        for i in range(M):
            f_sample = sample_f()
            g_sample = sample_g()
            grad_x = grad_x + (
                grad_x_F(x, y, f_sample)
                + lam * (grad_x_G(x, y, g_sample)
                         - grad_x_G(x, z, g_sample))
            ) / M
        # Normalized outer step  [alg. 1, l. 334]
        x = x - eta_x * grad_x / norm(grad_x)
        iterates.append(x)
    return uniform_choice(iterates)

Note on the algorithm's line 333: the same zeta_{t,i} sample appears in both x-gradient terms. The code above keeps that.

Algorithm 2, N2SGDA for minimax alg. 2. Setting G = -F and lambda = 0 reduces Algorithm 1 to a single ascent loop on y and a normalized outer step on x §4:

# N2SGDA: nested-loop normalized gradient descent ascent  [alg. 2]
# Inputs: x0, y0, eta_x, eta_y[t], clip_levels[t][k], T, K, M
# Output: x_hat, drawn uniformly from x_1 .. x_T  [alg. 2]

def n2sgda(x0, y0, eta_x, eta_y, clip_levels, T, K, M):
    # Inner ascent iterate carries over between outer steps  [alg. 2, l. 679]
    y_hat = y0
    x = x0
    iterates = []
    for t in range(T):
        for k in range(K):
            # Draw xi'_t and take a clipped ascent step  [alg. 2, l. 675-677]
            grad_y = grad_y_F(x, y_hat, sample_f())
            y_hat = y_hat + eta_y[t] * clip(grad_y, clip_levels[t][k])
        # Warm start for the next outer step  [alg. 2, l. 679]
        y = y_hat
        # Mini-batch x-gradient at (x, y)  [alg. 2, l. 680-681]
        grad_x = 0
        for i in range(M):
            grad_x = grad_x + grad_x_F(x, y, sample_f()) / M
        # Normalized outer step  [alg. 2, l. 682]
        x = x - eta_x * grad_x / norm(grad_x)
        iterates.append(x)
    return uniform_choice(iterates)

Key equations, if essential:

\nabla \mathcal{L}^*_\lambda(x) = \nabla_x f(x, y^*_\lambda(x)) + \lambda \big(\nabla_x g(x, y^*_\lambda(x)) - \nabla_x g(x, y^*(x))\big) §3.1

Lemma 3.1 bounds the average gradient norm of the surrogate by a descent term, a noise term of size (4 sigma_f + 8 lambda sigma_g) / M^((p-1)/p), and two inner-solver error terms weighted by lambda L_g §3.2, lemma 3.1. Lemma 3.2 gives the clipped-SGD guarantee for a strongly convex subproblem with heavy-tailed noise. Its step size is constant, and its clip levels decay in k §3.2, lemma 3.2.

Hyperparameters

Theorems 3.1 and 4.1 give parameter settings only up to order (tilde-Theta, tilde-O, O). The paper does not give numeric constants for them.

Not stated anywhere in the paper: the constants hidden by tilde-Theta, the values chosen for step sizes and lambda in the experiments, the number of outer iterations for hyper-cleaning, the per-sample scale of the Pareto noise, and the train/validation sizes for learnable regularization.

Results

Author-reported. Theoretical upper bounds on SFO calls (no experiment reports these numbers) and qualitative experimental results. Baselines and setups are as stated in the paper.

Complexity, Table 1 (bilevel, nonconvex-strongly-convex) tab. 1 §3.2:

Complexity, Table 2 (minimax, nonconvex-strongly-concave) tab. 2 §4:

Experiments (figures only; no numeric tables) §5:

Paper vs code

No official code.

Searched (in order): the paper body and LaTeX source for code or project links (none); GitHub repository search for "N2SBA", "N2SGDA" and "bilevel heavy-tailed noise" (no results); GitHub user search for "Zhuanghua Liu" (no results). A user search for "Luo Luo" returned only unrelated accounts, none identifiable as an author. No repository was found, so there is no code field.

Pitfalls

• Constants are hidden. Theorems 3.1 to 4.2 give step sizes, clip levels, K, M and lambda only up to tilde-Theta, tilde-O and O. The compiled paper gives no numeric recipe. The LaTeX source contains commented-out explicit constants (for example 400, 120 and 5400 in the Lemma 3.2 and Theorem 3.2 drafts, main.tex lines 482-487 and 563-573). They are not in the paper, so do not treat them as its values.
• R is undefined in Theorem 3.2. The penalty is lambda = max{kappa/sqrt(R), ...} and the text says R_y, R_z, R_0 and Delta follow Theorem 3.1, but no R is defined. Theorem 3.1 uses kappa/R_0, so reading R as R_0^2 makes the two theorems agree. This is inferred, not stated.
• Lemma 3.1 needs lambda >= 2 L_f / mu. The theorems set lambda through a max that is not checked against this condition. Verify it numerically.
• Output set is off by one. Algorithm 1 draws x_hat from x_1 .. x_T. The loop in the algorithm starts from x_0, and Lemma 3.1 and the high-probability theorem average over t = 0 .. T-1. Pick one convention and state it.
• Theorems differ in guarantee. Theorem 3.1 bounds E||grad phi(x_hat)|| for a random output. Theorem 3.2 bounds the average of ||grad phi(x_t)|| over iterates with probability 1 - delta. The step sizes change too: eps/(ell kappa^3) versus sqrt(Delta / (ell kappa^3 T)). The high-probability K drops the second term of Theorem 3.1, and the high-probability tau' has an extra factor 120 that the expectation version lacks. Use the theorem that matches the guarantee you need.
• Minimax R_y. The expectation theorem gives R_y^2 = O(R_0^2 + eps^2/ell^2) without constants; the high-probability theorem gives R_0^2 + 10 eps^2/ell^2.
• Step size scales with 1/lambda. The inner caps 1/(L_f + lambda L_g) and 1/(lambda mu K) shrink as lambda grows, and lambda ~ kappa^3 ell / eps grows as eps shrinks. Expect inner loops to need small steps for small eps. This is a consequence of the stated formulas, not an experiment.
• Algorithm numerics. The normalized update divides by ||g_x||, which the algorithm leaves undefined at g_x = 0. Guard it in code. The clip operator has an explicit zero case. The main text does not discuss the bias that clipping introduces.
• SFO counting. The complexity bounds count SFO calls, but the paper does not give a per-iteration count or say whether a gradient of G counts as one call or two. Algorithm 1 draws three samples per inner step (xi', zeta', zeta) and two per outer sample (xi, zeta). Work out the total yourself before comparing against baselines.
• Symbol collision. p is the moment order in the theory, but in Section 5 p also appears as the feature dimension (p = 101,631) and as the parameter vector size in hyper-cleaning. Keep the two apart.
• Synthetic setup is inconsistent. The text takes y as a scalar in R but defines K in R^{d x d}, while x^T K y is scalar. d_x = 30 is stated, but d is not explicitly set equal to d_x. Choose K so that x^T K y is scalar.
• Clipped baseline ambiguity. Section 5.1 names both "SGDmax" and "Clipped-SGDmax" and states T for each, but only SGDmax is defined as a baseline.
• Figure label mismatch. Figure 3's subpanel labels use p for the corruption ratio, while the caption and text use q.
• F2BSA constant. Table 1 lists kappa^11 for F2BSA. That is the authors' improvement of the published kappa^12 (footnote, appendix E), not the original number.
• Lower bound. Theorem 3.3 holds only for quadratic lower-level functions and requires C_f >= sqrt(13 L_f Delta) and eps <= min(sigma kappa_g/92, kappa_g sqrt(L_f Delta)/312). Its kappa_g exponent (3p-2)/(p-1) differs from the upper bound's kappa exponent (7p-3)/(p-1), and the paper leaves this gap open.
• Experiments are small and unreported numerically. Results appear only as plots of SFO calls against a metric. No table of final values and no code make them reproducible.

Lineage

• Builds on: the penalized reformulation and fully first-order approach of Kwon et al. 2023 §3.1; the p-BCM noise assumption and clipping analysis of Zhang et al. 2020 §2 §3.1; the normalized SGD analysis of heavy-tailed nonconvex optimization referenced as liunonconvex2025, used for the minimax lower bound §4.
• Compared against: F2BSA (Chen et al. 2025) in the bounded-variance table tab. 1 and in experiments §5.2 §5.3; stocBiO (Ji et al. 2021) and VRBO (Yang et al. 2021) in experiments §5.2; SGDA (Lin et al. 2020) and SGDmax (Jin et al. 2020) for minimax tab. 2 §5.1.
• Superseded by: not known.