Year: 2020

Abstract:
Bilevel optimization has arisen as a powerful tool for many machine learning problems such as meta-learning, hyperparameter optimization, and reinforcement learning. In this paper, we investigate the nonconvex-strongly-convex bilevel optimization problem. For deterministic bilevel optimization, we provide a comprehensive finite-time convergence analysis for two popular algorithms respectively based on approximate implicit differentiation (AID) and iterative differentiation (ITD). For the AID-based method, we orderwisely improve the previous finitetime convergence analysis due to a more practical parameter selection as well as a warm start strategy, and for the ITD-based method we establish the first theoretical convergence rate. Our analysis also provides a quantitative comparison between ITD and AID based approaches. For stochastic bilevel optimization, we propose a novel algorithm named stocBiO, which features a sample-efficient hypergradient estimator using efficient Jacobianand Hessian-vector product computations. We provide the finite-time convergence guarantee for stocBiO, and show that stocBiO outperforms the best known computational complexities orderwisely with respect to the condition number κ and the target accuracy . We further validate our theoretical results and demonstrate the efficiency of bilevel optimization algorithms by the experiments on meta-learning and hyperparameter optimization.

heading: 1 INTRODUCTION
text: Bilevel optimization has received significant attention recently and become an influential framework in various machine learning applications including meta-learning (Franceschi et al., 2018; Bertinetto et al., 2018; Rajeswaran et al., 2019; Ji et al., 2020a), hyperparameter optimization (Franceschi et al., 2018; Shaban et al., 2019; Feurer & Hutter, 2019), reinforcement learning (Konda & Tsitsiklis, 2000; Hong et al., 2020), and signal processing (Kunapuli et al., 2008; Flamary et al., 2014). A general bilevel optimization takes the following formulation.
min x∈Rp Φ(x) := f(x, y∗(x)) s.t. y∗(x) = arg min y∈Rq g(x, y), (1)
where the upper- and inner-level functions f and g are both jointly continuously differentiable. The goal of eq. (1) is to minimize the objective function Φ(x) w.r.t. x, where y∗(x) is obtained by solving the lower-level minimization problem. In this paper, we focus on the setting where the lower-level function g is strongly convex with respect to (w.r.t.) y, and the upper-level objective function Φ(x) is nonconvex w.r.t. x. Such types of geometrics commonly exist in many applications including meta-learning and hyperparameter optimization, where g corresponds to an empirical loss with a strongly-convex regularizer and x are parameters of neural networks.
A broad collection of algorithms have been proposed to solve such types of bilevel optimization problems. For example, Hansen et al. (1992); Shi et al. (2005); Moore (2010) reformulated the bilevel problem in eq. (1) into a single-level constrained problem based on the optimality conditions of the lower-level problem. However, such type of methods often involve a large number of constraints, and are hard to implement in machine learning applications. Recently, more efficient gradient-based bilevel optimization algorithms have been proposed, which can be generally categorized into the approximate implicit differentiation (AID) based approach (Domke, 2012; Pedregosa, 2016; Gould et al., 2016; Liao et al., 2018; Ghadimi & Wang, 2018; Grazzi et al., 2020; Lorraine
et al., 2020) and the iterative differentiation (ITD) based approach (Domke, 2012; Maclaurin et al., 2015; Franceschi et al., 2017; 2018; Shaban et al., 2019; Grazzi et al., 2020). However, most of these studies have focused on the asymptotic convergence analysis, and the finite-time analysis (that characterizes how fast an algorithm converges) has not been well explored except a few attempts recently. Ghadimi & Wang (2018) provided the finite-time analysis for the ITD-based approach. Grazzi et al. (2020) provided the iteration complexity for the hypergradient computation via ITD and AID, but did not characterize the finite-time convergence for the entire execution of algorithms.
• Thus, the first focus of this paper is to develop a comprehensive and enhanced theory, which covers a broader class of bilevel optimizers via ITD and AID based techniques, and more importantly, to improve the exiting analysis with a more practical parameter selection and order-level lower computational complexity.
The stochastic bilevel optimization often occurs in applications where fresh data need to be sampled as the algorithms run (e.g., reinforcement learning (Hong et al., 2020)) or the sample size of training data is large (e.g., hyperparameter optimization (Franceschi et al., 2018), Stackelberg game (Roth et al., 2016)). Typically, the corresponding objective function is given by
min x∈Rp
Φ(x) = f(x, y∗(x)) := { Eξ [F (x, y∗(x); ξ)] 1 n ∑n i=1 F (x, y ∗(x); ξi)
s.t. y∗(x) = arg min y∈Rq g(x, y) := { Eζ [G(x, y∗(x); ζ)] 1 m ∑m i=1G(x, y ∗(x); ζi), (2)
where f(x, y) and g(x, y) take either the expectation form w.r.t. the random variables ξ and ζ or the finite-sum form over given data Dn,m = {ξi, ζj , i = 1, ..., n; j = 1, ...,m} often with large sizes n and m. During the optimization process, the algorithms sample data batch via the distributions of ξ and ζ or from the setDn,m. For such a stochastic setting, Ghadimi & Wang (2018) proposed a bilevel stochastic approximation (BSA) method via single-sample gradient and Hessian estimates. Based on such a method, Hong et al. (2020) further proposed a two-timescale stochastic approximation (TTSA), and showed that TTSA achieves a better trade-off between the complexities of inner- and outer-loop optimization stages than BSA.
• The second focus of this paper is to design a more sample-efficient algorithm for bilevel stochastic optimization, which achieves an order-level lower computational complexity over BSA and TTSA.

heading: 1.1 MAIN CONTRIBUTIONS
text: Our main contributions lie in developing enhanced theory and provably faster algorithms for the nonconvex-strongly-convex bilevel deterministic and stochastic optimization problems, respectively. Our analysis involves several new developments, which can be of independent interest.
We first provide a unified finite-time convergence and complexity analysis for both ITD and AID based bilevel optimizers, which we call as ITD-BiO and AID-BiO. Compared to existing analysis in Ghadimi & Wang (2018) for AID-BiO that requires a continuously increasing number of innerloop steps to achieve the guarantee, our analysis allows a constant number of inner-loop steps as often used in practice. In addition, we introduce a warm start initialization for the inner-loop updates and the outer-loop hypergradient estimation, which allows us to backpropagate the tracking errors to previous loops, and results in an improved computational complexity. As shown in Table 1, the gradient complexities Gc(f, ), Gc(g, ), and Jacobian- and Hessian-vector product complexities JV(g, ) and HV(g, ) of AID-BiO to attain an -accurate stationary point improve those of Ghadimi & Wang (2018) by the order of κ, κ −1/4, κ, and κ, respectively, where κ is the condition number. In addition, our analysis shows that AID-BiO requires less computations of Jacobian- and Hessianvector products than ITD-BiO by an order of κ and κ1/2, which provides a justification for the observation in Grazzi et al. (2020) that ITD often has a larger memory cost than AID.
We then propose a stochastic bilevel optimizer (stocBiO) to solve the stochastic bilevel optimization problem in eq. (2). Our algorithm features a mini-batch hyper-gradient estimation via implicit differentiation, where the core design involves a sample-efficient Hypergradient estimator via the Neumann series. As shown in Table 2, the gradient complexities of our proposed algorithm w.r.t. F
and G improve upon those of BSA (Ghadimi & Wang, 2018) by an order of κ and −1, respectively. In addition, the Jacobian-vector product complexity JV(G, ) of our algorithm improves that of BSA by an order of κ. In terms of the target accuracy , our computational complexities improve those of TTSA (Hong et al., 2020) by an order of −1/2.
We further provide the theoretical complexity guarantee of ITD-BiO, AID-BiO and stocBiO in metalearning and hyperparameter optimization. The experiments validate our theoretical results for determinisitic bilevel optimization, and demonstrate the superior efficiency of stocBiO for stochastic bilevel optimization. Due to the space limitations, we present all theoretical and empirical results on hyperparameter optimization in the supplementary materials.

heading: 1.2 RELATED WORK
text: Bilevel optimization approaches: Bilevel optimization was first introduced by Bracken & McGill (1973). Since then, a number of bilevel optimization algorithms have been proposed, which include but not limited to constraint-based methods (Shi et al., 2005; Moore, 2010) and gradient-based methods (Domke, 2012; Pedregosa, 2016; Gould et al., 2016; Maclaurin et al., 2015; Franceschi et al., 2018; Ghadimi & Wang, 2018; Liao et al., 2018; Shaban et al., 2019; Hong et al., 2020; Liu et al., 2020; Li et al., 2020; Grazzi et al., 2020; Lorraine et al., 2020). Among them, Ghadimi & Wang (2018); Hong et al. (2020) provided the finite-time complexity analysis for their proposed methods for the nonconvex-strongly-convex bilevel optimization problem. For such a problem, this paper develops a general and enhanced finite-time analysis for gradient-based bilevel optimizers for the deterministic setting, and proposes a novel algorithm for the stochastic setting with order-level lower computational complexity than the existing results.
Some works have studied other types of loss geometries. For example, Liu et al. (2020); Li et al. (2020) assumed that the lower- and upper-level functions g(x, ·) and f(x, ·) are convex and stronglyconvex, and provided an asymptotic analysis for their methods. Ghadimi & Wang (2018); Hong et al. (2020) studied the setting where Φ(·) is strongly-convex or convex, and g(x, ·) is strongly-convex. Bilevel optimization in meta-learning: Bilevel optimization framework has been successfully employed in meta-learning recently (Snell et al., 2017; Franceschi et al., 2018; Rajeswaran et al., 2019; Zügner & Günnemann, 2019; Ji et al., 2020a;b). For example, Snell et al. (2017) proposed a bilevel optimization procedure for meta-learning to learn a common embedding model for all tasks. Rajeswaran et al. (2019) reformulated the model-agnostic meta-learning (MAML) (Finn et al., 2017) as a bilevel optimization problem, and proposed iMAML via implicit gradient. The paper provides a theoretical guarantee for two popular types of bilevel optimization algorithms, i.e., AID-BiO and ITD-BiO, for meta-learning.
Bilevel optimization in hyperparameter optimization: Hyperparameter optimization has become increasingly important as a powerful tool in the automatic machine learning (autoML) (Okuno et al.,
Algorithm 1 Deterministic bilevel optimization via AID or ITD 1: Input: Stepsizes α, β > 0, initializations x0, y0, v0. 2: for k = 0, 1, 2, ...,K do 3: Set y0k = y T k−1 if k > 0 and y0 otherwise
4: for t = 1, ...., T do 5: Update ytk = y t−1 k − α∇yg(xk, y t−1 k ) 6: end for 7: Hypergradient estimation via
• AID: 1) set v0k = vNk−1 if k > 0 and v0 otherwise 2) solve vNk from∇2yg(xk, yTk )v = ∇yf(xk, yTk ) via N steps of CG starting from v0k 3) compute Jacobian-vector product∇x∇yg(xk, yTk )vNk via automatic differentiation 4) compute ∇̂Φ(xk) = ∇xf(xk, yTk )−∇x∇yg(xk, yTk )vNk • ITD: compute ∇̂Φ(xk) = ∂f(xk,y T k )
xk via backpropagation w.r.t. xk
8: Update xk+1 = xk − β ∂f(xk,y T k ) ∂xk 9: end for
2018; Yu & Zhu, 2020). Recently, various bilevel optimization algorithms have been proposed in the context of hyperparameter optimization, which include implicit differentiation based methods (Pedregosa, 2016), dynamical system based methods via reverse or forward gradient computation (Franceschi et al., 2017; 2018; Shaban et al., 2019), etc. This paper demonstrates the superior efficiency of the proposed stocBiO algorithm in hyperparameter optimization.

heading: 2 ALGORITHMS
text: In this section, we describe two popular types of deterministic bilevel optimization algorithms, and propose a new algorithms for stochastic bilevel optimization.

