ES-BASED JACOBIAN ENABLES FASTER BILEVEL OP-
Year: 2021

Abstract:
Bilevel optimization (BO) has arisen as a powerful tool for solving many modern machine learning problems. However, due to the nested structure of BO, existing gradient-based methods require second-order derivative approximations via Jacobianor/and Hessian-vector computations, which can be very costly in practice, especially with large neural network models. In this work, we propose a novel BO algorithm, which adopts Evolution Strategies (ES) based method to approximate the response Jacobian matrix in the hypergradient of BO, and hence fully eliminates all second-order computations. We call our algorithm as ESJ (which stands for the ES-based Jacobian method) and further extend it to the stochastic setting as ESJ-S. Theoretically, we characterize the convergence guarantee and computational complexity for our algorithms. Experimentally, we demonstrate the superiority of our proposed algorithms compared to the state of the art methods on various bilevel problems. Particularly, in our experiment in the few-shot metalearning problem, we meta-learn the twelve millions parameters of a ResNet-12 network over the miniImageNet dataset, which evidently demonstrates the scalability of our ES-based bilevel approach and its feasibility in the large-scale setting.

heading: 1 INTRODUCTION
text: Bilevel optimization has recently arisen as a powerful tool to capture various modern machine learning problems, including meta-learning (Bertinetto et al., 2018; Franceschi et al., 2018; Rajeswaran et al., 2019; Ji et al., 2020a; Liu et al., 2021a), hyperparamater optimization (Franceschi et al., 2018; Shaban et al., 2019), neural architecture search (Liu et al., 2018a; Zhang et al., 2021), etc. Bilevel optimization generally takes the following mathematical form:
min x∈Rp Φ(x) := f(x, y∗(x)) s.t. y∗(x) = arg min y∈Rd g(x, y), (1)
where the outer and inner objectives f : Rp×Rd → R and g : Rp×Rd → R are both continuously differentiable with respect to (w.r.t.) the inner and outer variables x ∈ Rp and y ∈ Rd. In many machine learning problems, the loss functions f and g take a finite-sum form over a set of given data Dn,m = {ξi, ζj , i = 1, ..., n, j = 1, ...,m} as follows:
f(x, y) = 1n ∑n i=1 F (x, y; ξi), g(x, y) = 1 m ∑m i=1G(x, y; ζi) (2)
where the sample sizes n and m are typically very large. In this work, we consider a popular setting where the total objective function Φ(x) in eq. (1) is possibly nonconvex w.r.t. x and the inner function g is strongly convex w.r.t. y.
Gradient-based methods have served as a popular tool for solving bilevel optimization problems. Two types of approaches have been widely used: the iterative differentiation (ITD) method (Domke, 2012; Maclaurin et al., 2015; Finn et al., 2017; Franceschi et al., 2017; Shaban et al., 2019; Rajeswaran et al., 2019; Liu et al., 2020a) and the approximate iterative differentiation (AID) method (Domke, 2012; Pedregosa, 2016; Gould et al., 2016; Liao et al., 2018; Lorraine et al., 2020). In the ITD method, the objective value f(x, ŷ(x)) is first computed using an approximation point ŷ(x) close to y∗(x), and the hypergradient∇Φ(x) = ∂f(x,y
∗(x)) ∂x is then approximated by ∂f(x,ŷ(x)) ∂x using
either reverse or forward mode automatic differentiation. Alternatively, the AID method leverages
the Implicit Function Theorem (IFT) to establish an implicit expression for the hypergradient. Due to the bilevel structure of the problem, these gradient-based ITD and AID approaches typically involve second-order matrix computations. Existing efficient implementations of these methods often adopt Jacobian- or/and Hessian-vector computations, which can still be very costly in practice for high-dimensional problems with deep neural networks.
To overcome the computational challenge of current gradient-based methods, Evolution Strategies (ES) can be a good candidate to eliminate the second-order computations, which uses function values for estimating its gradient in blackbox optimization (Nesterov & Spokoiny, 2017). An ES-based bilevel optimizer has recently been proposed in Gu et al. (2021) for hyperparameter optimization, which uses the ES method to approximate the hypergradient of Φ(x), treating Φ(x) fully as a blackbox objective. However, the hypergradient in bilevel optimization has a much more involved structure (see eq. (3)) than the gradient in standard minimization problems, and hence direct use of the ES method without exploiting the structure of the hypergradient can encounter a large estimation bias and result in poor performance in practice, especially in deep learning as demonstrated in our experiments in Section 4. This hence motivates the following intriguing question:
Can we design a better hypergradient estimator based on the ES method by leveraging the analytical structure of the hypergradient, and then use such an estimator to construct more efficient ES-based bilevel optimizers, which are scalable for high-dimensional problems?
Furthermore, since the hypergradient in bilevel optimization is much more challenging to estimate than the gradient in standard minimization problems, it is unclear whether the ES-based estimator can yield a provably convergent algorithm for bilevel optimization. So far the ES-based bilevel optimizer in Gu et al. (2021) did not come with a polynomial-time complexity guarantee. Our empirical experiments in Section 4 show that such an algorithm fail to converge even with shallow neural networks. This thus motivates the second question we ask:
Can the ES method lead to bilevel optimizers with theoretical finite-time convergence guarantee?
This paper provides the affirmative answers to both questions.

heading: 1.1 MAIN CONTRIBUTIONS
text: Novel ES-based Jacobian bilevel optimizer. We propose a novel bilevel optimizer with the ESbased Jacobian estimation, which we call as ESJ. In contrast to the existing ES-based optimizer in Gu et al. (2021), which estimates the hypergradient by treating the outer-level objective Φ(x) fully as a blackbox, ESJ estimates only the respone Jacobian matrix (i.e., gradient of y∗(x) with respect to x), which is the major computational bottleneck, and then leverages the analytical structure of the hypergradient to construct a much more accurate estimator. Further, ESJ has only gradient computations, and is computationally much more efficient than existing AID and ITD bilevel optimizers that require Hessian- and/or Jacobian-vector product computations. We further propose a stochastic bilevel optimizer ESJ-S to allow the algorithm scalable over large dataset.
Convergence guarantee. Theoretically, we characterize the convergence rate and their computational complexity for both ESJ and ESJ-S to achieve an -accurate solution. Technically, in contrast to the standard analysis of ES-based methods on the smoothed blackbox function values, we develop tools to analyze the ES estimator on the output of the execution of an inner optimizer. Such an analysis can be of independent interest for ES-based bilevel optimization.
Superior performance. Experimentally, our algorithms ESJ and ESJ-S achieve superior performance compared to the current state of the art bilevel optimizers. Our approach is not only efficient and scalable, but also stable and robust across various bilevel problems. In particular, on the few-shot meta-learning problem, we meta-learn the twelve millions weights of a ResNet-12 network. Previous applications of other bilevel algorithms to meta-learning were limited only to settings where the backbone network has only up to few hundreds of thousands of parameters.

heading: 1.2 RELATED WORK
text: Bilevel optimization. Bilevel optimization has been studied for decades since Bracken & McGill (1973). A variety of bilevel optimization algorithms were then developed, including constraintbased approaches (Hansen et al., 1992; Shi et al., 2005; Moore, 2010), approximate implicit differ-
entiation (AID) approaches (Domke, 2012; Pedregosa, 2016; Gould et al., 2016; Liao et al., 2018; Lorraine et al., 2020) and iterative differentiation (ITD) approaches (Domke, 2012; Maclaurin et al., 2015; Franceschi et al., 2017; Finn et al., 2017; Shaban et al., 2019; Rajeswaran et al., 2019; Liu et al., 2020a). These methods often suffer from expensive computations of second-order information (e.g., Hessian-vector products). Hessian-free algorithms were recently proposed based on interiorpoint method (Liu et al., 2021b) and Evolution Strategies (ES) (Gu et al., 2021). This paper proposes a more efficient ES-based approach via exploiting the benign structure of the hypergradient.
Recently, the convergence rate has been established for gradient-based (Hessian involved) bilevel algorithms (Grazzi et al., 2020; Ji et al., 2021; Rajeswaran et al., 2019; Ji & Liang, 2021). This paper provides the convergence analysis for the proposed Hessian-free ES-based approach.
Stochastic bilevel optimization. Ghadimi & Wang (2018); Ji et al. (2021); Hong et al. (2020) proposed stochastic gradient descent (SGD) type bilevel optimization algorithms by employing Neumann series for Hessian-inverse-vector approximation. Recent works (Khanduri et al., 2021a;b; Chen et al., 2021; Guo et al., 2021; Guo & Yang, 2021; Yang et al., 2021) then leveraged momentumbased variance reduction to further reduce the computational complexity of existing SGD-type bilevel optimizers. In this paper, we propose a stochastic ES-based method, which eliminate the computation of second-order information required by all aforementioned stochastic methods.
Please see Appendix B for more works about applications of bilevel optimization and ES methods.

heading: 2 PROPOSED ALGORITHMS
text: We describe our proposed deterministic and stochastic algorithms for bilevel optimization.

heading: 2.1 HYPERGRADIENTS
text: The key step in popular gradient-based bilevel optimizers is the estimation of the hypergradient (i.e., the gradient of the objective with respect to the outer variable x), which takes the following form:
∇Φ(x) = ∇xf(x, y∗(x)) + J∗(x)>∇yf(x, y∗(x)) (3)
where the Jacobian matrix J∗(x) = ∂y ∗(x) ∂x ∈ R d×p. Following Lorraine et al. (2020), it can be seen that∇Φ(x) contains two components: the direct gradient∇xf(x, y∗(x)) and the indirect gradient J∗(x)>∇yf(x, y∗(x)). The direct component can be efficiently computed using the existing automatic differentiation techniques. The indirect component, however, is computationally much more complex, because J∗(x) takes the form of J∗(x) = − [ ∇2yg (x, y∗(x))
]−1∇x∇yg (x, y∗(x)) (if ∇2yg (x, y∗(x)) is invertible), which contains the Hessian inverse and the second-order mixed derivative. Some approaches mitigate the issue by designing Jacobian-vector and Hessian-vector products (Pedregosa, 2016; Franceschi et al., 2017; Grazzi et al., 2020) to replace second-order computations. But these algorithms still have poor scalability to modern large-scale bilevel problems with high-dimensional variables such as neural network parameters. We next introduce the ES approach which is at the core of our idea for designing scalable Hessian-free bilevel optimizers.

heading: 2.2 EVOLUTION STRATEGIES (ES) METHOD
text: Evolution Strategies (ES) is a powerful technique to estimate the gradient of a function based on function values, when it is not feasible (such as in black-box problems) or computationally costly to evaluate the gradient. The idea of the ES method in Nesterov & Spokoiny (2017) is to approximate the gradient of a general black-box function h : Rn → R using the following oracle based only on the function values (i.e., the zeroth-order information)
∇̂h(x;u) = h(x+ µu)− h(x) µ u (4)
where u ∈ Rn is a Gaussian random vector and µ > 0 is the smoothing parameter. Such an oracle can be shown to be an unbiased estimator of the gradient of the smoothed function Eu [h(x+ µu)].

