Provably Faster Algorithms for Bilevel Optimization
Year: 2022

Abstract:
Bilevel optimization has been widely applied in many important machine learning applications such as hyperparameter optimization and meta-learning. Recently, several momentum-based algorithms have been proposed to solve bilevel optimization problems faster. However, those momentum-based algorithms do not achieve provably better computational complexity than Õ( −2) of the SGD-based algorithm. In this paper, we propose two new algorithms for bilevel optimization, where the first algorithm adopts momentum-based recursive iterations, and the second algorithm adopts recursive gradient estimations in nested loops to decrease the variance. We show that both algorithms achieve the complexity of Õ( −1.5), which outperforms all existing algorithms by the order of magnitude. Our experiments validate our theoretical results and demonstrate the superior empirical performance of our algorithms in hyperparameter applications.

heading: 1 Introduction
text: Bilevel optimization has become a timely and important topic recently due to its great effectiveness in a wide range of applications including hyperparameter optimization [7, 5], meta-learning [33, 16, 1], reinforcement learning [14, 24]. Bilevel optimization can be generally formulated as the following minimization problem:
min x∈Rp Φ(x) := f(x, y∗(x)) s.t. y∗(x) = arg min y∈Rq g(x, y). (1)
Since the outer function Φ(x) := f(x, y∗(x)) depends on the variable x also via the optimizer y∗(x) of the inner-loop function g(x, y), the algorithm design for bilevel optimization is much more complicated and challenging than minimization and minimax optimization. For example, if the gradient-based approach is applied, then the gradient of the outer-loop function (also called hypergradient) will necessarily involve Jacobian and Hessian matrices of the inner-loop function g(x, y), which require more careful design to avoid high computational complexity.
This paper focuses on the nonconvex-strongly-convex setting, where the outer function f(x, y∗(x)) is nonconvex with respect to (w.r.t.) x and the inner function g(x, y) is strongly convex w.r.t. y for any x. Such a case often occurs in practical applications. For example, in hyperparameter optimization [7], f(x, y∗(x)) is often nonconvex with x representing neural network hyperparameters, but the inner function g(x, ·) can be strongly convex w.r.t. y by including a strongly-convex regularizer on y. In few-shot meta-learning [1], the inner function g(x, ·) often takes a quadratic form together with a strongly-convex regularizer. To efficiently solve the deterministic problem in eq. (1), various bilevel optimization algorithms have been proposed, which include two popular classes of deterministic gradient-based methods respectively based on approximate implicit differentiation (AID) [31, 9, 8] and iterative differentiation (ITD) [28, 6, 7].
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
Recently, stochastic bilevel opitimizers [8, 20] have been proposed, in order to achieve better efficiency than deterministic methods for large-scale scenarios where the data size is large or vast fresh data needs to be sampled as the algorithm runs.
In particular, such a class of problems adopt functions by:
Φ(x) := f(x, y∗(x)) := Eξ[F (x, y∗(x); ξ)], g(x, y) := Eζ [G(x, y; ζ)]
where the outer and inner functions take the expected values w.r.t. samples ξ and ζ, respectively.
Along this direction, [20] proposed a stochastic gradient descent (SGD) type optimizer (stocBiO), and showed that stocBiO attains a computational complexity of Õ( −2) in order to reach an - accurate stationary point. More recently, several studies [2, 11, 22] have tried to accelerate SGD-type bilevel optimizers via momentum-based techniques, e.g., by introducing a momentum (historical information) term into the gradient estimation. All of these optimizers follow a single-loop design, i.e., updating x and y simultaneously. Specifically, [22] proposed an algorithm MSTSA by updating x via a momentum-based recursive technique introduced by [3, 36]. [11] proposed an optimizer SEMA similarly to MSTSA but using the momentum recursive technique for updating both x and y. [2] proposed an algorithm STABLE, which applies the momentum strategy for updating the Hessian matrix, but the algorithm involves expensive Hessian inverse computation rather than hypergradient approximation loop. However, as shown in Table 1, SEMA, MSTSA and STABLE achieve the same complexity order of Õ( −2) as the SGD-type stocBiO algorithm, where the momentum technique in these algorithms does not exhibit the theoretical advantage. Such a comparison is not consistent with those in minimization [3] and minimax optimization [15], where the single-loop momentum-based recursive technique achieves provable performance improvements over SGD-type methods. This motivates the following natural but important question:
• Can we design a faster single-loop momentum-based recursive bilevel optimizer, which achieves order-wisely lower computational complexity than SGD-type stocBiO (and all other momentumbased algorithms), and is also easy to implement with efficient matrix-vector products?
Although the existing theoretical efforts on accelerating bilevel optimization algorithms have been exclusively focused on single-loop design1, empirical studies in [20] suggested that double-loop bilevel algorithms such as BSA [8] and stocBiO [20] achieve much better performances than single-loop algorithms such as TTSA [14]. A good candidate suitable for accelerating double-loop algorithms can be the popular variance reduction method, such as SVRG [21], SARAH [30] and SPIDER [4], which typically yield provably lower complexity. The basic idea is to construct low-variance gradient estimators using periodic high-accurate large-batch gradient evaluations. So far, there has not been any study on using variance reduction to accelerate double-loop bilevel optimization algorithms. This motivates the second question that we address in this paper:
• Can we develop a double-loop variance-reduced bilevel optimizer with improved computational complexity over SGD-type stocBiO (and all other existing algorithms)? If so, when such a double-loop algorithm holds advantage over the single-loop algorithms in bilevel optimization?