heading: 2.1 ALGORITHMS FOR DETERMINISTIC BILEVEL OPTIMIZATION
text: As shown in Algorithm 1, we describe two popular types of deterministic bilevel optimizers respectively based on AID and ITD (referred to as AID-BiO and ITD-BiO) for solving the problem eq. (1).
Both AID-BiO and ITD-BiO update in a nested-loop manner. In the inner loop, both of them run T steps of gradient decent (GD) to find an approximation point yTk close to y
∗(xk). Note that we choose the initialization y0k of each inner loop as the output y T k−1 of the preceding inner loop rather than a random start. Such a warm start allows us to backpropagate the tracking error ‖yTk −y∗(xk)‖ to previous loops, and yields an improved computational complexity.
At the outer loop, AID-BiO first solves vNk from a linear system ∇2yg(xk, yTk )v = ∇yf(xk, yTk )1 usingN steps of conjugate-gradient (CG) starting from v0k (where we also adopt a warm start scheme here by setting v0k = v N k−1), and then constructs
∇̂Φ(xk) = ∇xf(xk, yTk )−∇x∇yg(xk, yTk )vNk (3) as an estimate of the true hypergradient∇Φ(xk), whose form is given by the following proposition. Proposition 1. Recalling the definition Φ(x) := f(x, y∗(x)), it holds that
∇Φ(xk) =∇xf(xk, y∗(xk))−∇x∇yg(xk, y∗(xk))v∗k, (4)
where v∗k is the solution of the linear system∇2yg(xk, y∗(xk))v = ∇yf(xk, y∗(xk)).
As shown in Domke (2012); Grazzi et al. (2020), the construction of eq. (3) involves only Hessianvector products in solving vN via CG and Jacobian-vector product ∇x∇yg(xk, yTk )vNk , which can be efficiently computed and stored via existing automatic differentiation packages.
As a comparison, the outer loop of ITD-BiO computes the gradient ∂f(xk,y T k (xk))
∂xk as an approxima-
tion of the hyper-gradient ∇Φ(xk) = ∂f(xk,y ∗(xk))
∂xk via backpropagation, where we write yTk (xk)
1This is equivalent to solve a quadratic programing minv 12v T∇2yg(xk, yTk )v − vT∇yf(xk, yTk ).
Algorithm 2 Stochastic bilevel optimizer (stocBiO) 1: Input: Inner- and outer-loop stepsizes α, β > 0, initializations x0 and y0. 2: for k = 0, 1, 2, ...,K do 3: Set y0k = y T k−1 if k > 0 and y0 otherwise
4: for t = 1, ...., T do 5: Draw a sample batch St−1 6: Update ytk = y t−1 k − α∇yG(xk, y t−1 k ;St−1) 7: end for 8: Draw sample batch DF , and compute v0 = ∇yF (xk, yTk ;DF ) 9: Draw sample batch DH , and construct vQ via Algorithm 3
10: Draw sample batch DG, and compute Jacobian-vector product∇x∇yG(xk, yTk ;DG)vQ 11: Compute gradient estimate ∇̂Φ(xk) via eq. (6) 12: Update xk+1 = xk − β∇̂Φ(xk) 13: end for
Algorithm 3 Construct vQ given v0 1: Input: An integer Q, data samples DH = {Bj}Qj=1 and a constant η > 0. 2: for j = 1, 2, ..., Q do 3: Sample Bj and compute gradient Gj(y) = y − η∇yG(x, y;Bj) 4: end for 5: Set rQ = v0 6: for i = Q, ..., 1 do 7: ri−1 = ∂ ( Gi(y)ri ) /∂y = ri − η∇2yG(x, y;Bi)ri via automatic differentiation
8: end for 9: Return vQ = η ∑Q i=0 ri
because the output yTk of the inner loop has a dependence on xk through the inner-loop iterative GD updates. The explicit form of the estimate ∂f(xk,y T k (xk))
∂xk is given by the following proposition via the chain rule. For notation simplification, let ∏T−1 j=T (·) = I .
Proposition 2. The gradient ∂f(xk,y T k (xk))
∂xk takes the following analytical form:
∂f(xk, y T k )
∂xk = ∇xf(xk, yTk )− α T−1∑ t=0 ∇x∇yg(xk, ytk) T−1∏ j=t+1 (I − α∇2yg(xk, y j k))∇yf(xk, y T k ).
Proposition 2 shows that the differentiation involves the computations of second-order derivatives such as Hessian ∇2yg(·, ·). Since efficient Hessian-free methods such as CG have been successfully deployed in the existing automatic differentiation tools, computing these second-order derivatives reduces to more efficient computations of Jacobian- and Hessian-vector products.

heading: 2.2 ALGORITHM FOR STOCHASTIC BILEVEL OPTIMIZATION
text: We propose a new stochastic bilevel optimizer (stocBiO) in Algorithm 2 to solve the problem eq. (2). It has a double-loop structure similar to Algorithm 1, but runs T steps of stochastic gradient decent (SGD) at the inner loop to obtain an approximated solution yTk . Based on the output y T k of the inner loop, stocBiO first computes a gradient ∇yF (xk, yTk ;DF ) over a sample batch DF , and then computes a vector vQ via Algorithm 3, which takes a form of
vQ = η Q−1∑ q=−1 Q∏ j=Q−q (I − η∇2yG(xk, yTk ;Bj))∇yF (xk, yTk ;DF ), (5)
where {Bj , j = 1, ..., Q} are mutually-independent sample sets, Q and η are constants, and we let∏Q Q+1(·) = I for notational simplification. Note that our construction of vQ, i.e., Algorithm 3, is
motived by the Neumann series ∑∞ i=0 U
k = (I −U)−1, and involves only Hessian-vector products rather than Hessians, and hence is computationally and memory efficient. Then, we construct
∇̂Φ(xk) =∇xF (xk, yTk ;DF )−∇x∇yG(xk, yTk ;DG)vQ (6)
as an estimate of hypergradient ∇Φ(xk) given by Proposition 1. An important component of our algorithm is vQ, which serves as an estimate of v∗k in eq. (4) . Compared to the deterministic case, designing a sample-efficient Hypergradient estimator in the stochastic case is more challenging. For example, instead of choosing the same batch sizes for all Bj , j = 1, ..., Q in eq. (5), our analysis captures the different impact of components ∇2yG(xk, yTk ;Bj), j = 1, ..., Q on the Hypergradient estimation variance, and inspires an adaptive and more efficient choice by setting |BQ−j | to decay exponentially with j from 0 to Q− 1. By doing so, we achieve an improved complexity.

heading: 3 DEFINITIONS AND ASSUMPTIONS
text: Let z = (x, y) denote all parameters. For simplicity, suppose sample sets St for all t = 0, ..., T − 1, DG and DF have the sizes of S, Dg and Df , respectively. In this paper, we focus on the following types of loss functions for both the deterministic and stochastic cases. Assumption 1. The lower-level function g(x, y) is µ-strongly-convex w.r.t. y and the total objective function Φ(x) = f(x, y∗(x)) is nonconvex w.r.t. x. For the stochastic setting, the same assumptions hold for G(x, y; ζ) and Φ(x), respectively. Since the objective function Φ(x) is nonconvex, algorithms are expected to find an -accurate stationary point defined as follows. Definition 1. We say x̄ is an -accurate stationary point for the objective function Φ(x) in eq. (2) if E‖∇Φ(x̄)‖2 ≤ , where x̄ is the output of an algorithm.
In order to compare the performance of different bilevel algorithms, we adopt the following metrics of computational complexity. Definition 2. For a function f(x, y) and a vector v, let Gc(f, ) be the number of the partial gradient ∇xf or ∇yf , and let JV(g, ) and HV(g, ) be the number of Jacobian-vector products ∇x∇yg(x, y)v. and Hessian-vector products ∇2yg(x, y)v. For the stochastic case, similar metrics are adopted but w.r.t. the stochastic function F (x, y; ξ).
We take the following standard assumptions on the loss functions in eq. (2), which have been widely adopted in bilevel optimization (Ghadimi & Wang, 2018; Ji et al., 2020a). Assumption 2. The loss function f(z) and g(z) satisfy
• f(z) is M -Lipschitz, i.e., for any z, z′, |f(z)− f(z′)| ≤M‖z − z′‖.
• Gradients∇f(z) and ∇f(z) are L-Lipschitz, i.e., for any z, z′, ‖∇f(z)−∇f(z′)‖ ≤ L‖z − z′‖, ‖∇g(z)−∇g(z′)‖ ≤ L‖z − z′‖.
For the stochastic case, the same assumptions hold for F (z; ξ) and G(z; ζ) for any given ξ and ζ.
As shown in Proposition 1, the gradient of the objective function Φ(x) involves the second-order derivatives ∇x∇yg(z) and ∇2yg(z). The following assumption imposes the Lipschitz conditions on such high-order derivatives, as also made in Ghadimi & Wang (2018). Assumption 3. Suppose the derivatives∇x∇yg(z) and ∇2yg(z) are τ - and ρ- Lipschitz, i.e.,
• For any z, z′, ‖∇x∇yg(z)−∇x∇yg(z′)‖ ≤ τ‖z − z′‖.
• For any z, z′, ‖∇2yg(z)−∇2yg(z′)‖ ≤ ρ‖z − z′‖.
For the stochastic case, the same assumptions hold for∇x∇yG(z; ζ) and∇2yG(z; ζ) for any ζ.
As typically adopted in the analysis for stochastic optimization, we make the following boundedvariance assumption for the lower-level stochastic function G(z; ζ). Assumption 4. ∇G(z; ζ) has a bounded variance, i.e., Eξ‖∇G(z; ζ)−∇g(z)‖2 ≤ σ2 for some σ.

heading: 4 MAIN RESULTS FOR BILEVEL OPTIMIZATION

heading: 4.1 DETERMINISTIC BILEVEL OPTIMIZATION
text: We first characterize the convergence and complexity performance of the AID-BiO algorithm. Let κ = Lµ denote the condition number.
Theorem 1 (AID-BiO). Suppose Assumptions 1, 2, 3 hold. Define a smoothness parameter LΦ = L+ 2L 2+τM2
µ + ρLM+L3+τML µ2 + ρL2M µ3 = Θ(κ 3), choose the stepsizes α ≤ 1L , β = 1 8LΦ , and set
the inner-loop iteration number T ≥ Θ(κ) and the CG iteration number N ≥ Θ( √ κ), where the detailed forms of T,N can be found in Appendix E. Then, the outputs of AID-BiO satisfy
1
K K−1∑ k=0 ‖∇Φ(xk)‖2 ≤ 64LΦ(Φ(x0)− infx Φ(x)) + 5∆0 K , (7)
where ∆0 = ‖y0 − y∗(x0)‖2 + ‖v∗0 − v0‖2 > 0. In order to achieve an -accurate stationary point, we have
• Gradient complexity: Gc(f, ) = O(κ3 −1),Gc(g, ) = O(κ4 −1). • Jacobian- and Hessian-vector product: JV(g, ) = O ( κ3 −1 ) ,HV(g, ) = O ( κ3.5 −1 ) .
It can be seen from Table 1 that the complexities Gc(f, ),Gc(g, ), JV(g, ) and HV(g, ) of our analysis improves that of Ghadimi & Wang (2018) (eq. (2.30) therein) by the order of κ, κ −1/4, κ and κ. Such an improvement is achieved by a refined analysis with a constant number of innerloop steps, and by a warm start strategy to backpropagate the tracking errors ‖yTk − y∗(xk)‖ and ‖vNk − v∗k‖ to previous loops, as also demonstrated by our meta-learning experiments. We next characterize the convergence and complexity performance of the ITD-BiO algorithm. Theorem 2 (ITD-BiO). Suppose Assumptions 1, 2, and 3 hold. Define the parameter LΦ as in Theorem 1, and choose α ≤ 1L , β = 1 4LΦ
and T ≥ Θ(κ log 1 ), where the detailed form of T can be found in Appendix F. Then, the outputs of ITD-BiO satisfy
1
K K−1∑ k=0 ‖∇Φ(xk)‖2 ≤ 16LΦ(Φ(x0)− infx Φ(x)) K + 2 3 .
In order to achieve an -accurate stationary point, we have
• Gradient complexity: Gc(f, ) = O(κ3 −1),Gc(g, ) = O(κ4 −1 log( 1 )).
• Jacobian- and Hessian-vector product complexity: JV(g, ) = O ( κ4 −1 log −1 ) ,HV(g, ) = O ( κ4 −1 log −1 ) .
By comparing Theorem 1 and Theorem 2, it can be seen that the complexities Gc(g, ), JV(g, ), and HV(g, ) of AID-BiO are better than those of ITD-BiO by the order of log( 1 ), κ log( 1 ) and κ0.5 log( 1 ). This is in consistence with the comparison in Grazzi et al. (2020) that AID-BiO often has a lower memory cost than ITD-BiO.