heading: 2.3 ES-BASED JACOBIAN AND ITS ENABLED BILEVEL OPTIMIZERS
text: Drawback of an existing ES-based bilevel optimizer. In Gu et al. (2021), the ES method was directly applied to estimate the hypergradient of Φ(x), treating Φ(x) as a blackbox objective. Since
the hypergradient has a complex form as in eq. (3) and can be very sensitive to both inner and outer functions, such an estimation will likely have a large bias error. As our experiments in Section 4 demonstrate, such a method is not robust and even fails to converge in complex bilevel optimization problems (such as with neural network parameters).
Our key idea is to exploit the analytical structure of the hypergradient in eq. (3), where the derivatives ∇xf(x, y∗(x)) and ∇yf(x, y∗(x)) can be computed efficiently and accurately, and then use the ES method only to estimate the JacobianJ∗(x), which is the major term posing computational difficulty. In this way, our estimation of the hypergradient can be much more accurate and reliable.
We propose the following novel ES-based Jacobian estimator, which contains two ingredients: (i) for a given x, apply an algorithm to solve the inner optimization problem and use the output as an approximation of y∗(x); for example, the output yN (x) of N gradient descent steps of the inner problem can serve as an estimate for y∗(x). Then JN (x) = ∂y N (x) ∂x serves as an estimate of J∗(x); and (ii) construct an ES-based Jacobian estimator ĴN (x;u) ∈ Rd×p for JN (x) as
ĴN (x;u) = yN (x+ µu)− yN (x)
µ u> (5)
where u ∈ Rp is a Gaussian vector with independent and identically distributed (i.i.d.) entries. Then for any vector v ∈ Rd, the Jacobian-vector product can be efficiently computed using only vector-vector dot product ĴN (x;u)>v = 〈δ(x;u), v〉u, where δ(x;u) = y N (x+µu)−yN (x) µ ∈ R d.
ESJ: A novel bilevel optimizer with ES Jacobian oracles. We design a bilevel optimizer (see Algorithm 1) using the ES-based Jacobian oracle in eq. (5), which we call as the ESJ algorithm. At each step k of the algorithm, ESJ runs an N -step full GD to approximate yNk (xk). ESJ then samples Q Gaussian vectors {uk,j ∈ N (0, I), j = 1, ..., Q}, and for each sample uk,j , runs an N -step full GD to approximate yNk (xk+µuk,j), and then computes the Jacobian estimator ĴN (x;uk,j) as in eq. (5). Then the sample average over the Q estimators is used for constructing the following hypergradient estimator for updating the outer variable x.
∇̂Φ(xk) = ∇xf(xk, yNk ) + 1Q ∑Q j=1 Ĵ TN (xk;uk,j)∇yf(xk, yNk )
= ∇xf(xk, yNk ) + 1Q ∑Q j=1 〈 δ(xk;uk,j),∇yf(xk, yNk ) 〉 uk,j . (6)
Computationally, in contrast to the existing AID and ITD based bilevel optimizers (Pedregosa (2016), Franceschi et al. (2018), Grazzi et al. (2020)) that contains the complex Hessian- and/or Jacobian-vector product computations, ESJ has only gradient computations, and hence is computationally much more efficient as demonstrated in our experiments.
ESJ-S: A stochastic bilevel optimizer with ES Jacobian oracles. For the finite-sum problem with the objective functions given in eq. (2), we design a stochastic bilevel optimizer (see Algorithm 2 in Appendix A) based on ES Jacobian oracles, which we call as ESJ-S.
Differently from Algorithm 1, which applies GD updates to find yN (xk), ESJ-S adopts N stochastic gradient descent (SGD) steps to find {Y Nk , Y Nk,1, ..., Y Nk,Q} to the inner problem, each with the outer variable set to be xk + µuk,j . Note that all SGD runs follow the same batch sampling path {S0, ...,SN−1}. The Jacobian estimator ĴN (xk;uk,j) can then be computed as in eq. (5). At the outer level, ESJ-S samples a new batch DF independently from the inner batches {S0, ...,SN−1} to evaluate the stochastic gradients ∇xF (xk, Y Nk ;DF ) and ∇yF (xk, Y Nk ;DF ). The hypergradient is then estimated as follows.
∇̂Φ(xk) = ∇xF (xk, Y Nk ;DF ) + 1Q ∑Q j=1 〈 δ(xk;uk,j),∇yF (xk, Y Nk ;DF ) 〉 uk,j . (7)

heading: 3 CONVERGENCE ANALYSIS

heading: 3.1 TECHNICAL ASSUMPTIONS
text: In this work, we focus on the following types of loss functions. Assumption 1. The inner function g(x, y) is µg-strongly convex with respect to the variable y. For the finite-sum case, the same assumption holds for∇G(x, y; ζ).
Algorithm 1 Bilevel optimizer via ES Jacobians (ESJ) 1: Input: lower- and upper-level stepsizes α, β > 0, initializations x0 ∈ Rp and y0 ∈ Rd, inner and outer
iterations numbers K and N , and number of Gaussian vectors Q. 2: for k = 0, 1, 2, ...,K do 3: Set y0k = y0, y 0 k,j = y0, j = 1, ..., Q 4: for t = 1, 2, ..., N do 5: Update ytk = y t−1 k − α∇yg(xk, y t−1 k ) 6: end for 7: for j = 1, ..., Q do 8: Generate uk,j = N (0, I) ∈ Rp 9: for t = 1, 2, ..., N do
10: Update ytk,j = y t−1 k,j − α∇yg ( xk + µuk,j , y t−1 k,j ) 11: end for 12: Compute δj = yNk,j−y N k
µ
13: end for 14: Compute ∇̂Φ(xk) = ∇xf(xk, yNk ) + 1Q ∑Q j=1 〈 δj ,∇yf(xk, yNk ) 〉 uk,j
15: Update xk+1 = xk − β∇̂Φ(xk) 16: end for
We take the following assumptions on the inner and outer loss functions g(x, y) and f(x, y), as also adopted in Ghadimi & Wang (2018); Ji et al. (2021); Yang et al. (2021). Assumption 2. Let w = (x, y) ∈ Rd+p. The gradient ∇g(w) is Lg-Lipschitz continuous, i.e., for any w1, w2 ∈ Rd+p,
∥∥∇g(w1) −∇g(w2)∥∥ ≤ Lg∥∥w1 − w2∥∥; further, the derivatives ∇2yg(w) and ∇x∇yg(w) are ρ- and τ -Lipschitz continuous, i.e,
∥∥∇2yg(w1) − ∇2yg(w2)∥∥F ≤ ρ∥∥w1 − w2∥∥ and∥∥∇x∇yg(w1) − ∇x∇yg(w2)∥∥F ≤ τ∥∥w1 − w2∥∥. For the finite-sum case, the same assumptions hold for G(w; ζ). Assumption 3. Let w = (x, y) ∈ Rd+p. The objective f(w) and its gradient ∇f(w) are M - and Lf -Lipschitz continuous, i.e., for any w1, w2 ∈ Rd+p,∣∣f(w1)− f(w2)∣∣ ≤M∥∥w1 − w2∥∥, ∥∥∇f(w1)−∇f(w2)∥∥ ≤ Lf∥∥w1 − w2∥∥. For the finite-sum case, the same assumptions hold for F (w; ξ).
For the finite-sum setting, we take the standard bounded-variance assumption on∇G(w; ζ). Assumption 4. Gradient ∇G(w; ζ) has a bounded variance, i.e., Eζ‖∇G(w; ζ)−∇g(w)‖2 ≤ σ2 for some constant σ ≥ 0.

heading: 3.2 CONVERGENCE ANALYSIS FOR ESJ
text: Different from the standard zeroth-order analysis for a blackbox function, here we develop new techniques to analyze the ES-based estimator that depend on the entire inner optimization trajectory, which is specific to bilevel optimization. We first establish the following essential proposition which characterizes the Lipshitzness property of the approximate Jacobian matrix JN (x) = ∂y N (x) ∂x .
Proposition 1. Suppose that Assumptions 1 and 2 hold. Let LJ = (
1 + Lµg )( τ µg + ρLµ2g ) , with
L = max{Lf , Lg}. Then, the Jacobian JN (x) is Lipschitz continuous with constant LJ :∥∥JN (x1)− JN (x2)∥∥F ≤ LJ ∥∥x1 − x2∥∥ ∀x1, x2 ∈ Rp. (8) We next provide an upper-bound on the hypergradient estimation error for ESJ. Proposition 2. Suppose that Assumptions 1, 2, and 3 hold. Consider the ESJ algorithm. Its hypergradient estimation error can be upper-bounded as:
E ∥∥∇̂Φ(xk)−∇Φ(xk)∥∥2 ≤ O ((1− αµg)N + pQ + µ2dp3 + µ2dp4Q ) (9)
where E[·] is conditioned on xk and taken over the Gaussian vectors {uk,j : j = 1, ..., Q}.
By using the smoothness property in Proposition 1 and the upper-bound in Proposition 2, we provide the following characterization of the convergence rate for ESJ. Theorem 1 (Convergence of ESJ). Suppose that Assumptions 1, 2, and 3 hold. Choose the innerand outer-loop stepsizes respectively as α ≤ 1L and β = 1 4LΦ , where LΦ = L + 2L 2+τM2 µg + ρLM+L3+τML µ2g + ρL 2M µ3g . Then, the iterates xk for k = 0, ...,K − 1 of ESJ in Algorithm 1 satisfy:
1 K ∑K−1 k=0 E ∥∥∇Φ(xk)∥∥2 ≤ 16Lφ(Φ(x0)−Φ∗)K + ∆1 (10) with Φ∗ = infx Φ(x) and ∆1 = O ( (1− αµg)N + pQ + µ2dp4 Q + µ 2dp3 ) .
Further, choose K, N , µ, and Q at the order of O( 1 ), O(log 1 ), O(min{ 1 p √ dp , 1 p √ 1 dp}) and
O(p ). Then ESJ achieves an -accurate stationary point with O( p 2 log 1 ) computations of the gradient∇yg(x, y) and O( 1 ) computations of the gradients∇xf(x, y) and ∇yf(x, y).
Theorem 1 characterizes the sublinear convergence of ESJ with respect to the number K of outer iterations due to the nonconvexity of the outer objective. The convergence error in ∆1 is due to the estimation error of the Jacobian matrix J∗ via our ES oracle, which captures three types of errors: (a) the approximation error between JN and J∗ via inner-loop gradient descent, which decreases exponentially w.r.t. the number N of inner iterations due to the strong convexity of the inner objective; (b) the error between our ES oracle and the Jacobian Jµ of the smoothed output EuyN (xk + µu), which decreases sublinearly w.r.t. the batch size Q, and (c) the error between the Jacobians JN and Jµ, which can be controlled by the smoothness parameter µ. Theorem 1 also indicates that ESJ requires only gradient computations to converge, and eliminates Hessian- and Jacobian-vector products required by the existing AID and ITD based bilevel optimizers (Pedregosa (2016), Franceschi et al. (2018), Grazzi et al. (2020)). Thus, ESJ is computationally much more efficient particularly for high-dimensional problems.