heading: 1.1 Main Contributions
text: This paper proposes two algorithms for bilevel optimization, both outperforming all existing algorithms in terms of complexity order.
We first propose a single-loop momentum-based recursive bilevel optimizer (MRBO). MRBO updates variables x and y simultaneously, and uses the momentum recursive technique for constructing lowvariance mini-batch estimators for both the gradient ∇g(x, ·) and the hypergradient ∇Φ(·); in contrast to previous momentum-based algorithms that accelerate only one gradient or neither. Further, MRBO is easy to implement, and allows efficient computations of Jacobian- and Hessian-vector products via automatic differentiation. Theoretically, we show that MRBO achieves a computational complexity (w.r.t. computations of gradient, Jacobian- and Hessian-vector product) of Õ( −1.5),
1In the literature of bilevel optimization, although many hypergradient-based algorithms include an iteration loop of Hessian inverse estimation, such a loop is typically not counted when these algorithms are classified by the number of loops. This paper follows such a convention to be consistent with the existing literature. Namely, the single- and double-loop algorithms mentioned here can include an additional loop of Hessian inverse estimation in the hypergradient approximation.
which outperforms all existing algorithms by an order of −0.5. Technically, our analysis needs to first characterize the estimation property for the momentum-based recursive estimator for the Hessian-vector type hypergradient and then uses such a property to further bound the per-iteration error due to momentum updates for both inner and outer loops.
We then propose a double-loop variance-reduced bilevel optimizer (VRBO), which is the first algorithm that adopts the recursive variance reduction for bilevel optimization. In VRBO, each inner loop constructs a variance-reduced gradient (w.r.t. y) and hypergradient (w.r.t. x) estimators through the use of large-batch gradient estimations computed periodically at each outer loop. Similarly to MRBO, VRBO involves the computations of Jacobian- and Hessian-vector products rather than Hessians or Hessian inverse. Theoretically, we show that VRBO achieves the same near-optimal complexity of Õ( −1.5) as MRBO and outperforms all existing algorithms. Technically, differently from the use of variance reduction in minimization and minimax optimization, our analysis for VRBO needs to characterize the variance reduction property for the Hessian-vector type of hypergradient estimators, which only involves Hessian vector computation rather than Hessian. Such estimator introduces additional errors to handle in the telescoping and convergence analysis.
Our experiments2 show that VRBO achieves the highest accuracy among all comparison algorithms, and MRBO converges fastest among its same type of single-loop momentum-based algorithms. In particular, we find that our double-loop VRBO algorithm converges much faster than other singlrloop algorithms including our MRBO, which is in contrast to the existing efforts exclusively on accelerating the single-loop algorithms [2, 11, 22]. Such a result also differs from those phenomenons observed in minimization and minimax optimization, where single-loop algorithms often outperform double-loop algorithms.

heading: 1.2 Related Works
text: Bilevel optimization approaches: At the early stage of bilevel optimization studies, a class of constraint-based algorithms [13, 35, 29] were proposed, which tried to penalize the outer function with the optimality conditions of the inner problem. To further simplify the implementation of constraint-based bilevel methods, gradient-based bilevel algorithms were then proposed, which include but not limited to AID-based [33, 7, 34, 17], ITD-based [9, 31, 8, 16, 19] methods, and stochastic bilevel optimizers such as BSA [8], stocBiO [20], and TTSA [14]. The finite-time (i.e., non-asymptotic) convergence analysis for bilevel optimization has been recently studied in several works [8, 20, 14]. In this paper, we propose two novel stochastic bilevel algorithms using momentum recursive and variance reduction techniques, and show that they order-wise improve the computational complexity over existing stochastic bilevel optimizers.
Momentum-based recursive approaches: The momentum recursive technique was first introduced by [3, 36] for minimization problems, and has been shown to achieve improved computational complexity over SGD-based updates in theory and in practice. Several works [22, 2, 11] applied the
2Our codes are available online at https://github.com/JunjieYang97/MRVRBO
similar single-loop momentum-based strategy to bilevel optimization to accelerate the SGD-based bilevel algorithms such as BSA [8] and stocBiO [20]. However, the computational complexities of these momentum-based algorithms are not shown to outperform that of stocBiO. In this paper, we propose a new single-loop momentum-based recursive bilevel optimizer (MRBO), which we show achieves order-wisely lower complexity than existing stochastic bilevel optimizers.
Variance reduction approaches: Variance reduction has been studied extensively for conventional minimization problems, and many algorithms have been designed along this line, including but not limited to SVRG [21, 26], SARAH [30], SPIDER [4], SpiderBoost [37, 38, 18] and SNVRG [41]. Several works [27, 39, 40, 32] recently employed such techniques for minimax optimization to achieve better complexities. In this paper, we propose the first-known variance reduction-based bilevel optimizer (VRBO), which achieves a near-optimal computational complexity and outperforms existing stochastic bilevel algorithms.
Two concurrent works: As we were finalizing this submission, two concurrent studies were posted on arXiv recently ([23] was posted on May 8 and [12] was posted on May 5). Both studies overlap only with our MRBO algorithm, nothing similar to our VRBO. Specifically, [23] and [12] respectively proposed the SUSTAIN and RSVRB algorithms for bilevel optimization, both using momentumbased design as our MRBO. Although SUSTAIN and RSVRB have been shown to achieve the same theoretical complexity of O( −1.5) as our MRBO (and VRBO), both algorithms have major drawbacks in their design, so that their empirical performance (as we demonstrate in our experiments) is much worse that our MRBO (and even worse than our VRBO). SUSTAIN adopts only singlesample for each update (whereas MRBO uses minibatch for stability); and RSVRB requires to compute Hessian inverse at each iteration (whereas MRBO uses Hessian-vector products for fast computation). As an additional note, our experiments demonstrate that our VRBO significantly outperforms all these single-loop momentum-based algorithms SUSTAIN and RSVRB as well as our MRBO.