heading: 4.2 STOCHASTIC BILEVEL OPTIMIZATION
text: We first characterize the bias and variance of an important component vQ in eq. (5).
Proposition 3. Suppose Assumptions 1, 2 and 3 hold. Let the constant η ≤ 1L and choose the batch sizes |BQ+1−j | = BQ(1−ηµ)j−1 for j = 1, ..., Q, whereB ≥ 1Q(1−ηµ)Q−1 . Then, the bias satisfies∥∥EvQ − [∇2yg(xk, yTk )]−1∇yf(xk, yTk )∥∥ ≤ µ−1(1− ηµ)Q+1M. (8) Furthermore, the estimation variance is given by
E‖vQ − [∇2yg(xk, yTk )]−1∇yf(xk, yTk )‖2 ≤ 4η2L2M2 µ2 1 B + 4(1− ηµ)2Q+2M2 µ2 + 2M2 µ2Df . (9)
Proposition 3 shows that if we choose Q and B at the order level of O(log 1 ) and O(1/ ), the bias and variance are smaller than O( ), and the required number of samples is ∑Q j=1BQ(1 −
ηµ)j−1 = O ( −1 log 1 ) . Note that the chosen batch size |BQ+1−j | exponentially decays w.r.t. j.
In comparison, the uniform choice of all |Bj | would yield a worse complexity of O ( −1(log 1 ) 2 ) .
We next analyze stocBiO when the objective function Φ(x) := f(x, y∗(x)) is nonconvex.
Theorem 3. Suppose Assumptions 1, 2, 3 and 4 hold. Define parameter LΦ = L + 2L 2+τM2
µ + ρLM+L3+τML
µ2 + ρL2M µ3 = O ( κ3 ) , choose stepsize β = 14LΦ , and set η < 1 L in Algorithm 3. Set
T ≥ max { log ( 12+ 48β 2L2 µ2 (L+L 2 µ + Mτ µ + LMρ µ2 )2 )
2 log(L+µL−µ ) ,
log (√ β(L+L 2
µ + Mτ µ +
LMρ µ2 ) )
log(L+µL−µ )
} . Then, we have
1
K K−1∑ k=0 E‖∇Φ(xk)‖2 ≤ 32LΦ(Φ(x0)− infx Φ(x) + 52‖y0 − y ∗(x0)‖2) K + 72κ2M2(1− ηµ)2Q
+ 40 ( L+ L 2 µ + Mτ µ + LMρ µ2 )2 Lµ σ2 S + 16κ2M2 Dg + (8 + 32κ2)M2 Df + 64κ2M2 B . (10)
In order to achieve an -accurate stationary point, we have
• Gradient complexity: Gc(F, ) = O(κ5 −2),Gc(G, ) = O(κ9 −2).
• Jacobian- and Hessian-vector product: JV(G, ) = O(κ5 −2),HV(G, ) = Õ(κ6 −2).
Theorem 3 shows that stocBiO converges sublinearly with the convergence error decaying exponentially w.r.t. Q and sublinearly w.r.t. the batch sizes S,Dg, Df for gradient estimation and B for Hessian inverse estimation. In addition, it can be seen that the total number T of the inner-loop steps is chosen at nearly a constant level, rather than a typical choice of Θ(log( 1 )).
As shown in Table 2, the gradient complexities of our proposed algorithm in terms of F and G improve those of BSA in Ghadimi & Wang (2018) by an order of κ and −1, respectively. In addition, the Jacobian-vector product complexity JV(G, ) of our algorithm improves that of BSA by the order of κ. In terms of the accuracy , our gradient, Jacobian- and Hessian-vector product complexities improve those of TTSA in Hong et al. (2020) all by an order of −0.5.

heading: 5 APPLICATIONS TO META-LEARNING

heading: 5.1 META-LEARNING WITH COMMON EMBEDDING MODEL
text: Consider the few-shot meta-learning problem with m tasks {Ti, i = 1, ...,m} sampled from distribution PT . Each task Ti has a loss function L(φ,wi; ξ) over each data sample ξ, where φ are the parameters of an embedding model shared by all tasks, and wi are the task-specific parameters. The goal of this framework is to find good parameters φ for all tasks, and building on the embedded features, each task then adapts its own parameters wi by minimizing its loss.
The model training takes a bilevel procedure. In the lower-level stage, building on the embedded features, the base learner of task Ti searches w∗i as the minimizer of its loss function over a training set Si. In the upper-level stage, the meta-learner evaluates the minimizers w∗i , i = 1, ...,m on held-out test sets, and optimizes φ of the embedding model over all tasks. Specifically, let w̃ = (w1, ..., wm) denote all task-specific parameters. Then, the objective function is given by
min φ LD(φ, w̃∗) :=
1
m m∑ i=1
1 |Di| ∑ ξ∈Di
L(φ,w∗i ; ξ)︸ ︷︷ ︸ LDi (φ,w ∗ i ): task-specific upper-level loss
s.t. w̃∗ = arg min w̃ LS(φ, w̃) = arg min (w1,...,wm)
1
m m∑ i=1 ( 1 |Si| ∑ ξ∈Si
L(φ,wi; ξ) +R(wi)︸ ︷︷ ︸ LSi (φ,wi): task-specific lower-level loss
) , (11)
where Si and Di are the training and test datasets of task Ti, and R(wi) is a strongly-convex regularizer, e.g., L2. Note that the lower-level problem is equivalent to solving each w∗i as a minimizer of the task-specific loss LSi(φ,wi) for i = 1, ...,m. In practice, wi often corresponds to the parameters of the last linear layer of a neural network and φ are the parameters of the remaining layers (e.g., 4 convolutional layers in Bertinetto et al. (2018); Ji et al. (2020a)), and hence the lower-level function is strongly-convex w.r.t. w̃ and the upper-level function LD(φ, w̃∗(φ)) is generally nonconvex w.r.t. φ. In addition, due to the small sizes of datasets Di and Si in few-shot learning, all
updates for each task Ti use full gradient descent without data resampling. As a result, AID-BiO and ITD-BiO in Algorithm 1 can be applied here. In some applications where the number m of tasks is large, it is more efficient to sample a batch B of i.i.d. tasks from {Ti, i = 1, ...,m} at each meta (outer) iteration, and optimizes the mini-batch versions LD(φ, w̃;B) = 1|B| ∑ i∈B LDi(φ,wi)
and LS(φ, w̃;B) = 1|B| ∑ i∈B LSi(φ,wi) instead. The following theorem provides the convergence analysis of ITD-BiO for this case. Theorem 4. Suppose Assumptions 1, 2 and 3 hold and suppose each task loss LSi(φ,wi) is µstrongly-convex w.r.t. wi. Choose the same parameters β, T as in Theorem 2. Then, we have
1
K K−1∑ k=0 E‖∇Φ(φk)‖2 ≤ 16LΦ(Φ(φ0)− infφ Φ(φ)) K + 2 3 + ( 1 + L µ )2 M2 8|B| .
Theorem 4 shows that compared to the full batch (i.e., without task sampling) case in eq. (11), the task sampling introduces a variance termO( 1|B| ) due to the stochastic nature of the algorithm. Using an approach similar to Theorem 4, we can derive a similar result for AID-BiO.

heading: 5.2 EXPERIMENTS
text: To validate our theoretical results for deterministic bilevel optimization, we compare the performance among the following four algorithms: ITD-BiO, AID-BiO-constant (AID-BiO with a constant number of inner-loop steps as in our analysis), AID-BiO-increasing (AID-BiO with an increasing number of inner-loop steps under analysis in Ghadimi & Wang (2018)), and two popular meta-learning algorithms MAML2 (Finn et al., 2017) and ANIL3 (Raghu et al., 2019). We conduct experiments over a 5-way 5-shot task on two benchmark datasets: FC100 and miniImageNet, and the results are averaged over 10 trials with different random seeds. Due to the space limitations, we provide the model architectures, hyperparameter settings and additional experiments in Appendix B.
It can be seen from Figure 1 that for both the miniImageNet and FC100 datasets, AID-BiO-constant converges faster than AID-BiO-increasing in terms of both the training accuracy and test accuracy, and achieves a better final test accuracy than ANIL and MAML. This demonstrates the superior improvement of our developed analysis over existing analysis in Ghadimi & Wang (2018) for AIDBiO algorithm. Moreover, it can be observed that AID-BiO is slightly faster than ITD-BiO in terms of the training accuracy and test accuracy. This is also in consistence with our theoretical results.

heading: 6 CONCLUSION
text: In this paper, we develop a general and enhanced finite-time analysis for the nonconvex-stronglyconvex bilevel deterministic optimization, and propose a novel algorithm for the stochastic setting whose computational complexity outperforms the best known results order-wisely. We also provide the theoretical guarantee of various bilevel optimizers in meta-learning and hyperparameter optimization. The experiments validate our theoretical results and demonstrate the effectiveness of the proposed algorithm. We anticipate that the finite-time analysis that we develop will be useful for analyzing other bilevel optimization problems with different loss geometries, and the proposed algorithms will be useful for other applications such as reinforcement learning and Stackelberg game.
2MAML consists of an inner loop for task adaptation and an outer loop for meta initialization training. 3ANIL refers to almost no inner loop, which is an efficient MAML variant with task-specific adaption on
the last-layer of parameters.

heading: Supplementary Materials

heading: A APPLICATION TO HYPERPARAMETER OPTIMIZATION

heading: A.1 HYPERPARAMETER OPTIMIZATION
text: The goal of hyperparameter optimization (Franceschi et al., 2018; Feurer & Hutter, 2019) is to search for representation or regularization parameters λ to minimize the validation error evaluated over the learner’s parameters w∗, where w∗ is the minimizer of the inner-loop regularized training error. Mathematically, the objective function is given by
min λ LDval(λ) =
1 |Dval| ∑ ξ∈Dval L(w∗(λ); ξ)
s.t. w∗(λ) = arg min w
LDtr(w, λ) := 1 |Dtr| ∑ ξ∈Dtr ( L(w, λ; ξ) +R(w, λ) ) , (12)
where Dval and Dtr are validation and training data, L is the loss, andR(w, λ) is a regularizer. In practice, the lower-level function LDtr(w, λ) is often strongly-convex w.r.t. w. For example, for the data hyper-cleaning application proposed by Franceschi et al. (2018); Shaban et al. (2019), the predictor is modeled by a linear classifier, and the loss function L(w; ξ) is convex w.r.t. w and R(w, λ) is a strongly-convex regularizer, e.g., L2 regularization. In addition, the sample sizes of Dval and Dtr are often large, and stochastic algorithms are preferred for achieving better efficiency. As a result, the above hyperparameter optimization falls into the stochastic bilevel optimization we study in eq. (2), and we can apply the proposed stocBiO algorithm here and Theorem 3 establishes its finite-time performance guarantee.

heading: A.2 EXPERIMENTS
text: We compare our proposed stocBiO with the following baseline bilevel optimization algorithms.
• BSA (Ghadimi & Wang, 2018): implicit gradient based stochastic bilevel optimizer via singlesample data sampling.
• TTSA (Hong et al., 2020): two-time-scale stochastic optimizer via single-sample data sampling. • HOAG (Pedregosa, 2016): a hyperparameter optimization algorithm with approximate gradient.
We use the implementation in the repository https://github.com/fabianp/hoag. • reverse (Franceschi et al., 2017): an iterative differentiation based method that approximates
the hypergradient via backpropagation. We use its implementation in https://github.com/ prolearner/hypertorch.
• AID-FP (Grazzi et al., 2020): AID with the fixed-point method. We use its implementation in https://github.com/prolearner/hypertorch
• AID-CG (Grazzi et al., 2020): AID with the conjugate gradient method. We use its implementation in https://github.com/prolearner/hypertorch.
We demonstrate the effectiveness of the proposed stocBiO algorithm on two experiments: data hyper-cleaning and logistic regression.
Logistic Regression on 20 Newsgroup: We compare the performance of our algorithm stocBiO with the existing baseline algorithms reverse, AID-FP, AID-CG and HOAG over a logistic regression problem on 20 Newsgroup dataset Grazzi et al. (2020). The objective function of such a problem is given by
min λ E(λ,w∗) =
1 |Dval| ∑
(xi,yi)∈Dval
L(xiw ∗, yi)
s.t. w∗ = arg min w∈Rp×c ( 1 |Dtr| ∑ (xi,yi)∈Dtr L(xiw, yi) + 1 cp c∑ i=1 p∑ j=1 exp(λj)w 2 ij ) ,
where L is the cross-entropy loss, c = 20 is the number of topics, and p = 101631 is the feature dimension. Following Grazzi et al. (2020), we use SGD as the optimizer for the outer-loop update for all algorithms. For reverse, AID-FP, AID-CG, we use the suggested and well-tuned hyperparameter setting in their implementations https://github.com/prolearner/hypertorch on this application. In specific, they choose the inner- and outer-loop stepsizes as 100, the number of inner loops as 10, the number of CG steps as 10. For HOAG, we use the same parameters as reverse, AIDFP, AID-CG. For stocBiO, we use the same parameters as reverse, AID-FP, AID-CG, and choose η = 0.5, Q = 10. We use stocBiO-B as a shorthand of stocBiO with a batch size of B.
As shown in Figure 2, the proposed stocBiO achieves the fastest convergence rate as well as the best test accuracy among all comparison algorithms. This demonstrates the practical advantage of our proposed algorithm stocBiO. Note that we do not include BSA and TTSA in the comparison, because they converge too slowly with a large variance, and are much worse than the other competing algorithms. In addition, we investigate the impact of the batch size on the performance of our stocBiO in Figure 3. It can be seen that stocBiO outperforms HOAG under the batch sizes of 100, 500, 1000, 2000. This shows that the performance of stocBiO is not very sensitive to the batch size, and hence the tuning of the batch size is easy to handle in practice.
Data Hyper-Cleaning on MNIST. We first compare the performance of our proposed algorithm stocBiO with other baseline algorithms BSA, TTSA, HOAG4 on a hyperparameter optimization problem: data hyper-cleaning (Shaban et al., 2019) on a dataset derived from MNIST (LeCun et al., 1998), which consists of 20000 images for training, 5000 images for validation, and 10000 images for testing. Data hyper-cleaning is to train a classifier in a corrupted setting where each label of training data is replaced by a random class number with a probability p (i.e., the corruption rate).
4We do not include reverse, AID-CG and AID-FG because they perform similarly to HOAG.
The objective function is given by
min λ E(λ,w∗) =
1 |Dval| ∑
(xi,yi)∈Dval
L(w∗xi, yi)
s.t. w∗ = arg min w L(w, λ) := 1 |Dtr| ∑ (xi,yi)∈Dtr σ(λi)L(wxi, yi) + Cr‖w‖2,
where L is the cross-entropy loss, σ(·) is the sigmoid function, Cr is a regularization parameter. Following Shaban et al. (2019), we choose Cr = 0.001. All results are averaged over 10 trials with different random seeds. We adopt Adam (Kingma & Ba, 2014) as the optimizer for the outer-loop update for all algorithms. For stochastic algorithms, we set the batch size as 50 for stocBiO, and 1 for BSA and TTSA because they use the single-sample data sampling. For all algorithms, we use a grid search to choose the inner-loop stepsize from {0.01, 0.1, 1, 10}, the outer-loop stepsize from {10i, i = −4,−3,−2,−1, 0, 1, 2, 3, 4}, and the number T of inner-loop steps from {1, 10, 50, 100, 200, 1000}, where values that achieve the lowest loss after a fixed running time are selected. For stocBiO, BSA, and TTSA, we choose η from {0.5 × 2i, i = −3,−2,−1, 0, 1, 2, 3}, and Q from {3× 2i, i = 0, 1, 2, 3}.
It can be seen from Figure 4 that our proposed stocBiO algorithm achieves the fastest convergence rate among all competing algorithms in terms of both the training loss and the test loss. In addition, it is observed that such an improvement is more significant when the corruption rate p is smaller. We note that the stochastic algorithm TTSA converges very slowly with a large variance. This is because TTSA updates the costly outer loop more frequently than other algorithms, and has a larger variance due to the single-sample data sampling. As a comparison, our stocBiO achieves a much lower variance for hypergradient estimation as well as a much faster convergence rate. This verifies our theoretical results in Theorem 3.