heading: 3.3 CONVERGENCE ANALYSIS FOR ESJ-S
text: In this section, we apply the stochastic algorithm ESJ-S to the finite-sum objective in eq. (2) and analyze its convergence rate. We first note that the Lipschitzness property established in Proposition 1 is general and can also be used here. The following proposition establishes an upper bound on the estimation error of Jacobian J∗ by JN = ∂Y N ∂x , where Y N is the output of N inner SGD updates. Proposition 3. Suppose that Assumptions 1, 2, and 4 hold. Consider the ESJ-S algorithm. Choose the inner-loop stepsize as α = 2L+µg , where L = max{Lf , Lg}. Then, we have:
E ∥∥JN − J∗∥∥2F ≤CNγ L2µ2g + λ(L+ µg) 2(1− αµg)CN−1γ (L+ µg)2(1− αµg)− (L− µg)2 + Γ 1− Cγ .
where λ, Γ, and Cγ < 1 are constants (see appendix G for their forms).
We next provide an upper-bound on the estimation error of the hypergradient by ∇̂Φ(xk) in eq. (7). Proposition 4. Suppose that Assumptions 1, 2, 3, and 4 hold. Consider the ESJ-S algorithm. Set the inner-loop stepsize as α = 2L+µg where L = max{Lf , Lg}. Then, we have:
E ∥∥∇̂Φ(xk)−∇Φ(xk)∥∥2 ≤ O((1− αµg)N + 1
S +
1 Df + p Q + µ2dp3 +
µ2dp4
Q ) where S and Df are the sizes of the inner and outer mini-batches, respectively.
Based on Propositions 1, 3, and 4, we characterize the convergence rate for ESJ-S. Theorem 2 (Convergence of ESJ-S). Suppose that Assumptions 1, 2, 3, and 4 hold. Set the innerand outer-loop stepsizes respectively as α = 2L+µg and β = 1 4LΦ
, where L = max{Lf , Lg} and the constant LΦ are defined in Theorem 1. Then, the iterates xk, k = 0, ...,K − 1 of ESJ-S satisfy:
1 K ∑K−1 k=0 E ∥∥∇Φ(xk)∥∥2 ≤ 16(Φ(x0)−Φ∗)LΦK + ∆2, (11)
where ∆2 = O (
(1− αµg)N + 1S + 1 Df + pQ + µ2dp4 Q + µ 2dp3
) .
For ESJ-S to attain an -accurate stationary point, it suffices to select K, S, and Df at the order
of O( 1 ), and N , µ, and Q respectively at the order of O(log 1 ), O(min{ 1 p √ dp , 1 p √ 1 dp}) and
O(p ), which thus amount to O( p 3 log 1 ) computations of the gradient ∇yG(x, y, ξ) and O( 1 2 ) computations of the gradients∇xF (x, y, ξ) and ∇yF (x, y, ξ).
Comparing to the convergence error ∆1 in Theorem 1 for the deterministic algorithm ESJ, Theorem 2 for the stochastic algorithm ESJ-S captures two more sublinearly decreasing error terms 1S and 1Df respectively due to the sampling of inner and outer batches to estimate the objectives.
To compare between ESJ-S and ESJ on the finite-sum objective in eq. (2), Theorem 1 suggests that ESJ attains an -accurate stationary point of eq. (2) with O(pm 2 log 1 ) computations of ∇G and O(n ) computations of∇F . As a comparison, Theorem 2 suggests that ESJ-S achieves O( p 3 log 1 ) computations of ∇G and O( 1 2 ) computations of ∇F , which outperforms ESJ in the large sample regime where the sample sizes n,m > 1 . This also reflects the advantage of SGD over GD.

heading: 4 EXPERIMENTS
text: We validate our algorithms over three bilevel problems: shallow hyper-representation (HR) with linear/2-layer net embedding model on synthetic data, deep HR with LeNet network (LeCun et al., 1998) on MNIST dataset, and few-shot meta-learning with ResNet-12 on miniImageNet dataset. We run all models using a single NVIDIA Tesla P100 GPU. Hyperparamter Optimization (HO) experiments can be found in Appendix D. All running time are in seconds.

heading: 4.1 SHALLOW HYPER-REPRESENTATION ON SYNTHETIC DATA
text: The hyper-representation (HR) problem (Franceschi et al., 2018; Grazzi et al., 2020) searches for a regression (or classification) model following a two-phased optimization process. The inner-level identifies the optimal linear regressor parameters w, and the outer level solves for the optimal embedding model (i.e., representation) parameters λ. Mathematically, the problem can be modeled by the following bilevel optimization:
min λ∈Rp
f(λ) = 1 2n1 ‖T (X1;λ)w∗(λ)− Y1‖2 s.t. w∗(λ) = argmin
w∈Rd
1
2n2 ‖T (X2;λ)w − Y2‖2 +
γ 2 ‖w‖2 (12)
where X2 ∈ Rn2×m and X1 ∈ Rn1×m are matrices of synthesized training and validation data, and Y2 ∈ Rn2 , Y1 ∈ Rn1 are the corresponding response vectors. In the case of shallow HR, the embedding function T (·;λ) is either a linear transformation or a two-layer network. We generate data matrices X1, X2 and labels Y1, Y1 following the same process in Grazzi et al. (2020).
We compare our ESJ algorithm with the baseline bilevel optimizers AID-FP, AID-CG, ITD-R, and HOZOG (see Appendix C.1 for details about the baseline algorithms and hyperparameters used).
Figure 1 show the performance comparison among the algorithms under linear and two-layer net embedding models. It can be observed that for both cases, our proposed method ESJ converges faster than all the other approaches, and the advantage of ESJ becomes more significant in Figure 1 (d), which is under a higher-dimensional model of a two-layer net. In particular, ESJ outperforms the existing ES-based algorithm HOZOG. This is because HOZOG uses the ES technique to approximate the entire hypergradient, which likely incurs a large estimation error. In contrast, our ESJ exploits the structure of the hypergradient and uses ES only to estimate the response Jacobian so that the estimation of hypergradient is more accurate. Such an advantage is more evident under a twolayer net model, where HOZOG does not converge as shown in Figure 1 (d). This can be explained by the flat hypergradient norm as shown in Figure 1 (e), which indicates that the hypergradient estimator in HOZOG fails to provide a good descent direction for the outer optimizer. Figure 1 (c) and (f) further show that the convergence of ESJ does not change substantially with the number N of inner GD steps, and hence tuning of N in practice is not costly.

heading: 4.2 DEEP HYPER-REPRESENTATION ON MNIST DATASET
text: In order to demonstrate the advantage of our proposed algorithms in large neural net models, we perform deep hyper-representation to classify MNIST images by learning an entire LeNet network. The corresponding bilevel problem is given by
min λ Lout(λ) := 1|Dout| ∑ (xi,yi)∈Dout L (w∗(λ)f(xi;λ), yi)
s.t. w∗(λ) = arg min w∈Rc×p
Lin(w, λ) := 1|Din| ∑
(xi,yi)∈Din (L(wf(xi;λ), yi) + β2 ‖w‖
2). (13)
where f(xi;λ) ∈ Rp corresponds to features extracted from data point xi, L(·, ·) is the crossentropy loss function, c = 10 is the number of categories, and Din and Dout are data used to compute respectively inner and outer loss functions. Since the sizes of Dout and Din are large in the case of MNIST dataset, we apply the more efficient stochastic algorithm ESJ-S in Algorithm 2 with a minibatch size B = 256 to estimate the inner and outer losses Lin and Lout.
Figure 2 compares the classification accuracy on the outer datasetDout among the different methods. Our algorithm ESJ-S converges with the fastest rate and attains the best accuracy with the lowest variance among all algorithms. Note that ESJ-S is able to attain the same accuracy of 0.98+ obtained by the standard training of all parameters with one-phased optimization on the MNIST dataset using the same backbone network. All other methods fail to recover such a level of performance, and instead saturate around an accuracy of 0.93. Further, Figure 2(c) indicates that the convergence of ESJ-S does not change substantially with the number N of inner SGD steps. This demonstrates the robustness of our method when applied to complex function geometries such as deep nets.

heading: 4.3 FEW-SHOT META-LEARNING OVER MINIIMAGENET
text: To study our algorithms over larger neural nets, we study the few-shot image recognition problem, where classification tasks Ti, i = 1, ...,m are sampled over a distribution PT . In particular, we consider the ANIL meta-learning method (Raghu et al., 2019; Ji et al., 2020a), where all tasks share common embedding features parameterized by φ, and each task Ti has its task-specific parameter wi for i = 1, ...,m. More specifically, we set φ to be the parameters of the convolutional part of a deep CNN model (e.g., ResNet-12 network) and w includes the parameters of the last classification layer. All model parameters (φ, w) are trained following a bilevel procedure. In the inner-loop, the base learner of each task Ti fixes φ and minimizes its loss function over a training set Si to
obtain its adapted parameters w∗i . At the outer stage, the meta-learner computes the test loss for each task Ti using the parameters (φ, w∗i ) over a test set Di, and optimizes the parameters φ of the common embedding function by minimizing the meta-objective Lmeta over all classification tasks. The problem can be expressed as the following bilevel optimization
min φ Lmeta(φ, w̃∗) := 1m
∑m i=1 LDi (φ,w∗i )
s.t. w̃∗ = arg min w̃
Ladapt(φ, w̃) := 1m ∑m i=1 LSi (φ,wi) ,
(14)
where we collect all task-specific parameters into w̃ = (w1, ..., wm) and the corresponding minimizers into w̃∗ = (w∗1 , ..., w ∗ m). The functions LSi (φ,wi) = 1|Si| ∑ ζ∈Si(L(φ,wi; ζ) + R(wi))
and LDi (φ,w∗i ) = 1|Di| ∑ ξ∈Di L (φ,w ∗ i ; ξ) correspond respectively to the training and test loss functions for task Ti, with R a strongly-convex regularizer and L a classification loss function. In our setting, since the task-specific parameters correspond to the weights of the last linear layer, the inner-level objective Ladapt(φ, w̃) is strongly convex with respect to w̃ = (w1, ..., wm). We note that the problem studied in Section 4.2 can be seen as single-task instances of the more general multi-task learning problem in eq. (14). However, in contrast to the problem in Section 4.2, the sizes of the datasets Di and Si are usually small in few-shot learning and full GD can be applied here. Hence, we use ESJ (Algorithm 1) here. Also since the number m of tasks in few-shot classification datasets is often very large, it is preferable to sample a minibatch of i.i.d. tasks by PT at each meta (i.e., outer) iteration and update the meta parameters based on these tasks.
We conduct few-shot meta-learning on the miniImageNet dataset (Vinyals et al., 2016) using two different backbone networks for feature extraction: ResNet-12 and CNN4 (Vinyals et al., 2016). The dataset description and hyperparameter details can be found in Appendix C.4. We compare our algorithm ESJ with four baseline methods for few-shot meta-learning MAML (Finn et al., 2017), ANIL (Raghu et al., 2019), MetaOptNet (Lee et al., 2019), and ProtoNet (Snell et al., 2017). We run their efficient Pytorch Lightning implementations available at the learn2learn repository (Arnold et al., 2019).
Figure 3(a) and (b) show that our algorithm ESJ converges faster than the other baseline methods. Also, Comparing Figure 3(a) and (b), the advantage of our method over the baselines MAML and ANIL becomes more significant as the size of the network increases. Further, Figure 3(c) shows that MetaOptNet did not reach 69% accuracy after 20 hours of training with ResNet-12 network. In comparison, our ESJ is able to attain 69% in less than three hours, which is about 1.5 times less than the time taken for ProtoNet to reach the same performance level. Both ESJ and ProtoNet saturate around 70% accuracy after 10 hours of training.

heading: 5 CONCLUSION
text: In this paper, we propose a novel ES-based approach for bilevel optimization, which eliminates all second-order computations in the gradient-based approaches. Compared to the existing ES-based algorithm, our approach explores the analytical structure of the hypergradient, and hence leads to much more efficient and accurate hypergradient estimation. Thus, our algorithm outperforms the existing algorithms in the experiments, particularly in the high-dimensional applications. We also characterize the convergence rate for our proposed algorithms and show that the polynomial-time complexity can be achieved. We anticipate that our approach will be useful for accelerating bilevel algorithms in various machine learning problems.

heading: 6 REPRODUCIBILITY CHECKLIST
text: To ensure reproducibility, we use the Machine Learning Reproducibility Checklist v2.0, Apr. 7 2020 (Pineau et al., 2021). An earlier version of this checklist (v1.2) was used for NeurIPS 2019 (Pineau et al., 2021).
• For all models and algorithms presented, – A clear description of the mathematical settings, algorithm, and/or model. We
clearly describe all of the settings, formulations, and algorithms in Section 2. – A clear explanation of any assumptions. All assumptions are stated in Section 3.1
and details are clearly explained in Section 3.1. – An analysis of the complexity (time, space, sample size) of any algorithm. We
provide the time/computational complexity analysis for our algorithms in Theorems 1 and 2.
• For any theoretical claim, – A clear statement of the claim. A clear statement of theoretical claims are made in
Section 3. – A complete proof of the claim. The complete proofs of all claims are available in
Appendix E, Appendix F, and Appendix G. • For all datasets used, check if you include:
– The relevant statistics, such as number of examples. We use widely adopted datasets MNIST and miniImageNet in Section 4. The related statistics can be seen at http://yann.lecun.com/exdb/mnist/.
– The details of train/validation/test splits. We give this information in the Supplementary Appendix C.
– An explanation of any data that were excluded, and all pre-processing step. We did not exclude any data or perform any pre-processing.
– For new data collected,a complete description of the data collection process, such as instructions to annotators and methods for quality control. We do not collect or release new datasets.
• For all shared code related to this work, check if you include: – Training code. The training code is available in our code in supplementary material. – Evaluation code. The evaluation code is available in our code in supplementary ma-
terial. – (Pre-)trained model(s). We do not release any pre-trained models.
• For all reported experimental results, check if you include: – The range of hyper-parameters considered, method to select the best hyper-
parameter configuration, and specification of all hyper-parameters used to generate results. We provide all details of the hyper-parameter tuning in Supplementary Appendix C.
– A clear definition of the specific measure or statistics used to report results. We use the classification accuracy on test-set and the loss on the train-set.
– A description of results with central tendency (e.g. mean) & variation (e.g. error bars). We do not report the mean and standard deviation for experiments.
– The average runtime for each result, or estimated energy cost. We report the running time of the algorithms in Section 4.
– A description of the computing infrastructure used. All detailed descriptions are presented in Section 4.