heading: 2 Two New Algorithms
text: In this section, we propose two new algorithms for bilevel optimization. Firstly, we introduce the hypergradient of the objective function Φ(xk), which is useful for designing stochastic algorithms. Property 1. The (hyper)gradient of Φ(x) = f(x, y∗(x)) in eq. (1) takes a form of
∇Φ(x) = ∇xf(x, y∗(x))−∇x∇yg(x, y∗(x))[∇2yg(x, y∗(x))]−1∇yf(x, y∗(x)). (2)
However, it is not necessary to compute y∗ for updating x at every iteration, and it is not time and memory efficient to compute Hessian inverse matrix in eq. (2) explicitly. Here, we estimate the hypergradient similarly to [20, 8], which takes a form of
∇Φ(x) = ∇xf(x, y)−∇x∇yg(x, y)η Q−1∑ q=−1 (I − η∇2yg(x, y))q+1∇yf(x, y), (3)
where the Neumann series η ∑∞ i=0(I − ηG)i = G−1 is applied to approximate the Hessian inverse.

heading: 2.1 Momentum-based Recursive Bilevel Optimizer (MRBO)
text: As shown in Algorithm 1, we propose a Momentum-based Recursive Bilevel Optimizer (MRBO) for solving the bilevel problem in eq. (1).
MRBO updates in a single-loop manner, where the momentum recursive technique STORM [3] is employed for updating both x and y at each iteration simultaneously. To update y, at step k, MRBO first constructs the momentum-based gradient estimator uk based on the current ∇yG(xk, yk;By) and the previous∇yG(xk−1, yk−1;By) using a minibatch By of samples (see line 8 in Algorithm 1). Note that the hyperparameter βk decreases at each iteration, so that the gradient estimator uk is more determined by the previous uk−1, which improves the stability of gradient estimation, especially when yk is close to the optimal point. Then MRBO uses the gradient estimator for updating yk (see line 11). The stepsize ηk decreases at each iteration to reduce the convergence error.
To update x, at step k, MRBO first constructs the momentum-based recursive hypergradient estimator vk based on the current ∇̂Φ(xk;Bx) and the previous ∇̂Φ(xk−1;Bx) computed using several
Algorithm 1 Momentum-based Recursive Bilevel Optimizer (MRBO) 1: Input: Stepsize λ, γ > 0, Coefficients α0, β0, Initializers x0, y0, Hessian Estimation Number Q,
Batch Size S, Constant c1, c2,m, d > 0 2: for k = 0, 1, . . . ,K do 3: Draw Samples By,Bx = {Bj(j = 1, . . . , Q),BF ,BG} with batch size S for each component 4: if k = 0: then 5: vk = ∇̂Φ(xk;Bx), uk = ∇yG(xk, yk;By) 6: else 7: vk = ∇̂Φ(xk;Bx) + (1− αk)(vk−1 − ∇̂Φ(xk−1;Bx)) 8: uk = ∇yG(xk, yk;By) + (1− βk)(uk−1 −∇yG(xk−1, yk−1;By)) 9: end if
10: update: ηk = d3√m+k , αk+1 = c1η 2 k, βk+1 = c2η 2 k 11: xk+1 = xk − γηkvk, yk+1 = yk − ληkuk 12: end for
independent minibatches of samples Bx = {Bj(j = 1, . . . , Q),BF ,BG} (see line 7 in Algorithm 1). The hyperparameter αk decreases at each iteration, so that the new gradient estimation vk is more determined by the previous vk−1, which improves the stability of gradient estimation, especially when xk is around the optimal point. Specifically, the hypergradient estimator ∇̂Φ(xk;Bx) is designed based on the expected form in eq. (3), and takes a form of:
∇̂Φ(xk;Bx) = ∇xF (xk, yk;BF )
−∇x∇yG(xk, yk;BG)η Q−1∑ q=−1 Q∏ j=Q−q (I − η∇2yG(xk, yk;Bj))∇yF (xk, yk;BF ), (4)
Note that MRBO computes the above estimator recursively using only Hessian vectors rather than Hessians (see Appendix A) in order to reduce the memory and computational cost. Then MRBO uses the estimated gradient vk for updating xk (see line 11). The stepsize ηk decreases at each iteration to facilitate the convergence.