heading: B FURTHER SPECIFICATIONS ON META-LEARNING EXPERIMENTS

heading: B.1 DATASETS AND MODEL ARCHITECTURES
text: FC100 (Oreshkin et al., 2018) is a dataset derived from CIFAR-100 (Krizhevsky & Hinton, 2009), and contains 100 classes with each class consisting of 600 images of size 32. Following Oreshkin et al. (2018), these 100 classes are split into 60 classes for meta-training, 20 classes for metavalidation, and 20 classes for meta-testing. For all comparison algorithms, we use a 4-layer convolutional neural networks (CNN) with four convolutional blocks, in which each convolutional block
contains a 3× 3 convolution (padding = 1, stride = 2), batch normalization, ReLU activation, and 2× 2 max pooling. Each convolutional layer has 64 filters. The miniImageNet dataset (Vinyals et al., 2016) is generated from ImageNet Russakovsky et al. (2015), and consists of 100 classes with each class containing 600 images of size 84×84. Following the repository Arnold et al. (2019), we partition these classes into 64 classes for meta-training, 16 classes for meta-validation, and 20 classes for meta-testing. Following the repository (Arnold et al., 2019), we use a four-layer CNN with four convolutional blocks, where each block sequentially consists of a 3× 3 convolution, batch normalization, ReLU activation, and 2× 2 max pooling. Each convolutional layer has 32 filters.
B.2 IMPLEMENTATIONS AND HYPERPARAMETER SETTINGS
We adopt the existing implementations in the repository (Arnold et al., 2019) for ANIL and MAML. For all algorithms, we adopt Adam (Kingma & Ba, 2014) as the optimizer for the outer-loop update.
Parameter selection for the experiments in Figure 1(a): For ANIL and MAML, we adopt the suggested hyperparameter selection in the repository (Arnold et al., 2019). In specific, for ANIL, we choose the inner-loop stepsize as 0.1, the outer-loop (meta) stepsize as 0.002, the task sampling size as 32, and the number of inner-loop steps as 5L. For MAML, we choose the inner-loop stepsize as 0.5, the outer-loop stepsize as 0.003, the task sampling sizeas 32, and the number of inner-loop steps as 3. For ITD-BiO, AID-BiO-constant and AID-BiO-increasing, we use a grid search to choose the inner-loop stepsize from {0.01, 0.1, 1, 10}, the task sampling size from {32, 128, 256}, and the outer-loop stepsize from {10i, i = −3,−2,−1, 0, 1, 2, 3}, where values that achieve the lowest loss after a fixed running time are selected. For ITD-BiO and AID-BiO-constant, we choose the number of inner-loop steps from {5, 10, 15, 20, 50}, and for AID-BiO-increasing, we choose the number of inner-loop steps as dc(k + 1)1/4e as adopted by the analysis in Ghadimi & Wang (2018), where we choose c from {0.5, 2, 5, 10, 50}. For both AID-BiO-constant and AID-BiO-increasing, we choose the number N of CG steps for solving the linear system from {5, 10, 15}. Parameter selection for the experiments in Figure 1(b): For ANIL and MAML, we adopt the suggested hyperparameter selection in the repository (Arnold et al., 2019). Specifically, for ANIL, we choose the inner-loop stepsize as 0.1, the outer-loop (meta) stepsize as 0.001, the task sampling size as 32 and the number of inner-loop steps as 10. For MAML, we choose the inner-loop stepsize as 0.5, the outer-loop stepsize as 0.001, the task samling size as 32, and the number of inner-loop steps as 3. For ITD-BiO, AID-BiO-constant and AID-BiO-increasing, we adopt the same procedure as in the experiments in Figure 1(a).

heading: B.3 ADDITIONAL RESULTS FOR META LEARNING
text: In this subsection, we compare the robustness between bilevel optimizer ITD-BiO (AID-BiO performs similarly to ITD-BiO in terms of the convergence rate) and ANIL (ANIL outperforms MAML in general) to the number of inner-loop steps. For the experiments in Figure 5, we choose the innerloop stepsize as 0.05, the outer-loop (meta) stepsize as 0.002, the mini-batch size as 32, and the number T of inner-loop steps as 10 for both ANIL and ITD-BiO. For the experiments in Figure 6, we choose the inner-loop stepsize as 0.1, the outer-loop (meta) stepsize as 0.001, the mini-batch size as 32, and the number T of inner-loop steps as 20 for both ANIL and ITD-BiO.
It can be seen from Figure 5 and Figure 6 that when the number of inner-loop steps become larger, i.e., T = 10 for miniImageNet and T = 20 for FC100, the bilevel optimizer ITD-BiO converges stably with a small variance, whereas ANIL suffers from a sudden descent at 1500s on miniImageNet and even diverges after 2000s on FC100.

heading: C SUPPORTING LEMMAS
text: In this section, we provide some auxiliary lemmas used for proving the main convergence results.
First note that the Lipschitz properties in Assumption 2 imply the following lemma. Lemma 1. Suppose Assumption 2 holds. Then, the stochastic derivatives ∇F (z; ξ), ∇G(z; ξ), ∇x∇yG(z; ξ) and∇2yG(z; ξ) have bounded variances, i.e., for any z and ξ,
• Eξ ‖∇F (z; ξ)−∇f(z)‖2 ≤M2.
• Eξ ‖∇x∇yG(z; ξ)−∇x∇yg(z)‖2 ≤ L2.
• Eξ ∥∥∇2yG(z; ξ)−∇2yg(z)∥∥2 ≤ L2.
Recall that Φ(x) = f(x, y∗(x)) in eq. (2). Then, we use the following lemma to characterize the Lipschitz properties of∇Φ(x), which is adapted from Lemma 2.2 in Ghadimi & Wang (2018). Lemma 2. Suppose Assumptions 1, 2 and 3 hold. Then, we have, for any x, x′ ∈ Rp,
‖∇Φ(x)−∇Φ(x′)‖ ≤ LΦ‖x− x′‖, where the constant LΦ is given by
LΦ = L+ 2L2 + τM2 µ + ρLM + L3 + τML µ2 + ρL2M µ3 . (13)

heading: D PROOF OF PROPOSITIONS IN SECTION 2
text: In this section, we provide the proofs for Proposition 1 and Proposition 2 in Section 2.

heading: D.1 PROOF OF PROPOSITION 1
text: Using the chain rule over the gradient∇Φ(xk) = ∂f(xk,y ∗(xk))
∂xk , we have
∇Φ(xk) = ∇xf(xk, y∗(xk)) + ∂y∗(xk)
∂xk ∇yf(xk, y∗(xk)). (14)
Based on the optimality of y∗(xk), we have ∇yg(xk, y∗(xk)) = 0, which, using the implicit differentiation w.r.t. xk, yields
∇x∇yg(xk, y∗(xk)) + ∂y∗(xk)
∂xk ∇2yg(xk, y∗(xk)) = 0. (15)
Let v∗k be the solution of the linear system ∇2yg(xk, y∗(xk))v = ∇yf(xk, y∗(xk)). Then, multiplying v∗k at the both sides of eq. (15), yields
−∇x∇yg(xk, y∗(xk))v∗k = ∂y∗(xk)
∂xk ∇2yg(xk, y∗(xk))v∗k =
∂y∗(xk)
∂xk ∇yf(xk, y∗(xk)),
which, in conjunction with eq. (14) , yields the proof.

heading: D.2 PROOF OF PROPOSITION 2
text: Based on the iterative update of line 5 in Algorithm 1, we have yTk = y 0 k − α ∑T−1 t=0 ∇yg(xk, ytk), which, combined with the fact that ∇yg(xk, ytk) is differentiable w.r.t. xk, indicates that the inner output yTk is differentiable w.r.t. xk. Then, based on the chain rule, we have
∂f(xk, y T k )
∂xk = ∇xf(xk, yTk ) + ∂yTk ∂xk ∇yf(xk, yTk ). (16)
Based on the iterative updates that ytk = y t−1 k − α∇yg(xk, y t−1 k ) for t = 1, ..., T , we have
∂ytk ∂xk = ∂yt−1k ∂xk − α∇x∇yg(xk, yt−1k )− α ∂yt−1k ∂xk ∇2yg(xk, yt−1k )
= ∂yt−1k ∂xk (I − α∇2yg(xk, yt−1k ))− α∇x∇yg(xk, y t−1 k ).
Telescoping the above equality over t from 1 to T yields
∂yTk ∂xk = ∂y0k ∂xk T−1∏ t=0 (I − α∇2yg(xk, ytk))− α T−1∑ t=0 ∇x∇yg(xk, ytk) T−1∏ j=t+1 (I − α∇2yg(xk, y j k))
(i) = − α T−1∑ t=0 ∇x∇yg(xk, ytk) T−1∏ j=t+1 (I − α∇2yg(xk, y j k)). (17)
where (i) follows from the fact that ∂y 0 k
∂xk = 0. Combining eq. (16) and eq. (17) finishes the proof.