heading: Supplementary Material

heading: A STOCHASTIC BILEVEL OPTIMIZER BASED ON ES METHOD
text: In the following, we present the algorithm specification for our proposed stochastic bilevel optimizer based on ES method, which we call ESJ-S.
Algorithm 2 Stochastic bilevel optimizer via ES Jacobians (ESJ-S) 1: Input: lower- and upper-level stepsizes α, β > 0, initializations x0 ∈ Rp and y0 ∈ Rd, inner
and outer iterations numbers K and N , and number of Gaussian vectors Q. 2: for k = 0, 1, ...,K do 3: Set Y 0k = y0, Y 0 k,j = y0, j = 1, ..., Q 4: Generate uk,j = N (0, I) ∈ Rp, j = 1, ..., Q 5: for t = 1, 2, ..., N do 6: Draw a sample batch St−1 7: Update Y tk = Y t−1 k − α∇yG(xk, Y t−1 k ;St−1) 8: for j = 1, ..., Q do 9: Update Y tk,j = Y t−1 k,j − α∇yG ( xk + µuk,j , Y t−1 k,j ;St−1
) 10: end for 11: end for 12: Compute δj = Y Nk,j−Y N k
µ , j = 1, ..., Q
13: Draw a sample batch DF 14: Compute ∇̂Φ(xk) = ∇xF (xk, Y Nk ;DF ) + 1Q ∑Q j=1 〈 δj ,∇yF (xk, Y Nk ;DF ) 〉 uk,j
15: Update xk+1 = xk − β∇̂Φ(xk) 16: end for

heading: B ADDITIONAL RELATED WORK
text: In this section, we provide further related work on the applications of bilevel optimization and ES methods.
Bilevel optimization applications. Bilevel optimization has been employed in various applications such as few-shot meta-learning (Snell et al., 2017; Franceschi et al., 2018; Rajeswaran et al., 2019; Zügner & Günnemann, 2019; Ji et al., 2020a;b), hyperparameter optimization (Franceschi et al., 2017; Mackay et al., 2018; Shaban et al., 2019), neural architecture search (Liu et al., 2018a; Zhang et al., 2021), etc. This paper demonstrates the superior performance of the proposed ES-based bilevel optimizer in meta-learning and hyperparameter optimization.
ES applications. Evolution Strategies (ES) method has been studied for more than two decades (see review papers (Hansen, 2006; Varelas et al., 2018)). In particular, Nesterov & Spokoiny (2017) proposed a simple and effective ES-gradient via Gaussian smoothing, which was further extended to the stochastic setting by Ghadimi & Lan (2013). Such an ES technique has exhibited great effectiveness in various applications including meta-reinforcement learning (Song et al., 2019), hyperparameter optimization (Gu et al., 2021), adversarial machine learning (Ji et al., 2019; Liu et al., 2018b), minimax optimization (Liu et al., 2020b; Xu et al., 2020), etc. This paper proposes a novel ES-based Jacobian estimator for accelerating bilevel optimization. The ES method has also been used for studying hyperparameter optimization problems. For example, Mackay et al. (2018) models the response function itself as a neural network (where each layer involves an affine transformation of hyperparameters) using the Self-Tuning Networks (STNs). An improved and more stable version of STNs was further proposed in Bae & Grosse (2020), which focused on accurately approximating the response Jacobian rather than the response function itself.

heading: C FURTHER SPECIFICATIONS FOR EXPERIMENTS IN SECTION 4
text: We note that the smoothing parameter µ (in Algorithms 1 and 2) was easy to set and a value of 0.1 or 0.01 yields a good starting point across all our experiments. The batch size Q (in Algorithms 1 and 2) is fixed to 1 (i.e., we use one Jacobian oracle) in all our experiments.

heading: C.1 SPECIFICATIONS ON BASELINE BILEVEL APPROACHES IN SECTION 4.1
text: We compare our algorithm ESJ with the following baseline methods:
• HOZOG (Gu et al., 2021): a hyperparameter optimization algorithm that uses evolution strategies to estimate the entire hypergradient (both the direct and indirect component). We use our own implementation for this method.
• AID-CG (Grazzi et al., 2020; Rajeswaran et al., 2019): approximate implicit differentiation with conjugate gradient. We use its implementation provided at https://github.com/ prolearner/hypertorch
• AID-FP (Grazzi et al., 2020): approximate implicit differentiation with fixed-point. We experimented with its implementation at the repositery https://github.com/prolearner/ hypertorch
• ITD-R (REVERSE) (Franceschi et al., 2017): an iterative differentiation method that computes hypergradients using reverse mode automatic differention (RMAD). We use its implementation provided at https://github.com/prolearner/hypertorch.

heading: C.2 HYPERPARAMETERS DETAILS FOR SHALLOW HR EXPERIMENTS IN SECTION 4.1
text: For the linear embedding case, we set the smoothing parameter µ to be 0.01 for ESJ and HOZOG. We use the following hyperparameters for all compared methods. The number of inner GD steps is fixed to N = 20 with the learning rate of α = 0.001. For the outer optimizer, we use Adam (Kingma & Ba, 2014) with a learning rate of 0.05. The value of γ in eq. (12) is set to be 0.1. For the two-layer net case, we use µ = 0.1 for ESJ and HOZOG. For all methods, we set N = 10, α = 0.001, β = 0.001, and use Adam with a learning rate of 0.01 as the outer optimizer.

heading: C.3 SPECIFICATIONS ON BASELINE STOCHASTIC ALGORITHMS IN SECTION 4.2
text: We compare our stochastic algorithm ESJ-S with the following baseline stochastic bilevel algorithms.
• stocBiO: an approximate implicit differentiation method that uses Neumann Series to estimate the Hessian inverse. We use its implementation available at https://github.com/ JunjieYang97/StocBio.
• AID-CG-S and AID-FP-S: stochastic versions of AID-CG and AID-FP, respectively. We use their implementations in the repository https://github.com/prolearner/ hypertorch.

heading: C.4 SPECIFICATIONS FOR FEW-SHOT META-LEARNING IN SECTION 4.3
text: The miniImageNet dataset (Vinyals et al., 2016) is a large-scale benchmark for few-shot learning generated from ImageNet (Russakovsky et al., 2015) Russakovsky. The dataset consists of 100 classes with each class containing 600 images of size 84 × 84. Following (Arnold et al., 2019), we split the classes into 64 classes for meta-training, 16 classes for meta-validation, and 20 classes for meta-testing. More specfically, we use 20000 tasks for meta-training, 600 tasks for meta-validation, and 600 tasks for meta-testing. We normalize all image pixels by their means and standard deviations over RGB channels and do not perform any additional data augmentation. At each meta-iteration, we sample a batch of 16 training tasks and update the parameters based on these tasks. We set the smoothness parameter to be µ = 0.1 and use N = 30 inner steps. We use SGD with a learning rate of α = 0.01 as inner optimizer and Adam with a learning rate of β = 0.01 as outer (meta) optimizer.

heading: D EXPERIMENTS ON HYPERPARAMETER OPTIMIZATION
text: Hyperparameter optimization (HO) is the problem of finding the set of the best hyperparamters (either representational or regularization parameters) that yield the optimal value of some criterion of model quality (usually a validation loss on unseen data). HO can be posed as a bilevel optimization problem in which the inner problem corresponds to finding the model parameters by minimizing a training loss (usually regularized) for the given hyperparameters and then the outer problem minimizes over the hyperparameters. Hence, HO can be mathematically expressed as follows
min λ Lval(λ) := 1|Dval| ∑ ξ∈Dval L (w ∗(λ); ξ)
s.t. w∗(λ) = arg min w
Ltr(w, λ) := 1|Dtr| ∑ ζ∈Dtr(L(w, λ; ζ) +R(w, λ)),
(15)
where L is a loss function (e.g., logistic loss),R(w, λ) is a regularizer, andDtr andDval are respectively training and validation data. Note that the loss function used to identify hyperparameters must be different from the one used to find model parameters; otherwise models with higher complexities would be always favored. This is usually achieved in HO by using different data splits (here Dval and Dtr) to compute validation and training losses, and by adding a regularizer term on the training loss.
Following (Franceschi et al., 2017; Grazzi et al., 2020), we perform classification on the 20 Newsgroup dataset, where the classifier is modeled by an affine transformation, the cost function L is the cross-entrpy loss, and R(w, λ) is a strongly-convex regularizer. We set one l2-regularization hyperparameter for each weight in w, so that λ and w have the same size.
For ESJ and HOZOG, we use GD with a learning rate of 100 and a momentum of 0.9 to perform the inner updates. The outer learning rate is set to be 0.02. We set the smoothing parameter (µ in Algorithm 1) to be 0.01. For AID-FP, AID-CG, and REVERSE we use the suggested hyperparameters in their implementations accompanying the paper Grazzi et al. (2020).
It can be seen from Figure 4, our method ESJ slightly outperforms HOZOG and converges faster than the other AID and ITD based approaches. We note that the similar performance for ESJ and HOZOG can be explained by the fact that in HO, the hypergradient expession in eq. (6) contains only the second term (the first term is zero), which is very close to the approximation in HOZOG method. However, as we have seen in the other experiments (not for HO), ESJ is a more robust and stable bilevel optimizer than HOZOG, and it achieves good performance across many bilevel problems.