heading: 2.2 Variance Reduction Bilevel Optimizer (VRBO)
text: Although all of the existing momentum algorithms [2, 22, 11] (and two current studies [23, 12]) for bilevel optimization follow the single-loop design, empirical results in [20] suggest that double-loop bilevel algorithms can achieve much better performances than single-loop algorithms. Thus, as shown in Algorithm 2, we propose a double-loop algorithm called Variance Reduction Bilevel Optimizer (VRBO). VRBO adopts the variance reduction technique in SARAH [30]/SPIDER [4] for bilevel optimization, which is suitable for designing double-loop algorithms. Specifically, VRBO constructs the recursive variance-reduced gradient estimators for updating both x and y, where each update of x in the outer-loop is followed by (m+ 1) inner-loop updates of y. VRBO divides the outer-loop iterations into epochs, and at the beginning of each epoch computes the hypergradient estimator ∇̂Φ(xk, yk;S1) and the gradient ∇yG(xk, yk;S1) based on a relatively large batch S1 of samples for variance reduction, where ∇̂Φ(xk, yk;S1) takes a form of
∇̂Φ(xk, yk;S1) = 1
S1 S1∑ i=1 ( ∇xF (xk, yk; ξi)
−∇x∇yG(xk, yk; ζi)η Q−1∑ q=−1 Q∏ j=Q−q (I − η∇2yG(xk, yk; ζ j i ))∇yF (xk, yk; ξi) ) , (5)
where all samples in S1 = {ζji (j = 1, . . . , Q), ξi, ζi, i = 1, . . . , S1} are independent. Note that eq. (5) takes a different form from MRBO in eq. (4), but the Hessian-vector computation method for MRBO is still applicable here. Then, VRBO recursively updates the gradient estimators for ∇yG(x̃k,t, ỹk,t;S2) and ∇̂Φ(x̃k,t, ỹk,t;S2) (which takes the same form as eq. (5)) with a small sample batch S2 (see lines 11 to 16) during inner-loop iterations.
We remark that VRBO is the first algorithm that adopts the recursive variance reduction method for bilevel optimization. As we will shown in Section 3, VRBO achieves the same nearly-optimal computational complexity as MRBO (and outperforms all other existing algorithms). More interestingly, as a double-loop algorithm, VRBO empirically significantly outperforms all existing single-loop momentum algorithms including MRBO. More details and explanation are provided in Section 4.
Algorithm 2 Variance Reduction Bilevel Optimizer (VRBO) 1: Input: Stepsize β, α > 0, Initializer x0, y0, Hessian Q, Sample Size S1, S2, Periods q 2: for k = 0, 1, . . . ,K do 3: if mod(k, q) = 0: then 4: Draw a batch S1 of i.i.d. samples 5: uk = ∇yG(xk, yk;S1), vk = ∇̂Φ(xk, yk;S1) 6: else 7: uk = ũk−1,m+1, vk = ṽk−1,m+1 8: end if 9: xk+1 = xk − αvk 10: Set x̃k,−1 = xk, ỹk,−1 = yk, x̃k,0 = xk+1, ỹk,0 = yk, ṽk,−1 = vk, ũk,−1 = uk 11: for t = 0, 1, . . . ,m+ 1 do 12: Draw a batch S2 of i.i.d samples 13: ṽk,t = ṽk,t−1 + ∇̂Φ(x̃k,t, ỹk,t;S2)− ∇̂Φ(x̃k,t−1, ỹk,t−1;S2) 14: ũk,t = ũk,t−1 +∇yG(x̃k,t, ỹk,t;S2)−∇yG(x̃k,t−1, ỹk,t−1;S2) 15: x̃k,t+1 = x̃k,t, ỹk,t+1 = ỹk,t − βũk,t 16: end for 17: yk+1 = ỹk,m+1 18: end for

heading: 3 Main Results
text: In this section, we first introduce several standard assumptions for the analysis, and then present the convergence results for the proposed MRBO and VRBO algorithms.