heading: E CONVERGENCE PROOFS FOR AID-BIO IN SECTION 4.1
text: For notation simplification, we define the following quantities.
Γ =3L2 + 3τ2M2
µ2 + 6L2
( 1 + √ κ )2( κ+ ρM µ2 )2 , δT,N = Γ(1− αµ)T + 6L2κ (√κ− 1√ κ+ 1 )2N Ω =8 ( βκ2 + 2βML
µ2 +
2βLMκ
µ2
)2 , ∆0 = ‖y0 − y∗(x0)‖2 + ‖v∗0 − v0‖2. (18)
We first provide some supporting lemmas. The following lemma characterizes the Hypergradient estimation error ‖∇̂Φ(xk)−∇Φ(xk)‖, where ∇̂Φ(xk) is given by eq. (3) via implicit differentiation. Lemma 3. Suppose Assumptions 1, 2 and 3 hold. Then, we have
‖∇̂Φ(xk)−∇Φ(xk)‖2 ≤Γ(1− αµ)T ‖y∗(xk)− y0k‖2 + 6L2κ (√κ− 1√
κ+ 1
)2N ‖v∗k − v0k‖2.
where Γ is given by eq. (18).
Proof of Lemma 3. Based on the form of∇Φ(xk) given by Proposition 1, we have
‖∇̂Φ(xk)−∇Φ(xk)‖2 ≤3‖∇xf(xk, y∗(xk))−∇xf(xk, yTk )‖2 + 3‖∇x∇yg(xk, yTk )‖2‖v∗k − vNk ‖2
+ 3‖∇x∇yg(xk, y∗(xk))−∇x∇yg(xk, yTk )‖2‖v∗k‖2,
which, in conjunction with Assumptions 1, 2 and 3, yields
‖∇̂Φ(xk)−∇Φ(xk)‖2 ≤ 3L2‖y∗(xk)− yTk ‖2 + 3L2‖v∗k − vNk ‖2 + 3τ2‖v∗k‖2‖yTk − y∗(xk)‖2
(i) ≤3L2‖y∗(xk)− yTk ‖2 + 3L2‖v∗k − vNk ‖2 + 3τ2M2
µ2 ‖yTk − y∗(xk)‖2. (19)
where (i) follows from the fact that ‖v∗k‖ ≤ ‖(∇2yg(xk, y∗(xk)))−1‖‖∇yf(xk, y∗(xk))‖ ≤ Mµ .
For notation simplification, let v̂k = (∇2yg(xk, yTk ))−1∇yf(xk, yTk ). We next upper-bound ‖v∗k − vNk ‖ in eq. (19). Based on the convergence result of CG for the quadratic programing, e.g., eq. (17)
in Grazzi et al. (2020), we have ‖vNk − v̂k‖ ≤ √ κ (√ κ−1√ κ+1 )N ‖v0k − v̂k‖. Based on this inequality, we further have
‖v∗k − vNk ‖ ≤‖v∗k − v̂k‖+ ‖vNk − v̂k‖ ≤ ‖v∗k − v̂k‖+ √ κ (√κ− 1√
κ+ 1
)N ‖v0k − v̂k‖
≤ ( 1 + √ κ (√κ− 1√
κ+ 1
)N) ‖v∗k − v̂k‖+ √ κ (√κ− 1√
κ+ 1
)N ‖v∗k − v0k‖. (20)
Next, based on the definitions of v∗k and v̂k, we have
‖v∗k − v̂k‖ =‖(∇2yg(xk, yTk ))−1∇yf(xk, yTk )− (∇2yg(xk, y∗(xk))−1∇yf(xk, y∗(xk))‖ ≤ ( κ+ ρM
µ2
) ‖yTk − y∗(xk)‖. (21)
Combining eq. (19), eq. (20), eq. (21) yields ‖∇̂Φ(xk)−∇Φ(xk)‖2 ≤ ( 3L2 + 3τ2M2
µ2
) ‖y∗(xk)− yTk ‖2 + 6L2κ (√κ− 1√ κ+ 1 )2N ‖v∗k − v0k‖2
+ 6L2 ( 1 + √ κ (√κ− 1√
κ+ 1
)N)2( κ+ ρM
µ2
)2 ‖yTk − y∗(xk)‖2,
which, in conjunction with ‖yTk −y∗(xk)‖ ≤ (1−αµ) T 2 ‖y0k−y∗(xk)‖ and the notations in eq. (18), finishes the proof.
Lemma 4. Suppose Assumptions 1, 2 and 3 hold. Choose
T ≥ log (36κ(κ+ ρM µ2 )2 + 16(κ2 + 4LMκ µ2 )2β2Γ)/ log 1 1− α = Θ(κ)
N ≥1 2 log(8κ+ 48(κ2 + 2ML µ2 + 2LMκ µ2 )2β2L2κ)/ log √ κ+ 1√ κ− 1 = Θ( √ κ), (22)
where Γ is given by eq. (18). Then, we have
‖y0k − y∗(xk)‖2+‖v∗k − v0k‖2 ≤ (1
2
)k ∆0 + Ω k−1∑ j=0 (1 2 )k−1−j ‖∇Φ(xj)‖2, (23)
where Ω and ∆0 are given by eq. (18).
Proof of Lemma 4. Recall that y0k = yTk−1. Then, we have
‖y0k − y∗(xk)‖2 ≤2‖yTk−1 − y∗(xk−1)‖2 + 2‖y∗(xk)− y∗(xk−1)‖2
(i) ≤2(1− αµ)T ‖y0k−1 − y∗(xk−1)‖2 + 2κ2β2‖∇̂Φ(xk−1)‖2
≤2(1− αµ)T ‖y0k−1 − y∗(xk−1)‖2 + 4κ2β2‖∇Φ(xk−1)− ∇̂Φ(xk−1)‖2
+ 4κ2β2‖∇Φ(xk−1)‖2
(ii) ≤ ( 2(1− αµ)T + 4κ2β2Γ(1− αµ)T ) ‖y∗(xk−1)− y0k−1‖2
+ 24κ4L2β2 (√κ− 1√
κ+ 1
)2N ‖v∗k−1 − v0k−1‖2 + 4κ2β2‖∇Φ(xk−1)‖2. (24)
where (i) follows from Lemma 2.2 in Ghadimi & Wang (2018) and (ii) follows from Lemma 3. In addition, note that
‖v∗k − v0k‖2 =‖v∗k − vNk−1‖2 ≤ 2‖v∗k−1 − vNk−1‖2 + 2‖v∗k − v∗k−1‖2
(i) ≤4 ( 1 + √ κ )2( κ+ ρM
µ2
)2 (1− αµ)T ‖y0k−1 − y∗(xk−1)‖2
+ 4κ (√κ− 1√
κ+ 1
)2N ‖v∗k−1 − v0k−1‖2 + 2‖v∗k − v∗k−1‖2, (25)
where (i) follows from eq. (20). Combining eq. (25) with ‖v∗k−v∗k−1‖ ≤ (κ2+ 2MLµ2 + 2LMκ µ2 )‖xk− xk−1‖, we have
‖v∗k − v0k‖2 (i) ≤ ( 16κ ( κ+ ρM
µ2
)2 + 4 ( κ2 + 4LMκ
µ2
)2 β2Γ ) (1− αµ)T ‖y0k−1 − y∗(xk−1)‖2
+ ( 4κ+ 48 ( κ2 + 2ML
µ2 +
2LMκ
µ2
)2 β2L2κ )(√κ− 1√ κ+ 1 )2N ‖v∗k−1 − v0k−1‖2
+ 4 ( κ2 + 2ML
µ2 +
2LMκ
µ2
)2 β2‖∇Φ(xk−1)‖2, (26)
where (i) follows from Lemma 3. Combining eq. (24) and eq. (26) yields
‖y0k − y∗(xk)‖2 + ‖v∗k − v0k‖2 ≤ ( 18κ ( κ+ ρM
µ2
)2 + 8 ( κ2 + 4LMκ
µ2
)2 β2Γ ) (1− αµ)T ‖y0k−1 − y∗(xk−1)‖2
+ ( 4κ+ 24 ( κ2 + 2ML
µ2 +
2LMκ
µ2
)2 β2L2κ )(√κ− 1√ κ+ 1 )2N ‖v∗k−1 − v0k−1‖2
+ 8 ( κ2 + 2ML
µ2 +
2LMκ
µ2
)2 β2‖∇Φ(xk−1)‖2,
which, in conjunction with eq. (22),
‖y0k − y∗(xk)‖2 + ‖v∗k − v0k‖2 ≤ 1
2 (‖y0k−1 − y∗(xk−1)‖2 + ‖v∗k−1 − v0k−1‖2)
+ 8 ( βκ2 + 2βML
µ2 +
2βLMκ
µ2
)2 ‖∇Φ(xk−1)‖2. (27)
Telescoping eq. (27) over k and using the notations in eq. (18), we finish the proof.
Lemma 5. Under the same setting as in Lemma 4, we have
‖∇̂Φ(xk)−∇Φ(xk)‖2 ≤δT,N (1
2
)k ∆0 + δT,NΩ k−1∑ j=0 (1 2 )k−1−j ‖∇Φ(xj)‖2.
where δT,N , Ω and ∆0 are given by eq. (18).
Proof of Lemma 5. Based on Lemma 3, eq. (18) and using ab+cd ≤ (a+c)(b+d) for any positive a, b, c, d, we have
‖∇̂Φ(xk)−∇Φ(xk)‖2 ≤δT,N (‖y∗(xk)− y0k‖2 + ‖v∗k − v0k‖2),
which, in conjunction with Lemma 4, finishes the proof.

heading: E.1 PROOF OF THEOREM 1
text: In this subsection, provide the proof for Theorem 1 based on the supporting Lemma 5. Based on the smoothness of the function Φ(x) established in Lemma 2, we have
Φ(xk+1) ≤Φ(xk) + 〈∇Φ(xk), xk+1 − xk〉+ LΦ 2 ‖xk+1 − xk‖2
≤Φ(xk)− β〈∇Φ(xk), ∇̂Φ(xk)−∇Φ(xk)〉 − β‖∇Φ(xk)‖2 + β2LΦ‖∇Φ(xk)‖2
+ β2LΦ‖∇Φ(xk)− ∇̂Φ(xk)‖2 ≤Φ(xk)− (β
2 − β2LΦ
) ‖∇Φ(xk)‖2 + (β 2 + β2LΦ ) ‖∇Φ(xk)− ∇̂Φ(xk)‖2, (28)
which, combined with Lemma 5, yields Φ(xk+1) ≤Φ(xk)− (β
2 − β2LΦ
) ‖∇Φ(xk)‖2 + (β 2 + β2LΦ ) δT,N (1 2 )k ∆0
+ (β
2 + β2LΦ
) δT,NΩ k−1∑ j=0 (1 2 )k−1−j ‖∇Φ(xj)‖2. (29)
Telescoping eq. (29) over k from 0 to K − 1 yields(β 2 − β2LΦ )K−1∑ k=0 ‖∇Φ(xk)‖2 ≤ Φ(x0)− inf x Φ(x) + (β 2 + β2LΦ ) δT,N∆0
+ (β
2 + β2LΦ
) δT,NΩ K−1∑ k=1 k−1∑ j=0 (1 2 )k−1−j ‖∇Φ(xj)‖2,
which, using the fact that ∑K−1 k=1 ∑k−1 j=0 ( 1 2 )k−1−j ‖∇Φ(xj)‖2 ≤ ∑K−1 k=0 1 2k ∑K−1 k=0 ‖∇Φ(xk)‖
2 ≤ 2 ∑K−1 k=0 ‖∇Φ(xk)‖
2, yields(β 2 − β2LΦ − ( βΩ + 2Ωβ2LΦ ) δT,N )K−1∑ k=0 ‖∇Φ(xk)‖2
≤ Φ(x0)− inf x
Φ(x) + (β
2 + β2LΦ
) δT,N∆0. (30)
Choose N and T such that ( Ω + 2ΩβLΦ ) δT,N ≤ 1
4 , δT,N ≤ 1. (31)
Note that based on the definition of δT,N in eq. (18), it suffices to choose T ≥ Θ(κ) andN ≥ Θ( √ κ) to satisfy eq. (31). Then, substituting eq. (31) into eq. (30) yields(β 4 − β2LΦ )K−1∑ k=0 ‖∇Φ(xk)‖2 ≤ Φ(x0)− inf x Φ(x) + (β 2 + β2LΦ ) ∆0, (32)
which, in conjunction with β ≤ 18LΦ , yields
1
K K−1∑ k=0 ‖∇Φ(xk)‖2 ≤ 64LΦ(Φ(x0)− infx Φ(x)) + 5∆0 K . (33)
In order to achieve an -accurate stationary point, we obtain from eq. (33) that AID-BiO requires at most the total number K = O(κ3 −1) of outer iterations. Then, based on eq. (3), we have the following complexity results.
• Gradient complexity: Gc(f, ) = 2K = O(κ3 −1),Gc(g, ) = KT = O ( κ4 −1 ) .
• Jacobian- and Hessian-vector product complexities: JV(g, ) = K = O ( κ3 −1 ) ,HV(g, ) = KN = O ( κ3.5 −1 ) .
Then, the proof is complete.

heading: F CONVERGENCE PROOFS FOR ITD-BIO IN SECTION 4.1
text: We first characterize an important estimation property of the outer-loop gradient estimator ∂f(xk,y T k ) ∂xk in ITD-BiO for approximating the true gradient∇Φ(xk) based on Proposition 2. Lemma 6. Suppose Assumptions 1, 2 and 3 hold. Choose α ≤ 1L . Then, we have∥∥∥∂f(xk, yTk )
∂xk −∇Φ(xk) ∥∥∥ ≤(L(L+ µ)(1− αµ)T2 µ + 2M (τµ+ Lρ) µ2 (1− αµ) T−1 2 ) ‖y0k − y∗(xk)‖
+ LM(1− αµ)T
µ . Lemma 6 shows that the gradient estimation error ∥∥∂f(xk,yTk )
∂xk − ∇Φ(xk) ∥∥ decays exponentially w.r.t. the number T of the inner-loop steps. We note that Grazzi et al. (2020) proved a similar result via a fixed point based approach. As a comparison, our proof of Lemma 6 directly characterizes the rate of the sequence ( ∂ytk ∂xk , t = 0, ..., T ) converging to ∂y ∗(xk) ∂xk
via the differentiation over all corresponding points along the inner-loop GD path as well as the optimality of the point y∗(xk).
Proof of Lemma 6. Using∇Φ(xk) = ∇xf(xk, y∗(xk)) + ∂y ∗(xk) ∂xk ∇yf(xk, y∗(xk)) and eq. (16) , and using the triangle inequality, we have∥∥∥∂f(xk, yTk )
∂xk −∇Φ(xk) ∥∥∥ =‖∇xf(xk, yTk )−∇xf(xk, y∗(xk))‖+ ∥∥∥∥∂yTk∂xk − ∂y ∗(xk) ∂xk
∥∥∥∥ ‖∇yf(xk, yTk )‖ + ∥∥∥∂y∗(xk)
∂xk ∥∥∥∥∥∇yf(xk, yTk )−∇yf(xk, y∗(xk))∥∥ (i)
≤L‖yTk − y∗(xk)‖+M ∥∥∥∥∂yTk∂xk − ∂y ∗(xk) ∂xk ∥∥∥∥+ L∥∥∥∂y∗(xk)∂xk ∥∥∥‖yTk − y∗(xk)‖, (34)
where (i) follows from Assumption 2. Our next step is to upper-bound ∥∥∥∂yTk∂xk − ∂y∗(xk)∂xk ∥∥∥ in eq. (34).
Based on the updates ytk = y t−1 k −α∇yg(xk, y t−1 k ) for t = 1, ..., T in ITD-BiO and using the chain rule, we have ∂ytk ∂xk = ∂yt−1k ∂xk − α ( ∇x∇yg(xk, yt−1k ) + ∂yt−1k ∂xk ∇2yg(xk, yt−1k ) ) . (35) Based on the optimality of y∗(xk), we have ∇yg(xk, y∗(xk)) = 0, which, in conjunction with the implicit differentiation theorem, yields
∇x∇yg(xk, y∗(xk)) + ∂y∗(xk)
∂xk ∇2yg(xk, y∗(xk)) = 0. (36)
Substituting eq. (36) into eq. (35) yields ∂ytk ∂xk − ∂y ∗(xk) ∂xk
= ∂yt−1k ∂xk − ∂y ∗(xk) ∂xk − α
( ∇x∇yg(xk, yt−1k ) +
∂yt−1k ∂xk
∇2yg(xk, yt−1k ) )
+ α ( ∇x∇yg(xk, y∗(xk)) + ∂y∗(xk)
∂xk ∇2yg(xk, y∗(xk)) ) = ∂yt−1k ∂xk − ∂y ∗(xk) ∂xk − α ( ∇x∇yg(xk, yt−1k )−∇x∇yg(xk, y ∗(xk)) )
− α ( ∂yt−1k ∂xk − ∂y ∗(xk) ∂xk ) ∇2yg(xk, yt−1k ) + α ∂y∗(xk)
∂xk
( ∇2yg(xk, y∗(xk))−∇2yg(xk, yt−1k ) ) . (37)
Combining eq. (36) and Assumption 2 yields∥∥∥∥∂y∗(xk)∂xk ∥∥∥∥ = ∥∥∥∇x∇yg(xk, y∗(xk)) [∇2yg(xk, y∗(xk))]−1∥∥∥ ≤ Lµ . (38)
Then, combining eq. (37) and eq. (38) yields∥∥∥ ∂ytk ∂xk − ∂y ∗(xk) ∂xk ∥∥∥ (i)≤∥∥∥I − α∇2yg(xk, yt−1k )∥∥∥∥∥∥∂yt−1k∂xk − ∂y ∗(xk) ∂xk ∥∥∥ + α ( τ + Lρ
µ
) ‖yt−1k − y ∗(xk)‖
(ii) ≤ (1− αµ) ∥∥∥∂yt−1k ∂xk − ∂y ∗(xk) ∂xk ∥∥∥+ α(τ + Lρ µ ) ‖yt−1k − y ∗(xk)‖, (39)
where (i) follows from Assumption 3 and (ii) follows from the strong-convexity of g(x, ·). Based on the strong-convexity of the lower-level function g(x, ·), we have
‖yt−1k − y ∗(xk)‖ ≤ (1− αµ) t−1 2 ‖y0k − y∗(xk)‖. (40)
Substituting eq. (40) into eq. (39) and telecopting eq. (39) over t from 1 to T , we have∥∥∥∂yTk ∂xk − ∂y ∗(xk) ∂xk ∥∥∥ ≤(1− αµ)T∥∥∥ ∂y0k ∂xk − ∂y ∗(xk) ∂xk ∥∥∥ + α ( τ + Lρ
µ ) T−1∑ t=0 (1− αµ)T−1−t(1− αµ) t2 ‖y0k − y∗(xk)‖
=(1− αµ)T ∥∥∥ ∂y0k ∂xk − ∂y ∗(xk) ∂xk ∥∥∥+ 2 (τµ+ Lρ) µ2 (1− αµ) T−1 2 ‖y0k − y∗(xk)‖
≤L(1− αµ) T
µ +
2 (τµ+ Lρ)
µ2 (1− αµ)
T−1 2 ‖y0k − y∗(xk)‖, (41)
where the last inequality follows from ∂y 0 k
∂xk = 0 and eq. (38). Then, combining eq. (34), eq. (38),
eq. (40) and eq. (41) completes the proof.

heading: F.1 PROOF OF THEOREM 2
text: Based on the characterization on the estimation error of the gradient estimate ∂f(xk,y T k )
∂xk in Lemma 6,
we now prove Theorem 2.
Recall the notation that ∇̂Φ(xk) = ∂f(xk,y T k )
∂xk . Using an approach similar to eq. (28), we have Φ(xk+1) ≤Φ(xk)− (β
2 − β2LΦ
) ‖∇Φ(xk)‖2 + (β 2 + β2LΦ ) ‖∇Φ(xk)− ∇̂Φ(xk)‖2, (42)
which, in conjunction with Lemma 6 and use ‖y0k − y∗(xk)‖2 ≤ ∆, yields Φ(xk+1) ≤Φ(xk)− (β
2 − β2LΦ
) ‖∇Φ(xk)‖2
+ 3∆ (β
2 + β2LΦ )(L2(L+ µ)2 µ2 (1− αµ)T + 4M 2 (τµ+ Lρ) 2 µ4 (1− αµ)T−1 ) + 3 (β
2 + β2LΦ )L2M2(1− αµ)2T µ2 . (43)
Telescoping eq. (43) over k from 0 to K − 1 yields
1
K K−1∑ k=0 (1 2 − βLΦ ) ‖∇Φ(xk)‖2 ≤ Φ(x0)− infx Φ(x) βK + 3 (1 2 + βLΦ )L2M2(1− αµ)2T µ2
+3∆ (1
2 + βLΦ )(L2(L+ µ)2 µ2 (1− αµ)T + 4M 2 (τµ+ Lρ) 2 µ4 (1− αµ)T−1 ) . (44)
Substuting β = 14LΦ and T = log ( max { 3LM µ , 9∆L 2(1+ Lµ ) 2, 36∆M 2(τµ+Lρ)2 (1−αµ)µ4 } 9 2 ) / log 11−αµ = Θ(κ log 1 ) in eq. (44) yields
1
K K−1∑ k=0 ‖∇Φ(xk)‖2 ≤ 16LΦ(Φ(x0)− infx Φ(x)) K + 2 3 . (45)
In order to achieve an -accurate stationary point, we obtain from eq. (45) that ITD-BiO requires at most the total number K = O(κ3 −1) of outer iterations. Then, based on the gradient form given by Proposition 2, we have the following complexity results.
• Gradient complexity: Gc(f, ) = 2K = O(κ3 −1),Gc(g, ) = KT = O ( κ4 −1 log 1 ) .
• Jacobian- and Hessian-vector product complexities:
JV(g, ) = KT = O ( κ4 −1 log 1 ) ,HV(g, ) = KT = O ( κ4 −1 log 1 ) .
Then, the proof is complete.

heading: G PROOFS OF MAIN RESULTS FOR STOCHASTIC CASE IN SECTION 4.2
text: In this section, we provide proofs for the convergence and complexity results of the proposed algorithm for the stochastic case.

heading: G.1 PROOF OF PROPOSITION 3
text: Based on the definition of vQ in eq. (5) and conditioning on xk, yTk , we have
EvQ =Eη Q−1∑ q=−1 Q∏ j=Q−q (I − η∇2yG(xk, yTk ;Bj))∇yF (xk, yTk ;DF ),
= η Q∑ q=0 (I − η∇2yg(xk, yTk ))q∇yf(xk, yTk )
= η ∞∑ q=0 (I − η∇2yg(xk, yTk ))q∇yf(xk, yTk )− η ∞∑ q=Q+1 (I − η∇2yg(xk, yTk ))q∇yf(xk, yTk )
= η(η∇2yg(xk, yTk ))−1∇yf(xk, yTk )− η ∞∑
q=Q+1
(I − η∇2yg(xk, yTk ))q∇yf(xk, yTk ),
which, in conjunction with the strong-convexity of function g(x, ·), yields
∥∥EvQ − [∇2yg(xk, yTk )]−1∇yf(xk, yTk )∥∥ ≤ η ∞∑ q=Q+1 (1− ηµ)qM ≤ (1− ηµ) Q+1M µ . (46)
This finishes the proof for the estimation bias. We next prove the variance bound. Note that
E ∥∥∥∥η Q−1∑
q=−1 Q∏ j=Q−q (I − η∇2yG(xk, yTk ;Bj))∇yF (xk, yTk ;DF )− (∇2yg(xk, yTk ))−1∇yf(xk, yTk ) ∥∥∥∥2
(i) ≤2E ∥∥∥∥ η Q−1∑
q=−1 Q∏ j=Q−q (I − η∇2yG(xk, yTk ;Bj))− (∇2yg(xk, yTk ))−1 ∥∥∥∥2M2 + 2M2µ2Df
≤4E ∥∥∥∥η Q−1∑
q=−1 Q∏ j=Q−q (I − η∇2yG(xk, yTk ;Bj))− η Q∑ q=0 (I − η∇2yg(xk, yTk ))q ∥∥∥∥2M2
+ 4E ∥∥∥∥η Q∏
q=0
(I − η∇2yg(xk, yTk ))q)− (∇2yg(xk, yTk ))−1 ∥∥∥∥2M2 + 2M2µ2Df
(ii) ≤ 4η2E ∥∥∥∥ Q∑ q=0 Q∏ j=Q+1−q (I − η∇2yG(xk, yTk ;Bj))− Q∑ q=0 (I − η∇2yg(xk, yTk ))q ∥∥∥∥2M2 + 4(1− ηµ)2Q+2M2µ2 + 2M2µ2Df
(iii) ≤ 4η2M2QE Q∑ q=0 ∥∥∥∥ Q∏ j=Q+1−q (I − η∇2yG(xk, yTk ;Bj))− (I − η∇2yg(xk, yTk ))q ∥∥∥∥2︸ ︷︷ ︸
Mq
+ 4(1− ηµ)2Q+2M2
µ2 +
2M2 µ2Df (47)
where (i) follows from Lemma 1, (ii) follows from eq. (46), and (iii) follows from the CauchySchwarz inequality.
Our next step is to upper-bound Mq in eq. (47). For simplicity, we define a general quantity Mi for by replacing q in Mq with i. Then, we have
EMi =E ∥∥∥∥(I − η∇2yg(xk, yTk )) Q∏
j=Q+2−i (I − η∇2yG(xk, yTk ;Bj))− (I − η∇2yg(xk, yTk ))i
∥∥∥∥2
+ E ∥∥∥∥η(∇2yg(xk, yTk )−∇2yG(xk, yTk ;BQ+1−i)) Q∏
j=Q+2−i (I − η∇2yG(xk, yTk ;Bj))
∥∥∥∥2
+ 2E 〈 (I − η∇2yg(xk, yTk )) Q∏
j=Q+2−i (I − η∇2yG(xk, yTk ;Bj))− (I − η∇2yg(xk, yTk ))i,
η(∇2yg(xk, yTk )−∇2yG(xk, yTk ;BQ+1−i)) Q∏
j=Q+2−i (I − η∇2yG(xk, yTk ;Bj)) 〉 (i) =E ∥∥∥∥(I − η∇2yg(xk, yTk )) Q∏
j=Q+2−i (I − η∇2yG(xk, yTk ;Bj))− (I − η∇2yg(xk, yTk ))i
∥∥∥∥2
+ E ∥∥∥∥η(∇2yg(xk, yTk )−∇2yG(xk, yTk ;BQ+1−i)) Q∏
j=Q+2−i (I − η∇2yG(xk, yTk ;Bj)) ∥∥∥∥2 (ii)
≤ (1− ηµ)2EMi−1 + η2(1− ηµ)2i−2E‖∇2yg(xk, yTk )−∇2yG(xk, yTk ;BQ+1−i)‖2
(iii) ≤ (1− ηµ)2EMi−1 + η2(1− ηµ)2i−2 L2
|BQ+1−i| , (48)
where (i) follows from that fact that EBQ+1−i∇2yG(xk, yTk ;BQ+1−i) = ∇2yg(xk, yTk ), (ii) follows from the strong-convexity of function G(x, ·; ξ), and (iii) follows from Lemma 1.
Then, telescoping eq. (48) over i from 2 to q yields
EMq ≤ L2η2(1− ηµ)2q−2 q∑ j=1
1
|BQ+1−j | ,
which, in conjunction with the choice of |BQ+1−j | = BQ(1− ηµ)j−1 for j = 1, ..., Q, yields
EMq ≤η2(1− ηµ)2q−2 q∑ j=1 L2 BQ ( 1 1− ηµ )j−1
= η2L2
BQ (1− ηµ)2q−2
( 1
1−ηµ
)q−1 − 1
1 1−ηµ − 1
≤ ηL 2 (1− ηµ)µ 1 BQ (1− ηµ)q. (49)
Substituting eq. (49) into eq. (47) yields
E ∥∥∥∥η Q−1∑
q=−1 Q∏ j=Q−q (I − η∇2yG(xk, yTk ;Bj))∇yF (xk, yTk ;DF )− (∇2yg(xk, yTk ))−1∇yf(xk, yTk ) ∥∥∥∥2
≤4η2M2Q Q∑ q=0
ηL2
(1− ηµ)µ 1 BQ (1− ηµ)q + 4(1− ηµ) 2Q+2M2 µ2 + 2M2 µ2Df
≤4η 2L2M2 µ2 1 B + 4(1− ηµ)2Q+2M2 µ2 + 2M2 µ2Df , (50)
where the last inequality follows from the fact that ∑S q=0 x q ≤ 11−x . Then, the proof is complete.