heading: E SUPPORTING TECHNICAL LEMMAS
text: In this section, we provide auxiliary lemmas that are used for proving the convergence results for the algorithms ESJ and ESJ-S.
In the following proofs, we let L = max{Lg, Lf} and D be such that ‖y∗(x)‖ ≤ D.
First we recall that for any two matrices A ∈ Rm×r and B ∈ Rr×n, we have the following upperbound on the Frobenius norm of their product,∥∥AB∥∥ F ≤ ∥∥A∥∥∥∥B∥∥ F . (16)
The following lemma follows directly from the Lipschitz properties in Assumptions 2 and 3. Lemma 1. Suppose that Assumptions 2 and 3 hold. Then, the stochastic derivatives ∇F (x, y; ξ), ∇x∇yG(x, y; ξ), and ∇2yG(x, y; ξ) have bounded variances, i.e., for any (x, y) and ξ we have:
• Eξ ∥∥∇F (x, y; ξ)−∇f(x, y)∥∥2 ≤M2;
• Eξ ∥∥∇x∇yG(x, y; ξ)−∇x∇yg(x, y)∥∥2F ≤ L2;
• Eξ ∥∥∇2yG(x, y; ξ)−∇2yg(x, y)∥∥2F ≤ L2.
Using Lemma 2.2 in Ghadimi & Wang (2018), the following lemma characterizes the Lipschitz property of the gradient of the total objective Φ(x) = f(x, y∗(x)). Lemma 2. Suppose that Assumptions 1, 2, and 3 hold. Then, we have:∥∥∇Φ(x2)−∇Φ(x1)∥∥ ≤ LΦ∥∥x2 − x1∥∥ ∀x1 ∈ Rp, x2 ∈ Rp, where the constant LΦ = L+ 2L 2+τM2
µg + ρLM+L 3+τML µ2g + ρL 2M µ3g .
We next provide some essential properties of the ES-based (i.e., zeroth-order) gradient oracle in eq. (4), due to Nesterov & Spokoiny (2017). Lemma 3. Let h : Rn → R be a differentiable function with L-Lipsctitz gradient. Define its Gaussian smooth approximation hµ(x) = Eu [h(x+ µu)], where µ > 0 and u ∈ Rn is a standard Gaussian random vector. Then, hµ is differentiable and we have:
• The gradient of hµ takes the form
∇hµ(x) = Eu h(x+ µu)− h(x)
µ u.
• For any x ∈ Rn, ∥∥∇hµ(x)−∇h(x)∥∥ ≤ µ2 2 L(n+ 3)3/2.
• For any x ∈ Rn, Eu ∥∥∥h(x+ µu)− h(x) µ u ∥∥∥2 ≤ 4(n+ 4)∥∥∇hµ(x)∥∥2 + 3 2 µ2L2(n+ 5)3.
Note the first item in Lemma 3 implies that the oracle in eq. (4) is in indeed an unbiased estimator of the gradient of the smoothed function hµ. Lemma 4. Suppose that Assumptions 1 and 2 hold. The Jacobian J∗ = ∂y ∗(x) ∂x has bounded norm:∥∥J∗∥∥F ≤ Lµg . (17)
Proof of Lemma 4. From the first order optimality condition of y∗(x), we have∇yg (x, y∗(x)) = 0. Hence, the Implicit Function Theorem implies:
J∗ = − [ ∇2yg (x, y∗(x)) ]−1∇x∇yg (x, y∗(x)) . (18) Taking norms and applying eq. (16) together with Assumptions 1 and 2 yield the desired result∥∥J∗∥∥F ≤ ∥∥∇x∇yg (x, y∗(x)) ∥∥F∥∥ [∇2yg (x, y∗(x))]−1 ∥∥ ≤ Lµg . (19)
Lemma 5. Suppose that Assumptions 1 and 2 hold. The Jacobian JN = ∂y N k
∂xk has bounded norm:∥∥JN∥∥F ≤ Lµg . (20)
Proof of Lemma 5. The inner loop gradient descent updates writes
ytk = y t−1 k − α∇yg ( xk, y t−1 k ) , t = 1, . . . , N.
Taking derivatives w.r.t. xk yields Jt = Jt−1 − α∇x∇yg ( xk, y t−1 k ) − αJt−1∇2yg ( xk, y t−1 k ) = Jt−1 ( I − α∇2yg ( xk, y t−1 k )) − α∇x∇yg ( xk, y t−1 k ) .
Telescoping over t from 1 to N yields
JN = J0 N−1∏ t=0 ( I − α∇2yg ( xk, y t k )) − α N−1∑ t=0 ∇x∇yg ( xk, y t k ) N−1∏ m=t+1 ( I − α∇2yg (xk, ymk ) ) = −α
N−1∑ t=0 ∇x∇yg ( xk, y t k ) N−1∏ m=t+1 ( I − α∇2yg (xk, ymk ) ) . (21)
Hence, we have
∥∥JN∥∥F ≤ αN−1∑ t=0 ∥∥∇x∇yg (xk, ytk) ∥∥F∥∥ N−1∏ m=t+1 ( I − α∇2yg (xk, ymk ) ) ∥∥ (i)
≤ α N−1∑ t=0 L N−1∏ m=t+1 ∥∥I − α∇2yg (xk, ymk )∥∥ (ii)
≤ αL N−1∑ t=0 (1− αµg)N−1−t
= αL N−1∑ t=0 (1− αµg)t ≤ L µg
where (i) follows from Assumption 2 and (ii) applies the strong-convexity of function g(x, ·). This completes the proof.
Lemma 6. Suppose that Assumptions 1 and 2 hold. Then, the Jacobian in the stochastic algorithm ESJ-S JN = ∂Y N k
∂xk has bounded norm, as shown below.∥∥JN∥∥F ≤ Lµg . (22)
Proof of Lemma 6. The proof follows similarly to Lemma 5.

heading: E.1 PROOF OF PROPOSITION 1
text: Proposition 1 (Re-stated). Suppose that Assumptions 1 and 2 hold. Define the constant
LJ =
( 1 + L
µg
)( τ
µg + ρL µ2g
) . (23)
Then, the Jacobian JN (x) = ∂y N (x) ∂x is LJ -Lipschitz with respect to x under the Frobenious norm:∥∥JN (x1)− JN (x2)∥∥F ≤ LJ ∥∥x1 − x2∥∥ ∀x1 ∈ Rp, x2 ∈ Rp. (24)
Proof of Proposition 1. Using eq. (21), we have
JN (x) = −α N−1∑ t=0 ∇x∇yg ( x, yt(x) ) N−1∏ m=t+1 ( I − α∇2yg (x, ym(x)) ) , x ∈ Rp
Hence, for x1 ∈ Rp and x2 ∈ Rp, we have∥∥JN (x1)− JN (x2)∥∥F =α ∥∥∥N−1∑ t=0 ∇x∇yg ( x1, y t(x1) ) N−1∏ m=t+1 ( I − α∇2yg (x1, ym(x1))
) − N−1∑ t=0 ∇x∇yg ( x2, y t(x2) ) N−1∏ m=t+1 ( I − α∇2yg (x2, ym(x2)) ) ∥∥∥ F
≤α N−1∑ t=0 ∥∥∇x∇yg (x1, yt(x1)) N−1∏ m=t+1 ( I − α∇2yg (x1, ym(x1)) ) −∇x∇yg ( x2, y
t(x2) ) N−1∏ m=t+1 ( I − α∇2yg (x2, ym(x2)) ) ∥∥ F
(i) ≤α N−1∑ t=0 ∥∥∇x∇yg (x1, yt(x1)) ∥∥F∥∥At(x1)−At(x2)∥∥ + α
N−1∑ t=0 ∥∥At(x2)∥∥∥∥∇x∇yg (x1, yt(x1))−∇x∇yg (x2, yt(x2)) ∥∥F (25) where we define At(x) = ∏N−1 m=t+1 ( I − α∇2yg (x, ym(x)) ) and (i) follows from eq. (16).
Next we upper bound the quantity ∥∥At(x1)−At(x2)∥∥, as shown below.∥∥At(x1)−At(x2)∥∥ ≤∥∥α (∇2yg (x2, yt+1(x2))−∇2yg (x2, yt+1(x2)))At+1(x1)
+ ( I − α∇2yg ( x2, y t+1(x2) )) (At+1(x1)−At+1(x2)) ∥∥
≤α ∥∥At+1(x1)∥∥∥∥∇2yg (x1, yt+1(x1))−∇2yg (x2, yt+1(x2)) ∥∥
+ ∥∥I − α∇2yg (x2, yt+1(x2)) ∥∥∥∥At+1(x1)−At+1(x2)∥∥
≤(1− αµg) ∥∥At+1(x1)−At+1(x2)∥∥
+ αρ(1 + L
µg )(1− αµg)N−t−2 ∥∥x1 − x2∥∥, (26) where the last inequality follows from Lemma 5 and Assumptions 1 and 2.
Telescoping eq. (26) over t yields∥∥At(x1)−At(x2)∥∥ ≤(1− αµg)N−t−2∥∥AN−2(x1)−AN−2(x2)∥∥ +
N−t−3∑ m=0 αρ(1 + L µg )(1− αµg)N−t−2−m(1− αµg)m ∥∥x1 − x2∥∥ =(1− αµg)N−t−2
∥∥∇2yg (x1, yN−1(x1))−∇2yg (x2, yN−1(x2)) ∥∥ +
N−t−3∑ m=0 αρ(1 + L µg )(1− αµg)N−t−2 ∥∥x1 − x2∥∥ (i) ≤αρ(1 + L µg )(1− αµg)N−t−2 ∥∥x1 − x2∥∥
+ (N − t− 2)αρ(1 + L µg
)(1− αµg)N−t−2 ∥∥x1 − x2∥∥
=αρ(1 + L
µg )(N − t− 1)(1− αµg)N−t−2 ∥∥x1 − x2∥∥, (27) where (i) follows from Lemma 5 and Assumption 2. Replacing eq. (27) in eq. (25) and using Assumption 2, we have∥∥JN (x1)− JN (x2)∥∥F
≤ α N−1∑ t=0 Lαρ(1 + L µg )(N − t− 1)(1− αµg)N−t−2 ∥∥x1 − x2∥∥ + α
N−1∑ t=0 τ(1 + L µg )(1− αµg)N−t−1 ∥∥x1 − x2∥∥ ≤ α2Lρ(1 + L
µg ) ∥∥x1 − x2∥∥N−1∑
m=0
m(1− αµg)m−1 + τ
µg (1 +
L µg ) ∥∥x1 − x2∥∥
≤ ρL µ2g (1 + L µg ) ∥∥x1 − x2∥∥+ τ µg (1 + L µg ) ∥∥x1 − x2∥∥ (28)
where we use ∑N−1 m=0 mx m−1 ≤ 1α2µ2g in eq. (28), which can be obtained by taking derivatives for
the expression ∑N−1 m=0 x
m with respect to x. Hence, rearranging and using the definition of LJ in eq. (23) finishes the proof.
Lemma 7. Suppose that Assumptions 1 and 2 hold. Define the constant
LJ =
( 1 + L
µg
)( τ
µg + ρL µ2g
) .
Then, the Jacobian JN (x) = ∂Y N (x;·) ∂x is LJ -Lipschitz with respect to x under the Frobenius norm:∥∥JN (x1; ·)− JN (x2; ·)∥∥F ≤ LJ ∥∥x1 − x2∥∥ ∀x1 ∈ Rp, x2 ∈ Rp.
Proof of Lemma 7. The proof follows similarly to that for Proposition 1.

heading: F PROOFS FOR DETERMINISTIC BILEVEL OPTIMIZATION
text: For notation convenience, we define the following quantities:
ĴN,j = ĴN (xk, uj) =  yN1 (xk+µuj)−y N 1 (xk) µ u > j
... yNd (xk+µuj)−y N d (xk)
µ u > j
 , JN = ∂yNk∂xk , J∗ = ∂y ∗ k ∂xk , (29)
where uj ∈ Rp, j = 1, . . . , Q are standard Gaussian vectors. Let yNi,µ(xk) be the Gaussian smooth approximation of yNi (xk). We collect y N i,µ(xk) for i = 1, . . . , d together as a vector y N µ (xk), which is the Gaussian approximation of the vector yN (xk). If µ > 0, yNµ (xk) is differentiable and we let Jµ be the Jocobian given by
Jµ = ∂yNµ (xk)
∂xk . (30)
We approximate ∂y N k
∂xk using the average ES estimator given by ĴN = 1Q ∑Q j=1 ĴN,j . The hypergra-
dient is then approximated as
∇̂Φ(xk) = ∇xf(xk, yNk ) + Ĵ>N∇yf(xk, yNk )
= ∇xf(xk, yNk ) + 1
Q Q∑ j=1 Ĵ>N,j∇yf(xk, yNk ). (31)
Let δj = yN (xk+µuj)−yN (xk)
µ , and let δi,j be the i-th component of δj . Hence, we have
ĴN,j =  δ1,ju > j δ2,ju > j ...
δd,ju > j
 .
Ĵ>N,j∇yf(xk, yNk ) = ( δ1,juj δ2,juj . . . δd,juj ) ∇yf(xk, yNk )
= 〈 δj ,∇yf(xk, yNk ) 〉 uj . (32)
Using eq. (31) and eq. (32), the estimator for the hypergradient can thus be computed as
∇̂Φ(xk) = ∇xf(xk, yNk ) + 1
Q Q∑ j=1 〈 δj ,∇yf(xk, yNk ) 〉 uj .