heading: 3.1 Technical Assumptions and Definitions
text: Assumption 1. Assume that the inner function G(x, y; ζ) is µ-strongly-convex w.r.t. y for any ζ.
We then make the following assumptions on the Lipschitzness and bounded variance, as adopted by the existing studies [8, 20, 14] on stochastic bilevel optimization. Assumption 2. Let z := (x, y). Assume the functions F (z; ξ) and G(z; ζ) satisfy, for any ξ and ζ,
a) F (z; ξ) is M -Lipschitz, i.e., for any z, z′, |F (z; ξ)− F (z′; ξ)| ≤M‖z − z′‖.
b) ∇F (z; ξ) and ∇G(z; ζ) are L-Lipschitz, i.e., for any z, z′,
‖∇F (z; ξ)−∇F (z′; ξ)‖ ≤ L‖z − z′‖, ‖∇G(z; ζ)−∇G(z′; ζ)‖ ≤ L‖z − z′‖.
c) ∇x∇yG(z; ζ) is τ -Lipschitz, i.e., for any z, z′, ‖∇x∇yG(z; ζ)−∇x∇yG(z′; ζ)‖ ≤ τ‖z− z′‖.
d) ∇2yG(z; ζ) is ρ-Lipschitz, i.e., for any z, z′, ‖∇2yG(z; ζ)−∇2yG(z′; ζ)‖ ≤ ρ‖z − z′‖.
Note that Assumption 2 also implies that Eξ‖∇F (z; ξ) − ∇f(z)‖2 ≤ M2, Eζ‖∇x∇yG(z; ζ) − ∇x∇yg(z)‖2 ≤ L2 and Eζ‖∇2yG(z; ζ)−∇2yg(z)‖2 ≤ L2. Assumption 3. Assume that∇G(z; ξ) has bounded variance, i.e., Eξ‖∇G(z; ξ)−∇g(z)‖2 ≤ σ2.
Assumptions 2 and 3 require the Lipschitzness conditions to hold for the gradients and secondorder derivatives of the inner and outer objective functions, which further imply the gradient of the outer objective function is bounded. Such assumptions have also been adopted by the existing studies [20, 14, 23, 22, 12] for stochastic bilevel optimization. Furthermore, these assumptions are mild in practice as long as the iterates along practical training paths are bounded. All our experiments
indicate that these iterates are well located in a bounded regime. It is also possible to consider a bilevel problem over a convex compact set, which relaxes the boundedness assumption. By introducing a projection of iterative updates into such a set, our analysis for the unconstrained setting can be extended easily to such a constrained problem.
We next define the -stationary point for a nonconvex function as the convergence criterion. Definition 1. We call x̄ an -stationary point for a function Φ(x) if ‖∇Φ(x̄)‖2 ≤ .