heading: G.2 AUXILIARY LEMMAS FOR PROVING THEOREM 3
text: We first use the following lemma to characterize the first-moment error of the gradient estimate ∇̂Φ(xk), whose form is given by eq. (6). Lemma 7. Suppose Assumptions 1, 2 and 3 hold. Then, conditioning on xk and yTk , we have∥∥E∇̂Φ(xk)−∇Φ(xk)∥∥2 ≤ 2(L+ L2
µ + Mτ µ + LMρ µ2
)2 ‖yTk − y∗(xk)‖2 + 2L2M2(1− ηµ)2Q
µ2 .
Proof of Lemma 7. To simplify notations, we define
∇̃ΦT (xk) = ∇xf(xk, yTk )−∇x∇yg(xk, yTk ) [ ∇2yg(xk, yTk ) ]−1∇yf(xk, yTk ). (51) Based on the definition of ∇̂Φ(xk) in eq. (6) and conditioning on xk and yTk , we have
E∇̂Φ(xk) =∇xf(xk, yTk )−∇x∇yg(xk, yTk )EvQ =∇̃ΦT (xk)−∇x∇yg(xk, yTk )(EvQ − [∇2yg(xk, yTk )]−1∇yf(xk, yTk )),
which further implies that∥∥E∇̂Φ(xk)−∇Φ(xk)∥∥2 ≤2E‖∇̃ΦT (xk)−∇Φ(xk)‖2 + 2‖E∇̂Φ(xk)− ∇̃ΦT (xk)‖2
≤2E‖∇̃ΦT (xk)−∇Φ(xk)‖2 + 2L2‖EvQ − [∇2yg(xk, yTk )]−1∇yf(xk, yTk )‖2
≤2E‖∇̃ΦT (xk)−∇Φ(xk)‖2 + 2L2M2(1− ηµ)2Q+2
µ2 , (52)
where the last inequality follows from Proposition 3. Our next step is to upper-bound the first term at the right hand side of eq. (52). Using the fact that ∥∥∇2yg(x, y)−1∥∥ ≤ 1µ and based on Assumptions 2
and 3, we have
‖∇̃ΦT (xk)−∇Φ(xk)‖ ≤‖∇xf(xk, yTk )−∇xf(xk, y∗(xk))‖
+ L2
µ ‖yTk − y∗(xk)‖+
Mτ
µ ‖yTk − y∗(xk)‖
+ LM ∥∥∇2yg(xk, yTk )−1 −∇2yg(xk, y∗(xk))−1∥∥
≤ ( L+ L2
µ + Mτ µ + LMρ µ2
) ‖yTk − y∗(xk)‖, (53)
where the last inequality follows from the inequality ‖M−11 −M −1 2 ‖ ≤ ‖M −1 1 M −1 2 ‖‖M1 −M2‖ for any two matrices M1 and M2. Combining eq. (52) and eq. (53) yields∥∥E∇̂Φ(xk)−∇Φ(xk)∥∥2 ≤ 2(L+ L2 µ + Mτ µ + LMρ µ2 )2 ‖yTk − y∗(xk)‖2 + 2L2M2(1− ηµ)2Q µ2 ,
which completes the proof.
Then, we use the following lemma to characterize the variance of the estimator ∇̂Φ(xk). Lemma 8. Suppose Assumptions 1, 2 and 3 hold. Then, we have
E‖∇̂Φ(xk)−∇Φ(xk)‖2 ≤ 4L2M2 µ2Dg + (8L2 µ2 + 2 )M2 Df + 16η2L4M2 µ2 1 B + 16L2M2(1− ηµ)2Q µ2
+ ( L+ L2
µ + Mτ µ + LMρ µ2
)2 E‖yTk − y∗(xk)‖2.
Proof of Lemma 8. Based on the definitions of ∇Φ(xk) and ∇̃ΦT (xk) in eq. (4) and eq. (51) and conditioning on xk and yTk , we have
E‖∇̂Φ(xk)−∇Φ(xk)‖2
(i) =E‖∇̂Φ(xk)− ∇̃ΦT (xk)‖2 + ‖∇̃ΦT (xk)−∇Φ(xk)‖2
(ii) ≤ 2E ∥∥∇x∇yG(xk, yTk ;DG)vQ −∇x∇yg(xk, yTk )[∇2yg(xk, yTk )]−1∇yf(xk, yTk )∥∥2 + 2M2Df
+ ( L+ L2
µ + Mτ µ + LMρ µ2
)2 ‖yTk − y∗(xk)‖2
(iii) ≤ 4M 2
µ2 E‖∇x∇yG(xk, yTk ;DG)−∇x∇yg(xk, yTk )‖2 + 4L2E‖vQ −
[ ∇2yg(xk, yTk ) ]−1∇yf(xk, yTk )‖2 + ( L+ L2
µ + Mτ µ + LMρ µ2
)2 ‖yTk − y∗(xk)‖2 + 2M2
Df , (54)
where (i) follows from the fact that EDG,DH ,DF ∇̂Φ(xk) = ∇̃ΦT (xk), (ii) follows from Lemma 1 and eq. (53), and (iii) follows from the Young’s inequality and Assumption 2.
Using Lemma 1 and Proposition 3 in eq. (54), yields
E‖∇̂Φ(xk)−∇Φ(xk)‖2 ≤ 4L2M2
µ2Dg +
16η2L4M2
µ2 1 B +
16(1− ηµ)2QL2M2
µ2 +
8L2M2
µ2Df + ( L+ L2
µ + Mτ µ + LMρ µ2
)2 ‖yTk − y∗(xk)‖2 + 2M2
Df , (55)
which, unconditioning on xk and yTk , completes the proof.
It can be seen from Lemmas 7 and 8 that the upper bounds on both the estimation error and bias depend on the tracking error ‖yTk − y∗(xk)‖2. The following lemma provides an upper bound on such tracking error ‖yTk − y∗(xk)‖2.
Lemma 9. Suppose Assumptions 1, 2 and 4 hold. Define constants λ = (L− µ L+ µ )2T( 2 + 4β2L2 µ2 ( L+ L2 µ + Mτ µ + LMρ µ2 )2) ∆ = 4L2M2
µ2Dg + (8L2 µ2 + 2 )M2 Df + 16η2L4M2 µ2 1 B + 16L2M2(1− ηµ)2Q µ2
ω = 4β2L2
µ2 (L− µ L+ µ )2T . (56)
Choose T such that λ < 1 and set inner-loop stepsize α = 2L+µ . Then, we have
E‖yTk − y∗(xk)‖2 ≤λk ((
L− µ L+ µ
)2T ‖y0 − y∗(x0)‖2 + σ2
LµS
) + ω
k−1∑ j=0 λk−1−jE‖∇Φ(xj)‖2 + ω∆ + σ 2 LµS 1− λ .
Proof of Lemma 9. First note that for an integer t ≤ T ‖yt+1k − y ∗(xk)‖2
=‖yt+1k − y t k‖2 + 2〈yt+1k − y t k, y t k − y∗(xk)〉+ ‖ytk − y∗(xk)‖2
=α2‖∇yG(xk, ytk;St)‖2 − 2α〈∇yG(xk, ytk;St), ytk − y∗(xk)〉+ ‖ytk − y∗(xk)‖2. (57) Conditioning on ytk and taking expectation in eq. (57), we have
E‖yt+1k − y ∗(xk)‖2
(i) ≤α2 (σ2 S + ‖∇yg(xk, ytk)‖2 ) − 2α〈∇yg(xk, ytk), ytk − y∗(xk)〉
+ ‖ytk − y∗(xk)‖2
(ii) ≤ α 2σ2
S + α2‖∇yg(xk, ytk)‖2 − 2α
( Lµ
L+ µ ‖ytk − y∗(xk)‖2 +
‖∇yg(xk, ytk)‖2
L+ µ ) + ‖ytk − y∗(xk)‖2
= α2σ2
S − α
( 2
L+ µ − α
) ‖∇yg(xk, ytk)‖2 + ( 1− 2αLµ
L+ µ
) ‖ytk − y∗(xk)‖2 (58)
where (i) follows from the third item in Assumption 2, (ii) follows from the strong-convexity and smoothness of the function g. Since α = 2L+µ , we obtain from eq. (58) that
E‖yt+1k − y ∗(xk)‖2 ≤ ( L− µ L+ µ )2 ‖ytk − y∗(xk)‖2 +
4σ2
(L+ µ)2S . (59)
Unconditioning on ytk in eq. (59) and telescoping eq. (59) over t from 0 to T − 1 yields E‖yTk − y∗(xk)‖2 ≤ ( L− µ L+ µ )2T E‖y0k − y∗(xk)‖2 + σ2 LµS
= ( L− µ L+ µ )2T E‖yTk−1 − y∗(xk)‖2 + σ2 LµS , (60)
where the last inequality follows from Algorithm 2 that y0k = y T k−1. Note that
E‖yTk−1 − y∗(xk)‖2 ≤2E‖yTk−1 − y∗(xk−1)‖2 + 2E‖y∗(xk−1)− y∗(xk)‖2
(i) ≤2E‖yTk−1 − y∗(xk−1)‖2 + 2L2
µ2 E‖xk − xk−1‖2
≤2E‖yTk−1 − y∗(xk−1)‖2 + 2β2L2
µ2 E‖∇̂Φ(xk−1)‖2
≤2E‖yTk−1 − y∗(xk−1)‖2 + 4β2L2
µ2 E‖∇Φ(xk−1)‖2
+ 4β2L2
µ2 E‖∇̂Φ(xk−1)−∇Φ(xk−1)‖2, (61)
where (i) follows from Lemma 2.2 in Ghadimi & Wang (2018). Using Lemma 8 in eq. (61) yields
E‖yTk−1 − y∗(xk)‖2 ≤ ( 2 + 4β2L2
µ2
( L+ L2
µ + Mτ µ + LMρ µ2
)2) E‖yTk−1 − y∗(xk−1)‖2 + 4β2L2
µ2 E‖∇Φ(xk−1)‖2
+ 4β2L2
µ2
( 4L2M2 µ2Dg + (8L2 µ2 + 2 )M2 Df + 16η2L4M2 µ2 1 B + 16L2M2(1− ηµ)2Q µ2 ) . (62)
Combining eq. (60) and eq. (62) yields
E‖yTk − y∗(xk)‖2 ≤ (L− µ L+ µ )2T( 2 + 4β2L2 µ2 ( L+ L2 µ + Mτ µ + LMρ µ2 )2) E‖yTk−1 − y∗(xk−1)‖2
+ (L− µ L+ µ )2T 4β2L2 µ2 ( 4L2M2 µ2Dg + (8L2 µ2 + 2 )M2 Df + 16η2L4M2 µ2 1 B + 16L2M2(1− ηµ)2Q µ2 ) + 4β2L2
µ2 (L− µ L+ µ )2T E‖∇Φ(xk−1)‖2 + σ2 LµS . (63)
Based on the definitions of λ, ω,∆ in eq. (56), we obtain from eq. (63) that
E‖yTk − y∗(xk)‖2 ≤λE‖yTk−1 − y∗(xk−1)‖2 + ω∆ + σ2
LµS + ωE‖∇Φ(xk−1)‖2. (64)
Telescoping eq. (64) over k yields
E‖yTk − y∗(xk)‖2
≤λkE‖yT0 − y∗(x0)‖2 + ω k−1∑ j=0 λk−1−jE‖∇Φ(xj)‖2 + ω∆ + σ 2 LµS 1− λ
≤λk ((
L− µ L+ µ
)2T ‖y0 − y∗(x0)‖2 + σ2
LµS
) + ω
k−1∑ j=0 λk−1−jE‖∇Φ(xj)‖2 + ω∆ + σ 2 LµS 1− λ ,
which completes the proof.