heading: F.1 PROOF OF PROPOSITION 2
text: Proposition 2 (Re-stated). Suppose that Assumptions 1, 2, and 3 hold. Then, the expected estimation error can be upper-bounded as follows: E ∥∥∇̂Φ(xk)−∇Φ(xk)∥∥2 ≤2L2D2(1− αµg)N + 4L4
µ2g D2(1− αµg)N + 24(4p+ 15)
L2M2
Qµ2g
+ µ2
Q L2JM
2dP4(p) + 24L2M2(1− αµg)2N
µ2g + 6µ2L2JM 2d(p+ 3)3
+ 48M2(τµg + Lρ) 2
µ4g (1− αµg)N−1D2, (33)
where the expectation E[·] is conditioned on xk and yNk .
Proof of Proposition 2. Based on the definitions of ∇Φ(xk) and ∇̂Φ(xk) and conditioning on xk and yNk , we have
E ∥∥∇̂Φ(xk)−∇Φ(xk)∥∥2 ≤2 ∥∥∇xf(xk, yNk )−∇xf(xk, y∗k)∥∥2 + 2E∥∥Ĵ>N∇yf(xk, yNk )− J>∗ ∇yf(xk, y∗k)∥∥2
≤2L2 ∥∥yNk − y∗k∥∥2 + 4∥∥J∗∥∥2F∥∥∇yf(xk, yNk )−∇yf(xk, y∗k)∥∥2
+ 4E ∥∥ĴN − J∗∥∥2F∥∥∇yf(xk, yNk )∥∥2
(i) ≤2L2D2(1− αµg)N + 4 L4
µ2g ∥∥yNk − y∗k∥∥2 + 4M2E∥∥ĴN − J∗∥∥2F (ii)
≤ 2L2D2(1− αµg)N + 4 L4
µ2g D2(1− αµg)N + 4M2E ∥∥ĴN − J∗∥∥2F (34) where (i) follows from Lemma 4 and Assumption 3, and (i) and (ii) also use the following result for full GD (when applied to a strongly-convex function).∥∥yNk − y∗k∥∥2 ≤ (1− αµg)ND2. Next, we upper-bound the last term E
∥∥ĴN − J∗∥∥2F at the last line of eq. (34). First note that E ∥∥ĴN − J∗∥∥2F ≤ 3E∥∥ĴN − Jµ∥∥2F + 3∥∥JN − J∗∥∥2F + 3∥∥Jµ − JN∥∥2F . (35)
We then upper-bound each term of the right hand side of eq. (35). For the first term, we have
E ∥∥ĴN − Jµ∥∥2F =E∥∥ 1Q Q∑ j=1 ĴN,j − Jµ ∥∥2 F
= 1 Q2 E ∥∥ Q∑ j=1 ( ĴN,j − Jµ )∥∥2 F
= 1
Q2 E  Q∑ j=1 ∥∥ĴN,j − Jµ∥∥2F + 2∑ i<j 〈 ĴN,i − Jµ, ĴN,j − Jµ 〉 = 1
Q2 Q∑ j=1 E ∥∥ĴN,j − Jµ∥∥2F
= 1 Q E ∥∥ĴN,j − Jµ∥∥2F , j ∈ {1, . . . , Q}. (36)
We next upper-bound the term E ∥∥ĴN,j − Jµ∥∥2F in eq. (36).
E ∥∥ĴN,j − Jµ∥∥2F =E∥∥ĴN,j∥∥2F − ∥∥Jµ∥∥2F (37)
(i) ≤ d∑ i=1 ( 4(p+ 4) ∥∥∇yNi,µ∥∥2 + 32µ2L2J (p+ 5)3 ) − d∑ i=1 ∥∥∇yNi,µ∥∥2 ≤
d∑ i=1 ( (4p+ 15) ∥∥∇yNi,µ∥∥2 + 32µ2L2J (p+ 5)3 ) , (38)
where (i) follows by applying Lemma 3 to the components of vector yN (xk) which have Lipschitz gradients by Proposition 1. Then, noting that ∥∥∇yNi,µ∥∥2 ≤ 2∥∥∇yNi ∥∥2+ 12µ2L2J (p+3)3 and replacing in eq. (38), we have
E ∥∥ĴN,j − Jµ∥∥2F ≤ d∑
i=1
( 2(4p+ 15) ∥∥∇yNi ∥∥2 + µ2L2JP4(p)) ≤2(4p+ 15)
∥∥JN∥∥2F + µ2L2J dP4(p) (i)
≤2(4p+ 15)L 2
µ2g + µ2L2J dP4(p), (39)
where (i) follows from Lemma 5 and P4 is a polynomial of degree 4 in p. Combining eq. (36) and eq. (38) yields
E ∥∥ĴN − Jµ∥∥2F ≤2(4p+ 15) L2Qµ2g + µ 2 Q L2J dP4(p). (40)
We next upper-bound the second term at the right hand side of eq. (35), which can be upper-bounded using eq. (41) in Ji et al. (2021), as shown below.∥∥JN − J∗∥∥2 ≤ 2L2(1− αµg)2N
µ2g +
4(τµg + Lρ) 2
µ4g (1− αµg)N−1D2. (41)
We finally upper-bound the last term at the right hand side of eq. (35) using Lemma 3.
∥∥Jµ − JN∥∥2F = d∑ i=1 ∥∥∇yNi,µ −∇yNi ∥∥2 ≤µ 2
2 L2J d(p+ 3) 3. (42)
Substituting eq. (40), eq. (41) and eq. (42) into eq. (35) yields
E ∥∥ĴN − J∗∥∥2F ≤6(4p+ 15) L2Qµ2g + µ 2 Q L2J dP4(p) + 6L2(1− αµg)2N µ2g
+ 12(τµg + Lρ)
2
µ4g (1− αµg)N−1D2 +
3µ2
2 L2J d(p+ 3) 3. (43)
Finally, the bound for the expected estimation error in eq. (34) becomes E ∥∥∇̂Φ(xk)−∇Φ(xk)∥∥2 ≤2L2D2(1− αµg)N + 4L4
µ2g D2(1− αµg)N + 24(4p+ 15)
L2M2
Qµ2g
+ µ2
Q L2JM
2dP4(p) + 24L2M2(1− αµg)2N
µ2g + 6µ2L2JM 2d(p+ 3)3
+ 48M2(τµg + Lρ) 2
µ4g (1− αµg)N−1D2. (44)
This completes the proof.

heading: F.2 PROOF OF THEOREM 1
text: Theorem 1 (Re-stated). Suppose that Assumptions 1, 2, and 3 hold. Choose the inner- and outerloop stepsizes respectively as α ≤ 1L and β = 1 4LΦ , where LΦ = L+ 2L 2+τM2 µg + ρLM+L 3+τML µ2g + ρL2M µ3g . Then, the iterates xk for k = 0, ...,K − 1 of ESJ in Algorithm 1 satisfy:
1
K K−1∑ k=0 E ∥∥∇Φ(xk)∥∥2 ≤ 16Lφ(Φ(x0)− Φ∗) K + 3D (45)
with Φ∗ = infx Φ(x) and D is the upper-bound established in Proposition 2 and is given by
D =2L2D2(1− αµg)N + 4 L4
µ2g D2(1− αµg)N + 24(4p+ 15)
L2M2
Qµ2g
+ µ2
Q L2JM
2dP4(p) + 24L2M2(1− αµg)2N
µ2g + 6µ2L2JM 2d(p+ 3)3
+ 48M2(τµg + Lρ) 2
µ4g (1− αµg)N−1D2. (46)
Proof of Theorem 1. Using the Lipschitzness of function Φ(xk), we have
Φ(xk+1) ≤Φ(xk) + 〈∇Φ(xk), xk+1 − xk〉+ Lφ 2 ∥∥xk+1 − xk∥∥2 ≤Φ(xk)− β〈∇Φ(xk), ∇̂Φ(xk)〉+ Lφ 2 β2 ∥∥∇̂Φ(xk)∥∥2
≤Φ(xk)− β〈∇Φ(xk), ∇̂Φ(xk)−∇Φ(xk)〉 − β ∥∥∇Φ(xk)∥∥2
+ Lφβ 2 (∥∥∇Φ(xk)∥∥2 + ∥∥∇̂Φ(xk)−∇Φ(xk)∥∥2)
≤Φ(xk)− ( β
2 − Lφβ2) ∥∥∇Φ(xk)∥∥2 + (β 2 + Lφβ 2) ∥∥∇̂Φ(xk)−∇Φ(xk)∥∥2. (47)
Let Ek[·] = Euk,1:Q [·|xk, yNk ] be the expectation over the Gaussian vectors uk,1, . . . , uk,Q conditioned on xk and yNk . Applying the expectation Ek[·] to eq. (47) yields
EkΦ(xk+1) ≤Φ(xk)− ( β
2 − Lφβ2) ∥∥∇Φ(xk)∥∥2 + (β 2 + Lφβ 2)Ek ∥∥∇̂Φ(xk)−∇Φ(xk)∥∥2 ≤Φ(xk)− ( β
2 − Lφβ2) ∥∥∇Φ(xk)∥∥2 + (β 2 + Lφβ 2)D, (48)
whereD represent the upper-bound established for the estimation error in Proposition 2. Now taking total expectation over Uk = {u1,1:Q, . . . , uk,1:Q}, we have
Ek+1 ≤Ek − ( β
2 − Lφβ2)EUk ∥∥∇Φ(xk)∥∥2 + (β 2 + Lφβ 2)D (49)
where Ek = EUk−1Φ(xk). Summing up the inequalities in eq. (49) for k = 0, . . . ,K − 1 yields
Ek ≤E0 − ( β
2 − Lφβ2) K−1∑ k=0 EUk ∥∥∇Φ(xk)∥∥2 + (β 2 + Lφβ 2)KD. (50)
Setting β = 14Lφ , denoting by Φ ∗ = infx Φ(x), and rearranging eq. (50), we have
1
K K−1∑ k=0 EUk ∥∥∇Φ(xk)∥∥2 ≤ 16Lφ(Φ(x0)− Φ∗) K + 3D. (51)
Hence, the proof is finished.

heading: G PROOFS FOR STOCHASTIC BILEVEL OPTIMIZATION
text: Define the following quantities
ĴN,j = ĴN (xk, uj) =  Y N1 (xk+µuj ;S)−Y N 1 (xk;S) µ u > j
... Y Nd (xk+µuj ;S)−Y N d (xk;S)
µ u > j
 , JN = ∂Y Nk∂xk , J∗ = ∂y ∗ k ∂xk
where uj ∈ Rp, j = 1, . . . , Q are standard Gaussian vectors and Y Nk is the output of SGD obtained with the minibatches {S0, ...,SN−1}. Conditioning on xk and Y Nk and taking expectation over uj yields
Euj ĴN,j = Euj  Y N1 (xk+µuj ;S)−Y N 1 (xk;S) µ u > j
... Y Nd (xk+µuj ;S)−Y N d (xk;S)
µ u > j

=  ∇>x Y N1,µ(xk;S) ...
∇>x Y Nd,µ(xk;S)  = Jµ(S)
where Y Ni,µ(xk;S) is the i-th component of vector Y Nµ (xk;S), which is the entry-wise Gaussian smooth approximation of vector Y N (xk;S). Let Ek[·] = E[·|xk, Y Nk ] = EDF ,u1:q be the expectation over the Gaussian vectors and the sample minibatch DF conditioned on xk and Y Nk .