heading: 3.2 Convergence Analysis of MRBO Algorithm
text: To analyze the convergence of MRBO, bilevel optimization presents two major challenges due to the momentum recursive method in MRBO, beyond the previous studies of momentum in conventional minimization and minimax optimization. (a) Outer-loop updates of bilevel optimization use hypergradients, which involve both the first-order gradient and the Hessian-vector product. Thus, the analysis of the momentum recursive estimator for such a hypergradient is much more complicated than that for the vanilla gradient. (b) Since MRBO applies the momentum-based recursive method to both inner- and outer-loop iterations, the analysis needs to capture the interaction between the inner-loop gradient estimator and the outer-loop hypergradient estimator. Below, we will provide two major properties for MRBO, which develop new analysis for handling the above two challenges.
In the following proposition, we characterize the variance bound for the hypergradient estimator in bilevel optimization, and further use such a bound to characterize the variance of the momentum recursive estimator of the hypergradient. Proposition 1. Suppose Assumptions 1, 2 and 3 hold and η < 1L , the hypergradient estimator ∇̂Φ(xk;Bx) w.r.t. x based on a minibatch Bx has bounded variance
E‖∇̂Φ(xk;Bx)−∇Φ(xk)‖2 ≤ G2, (6)
where G2 = 2M 2 S + 12M2L2η2(Q+1)2 S + 4M2L2(Q+2)(Q+1)2η4σ2
S . Further, let ̄k = vk − ∇Φ(xk), where vk denotes the momentum recursive estimator for the hypergradient. Then the per-iteratioon variance bound of vk satisfies
E‖̄k‖2 ≤ E[2α2kG2 + 2(1− αk)2L2Q‖xk − xk−1‖2
+ 2(1− αk)2L2Q‖yk − yk−1‖2 + (1− αk)2‖̄k−1‖2], (7)
where L2Q = 2L 2 + 4τ2η2M2(Q+ 1)2 + 8L4η2(Q+ 1)2 + 2L2η4M2ρ2Q2(Q+ 1)2.
The variance bound G of the hypergradient in eq. (6) scales with the number Q of Neumann series terms (i.e., the number of Hessian vectors) and can be reduced by that minibatch size S.
Then the bound eq. (7) further captures how the variance ‖̄k‖ of momentum recursive hypergradient estimator changes after one iteration. Clearly, the term (1 − αk)2‖̄k−1‖2 indicates a variance reduction per iteration, and the remain three terms capture the impact of the randomness due to the update in step k, including the variance of the stochastic hypergradient estimator G2 (as captured in eq. (6)) and the stochastic update of both variables x and y. In particular, the variance reduction term plays a key role in the performance improvement for MRBO over other existing algorithms. Proposition 2. Suppose Assumptions 1, 2, 3 hold. Let η < 1L and γ ≤ 1 4LΦηk , where LΦ = L+ 2L 2+τM2
µ + ρLM+L3+τML µ2 + ρL2M µ3 . Then, we have
E[Φ(xk+1)] ≤ E[Φ(xk)] + 2ηkγ(L′ 2‖yk − y∗(xk)‖2 + ‖̄k‖2 + C2Q)−
1
2γηk ‖xk+1 − xk‖2,
where CQ = (1−ηµ)Q+1ML µ , L ′2 = max{(L+ L 2 µ + Mτ µ + LMρ µ2 ) 2, L2Q}.
Proposition 2 characterizes how the objective function value decreases (i.e., captured by E[Φ(xk+1)]− E[Φ(xk)]) due to one-iteration update ‖xk+1 − xk‖2 of variable x (last term in the bound). Such a value reduction is also affected by the tracking error ‖yk − y∗(xk)‖2 of the variable y (i.e., yk does not equal the desirable y∗(xk)), the variance ‖̄k‖2 of momentum recursive hypergradient estimator, and the Hessian inverse approximation error CQ w.r.t. hypergradient.
Based on Propositions 1 and 2, we next characterize the convergence of MRBO.
Theorem 1. Apply MRBO to solve the problem eq. (1). Suppose Assumptions 1, 2, and 3 hold. Let hyperparameters c1 ≥ 23d3 + 9λµ 4 , c2 ≥ 2 3d3 + 75L′2λ 2µ ,m ≥ max{2, d 3, (c1d) 3, (c2d)
3}, y1 = y∗(x1), η < 1 L , 0 ≤ λ ≤ 1 6L , 0 ≤ γ ≤ min{ 1 4LΦηK , λµ√ 150L′2L2/µ2+8λµ(L2Q+L 2) }. Then, we have
1
K K∑ k=1
( L′ 2
4 ‖y∗(xk)− yk‖2 +
1 4 ‖̄k‖2 + 1 4γ2η2k ‖xk+1 − xk‖2
) ≤ M ′
K (m+K)1/3, (8)
where L′2 is defined in Proposition 2, and M ′ = Φ(x1)−Φ ∗ γd + ( 2G2(c11+c 2 2)d 2 λµ + 2C2Qd 2
η2K
) log(m +
K) + 2G 2
Sλµdη0 .
Theorem 1 captures the simultaneous convergence of the variables xk, yk and ‖̄k‖: the tracking error ‖y∗(xk)− yk‖ converges to zero, and the variance ‖̄k‖ of the momentum recursive hypergradient estimator reduces to zero, both of which further facilitate the convergence of xk and the algorithm.
By properly choosing the hyperparameters in Algorithm 1 to satisfy the conditions in Theorem 1, we obtain the following computational complexity for MRBO. Corollary 1. Under the same conditions of Theorem 1 and choosingK = O( −1.5), Q = O(log( 1 )), MRBO in Algorithm 1 finds an -stationary point with the gradient complexity of O( −1.5) and the (Jacobian-) Hessian-vector complexity of Õ( −1.5).
As shown in Corollary 1, MRBO achieves the computational complexity of Õ( −1.5), which outperforms all existing stochastic bilevel algorithms by a factor of Õ( −0.5) (see Table 1). Further, this also achieves the best known complexity of Õ( −1.5) for vanilla nonconvex optimization via first-order stochastic algorithms. As far as we know, this is the first result to demonstrate the improved performance of single-loop recursive momentum over SGD-type updates for bilevel optimization.