heading: G.3 PROOF OF THEOREM 3
text: In this subsection, we provide the proof for Theorem 3, based on the supporting lemmas we develop in Appendix G.2.
Based on the smoothness of the function Φ(x) in Lemma 2, we have
Φ(xk+1) ≤Φ(xk) + 〈∇Φ(xk), xk+1 − xk〉+ LΦ 2 ‖xk+1 − xk‖2
≤Φ(xk)− β〈∇Φ(xk), ∇̂Φ(xk)〉+ β2LΦ‖∇Φ(xk)‖2 + β2LΦ‖∇Φ(xk)− ∇̂Φ(xk)‖2. For simplicity, let Ek = E(· |xk, yTk ). Note that we choose β = 14Lφ . Then, taking expectation over the above inequality, we have
EΦ(xk+1) ≤EΦ(xk)− βE〈∇Φ(xk),Ek∇̂Φ(xk)〉+ β2LΦE‖∇Φ(xk)‖2
+ β2LΦE‖∇Φ(xk)− ∇̂Φ(xk)‖2
(i) ≤EΦ(xk) + β
2 E‖Ek∇̂Φ(xk)−∇Φ(xk)‖2 −
β 4 E‖∇Φ(xk)‖2 + β 4 E‖∇Φ(xk)− ∇̂Φ(xk)‖2
(ii)
≤ EΦ(xk)− β
4 E‖∇Φ(xk)‖2 +
βL2M2(1− ηµ)2Q
µ2
+ β
4
( 4L2M2 µ2Dg + (8L2 µ2 + 2 )M2 Df + 16η2L4M2 µ2 1 B + 16L2M2(1− ηµ)2Q µ2 ) + 5β
4
( L+ L2
µ + Mτ µ + LMρ µ2
)2 E‖yTk − y∗(xk)‖2 (65)
where (i) follows from Cauchy-Schwarz inequality, and (ii) follows from Lemma 7 and Lemma 8. To simplify notations, Let
ν = 5
4
( L+ L2
µ + Mτ µ + LMρ µ2
)2 . (66)
Then, applying Lemma 9 in eq. (65) and using the definitions of ω,∆, λ in eq. (56), we have
EΦ(xk+1) ≤EΦ(xk)− β
4 E‖∇Φ(xk)‖2 +
βL2M2(1− ηµ)2Q
µ2
+ β
4 ∆ + βνλk (( L− µ L+ µ )2T ‖y0 − y∗(x0)‖2 + σ2 LµS )
+ βνω k−1∑ j=0 λk−1−jE‖∇Φ(xj)‖2 + βν(ω∆ + σ 2 LµS ) 1− λ ,
Telescoping the above inequality over k from 0 to K − 1 yields
EΦ(xK) ≤ Φ(x0)− β
4 K−1∑ k=0 E‖∇Φ(xk)‖2 + βνω K−1∑ k=1 k−1∑ j=0 λk−1−jE‖∇Φ(xj)‖2
+ Kβ∆ 4 + ((L− µ L+ µ )2T ‖y0 − y∗(x0)‖2 + σ2 LµS ) βν 1− λ
+ KβL2M2(1− ηµ)2Q µ2 + Kβν(ω∆ + σ
2
LµS )
1− λ ,
which, using the fact that
K−1∑ k=1 k−1∑ j=0 λk−1−jE‖∇Φ(xj)‖2 ≤ ( K−1∑ k=0 λk ) K−1∑ k=0 E‖∇Φ(xk)‖2 < 1 1− λ K−1∑ k=0 E‖∇Φ(xk)‖2,
yields
(1 4 − νω 1− λ ) 1 K K−1∑ k=0 E‖∇Φ(xk)‖2
≤Φ(x0)− infx Φ(x) βK
+ ν ( (L−µL+µ ) 2T ‖y0 − y∗(x0)‖2 + σ 2 LµS ) K(1− λ) + ∆ 4 + L2M2(1− ηµ)2Q µ2
+ ν(ω∆ + σ
2
LµS )
1− λ . (67)
Since β = 14LΦ and T ≥ log ( 12+ 48β 2L2 µ2 (L+L 2 µ + Mτ µ + LMρ µ2 )2 ) 2 log(L+µL−µ ) , we have λ ≤ 16 , and hence eq. (67) is further simplified to
(1 4 −6 5 νω ) 1 K K−1∑ k=0 E‖∇Φ(xk)‖2
≤Φ(x0)− infx Φ(x) βK
+ 2ν ( (L−µL+µ ) 2T ‖y0 − y∗(x0)‖2 + σ 2 LµS ) K + ∆ 4 + L2M2(1− ηµ)2Q µ2
+ 2ν ( ω∆ + σ2
LµS
) . (68)
By the definitions of ω in eq. (56) and ν in eq. (66) and T ≥ log ( 12+ 48β 2L2 µ2 (L+L 2 µ + Mτ µ + LMρ µ2 )2 )
2 log(L+µL−µ ) , we
have
νω = 5β2L2
µ2 (L− µ L+ µ )2T( L+ L2 µ + Mτ µ + LMρ µ2 )2 < 5β2L2 µ2 ( L+ L 2 µ + Mτ µ + LMρ µ2 )2 12 + 48β 2L2
µ2 (L+ L2 µ + Mτ µ + LMρ µ2 )
2 ≤ 5 48 . (69)
In addition, since T > log (√ β ( L+L 2 µ + Mτ µ + LMρ µ2 )) log(L+µL−µ ) , we have
ν (L− µ L+ µ )2T = 5 4 (L− µ L+ µ )2T( L+ L2 µ + Mτ µ + LMρ µ2 )2 < 5 4β . (70)
Substituting eq. (69) and eq. (70) in eq. (68) yields
1
K K−1∑ k=0 E‖∇Φ(xk)‖2 ≤ 8(Φ(x0)− infx Φ(x) + 52‖y0 − y ∗(x0)‖2) βK + ( 1 + 1 K )16νσ2 LµS
+ 11
3 ∆ +
8L2M2
µ2 (1− ηµ)2Q, (71)
which, in conjunction with eq. (56) and eq. (66), yields eq. (10) in Theorem 3.
Then, based on eq. (10), in order to achieve an -accurate stationary point, i.e., E‖∇Φ(x̄)‖2 ≤ with x̄ chosen from x0, ..., xK−1 uniformly at random, it suffices to choose
K = 32LΦ(Φ(x0)− infx Φ(x) + 52‖y0 − y ∗(x0)‖2) = O (κ3 ) , T = Θ(κ)
Q =κ log κ2 , S = O
(κ5 ) , Dg = O ( κ2 ) , Df = O ( κ2 ) , B = O ( κ2 ) .
Note that the above choices of Q and B satisfy the condition that B ≥ 1 Q(1−ηµ)Q−1 required in Proposition 3.
Then, the gradient complexity is given by Gc(F, ) = KDf = O(κ5 −2),Gc(G, ) = KTS = O(κ9 −2). In addition, the Jacobian- and Hessian-vector product complexities are given by JV(G, ) = KDg = O(κ5 −2) and
HV(G, ) = K Q∑ j=1 BQ(1− ηµ)j−1 = KBQ ηµ ≤ O ( κ6 2 log κ2 ) .
Then, the proof is complete.