heading: G.1 PROOF OF PROPOSITION 3
text: Proposition 3 (Re-stated). Suppose that Assumptions 1, 2, and 4 hold. Choose the inner-loop stepsize as α = 2L+µg . Define the constants
Cγ =(1− αµg) ( 1− αµg + α
γ + αL
γµg
) , Cxy = α ( α+ γ(1− αµg) + α L
µg
) , Cy = L
µg Cxy
Γ =2(τ2Cxy + ρ 2Cy)
σ2
µgLS + 2
L2
S (Cxy + Cy), λ = 2(τ
2Cxy + ρ 2Cy)D 2, (52)
where γ is such that γ ≥ L+µgµ2g . Then, we have:
E ∥∥JN − J∗∥∥2F ≤CNγ L2µ2g + λ(L+ µg) 2(1− αµg)CN−1γ (L+ µg)2(1− αµg)− (L− µg)2 + Γ 1− Cγ .
Proof of Proposition 3. Based on the SGD updates, we have
Y tk = Y t−1 k − α∇yG ( xk, Y t−1 k ;St−1 ) , t = 1, . . . , N
Taking the derivatives w.r.t. xk yields Jt =Jt−1 − α∇x∇yG ( xk, Y t−1 k ;St−1 ) − αJt−1∇2yG ( xk, Y t−1 k ;St−1 ) Jt − J∗ =Jt−1 − J∗ − α∇x∇yG ( xk, Y t−1 k ;St−1 ) − αJt−1∇2yG ( xk, Y t−1 k ;St−1
) + α ( ∇x∇yg (xk, y∗k) + J∗∇2yg (xk, y∗k)
) =Jt−1 − J∗ − α ( ∇x∇yG ( xk, Y t−1 k ;St−1 ) −∇x∇yg (xk, y∗k)
) − α (Jt−1 − J∗)∇2yG ( xk, Y t−1 k ;St−1
) + αJ∗ ( ∇2yg (xk, y∗k)−∇2yG ( xk, Y t−1 k ;St−1 )) .
Hence, using the triangle inequality, we have∥∥Jt − J∗∥∥F (i)≤∥∥ (Jt−1 − J∗) (I −∇2yG (xk, Y t−1k ;St−1)) ∥∥F + α
∥∥∇x∇yG (xk, Y t−1k ;St−1)−∇x∇yg (xk, y∗k)∥∥F + α
∥∥J∗ (∇2yG (xk, Y t−1k ;St−1)−∇2yg (xk, y∗k)) ∥∥F , where (i) follows from Assumption 1. We then further have∥∥Jt − J∗∥∥2F ≤
(1− αµg)2 ∥∥Jt−1 − J∗∥∥2F + α2∥∥∇x∇yG (xk, Y t−1k ;St−1)−∇x∇yg (xk, y∗k) ∥∥2F
+ α2 L2
µ2g ∥∥∇2yG (xk, Y t−1k ;St−1)−∇2yg (xk, y∗k) ∥∥2F + 2α(1− αµg)
∥∥Jt−1 − J∗∥∥F∥∥∇x∇yG (xk, Y t−1k ;St−1)−∇x∇yg (xk, y∗k) ∥∥F︸ ︷︷ ︸ P1
+ 2α(1− αµg) L
µg ∥∥Jt−1 − J∗∥∥F∥∥∇2yG (xk, Y t−1k ;St−1)−∇2yg (xk, y∗k) ∥∥F︸ ︷︷ ︸ P2
+2α2 L
µg ∥∥∇2yG (xk, Y t−1k ;St−1)−∇2yg (xk, y∗k) ∥∥F∥∥∇x∇yG (xk, Y t−1k ;St−1)−∇x∇yg (xk, y∗k) ∥∥F︸ ︷︷ ︸ P3 .
The terms P1, P2 and P3 in the above inequality can be transformed as follows using the Peter-Paul version of Young’s inequality.
P1 ≤ 1
2γ ∥∥Jt−1 − J∗∥∥2F + γ2∥∥∇x∇yG (xk, Y t−1k ;St−1)−∇x∇yg (xk, y∗k)∥∥2F , γ > 0 P2 ≤ 1
2γ ∥∥Jt−1 − J∗∥∥2F + γ2∥∥∇2yG (xk, Y t−1k ;St−1)−∇2yg (xk, y∗k)∥∥2F , γ > 0 P3 ≤ 1
2 ∥∥∇2yG (xk, Y t−1k ;St−1)−∇2yg (xk, y∗k)∥∥2F + 1
2 ∥∥∇x∇yG (xk, Y t−1k ;St−1)−∇x∇yg (xk, y∗k)∥∥2F Note that the trade-off constant γ controls the contraction coefficient (factor in front of
∥∥Jt−1 − J∗ ∥∥2 F
). Hence, we have∥∥Jt − J∗∥∥2F
≤ ( (1− αµg)2 + α
γ (1− αµg) +
αL γµg (1− αµg) )∥∥Jt−1 − J∗∥∥2F + ( α2 + αγ(1− αµg) + α2 L
µg )∥∥∇x∇yG (xk, Y t−1k ;St−1)−∇x∇yg (xk, y∗k)∥∥2F + ( α2 L2
µ2g + αγ
L µg (1− αµg) + α2 L µg )∥∥∇2yG (xk, Y t−1k ;St−1)−∇2yg (xk, y∗k)∥∥2F . Let Et−1[·] = E[·|xk, Y t−1k ]. Conditioning on xk and Y t−1 k and taking expectations yield
Et−1 ∥∥Jt − J∗∥∥2F ≤ Cγ ∥∥Jt−1 − J∗∥∥2F + CxyEt−1∥∥∇x∇yG (xk, Y t−1k ;St−1)−∇x∇yg (xk, y∗k)∥∥2F
+ CyEt−1 ∥∥∇2yG (xk, Y t−1k ;St−1)−∇2yg (xk, y∗k)∥∥2F , (53)
where Cγ , Cxy and Cy are defined as follows Cγ = (1− αµg) ( 1− αµg + α
γ + αL
γµg
) , Cxy = α ( α+ γ(1− αµg) + α L
µg
) , Cy = L
µg Cxy.
Conditioning on xk and Y t−1k , we have
Et−1 ∥∥∇x∇yG (xk, Y t−1k ;St−1)−∇x∇yg (xk, y∗k)∥∥2F
≤2Et−1 ∥∥∇x∇yg (xk, Y t−1k )−∇x∇yg (xk, y∗k)∥∥2F
+ 2Et−1 ∥∥∇x∇yG (xk, Y t−1k ;St−1)−∇x∇yg (xk, Y t−1k ) ∥∥2F
(i) ≤2L 2
S + 2τ2 ∥∥Y t−1k − y∗k∥∥2, (54) where (i) follows from Lemma 1 and Assumption 2. Similarly we can derive
Et−1 ∥∥∇2yG (xk, Y t−1k ;St−1)−∇2yg (xk, y∗k)∥∥2F ≤2L2S + 2ρ2∥∥Y t−1k − y∗k∥∥2. (55)
Combining eq. (53), eq. (54), and eq. (55) we obtain Et−1 ∥∥Jt − J∗∥∥2F ≤Cγ∥∥Jt−1 − J∗∥∥2F + 2(τ2Cxy + ρ2Cy)∥∥Y t−1k − y∗k∥∥2
+ 2 L2
S (Cxy + Cy). (56)
Unconditioning on xk and Y t−1k and taking total expectations of eq. (56) yield
E ∥∥Jt − J∗∥∥2F ≤CγE∥∥Jt−1 − J∗∥∥2F + 2(τ2Cxy + ρ2Cy)E∥∥Y t−1k − y∗k∥∥2 + 2L2S (Cxy + Cy)
(i) ≤CγE ∥∥Jt−1 − J∗∥∥2F + 2(τ2Cxy + ρ2Cy) (( L− µg L+ µg )2(t−1) D2 + σ2 µgLS )
+ 2 L2
S (Cxy + Cy),
where (i) follows from the analysis of SGD for a strongly-convex function. Let Γ = 2(τ2Cxy + ρ2Cy) σ2
µgLS + 2L
2
S (Cxy + Cy) and λ = 2(τ 2Cxy + ρ 2Cy)D 2. Then, we have
E ∥∥Jt − J∗∥∥2F ≤CγE∥∥Jt−1 − J∗∥∥2F + λ(L− µgL+ µg )2(t−1) + Γ. (57)
Telescoping eq. (57) over t from N down to 1 yields
E ∥∥JN − J∗∥∥2F ≤CNγ E∥∥J0 − J∗∥∥2F + λN−1∑
t=0
( L− µg L+ µg )2t CN−1−tγ + Γ N−1∑ t=0 Ctγ (58)
which, in conjunction with ( L−µg L+µg )2 ≤ 1− αµg and γ ≥ L+µgµ2g such that Cγ ≤ 1− αµg , yields
E ∥∥JN − J∗∥∥2F ≤CNγ L2µ2g + λCN−1γ N−1∑ t=0 ( (L− µg)2 (L+ µg)2(1− αµg) )t + Γ 1− Cγ
≤CNγ L2
µ2g +
λ(L+ µg) 2(1− αµg)CN−1γ
(L+ µg)2(1− αµg)− (L− µg)2 +
Γ
1− Cγ . (59)
The proof is then complete.
Lemma 8. Suppose that Assumptions 1, 2, 3, and 4 hold. Set the inner-loop stepsize as α = 2L+µg . Then, we have
E ∥∥Ek∇̂Φ(xk)−∇Φ(xk)∥∥2
≤8M2 ( CNγ L2
µ2g +
λ(L+ µg) 2(1− αµg)CN−1γ
(L+ µg)2(1− αµg)− (L− µg)2 +
Γ 1− Cγ + µ2 2 L2J d(p+ 3) 3
)
+ 2L2 ( 1 + 2 L2
µ2g )(( L− µg L+ µg )2N D2 + σ2 µgLS ) , (60)
where the expectation Ek[·] is conditioned on xk and Y Nk .
Proof of Lemma 8. Conditioning on xk and Y Nk , we have
Ek∇̂Φ(xk) = ∇xf(xk, Y Nk ) + J>µ ∇yf(xk, Y Nk ).
Recall ∇Φ(xk) = ∇xf(xk, y∗k) + J>∗ ∇yf(xk, y∗k). Thus, we have∥∥Ek∇̂Φ(xk)−∇Φ(xk)∥∥2 ≤2 ∥∥∇xf(xk, Y Nk )−∇xf(xk, y∗k)∥∥2 + 2∥∥J>µ ∇yf(xk, Y Nk )− J>∗ ∇yf(xk, y∗k)∥∥2
≤2L2 ∥∥Y Nk − y∗k∥∥2 + 4∥∥J>µ ∇yf(xk, Y Nk )− J>∗ ∇yf(xk, Y Nk )∥∥2
+ 4 ∥∥J>∗ ∇yf(xk, Y Nk )− J>∗ ∇yf(xk, y∗k)∥∥2
(i) ≤2L2 ∥∥Y Nk − y∗k∥∥2 + 4M2∥∥Jµ − J∗∥∥2F + 4L2µ2g ∥∥∇yf(xk, Y Nk )−∇yf(xk, y∗k)∥∥2
≤2L2 ∥∥Y Nk − y∗k∥∥2 + 8M2∥∥Jµ − JN∥∥2F + 8M2∥∥JN − J∗∥∥2F + 4L4µ2g ∥∥Y Nk − y∗k∥∥2.
where (i) applies Lemma 4 and Assumption 3. Taking expectation of the above inequality yields E ∥∥Ek∇̂Φ(xk)−∇Φ(xk)∥∥2 ≤2L2(1 + 2L2
µ2g
) E ∥∥Y Nk − y∗k∥∥2 + 8M2E∥∥JN − J∗∥∥2F
+ 8M2E ∥∥Jµ − JN∥∥2F . (61)
Using the fact that Y Ni (xk; ·) has Lipschitz gradient (Proposition 1) and Lemma 3, the last term at the right hand side of eq. (61) can be directly upper-bounded as
∥∥Jµ − JN∥∥2F = d∑ i=1 ∥∥∇Y Ni,µ(xk;S)−∇Y Ni (xk;S)∥∥2 ≤ µ22 L2J d(p+ 3)3, (62) whereLJ is the Lipschitz constant of the JacobianJN (and also of its rows∇Y Ni (xk;S)) as defined in Proposition 1. Combining eq. (61), eq. (62), Proposition 3, and SGD analysis (as in eq. (60) in Ji et al. (2021)) yields
E ∥∥Ek∇̂Φ(xk)−∇Φ(xk)∥∥2
≤8M2 ( CNγ L2
µ2g +
λ(L+ µg) 2(1− αµg)CN−1γ
(L+ µg)2(1− αµg)− (L− µg)2 +
Γ 1− Cγ + µ2 2 L2J d(p+ 3) 3 ) + 2L2 ( 1 + 2 L2
µ2g )(( L− µg L+ µg )2N D2 + σ2 µgLS ) . (63)
This finishes the proof.