heading: 3.3 Convergence Analysis of VRBO Algorithm
text: To analyze the convergence of VRBO, we need to first characterize the statistical properties of the hypergradient estimator, in which all the gradient, Jacobian-vector, and Hessian-vector have recursive variance reduction forms. We then need to characterize how the inner-loop tracking error affects the outer-loop hypergradient estimation error in order to establish the overall convergence. The complication in the analysis is mainly due to the hypergradient in bilevel optimization, which does not exist in the previous studies of variance reduction in conventional minimization and minimax optimization. Below, we provide two properties of VRBO for handling the aforementioned challenges.
In the following proposition, we characterize the variance of the hypergradient estimator, and further use such a bound to characterize the cumulative variances of both the hypergradient and inner-loop gradient estimators based on the recursive variance reduction technique over all iterations. Proposition 3. Suppose Assumptions 1, 2, 3 hold. Let η < 1L . Then the hypergradient estimator ∇̂Φ(xk, yk;S1) defined in eq. (5) w.r.t. x has bounded variance as
E‖∇̂Φ(xk, yk;S1)−∇Φ(xk)‖2 ≤ σ′2
S1 , (9)
where σ′2 = 2M2 + 28L2M2η2(Q+ 1)2. Let ∆k = E(‖vk −∇Φ(xk)‖2 + ‖uk −∇yg(xk, yk)‖2), where vk and uk denote the recursive variance reduction estimators for hypergradient and inner-loop gradient respectively. Then, the cumulative variance of vk and uk is bounded by
K−1∑ k=0 ∆k ≤ 4σ′2K S1 + 22α2L2Q K−2∑ k=0 E‖vk‖2 + 4 3 E‖∇yg(x0, y0)‖2. (10)
As shown in eq. (9), the variance bound of the hypergradient estimator increases with the number Q of Hessian-vector products for approximating the Hessian inverse and can be reduced by the batch size S1. Then eq. (10) further provides an upper bound on the cumulative variance ∑K−1 k=0 ∆k of the recursive hypergradient estimator and inner-loop gradient estimator.
Proposition 4. Suppose Assumptions 1, 2, 3 hold. Let η < 1L . Then, we have
E[Φ(xk+1)] ≤ E[Φ(xk)] + αL ′2 µ2 E‖∇yg(xk, yk)‖ 2 + αE‖∇̃Φ(xk)− vk‖2 − (α2 − α2 2 LΦ)E‖vk‖ 2,
where L′2 = (L+ L 2
µ + Mτ µ + LMρ µ2 ) 2 and ∇̃Φ(xk) takes a form of
∇̃Φ(xk) = ∇xf(xk, yk)−∇x∇yg(xk, yk)[∇2yg(xk, yk)]−1∇yf(xk, yk). (11)
Proposition 4 characterizes how the objective function value decreases (i.e., captured by E[Φ(xk+1)]− E[Φ(xk)]) due to one iteration update ‖vk‖2 of variable x (last term in the bound). Such a value reduction is also affected by the moments of gradient w.r.t. y and the variance of recursive hypergradient estimator.
Based on Propositions 3 and 4, we next characterize the convergence of VRBO. Theorem 2. Apply VRBO to solve the problem eq. (1). Suppose Assumptions 1, 2, 3 hold. Let α =
1 20L3m , β = 213LQ , η < 1 L , S2 ≥ 2( L µ + 1)Lβ,m = 16 µβ − 1, q = µLβS2 µ+L where Lm = max{LQ, LΦ}. Then, we have
1
K K−1∑ k=0 E‖∇Φ(xk)‖2 ≤ O( Q4 K + Q6 S1 +Q4(1− ηµ)2Q). (12)
Theorem 2 shows that VRBO converges sublinearly w.r.t. the number K of iterations with the convergence error consisting of two terms. The first error term Q 6
S1 is caused by the minibatch gradient
and hypergradient estimation at outer loops and can be reduced by increasing the batch size S1 (in fact, Q scales only logarithmically with S1). The second error term Q4(1 − ηµ)2Q is due to the approximation error of the Hessian-vector type of hypergradient estimation, which decreases exponentially fast w.r.t. Q. By properly choosing the hyperparameters in Algorithm 2, we obtain the following complexity result for VRBO. Corollary 2. Under the same conditions of Theorem 2, choose S1 = O( −1), S2 = O( −0.5), Q = O(log( 1 0.5 )),K = O(
−1). Then, VRBO finds an -stationary point with the gradient complexity of Õ( −1.5) and Hessian-vector complexity of Õ( −1.5).
Similarly to MRBO, Corollary 2 indicates that VRBO also outperforms all existing stochastic algorithms for bilevel optimization by a factor of Õ( −0.5) (see Table 1). Further, although MRBO and VRBO achieve the same theoretical computational complexity, VRBO empirically performs much better than MRBO (as well as other single-loop momentum-based algorithms MSTSA [22], STABLE [2], and SEMA [11]), as will be shown in Section 4.
We note that although our theory requires Q to scale as O(log( 1 0.5 )), a very small Q is sufficient to attain a fast convergence speed in experiments. For example, we choose Q = 3 in our hyper-cleaning experiments as other benchmark algorithms such as AID-FP, reverse, and stocBiO.