heading: H PROOF OF THEOREM 4 ON META-LEARNING
text: To prove Theorem 4, we first establish the following lemma to characterize the estimation variance EB ∥∥∂LD(φk,w̃Tk ;B)
∂φk − ∂LD(φk,w̃ T k ) ∂φk ∥∥2, where w̃Tk is the output of T inner-loop steps of gradient descent at the kth outer loop.
Lemma 10. Suppose Assumptions 2 and 3 are satisfied and suppose each task loss LSi(φ,wi) is µ-strongly-convex w.r.t. wi. Then, we have
EB ∥∥∥∂LD(φk, w̃Tk ;B)
∂φk − ∂LD(φk, w̃ T k ) ∂φk ∥∥∥2 ≤ (1 + L µ )2M2 |B| .
Proof. Let w̃Tk = (w T 1,k, ..., w T m,k) be the output of T inner-loop steps of gradient descent at the k th outer loop. Using Proposition 2, we have, for task Ti,∥∥∥∂LDi(φk, wTi,k) ∂φk
∥∥∥ ≤‖∇φLDi(φk, wTi,k)‖ + ∥∥∥α T−1∑
t=0 ∇φ∇wiLSi(φk, wti,k) T−1∏ j=t+1 (I − α∇2wiLSi(φk, w j i,k))∇wiLDi(φk, w T i.k) ∥∥∥
(i) ≤M + αLM T−1∑ t=0 (1− αµ)T−t−1 = M + LM µ , (72)
where (i) follows from assumptions 2 and strong-convexity of LSi(φ, ·). Then, using the definition of LD(φ, w̃;B) = 1|B| ∑ i∈B LDi(φ,wi), we have
EB ∥∥∥∂LD(φk, w̃Tk ;B)
∂φk − ∂LD(φk, w̃ T k ) ∂φk ∥∥∥2 = 1|B|Ei∥∥∥∂LDi(φk, wTi,k)∂φk − ∂LD(φk, w̃ T k ) ∂φk ∥∥∥2 (i)
≤ 1 |B|
Ei ∥∥∥∂LDi(φk, wTi,k)
∂φk ∥∥∥2 (ii)
≤ ( 1 + L
µ )2M2 |B| . (73)
where (i) follows from Ei ∂LDi (φk,w T i,k)
∂φk = ∂LD(φk,w̃Tk ) ∂φk
and (ii) follows from eq. (72). Then, the proof is complete.
Proof of Theorem 4. Recall Φ(φ) := LD(φ, w̃∗(φ)) be the objective function, and let ∇̂Φ(φk) = ∂LD(φk,w̃Tk )
∂φk . Using an approach similar to eq. (42), we have
Φ(φk+1) ≤Φ(φk) + 〈∇Φ(φk), φk+1 − φk〉+ LΦ 2 ‖φk+1 − φk‖2
≤Φ(φk)− β 〈 ∇Φ(φk),
∂LD(φk, w̃Tk ;B) ∂φk
〉 + β2LΦ
2 ∥∥∥∂LD(φk, w̃Tk ;B) ∂φk ∥∥∥2. (74) Taking the expectation of eq. (74) yields
EΦ(φk+1) (i) ≤EΦ(φk)− βE 〈 ∇Φ(φk), ∇̂Φ(φk) 〉 + β2LΦ
2 E‖∇̂Φ(φk)‖2
+ β2LΦ 2 E ∥∥∥∇̂Φ(φk)− ∂LD(φk, w̃Tk ;B)
∂φk ∥∥∥2 (ii)
≤ EΦ(φk)− βE 〈 ∇Φ(φk), ∇̂Φ(φk) 〉 + β2LΦ
2 E‖∇̂Φ(φk)‖2 + β2LΦ 2
( 1 + L
µ )2M2 |B|
≤EΦ(φk)− (β
2 − β2LΦ
) E‖∇Φ(φk)‖2 + (β 2 + β2LΦ ) E‖∇Φ(φk)− ∇̂Φ(φk)‖2
+ β2LΦ
2
( 1 + L
µ )2M2 |B| , (75)
where (i) follows from EBLD(φk, w̃Tk ;B) = LD(φk, w̃Tk ) and (ii) follows from Lemma 10. Using Lemma 6 in eq. (75) and rearranging the terms, we have
1
K K−1∑ k=0 (1 2 − βLΦ ) E‖∇Φ(φk)‖2
≤Φ(φ0)− infφ Φ(φ) βK
+ 3 (1
2 + βLΦ )L2M2(1− αµ)2T µ2 + βLΦ 2 ( 1 + L µ )2M2 |B|
+ 3∆ (1
2 + βLΦ )(L2(L+ µ)2 µ2 (1− αµ)T + 4M 2 (τµ+ Lρ) 2 µ4 (1− αµ)T−1 ) ,
where ∆ = maxk ‖w̃0k− w̃∗(φk)‖2 <∞. Choose the same parameters β, T as in Theorem 2. Then, we have
1
K K−1∑ k=0 E‖∇Φ(φk)‖2 ≤ 16LΦ(Φ(φ0)− infφ Φ(φ)) K + 2 3 + ( 1 + L µ )2 M2 8|B| .
Then, the proof is complete.

Review summaries:
The paper presents two algorithms - one for the deterministic and one for stochastic bilevel optimization. The paper claims the methods are lower cost in computational complexity for various terms and easy to implement. A finite-time convergence proof is provided for the algorithms.  Empirical results are presented for meta-learning, and (in the appendix) hyperparameter optimization.