heading: G.2 PROOF OF PROPOSITION 4
text: Proposition 4 (Re-stated). Suppose that Assumptions 1, 2, 3, and 4 hold. Set the inner-loop stepsize as α = 2L+µg . Then, we have:
E ∥∥∇̂Φ(xk)−∇Φ(xk)∥∥2 ≤ ∆ + B1 (64)
where ∆ = 8M2 ((
1 + 1 Df ) 4p+15 Q + 1 Df ) L2 µ2g +2M 2 Df + ( 1 + 1 Df ) 4M2 Q µ2dL2JP4(p)+ 4M 2 Df µ2dL2JP3(p)
and B1 respresents the upper bound established in Lemma 8.
Proof of Proposition 4. We have, conditioning on xk and Y Nk Ek ∥∥∇̂Φ(xk)−∇Φ(xk)∥∥2 =Ek∥∥∇̂Φ(xk)− Ek∇̂Φ(xk)∥∥2 + ∥∥Ek∇̂Φ(xk)−∇Φ(xk)∥∥2. (65)
Our next step is to upper-bound the first term in eq. (65). Ek ∥∥∇̂Φ(xk)−Ek∇̂Φ(xk)∥∥2
≤2Ek ∥∥∇xF (xk, Y Nk ;DF )−∇xf(xk, Y Nk )∥∥2
+ 2Ek ∥∥Ĵ>N∇yF (xk, Y Nk ;DF )− J>µ ∇yf(xk, Y Nk )∥∥2
≤ 2M 2
Df + 4Ek ∥∥∇yF (xk, Y Nk ;DF )∥∥2∥∥ĴN − Jµ∥∥2F + 4Ek
∥∥Jµ∥∥2F∥∥∇yF (xk, Y Nk ;DF )−∇yf(xk, Y Nk )∥∥2 ≤ 2M 2
Df + 4M2
( 1 + 1
Df ) Ek ∥∥ĴN − Jµ∥∥2F + 4M2Df ∥∥Jµ∥∥2F , (66)
where the last two steps follows from Lemma 1. Next, we upper-bound the term Ek ∥∥ĴN − Jµ∥∥2F .
Ek ∥∥ĴN − Jµ∥∥2F = 1QEk∥∥ĴN,j − Jµ∥∥2F , j ∈ {1, . . . , Q}
≤ 1 Q ( Ek ∥∥ĴN,j∥∥2F − ∥∥Jµ∥∥2F)
≤ 1 Q d∑ i=1 ( Ek ∥∥∥Y Ni (xk + µuj ;S)− Y Ni (xk;S) µ uj ∥∥∥2 − ∥∥∇Y Ni,µ(xk;S)∥∥2) . (67) Recall that for a function h with L-Lipschitz gradient, we have
Eu ∥∥∥h(x+ µu)− h(x) µ u ∥∥∥2 ≤ 4(p+ 4)∥∥∇hµ(x)∥∥2 + 3 2 µ2L2(p+ 5)3. (68)
Then, applying eq. (68) to function Y Ni (·;S) yields
Euj ∥∥∥Y Ni (xk + µuj ;S)− Y Ni (xk;S)
µ uj ∥∥∥2 ≤ 4(p+ 4)
∥∥∇Y Ni,µ(xk;S)∥∥2 + 32µ2L2J (p+ 5)3.
Hence, eq. (67) becomes
Ek ∥∥ĴN − Jµ∥∥2 ≤4p+ 15
Q d∑ i=1 ∥∥∇Y Ni,µ(xk;S)∥∥2 + 3µ2dL2J2Q (p+ 5)3 ≤2(4p+ 15)
Q
d∑ i=1 ∥∥∇Y Ni (xk;S)∥∥2 + µ2dL2JQ P4(p) ≤2(4p+ 15)
Q ∥∥JN∥∥2F + µ2dL2JQ P4(p), (69) where P4 is a polynomial of degree 4 in p. Combining eq. (66), eq. (69) and Lemma 3 yields
Ek ∥∥∇̂Φ(xk)− Ek∇̂Φ(xk)∥∥2 ≤2M2
Df + 4M2
( 1 + 1
Df
)(2(4p+ 15) Q ∥∥JN∥∥2F + µ2dL2JQ P4(p)) + 4M2
Df
( 2 ∥∥JN∥∥2F + µ2dL2JP3(p))
≤ 8M2 (( 1 + 1
Df
)4p+ 15 Q + 1
Df )L2 µ2g + 2 M2 Df
+ ( 1 + 1
Df
)4M2 Q µ2dL2JP4(p) + 4M2 Df µ2dL2JP3(p)
= ∆, (70) where ∆ = 8M2 ((
1 + 1 Df ) 4p+15 Q + 1 Df ) L2 µ2g + 2M 2 Df + ( 1 + 1 Df ) 4M2 Q µ2dL2JP4(p) +
4M2
Df µ2dL2JP3(p).
Taking total expectations of both eq. (65) and eq. (70) and combining, we have E ∥∥∇̂Φ(xk)−∇Φ(xk)∥∥2 =E∥∥∇̂Φ(xk)− Ek∇̂Φ(xk)∥∥2 + E∥∥Ek∇̂Φ(xk)−∇Φ(xk)∥∥2
≤∆ + B1 (71) where B1 representing the upper-bound established in eq. (63). The proof is then complete.

heading: G.3 PROOF OF THEOREM 2
text: Theorem 2 (Re-stated). Suppose that Assumptions 1, 2, 3, and 4 hold. Set the inner- and outer-loop stepsizes respectivelly as α = 2L+µg and β = 1 4LΦ
, where L = max{Lf , Lg} and the constant LΦ is defined as in Theorem 1. Then, the iterates xk, k = 0, ...,K − 1 of the ESJ-S algorithm satisfy
1 K ∑K−1 k=0 E ∥∥∇Φ(xk)∥∥2 ≤ 16(Φ(x0)−Φ∗)LΦK + 3B1 + ∆, (72) where ∆ and B1 have the following forms
∆ =8M2 (( 1 + 1
Df
) 4p+ 15
Q +
1
Df
) L2
µ2g + 2
M2 Df +
( 1 + 1
Df
) 4M2
Q µ2dL2JP4(p)
+ 4M2
Df µ2dL2JP3(p)
B1 =8M2 ( CNγ L2
µ2g +
λ(L+ µg) 2(1− αµg)CN−1γ
(L+ µg)2(1− αµg)− (L− µg)2 +
Γ 1− Cγ + µ2 2 L2J d(p+ 3) 3
)
+ 2L2 ( 1 + 2 L2
µ2g )(( L− µg L+ µg )2N D2 + σ2 µgLS ) and the constants Γ, λ, and Cγ are defined in Proposition 3.
Proof of Theorem 2. Using the Lipschitzness of function Φ(xk), we have
Φ(xk+1) ≤Φ(xk) + 〈∇Φ(xk), xk+1 − xk〉+ Lφ 2
∥∥xk+1 − xk∥∥2
≤Φ(xk)− β〈∇Φ(xk), ∇̂Φ(xk)〉+ Lφ 2 β2 ∥∥∇̂Φ(xk)−∇Φ(xk) +∇Φ(xk)∥∥2
≤Φ(xk)− β〈∇Φ(xk), ∇̂Φ(xk)〉+ Lφβ2 (∥∥∇Φ(xk)∥∥2 + ∥∥∇̂Φ(xk)−∇Φ(xk)∥∥2)
Hence, taking expectation over the above inequality yields EΦ(xk+1) ≤EΦ(xk)− βE〈∇Φ(xk), ∇̂Φ(xk)〉+ Lφβ2E ∥∥∇Φ(xk)∥∥2
+ Lφβ 2E ∥∥∇̂Φ(xk)−∇Φ(xk)∥∥2. (73)
Also, we have −E〈∇Φ(xk), ∇̂Φ(xk)〉 =− E〈∇Φ(xk), ∇̂Φ(xk)−∇Φ(xk)〉 − E ∥∥∇Φ(xk)∥∥2
=E [ −〈∇Φ(xk),Ek∇̂Φ(xk)−∇Φ(xk)〉 ] − E ∥∥∇Φ(xk)∥∥2 ≤E ( 1
2 ∥∥∇Φ(xk)∥∥2 + 1 2 ∥∥Ek∇̂Φ(xk)−∇Φ(xk)∥∥2)− E∥∥∇Φ(xk)∥∥2 = 1
2 E ∥∥Ek∇̂Φ(xk)−∇Φ(xk)∥∥2 − 1 2 E ∥∥∇Φ(xk)∥∥2,
which, in conjunction with eq. (73), yields
EΦ(xk+1) ≤EΦ(xk) + β 2 E ∥∥Ek∇̂Φ(xk)−∇Φ(xk)∥∥2 − (β 2 − Lφβ2 ) E ∥∥∇Φ(xk)∥∥2
+ Lφβ 2E ∥∥∇̂Φ(xk)−∇Φ(xk)∥∥2. (74)
Setting β = 14LΦ and using the bounds established in Lemma 8 and Proposition 4, we have
EΦ(xk+1) ≤EΦ(xk) + β
2 B1 −
β 4 E ∥∥∇Φ(xk)∥∥2 + β 4 (B1 + ∆).
Summing up the above inequality over k from k = 0 to k = K − 1 yields
EΦ(xK) ≤EΦ(x0) + β
2 KB1 −
β
4 K−1∑ k=0 E ∥∥∇Φ(xk)∥∥2 + β 4 K(B1 + ∆).
Rearranging the above inequality yields
1
K K−1∑ k=0 E ∥∥∇Φ(xk)∥∥2 ≤4 (EΦ(x0)− EΦ(xK)) βK + 3B1 + ∆
≤16 (Φ(x0)− Φ ∗)LΦ
K + 3B1 + ∆.
The proof is then complete.

Review summaries:
The paper introduce two new algorithms to tackle bi-level optimisation, which has wide-ranging application in machine learning. Crucially, these two algorithms forego second order gradient computation completely, favoring instead an evolution strategies approach.  While this approach has been proposed before in this context, the new algorithms are more efficient because they do not treat the whole problem as a black box, and leverage the structure of the hypergradient. Theoretical convergence analyses are then provided for both the full gradient and the stochastic gradient alternatives, as well as experiments on various bi-level optimisation problems, showcasing improved performance over previous approaches.

This work mainly focuses on proposing a computational efficient for approximating the response Jacobian matrix in the hypergradient of bilevel optimization with only first-order information of the loss functions. The proposed methods are applicable for developing both deterministic and stochastic bilevel optimization algorithms. The author further show the convergence guarantees and computational complexity of the proposed algorithms. Multiple numerical results showcase the efficiency of the proposed method in applications of machine learning.

The paper aims to reduce the computational cost in solving a bilevel problem. In bilevel optimization, the outer gradient can be decomposed into direct and indirect gradients. The indirect gradient requires computing the response Jacobian, which involves computing & inverting the Hessian and computing mixed derivative. To avoid expensive second-order approximations, the authors propose to use the evolution strategies (ES) based method to approximate the response Jacobian. The authors provide the convergence guarantee of the algorithm and empirically demonstrate that the proposed algorithm performs well on large-scale bilevel problems. 

This paper studies the evolution strategies (ES) based method to estimate the Jacobian matrix used in bilevel optimization. The main idea is to use first-order information to estimate the second-order information in the Jacobian matrix. Some convergence results have been provided, companied with some simulations.