heading: 4 Experiments
text: In this section, we compare the performances of our proposed VRBO and MRBO algorithms with the following bilevel optimization algorithms: AID-FP [10], reverse [6] (both are double-loop deterministic algorithms), BSA [8] (double-loop stochastic algorithm), MSTSA [22] and SUSTAIN [23] (single-loop stochastic algorithms), STABLE [2] (single-loop stochastic algorithm with Hessian inverse computations), and stocBiO [20] (double-loop stochastic algorithm). SEMA [11] is not included in the list because it performs similarly to SUSTAIN. RSVRB [12] is not included since it performs similarly to STABLE. Our experiments are run over a hyper-cleaning application on MNIST3. We provide the detailed experiment specifications in Appendix B.
As shown in Figure 1 (a) and (b), the convergence rate (w.r.t. running time) of our VRBO and the SGD-type stocBiO converge much faster than other algorithms in comparison. Between VRBO and stocBiO, they have comparable performance, but our VRBO achieves a lower training loss as well as
3The experiments on CIFAR10 are still ongoing.
a more stable convergence. Further, our VRBO converges significantly faster than all single-loop momentum-based methods. This provides some evidence on the advantage of double-loop algorithms over single-loop algorithms for bilevel optimization. Moreover, our MRBO achieves the fastest convergence rate among all single-loop momentum-based algorithms, which is in consistent with our theoretical results. In Figure 1 (c), we compare our algorithms MRBO and VRBO with three momentum-based algorithms, i.e., MSTSA, STABLE, and SUSTAIN, where SUSTAIN (proposed in the concurrent work [23]) achieves the same theoretical complexity as our MRBO and VRBO. However, it can be seen that MRBO and VRBO are significantly faster than the other three algorithms.
All three plots suggest an interesting observation that double-loop algorithms tend to converge faster than single-loop algorithms as demonstrated by (i) double-loop VRBO performs the best among all algorithms; and (ii) double-loop SGD-type stocBiO, GD-type reverse and AID-FP perform even better than single-loop momentum-accelerated stochastic algorithm MRBO; and (iii) double-loop SGDtype BSA (with single-sample updates) converges faster than single-loop momentum-accelerated stochastic MSTSA, STABLE and SUSTAIN (with single-sample updates). Such a phenomenon has been observed only in bilevel optimization (to our best knowledge), and occurs oppositely in minimization and minimax problems, where single-loop algorithms substantially outperform double-loop algorithms. The reason for this can be that the hypergradient estimation at the outerloop in bilevel optimization is very sensitive to the inner-loop output. Thus, for each outer-loop iteration, sufficient inner-loop iterations in the double-loop structure provides a much more accurate output close to y∗(x) than a single inner-loop iteration, and thus helps to estimate a more accurate hypergradient in the outer loop. This further facilitates better outer-loop iterations and yields faster overall convergence.

heading: 5 Conclusion
text: In this paper, we proposed two novel algorithms MRBO and VRBO for the nonconvex-stronglyconvex bilevel stochastic optimization problem, and showed that their computational complexities outperform all existing algorithms order-wise. In particular, MRBO is the first momentum algorithm that exhibits the order-wise improvement over SGD-type algorithms for bilevel optimization, and VRBO is the first that adopts the recursive variance reduction technique to accelerate bilevel optimization. Our experiments demonstrate the superior performance of these algorithms, and further suggest that the double-loop design may be more suitable for bilevel optimization than the singleloop structure. We anticipate that our analysis can be applied to studying bilevel problems under various other loss geometries. We also hope that our study can motivate further comparison between double-loop and single-loop algorithms in bilevel optimization.

heading: Acknowledgements
text: The work was supported in part by the U.S. National Science Foundation under the grants ECCS2113860, DMS-2134145 and CNS-2112471.

Review summaries:
This paper presents two algorithms (one single-loop and one double-loop) which exploit variance reduction techniques to improve the best known sample complexity for stochastic bilevel optimization problems with strongly convex lower-level from $O(\epsilon^{-2})$ to  $O(\epsilon^{-1.5})$.    Bilevel optimization consists in two nested optimization problems: the upper-level aims at minimizing an objective function where some of the variables are a solution of the lower-level problem. In the stochastic setting, upper and lower-level objectives are not accessible and are instead replaced by unbiased estimators. In addition to several smoothness assumptions, this work and several recent others consider the case when the lower-level problem is strongly-convex, which guarantees that the lower-level solution is unique, the existence of the upper-level gradient, and allows to provide convergence guarantees to stationary points even when the lower-level solution is not computed exactly.   This work studies two algorithms (MRBO and VRBO) which both achieve a sample complexity of $O(\epsilon^{-1.5})$, which improves upon the previous optimal complexity of $O(\epsilon^{-2})$. Both algorithms exploit variance reduction  to achieve the desired complexity. MRBO is single loop and combines variance reduction and momentum while VRBO is double loop and uses large-batch gradients to reduce the variance. Recent works have already applied momentum and variance reduction to solve stochastic bilevel problems but this one is the first (together with two concurrent works) to achieve a sample complexity this low. The authors present experiments on the data hyperclearner setting which show that VRBO performs the best and double-loop algorithms seem to outperform single-loop ones. 

The authors developed two algorithms MRBO and VRBO for  bilevel stochastic optimization problem, and showed that their computational complexities ($\mathcal{O}(\epsilon^{-1.5})$) outperform all existing algorithms ($\mathcal{O}(\epsilon^{-2})$) orderwisely. The authors claim two main contributions. First, MRBO is the first momentum algorithm that exhibits the orderwise improvement over SGD-type algorithms for bilevel optimization. Second, VRBO is the first that adopts the recursive variance reduction technique to accelerate bilevel optimization. Their experiments demonstrate that these algorithms  outperform existing algorithms, and suggest that the double-loop design may be more suitable for bilevel optimization than the single loop structure.