ES-Based Jacobian Enables Faster Bilevel Optimization
Authors: Daouda Sow
Kaiyi Ji
Yingbin Liang
Venue: ICLR 2022 Submitted (2022)
Keywords: Evolution Strategies
Computational Complexity
Jacobian Matrix
Stochastic Algorithm
Bilevel Optimization

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 Jacobian- or/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 meta-learning 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. 

Under review as a conference paper at ICLR 2022
ES-BASED JACOBIAN ENABLES FASTER BILEVEL OP- TIMIZATION
Anonymous authors Paper under double-blind review
ABSTRACT
Bilevel optimization (BO) has arisen as powerful tool for solving many modern machine learning problems. However, due the nested structure of BO, exist- ing gradient- based methods require second-orde derivative approximations via Jacobian- or/and Hessian-vector computations, which can be very costly in prac- tice, especially with large neural network models. In this work, we propose novel BO algorithm, which adopts Evolution Strategies (ES) based method approx- imate the response Jacobian matrix in the hypergradient of BO, and hence fully eliminates all second-order computations. stands for the ES-basei Jacobian method) and further extend it to the stochas- tic 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 various bilevel problems. Particularly, in experiment in the few-shot meta- learning problem, we meta-learn the twelve millions parameters of a ResNet- 12 network over the minilmageNet dataset, which evidently demonstrates the scala- bility of our ES- based bilevel approach and its feasibility in the large-scale setting.
call our algorithm as ESJ (which
INTRODUCTION
Bilevel optimization has recently arisen as a powerful tool to capture various modern machine learn- ing problems, including meta-learning (Bertinetto et al.) 2018 Franceschi et al. [0018] Rajeswaran et al. 2019 Ji et al. 2020a Liu et al. 2021a), hyperparamater optimization (Franceschi al. 2018
Shaban et al. 2019 neural architecture search (Liu et al., 2018a| Zhang et al. (022)) etc. Bilevel optimization generally takes the following mathematical form:
min (xx) = f(x, (e.vv((r) (I:) s.t. y²(x) arg min g(r, VER*
(1)
ZER
where the outer and inner objectives f RP x Rd + R and g RP differentiable with respect to (w..t.) the inner and outer variables x E RP and y E Rd. In many machine learning problems, the loss functions f and g take a finite-sum form over set of given data Dn,m {Ei, Sj 1, I,...,'.,, n, = I,...,m) as follows:
X xd Rd + R are both continuously
f(x, v) 1 E F(x,
g(z,y)
E G(z,
where the sample sizes and m are typically very large. In this work, we consider popular setting where the total objective function (xx) in eq. [ is possibly nonconvex w.r.. 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. [017 Franceschi et al.) 2017 Shaban et al.] 2019 Ra
jeswaran et al. 2019 Liu al. 2020a) and the approximate iterative differentiation (AID) method
(Domke 2012 Pedregosa 2016 Gould et al., 2016 [iio et af.] 2018 Lorraine et al.; 2020). In
the ITD method, the objective value f(x, ((x)) first computed using an approximation point (x)
close to y (z), and the hypergradient V = afiz.y' (a)) is then approximated by 2f(x.0(x))
dy using
dgg
either reverse forward mode automatic differentiation. Alternatively, the AID method leverages

---

Under review as a conference paper at ICLR 2022
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 in- volve 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- methods, Evolution Strategies (ES) can be a good candidate eliminate the second-order computations, which uses function val- ues for estimating its gradient in blackbox optimization (Nesterov & Spokoiny (017). An ES-based bilevel optimizer has recently been proposed in Gu al.) (2021) for hyperparameter optimization, which uses the ES method to approximate the hypergradient of ((x), treating D (z) fully as a black- box objective. However, the hypergradient bilevel optimization has a much more involved struc- ture (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 design 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-basel bilevel optimizers, which are scalable for high- sigh-dimensionll 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- estimator can yield provably convergent algorithm for bilevel optimization. So far the ES-based bilevel optimizer in Gu et al. (2021) did not come with polynomial-tine 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 both questions.
1.1 MAIN CONTRIBUTIONS
Novel ES-based Jacobian bilevel optimizer. We propose a novel bilevel optimizer with the ES- based Jacobian estimation, which we call as ESJ. In contrast to the existing ES-based optimizer in Gu et al. 22021), 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 (2) with respect 2)), 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 compu- tations, and computationally much more efficient than existing AID and ITD bilevel optimizers that require Hessian- and/or Jacobian-vector product computations. We further propose stochastic bilevel optimizer ESJ-S to allow the algorithm scalable over large dataset.
Convergence guarantee. Theoretically, we characterize the convergence rate and their computa- tional complexity for both ESJ and ESJ-S to achieve an e-accurate solution. Technically, in contrast to the standard analysis of ES-based methods on the smoothed blackbox function values, we de- velop 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 perfor- mance 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,
meta-leam the twelve millions weights of ResNet- 12 network. Previ- ous 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.
1.2 RELATED WORK
Bilevel optimization. Bilevel optimization has been studied for decades since Bracken & McGill (1973). A variety of bilevel optimization algorithms were then developed, including constraint- based approaches (Hansen et al. 1992 Shi et al. 2005 Moore 2010), approximate implicit differ-
2

---

Under review as a conference paper at ICLR 2022
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.
2013
Franceschi et al. 2017, Finn al. 2017 Shaban et al., 2019; Rajeswaran al. [019 Liu
et al.
2020a). These methods often suffer from expensive computations of second-order information
(e.g.- Hessian-vector products). Hessian-free algorithms
recently proposed based interior-
point method (Liu et al. (2221b) and Evolution Strategies (ES) (Gu et al. (022). 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 al.) 2021 Rujeswaran al., 2019 Ji & Liang 202T). This paper provides the convergence analysis for the proposed Hessian- free ES- based approach.
Stochastic bilevel optimization. Ghadimi & Wang (018); Ji et al. (2021), Hong et al. (2020) pro- posed stochastic gradient descent (SGD) type bilevel optimization algorithms by employing Neu- mann series for flessian-inversee inverse- vector approximation. Recent works (Khanduri al., [2021a bj Chen et al. 2022 Guo et al. [021 Guo & Yang Yang et al., (202)) then leveraged momentum- based variance reduction to further reduce the computational complexity of existing SGD-type bilevel optimizers. In this paper, we propose a stochastic ESSbased method, which eliminate the computation of second-order information required by all aforementioned stochastic methods.
Please Appendix B for more works about applications of bilevel optimization and ES methods.
2
PROPOSED ALGORITTM
We describe our proposed deterministic and stochastic algorithms for bilevel optimization.
2.1
HYPERGRADIENTS
The key step in popular gradient- 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(x) V y (x)) + J=(x) (T))T))))) V y (xx)
(3)
where the Jacobian matrix [-(x) = (x)(z) Rd Following Lorraine et al. 2020 it can be thr seen that (4b(x) contains two components: the direct gradient Vzf(z, y (x)) and the indirect gra- dient J. (x) Vy vf y" (z)). The direct component can be efficiently computed using the existing automatic differentiation techniques. The indirect component, however, is computationally much more complex, because (x) takes the form of J= = [V (I. y" y=))] V.V,g y (z, y (z)) (if V (x, y" (x:) invertible). which contains the Hessian inverse and the second-order mixed derivative. Some approaches mitigate the issue by designing iacobian-vecton and Hessian-vector products (Pedregosal 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 prob- lems with high-dimg sigh-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.
2.2 EVOLUTION STRATEGIES (ES) METHOD
Evolution Strategies (ES) is a powerful technique to estimate the gradient of function based on function values, when 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-bos box function h R" + R using the following oracle based only on the function values (i.c., the zeroth-order information)
ha( h(x + um) h(x)
(4)
where u E I is a Gaussian random vector and μ > is the smoothing parameter. Such an oracle can be shown to be unbiased estimator of the gradient of the smoothed function E,, [h(x μou)).
2.3 ES ES-BASED JACOBIAN AND ITS ENABLED BILEVEL OPTIMIZERS
Drawback of an existing ES-based bilevel optimizer. In Gu et al. (002I). the ES method was directly applied to estimate the hypergradient of ((:),, treating ((x) blackbox objective. Since
3

---

Under review as a conference paper at ICLR 2022
the hypergradient has complex form in eq. ( and can be very sensitive both inner and outer functions, such an estimation will likely have a large bias error. As our experiments in Section 4 demonstrate, such method is not robust and even fails to converge in complex bilevel optimization problems (such as with neural network parameters).
Our key idea is exploit the analytical structure of the hypergradient in eq. ( where the derivatives [ y" (x)) and f y (x)) can be computed efficiently and accurately, and then use the ES method only estimate the Jacobian (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 ingredients: (i) for a given I, apply an algorithm to solve the inner optimization problem and use the output as an approximation of y (a); for example, the output y (z) of N gradient descent steps of the inner problem can serve as estimate for y" (z). Then IN (z) (x) serves an estimate of [x);;;
8m
and (ii) construct an ES-based Jacobian estimator IN (z;; u) E R for IN (x) as
IN (x; u) (x + pm) y (x)
(5)
where и E R2 a Gaussian vector with independent and identically distributed (i.i.d.) entries. Then for any vector E Rdd, the Jacobian-vector product can be efficiently computed using only vector-vector dot product [x( u) v = (6(x; u), v) u, where 6(x; u) = (w+ww)- (w) E Rd.
ESJ: A novel bilevel optimizer with ES Jacobian oracles. We design a bilevel optimizer (see Algorithm D using the ES-based Jacobian oracle in eq. (5), which we call the ESJ algorithm. At each step k of the algorithm, ESJ runs an -tep full GD to approximate y (Ix). ESJ then samples Q Gaussian vectors {uR, E (O, [) 1,... Q), and for each sample ukyj runs N-step full GD to approximate V (xx + s) .) and then computes the Jacobian estimator IN (z; uR.;) 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.
( V V.I(ua
23=1
Q Ha V (
V VZf(Ze, v
2 23=1 (d(Iz; ((s)) ( (x unj-
(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 computa- tionally much more efficient as demonstrated in our experiments.
ESJ-S: A stochastic bilevel optimizer with ES Jacobian oracles. For the finite-sum sum problem with the objective functions given in eq. (2) we design stochastic bilevel optimizer (see Algorithm in Appendix A based on ES Jacobian oracles, which we call as ESJ-S.
Differently from Algorithm 1 which applies GD updates to find (xx), ESJ-S adopts N stochas- tic gradient descent (SGD) steps to find (Yg", ye ) to the inner problem, each with the outer variable set to be I& + pukej. Note that all SGD runs follow the same batch sampling path {So. .N., The Jacobian estimator IN (1z3 uk,j) can then be computed as in eq. ( At the outer level, ESJ-S samples a new batch Dr independently from the inner batches (So,..,, SN-1} to evaluate the stochastic gradients V., F(*,, ; Dr) and F,F( V Dr). The hypergradient then estimated as follows.
(xx))) = V VZF(EN, YN; ; ④ EJ (5(0)) uns), V ((K,, ;D))) un.j.
(7)
3 CONVERGENCE ANALYSIS
3.1 TECHNICAL. ASSUMPTIONS
In this work, we focus on the following types of loss functions. Assumption 1. The inner function (I y) is pyy strongly convex with respect to the variable y. For the sumi case, the same assumption holds for VG(x, ().
4

---

Under review as a conference paper at ICLR 2022
Algorithm Bilevel optimizer via ES Jacobians (ESJ)
1: Input: lower- and upper-level stepsizes a, 8 > 0, initializations Io E R and yo E R", inner and outer iterations numbers K and N, and number of Gaussian vectors Q.
2: for 1,,,,,.... 2,.,,....
. K do
3: Sct yl y0, yk.5 yo 1,.,.
4: for = 1, 2, N do
5: Update yi = vi 1-1 aV Vog(en. 11
end for
for 1= Q
8: Generate uk. = (0, 1) E R
9:
for 1,2, N do
10:
Update Uh.s = Hki 1-1
aV aa (It + peanss s
11:
end for
Compute 63 -
12:
1
13:
end for
14:
Compute (xx) = Vaf(on, + Egnn (6,, Var(zn, uns
15:
Update 2k+1 = - (z4(
16:
end for
We take the following assumptions on the inner and outer loss functions g(,, y) and f(x,, v), as also adopted in Ghadimi & Wang (2018): al. (2021); ); Yang al. (202T).
Assumption 2. Lef w (x, v) E Rd+p. The gradient Vg(w) is L-g-iipschitz continuous, i.e., for any w1, W2 E TRd+P, [ Vg(w2))| I L.g [1] wa]l; further, the derivatives Vzg(w) and V.V. vug(w) are p- and r-Lipschitz continuous, i.e, (2z((n)) Vig(w2)l1 pl[w1 w2ll and VyVgg(u - V v 1 r||ur - ww2l- For the finite-sum case, the same assumptions hold for G(w; C).
Assumption 3. Let w (x. Rd+P. The objective f(w) and its gradient V f(w) are M- and Lipsschht continuous, i.e., for any W1, w2 E Rdtp.
([w1))f(w>) -
S M || wall, [[Vf(w)) V f((22)] s L w2/l
For the finite-sum sum case, the same assumptions hold for F(w: E).
For the finite-sum setting, take the standard bounded-variance assumption VG(w; Assumption 4. Gradient (G(;; () has a bounded variance, i.e., [el|VG((w; [)(w)]| - o² for some constant o 2
3.2 CONVERGENCE ANALYSIS FOR ESJ
Different from the standard zeroth-order analysis for blackbox function, here we develop new techniques to analyze the ES-based estimator that depend on the entire inner optimization trajectory, which is specific bilevel optimization. We first establish the following essential proposition which
8y "
characterizes the Lipshitzness property of the approximate Jacobian matrix IN
Proposition Suppose that Assumptions 1 and hold. Let LJ
(1+ (
I el with
I. max [1 L.g). Then, the Jacobian [N(x) is Lipschitz continuous with constant l.g: [N((F1) (x2) s Ly || 2211 Vzl, 222 E RP.
(8)
We next provide an upper-bound on the hypergradient estimation error for ESJ. Proposition 2. Suppose that Assumptions 2 and
3 hold. Consider the ESJ algorithm. Its hyper-
gradient estimation error can be upper-bounded as:
E [)()))=> (=(x))]| s (11-app)))
+ p2 + pd
(9)
Q
where [[|] is conditioned 3g and taken over the Gaussian vectors {uR. ; 1,...QQ
5

---

Under review as a conference paper at ICLR 2022
By using the smoothness property in Proposition and the upper-bound in Proposition 2 we provide the following characterization of the convergence rate for ESJ.
Theorem (Convergence of ESJ). Suppose that Assumptions
7 and hold. Choose the inner-
and outer-loop stepsizes respectively as a S T, and B A1.o" where L. = L. 21 +M +
pLM4 + +rMM + pt Then, the iterates Ikk Fx for 1 O,.... K of ESJ in Algorithm satisfy:
K
1 Ek= E. (
(0/00 4-1 + A
K
(10)
with x- = =inf, (xx) and A O
+ Edp + p ppdp')
Further, choose K, N. Jt, and Q at the order of O Ollog 1) Ominn P V dp a 1) and (E)). Then ESJ achieves an e-accurate stationary point with O log }) computations of the gradient V Vog(1,) and (2) computations of the gradients (a y) and V (x.
Theorem I characteriizss 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 A 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 IN 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 Jn of the smoothed output E.y² (xrx + nu), which decreases sublinearly w.r.t. the batch size Q, and (c) the error between the Jacobians IN and Ju. which can be controlled by the smoothness parameter μ.
Theorem D also indicates that ESJ requires only gradient computations to converge, and eliminates Hessian- and Jacobian-vector vector products required by the existing AID and ITD based bilevel optimiz- (Pedregosa (2016). Franceschi et al. (2018), , Grazzi et al. (2020)) Thus, ESJ is computationally much more efficintt particularly for high- -igh--imensional problems.
3.3 CONVERGENCE ANALYSIS FOR ESJ-S
In this section, we apply the stochastic algorithm ESJ-S to the finite-sum objective in eq. (2) and an- alyze its convergence rate. We first note that the Lipschitzness property established in Proposition is general and can also be used here. The following proposition establishes an upper bound on the
estimation error of Jacobian J. by IN a where Y is the output of N inner SGD updates.
8m
Proposition 3. Suppose that Assumptions D and 4 hold. Consider the ESJ-S algorithm. Choose the inner-loop stepsize as a Ltm.
2 where I. max (T. L.g. Then, we have:
L2
X(T + (1 C
r
J J. ICN
+ 1- Cry
(L+ 1,) axpl)) - p3)
where X, I, and C2 < are constants (see appendix G for their forms).
We next provide an upper- upper-bound bound on the estimation error of the hypergradient by (4 in eq. 7
Proposition 4. Suppose that Assumptions 7 3 and 4 hold. Consider the ESJ-S algorithm. Ser 2 the inner-loop stepsize as L+p where L max L,, L,L)) Then, we have:
E ()(( vo(aw)| s (1-
1
1
P +
dd
+
# pd
S
D
Q
Q
where S and D, are the sizes of the inner and outer mini-batchen respectively.
Based on Propositions
and 4 we characterize the convergence rate for ESJ-S
Theorem 2 (Convergence of ESJ-S), Suppose that Assumptions and outer-lopp stepsizes respectively as a 2_ and TL. where
2 3 and f hold. Set the inner-
max max(LFf Lg} and
Ltpy Then, the iterates Ik. k = 0, .... K - 1 of ESJ-S satisfy:
the constant L are defined in Theorem
R Σ (o
16(4(zz) -- LL + A2;
(11)

---

Under review as a conference paper at ICLR 2022
where A2 = O (11 a1)))))))* x4 + Dj + s'de" Q + upppp)
For ESJ-S to attain e-accurate stationary point, it suffices select K, S, and D, at the order of ) and N. μ., and Q respectively at the order of Ollog ) O fe p dp p dp)) and ( which thus amount to O & log ( computations of the gradient VyG(r, y: E) and ( computations of the gradients Fx E.E) and V Fa, V.E).
Comparing to the convergence emror A in Theorem 1 for the deterministic algorithm ESJ, Theo- rem for the stochastic algorithm ESJ-S captures two more sublinearly decreasing error terms $ and Dj 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. ( Theorem D suggests that
ESJ attains an e-accurate stationary point of eq. (2) with O pm ( computations of F. As a comparison, Theorem suggests that ESJ-S achieves O fs log #) computations of VG and (+)) computations of V., which outperforms ESJ in the large sample regime where the sample sizes n, m > 1 This also reflects the advantage of SGD over GD.
e log = computations of VG and
EXPERIMENTS
We validate our algorithms over three bilevel problems: shallow hyper-representation (HR) with linear/2-layer net embedding model synthetic data, deep HR with LeNet network (LeCun al. 1998) on MNIST dataset, and few-shot meta-learning with ResNet-11 12 on minilmageNet 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.
(a) d = 128
(b) d = 256
(c) ESJ-N: N inner GD steps
12
AB.c0
AD-02
10:00
455
135
130
(d) Outer loss
(e) Hypergradient norm
(f) ESJ-N: N inner GD steps
10200
AB.05
KD1F
mBn
tg
(92)
Sahaapone
Figure 1: First row: HR with linear embedding model. Second row: HR with two-layer net.
4.1 SHALLOW SEPPER-REPRESENTATUUE ON SYNTHETIC DATA
The hyper-representation (HR) problem (Franceschi et al. 2018 Grazzi et al., 2020 searches for a regression (or classification) model following two-phased optimization process. The inner-level identifies the optimal linear regressor parameters and the outer level solves for the optimal em- bedding model (i.e., representation) parameters X. Mathematically, the problem can be modeled by the following bilevel optimization:
I[T(X;; - Yill s.t. w ((()) argmin 1 [[ A)u YaIl [| (12)
min f(A) RE
2m
wekd 2nz
where X2 E Xm and X1 E x are matrices of synthesized training and validation data, and Y, E "2 Y, E Td 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 X, X2 and labels Y,, Y, 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).
7

---

Under review as a conference paper at ICLR 2022
Figure 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 (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 ap- proximate the entire hypergradient, which likely incurs large estimation error. In contrast, our ESJ exploits the structure of the hypergradient and uses ES only to estimate the response Jacobian that the estimation of hypergradient is more accurate. Such an advantage is more evident under a two- layer net model, where HOZOG does not converge as shown in [igure[[[d) This can be explained by the flat hypergradient norm shown in Figure which indicates that the hypergradient esti- mator in HOZOG fails to provide a good descent direction for outer optimizer. Figure (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 N in practice is not costly.
4.2 DEEP HYPER- OOPEERRRERESENTTTION ON MNIST DATASET
In order to demonstrate the advantage of our proposed algorithms in large neural net models, we perform deep repeeetepresnnttaiion to classify MNIST images by learning entire LeNet network. The corresponding bilevel problem given by
min Lout (A) Dout A
Σ (A)((xA)) y.)
(Z..B))EDom.
(13)
s.t. (²(A)= = arg min Lin (w. A) Din wERexp
Σ (L( A). + {
(E
where f(x;; A) E I corresponds to features extracted from data point Vi, L(:,)) is cross- entropy loss function, = 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, apply the more efficient stochastic algorithm ESJ-S in Algorithm With with a minibatch size B = 256 to estimate the inner and outer losses Lin and Lout-
Eq
85-5
AD-CS
AB-9P
S0.50
5555
AD-G
05,5-10
Naveeng time
(a) Accuracy on outer data
(b) Outer loss value
(c) ESJ-S-N: N inner SGD steps
Figure 2: Deep HR the MNIST dataset.
Figure 2 compares the classification accuracy on the outer dataset Dout 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 level of performance, and instead saturate around an accuracy of 0.93. Further, Figure (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.
4.3 FEW-SHOT
OVER MINIIMAGENET
To study our algorithms over larger neural nets, we study the few-shot image recognition problem, where classification tasks T, 1,..., are sampled over distribution Py. In particular, we consider the ANIL meta-learning method (Raghu et al. 2019:1 Ji et al., 2020a), where all tasks share common embedding features parameterized by Ø, and each task T: has its task-specific parameter w, 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 T. fixes d and minimizes its loss function over a training set S, to
8

---

Under review as a conference paper at ICLR 2022
obtain its adapted parameters w.- At the outer stage, the meta-learner computes the test loss for each task Ti using the parameters (Ø, w() over a test set D,, and optimizes the parameters ø of the common embedding function by minimizing the meta-objective meta over all classification tasks. The problem can be expressed as the following bilevel optimization
min L (6, ) E CD. (Ø, w()
+
(14)
s.t. = arg min Ladant (0, (1) == 1 m Ei= Ls, (0. wi);
where we collect all task-specific parameters into (Wi, ... Wm and the corresponding mini- mizers into w" = (wi, ..., w The functions Ls. (Ø, wx) = |s,1 Eces, (C(o, w;: () + R(w;)) and LD. (0,w;) w;) = D. ΣEED, C ( Efi correspond respectively to the training and test loss functions for task Tr with R a strongly-convex regularizer and L classification loss function. In our setting, since the task-specific parameters correspond the weights of the last linear layer, the inner-level objective Cadapt (Ø, 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 the sizes of the datasets D, and S, are usually small in few-shot learning and full GD can be applied here. Hence, we use ESJ (Algorithm ① here. Also since the number of tasks in few-shot classification datasets often very large, 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 these tasks.
1
Algorithm
Time (hours)
ESJ
2.8
E3
Prstatit
MetaOptNet
20+
VetsOppiet
AKL
WL
ProtNet
4.4
MAN.
Saniingdte 2ite
Raniing 201e
(a) Test accuracy (ResNet-12 12)
(b) Test accuracy (CNN4)
(c) Time to reach 69% (ResNet-12)
Figure 3: 5way-5shot image classification on the minilmageNet dataset on single GPU.
We conduct few-shot meta-learning on the minilmageNet dataset (Vinyals et al. (016) 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),
ANII (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
ct al. 2019).
Figure a)) and (b) show that our algorithm ESJ converges faster than the other baseline methods. Also, Comparing FigureB] Fa) 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 3c) shows that MetaOptNet did not reach 69% accuracy after 20 hours 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.
CONCLUSION
In this paper, we propose a novel ES-based approach for bilevel optimization, which eliminates all second-order computations in the gradient-bard baadee 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-tim complexity can be achieved. We anticipate that our approach will be useful for accelerating bilevel algorithms in various machine learning problems.

---

Under review as a conference paper at ICLR 2022
REPRODUCIBILITY CHECKLIST
To ensure reproducibility, we use the Machine Learning Reproducibility Checklist v2.0, Apr. 2020 (Pineau 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
clear explanation of any assumptions. All assumptions are stated 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 tine/computational complexity analysis for our algorithms in Theorems 2
For any theoretical claim,
A clear statement of the claim. A clear statement of theoretical claims are made in Section
complete proof of the claim. The complete proofs of all claims are available in Appendix E Appendix[F] and Appendix
For all datasets used, check if you include:
The relevant statistics, such as number of examples. We use widely adopted datasets MNIST and minilmageNee Section 4 The related statistics can be seen at htt Lecun. cww/excb///eeooootttcttttpt....c hotos//ni////
The details of train/validation//es splits. We give this information in the Supple- mentary Appendix
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 this work, check you include:
Training code. The training code is available in code supplementary material. Evaluation code. The evaluation code available in our code supplementary ma- terial.
(Pre-)trained model(s). We do not release any pre-trdined 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 erate results. We provide all details of the hyper-paramete tuning in Supplemen- tary Appendix
hyper-parameters used to gen-
clear definition of the specific measure or statistics used to report results. We the classification accuracy test-set and the loss the train-set.
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, estimated energy cost. We report the running time of the algorithms in Section
- A description of the computing infrastructure used. All detailed descriptions are presented in Section
10

---

Under review as a conference paper at ICLR 2022
REFERENCES
Sebastien M.R. Amold, Praateek Mahajan, Debajyoti Datta, and Ian Bunner. learn2learn, 2019.
Juhan Bae and Roger B. Grosse. Delta-stn: Efficient bilevel optimization for neural networks using structured response jacobians. CoRR, abs/2010. 2020.
Luca Bertinetto, Joao F Henriques, Philip Torr, and Andrea Vedaldi. Meta-learning with differ- entiable closed-form solvers. In International Conference on Learning Representations (ICLR), 2018.
Jerome Bracken and James T McGill. Mathematical programs with optimization problems in the constraints. Operations Research, 21(1):37-44 44, 1973.
Tianyi Chen, Yuejiao Sun, and Wotao Yin. A timescale stochastic bilevel optimization method. arXiv preprint arXiv:2102.04400 02.01.00..,..... 2021.
Justin Domke. Generic methods for optimization-based modeling. International Conference Artificial Intelligence and Statistics (AISTATS), pp. 318-326, 2012.
Chelsea Finn, Pieter Abbeel, and Sergey Levine. Model-agnostic meta-learning for fast adaptation of deep networks. In Proc. International Conference on Machine Learning (ICML), pp. 1126-1 135, 2017.
Luca Franceschi, Michele Donini, Paolo Frasconi, and Massimiliano Pontil. Forward and reverse gradient- hyperparameter optimization. International Conference on Machine Learning (ICML), pp. 1165-1173, 2017.
Luca Franceschi, Paolo Frasconi, Saverio Salzo, Riccardo Grazzi, and Massimiliano Pontil. Bilevel programming for hyperparameter optimization and meta-learning. In International Conference on Machine Learning (ICML), pp. 1568-1577, 2018.
Saeed Ghadimi and Guanghui Lan. Stochastic first-and zeroth-order methods for nonconvex stochas- tic programming. SIAM Journal on Optimization, 23(4):2341-255 2013.
Saeed Ghadimi and Mengdi Wang. Approximation methods for bilevel programming. arXiv preprint 10X2:0/0000022466 2018.
Stephen Gould, Basura Fernando, Anoop Cherian, Peter Anderson, Rodrigo Santa Cruz, and Edison Guo. On differentiating parameterized argmin and argmax problems with application to bi-level optimization. arXiv preprint a607v: 00547,, 2016.
Riccardo Grazzi, Luca Franceschi, Massimiliano Pontil, and Saverio Salzo. On the iteration com- plexity of hypergradient computation. International Conference on Machine Learning (ICML)), 2020.
Bin Gu, Guodong Liu, Yanfu Zhang, Xiang Geng, and Heng Huang. Optimizing large-scale hy- perparameters via automated learning algorithm. CoRR, 2021. URL httptt://///xxi...o.. abs/2102. 09026
Zhishuai Guo and Tianbao Yang. Randomized stochastic variance-reduced methods for stochastic bilevel optimization. arXiv preprint arXiv:2105.02266, 2021.
Zhishuai Guo, Yi Xu, Wotao Yin, Rong Jin, and Tianbao Yang. stochastic moving-average
estimators for non-convex optimization. arXiv preprint arXiv:2104. 14840, 2021.
Nikolaus Hansen. The ema evolution strategy: a comparing review. Towards a new evolutionary computation, pp. 75- 102, 2006.
Pierre Hansen, Brigitte Jaumard, and Gilles Savard. New branch-and bound rules for linear bilevel programming. SIAM Journal Scientific and Statistical Computing,
1992.

---

Under review as a conference paper at ICLR 2022
Mingyi Hong, Hoi-To Wai, Zhaoran Wang, and Zhuoran Yang. A two-timescale framework for bilevel optimization: Complexity analysis and application to actor-critic. arXiv preprint arXiv:2007.05170, 2020.
Kaiyi Ji and Yingbin Liang. Lower bounds and accelerated algorithms for bilevel optimization. arXiv preprint Xiv:22002.0392. 2021.
Kaiyi Ji, Zhe Wang. Yi Zhou, and Yingbin Liang. Improved zeroth-order variance reduced algo- rithms and analysis for nonconvex optimization. In International Conference on Machine Learn- ing (ICML), pp. 3100-3109, 2019.
Kaiyi Ji, Jason D Lee, Yingbin Liang, and H Vincent Poor. Convergence of meta-learning with task-specific adaptation over partial parameter. In Advances Neural Information Processing Systems (NeurIPS), 2020a.
Kaiyi Ji, Junjie Yang, and Yingbin Liang. Multi-step madel-agn meta-learning: Convergence
and improved algorithms. arXiv preprint arXtv:2002.07836, 2020b.
Kayi Ji, Junjie Yang, and Yingbin Liang. Bilevel optimization: Convergence analysis and enhanced design. International Conference on Machine Learning (ICML)), 2021.
Prashant Khanduri, Siliang Zeng, Mingyi Hong, Hoi-To Wai, Zhaoran Wang, and Zhuoran Yang. A momentum momentum-assisted assisted single-time timescale stochastic approximation algorithm for bilevel optimization. arXiv e-prints, pp. arXiv-2108 2021a.
Prashant Khanduri, Siliang Zeng, Mingyi Hong, Hoi-To Wai, Zhaoran Wang, and Zhuoran Yang. A near-optimal algorithm for stochastic bilevel optimization via double-momentum. arXiv preprint arXiv:2102.07367, 2021b.
Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. International Conference Learning Representations (ICLR), 2014.
Yann LeCun, Léon Bottou, Yoshua Bengio, and Patrick Haffner. Gradient-based learning applied to document recognition. Proceedings of the IEFE, (6(11):2278- -2224, 1998.
Kwonjoon Lee, Subhransu Maji, Avinash Ravichandran, and Stefano Soatto. Meta-learning with differentiable convex optimization. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 10657-10665, 2019.
Renjie Liao, Yuwen Xiong, Ethan Fetaya, Lisa Zhang, Xaq Yoon, KiJung Pitkow, Raquel Urtasun, and Richard Zemel. Reviving and improving recurrent puekppropagttion. International Confer- ence on Machine Learning (ICML)), 2018.
Hanxiao Liu, Karen Simonyan, and Yiming Yang. Darts: Differentiable architecture search. arXiv 2018a.
preprint 17Xiv:1806.040
Risheng Liu, Pan Mu, Xiaoming Yuan, Shangzhi Zeng, and Jin Zhang. A generic first-order al- gorithmic framework for bi-level programming beyond lower-level singleton. In International Conference on Machine Learning (ICML), 2020a.
Risheng Liu, Jiaxin Gao, Jin Zhang, Deyu Meng, and Zhouchen Lin. Investigating bi-level opti- mization for learning and vision from unified perspective: A survey and beyond. arXiv preprint 101. 2021a
Risheng Liu, Xuan Liu, Xiaoming Yuan, Shangzhi Zeng, and Jin Zhang. A foilee-funniion-baee interior-poine method for non-convex bi-level optimization. In International Conference on Ma- chine Learning (ICML), 2021b.
Sijia Liu, Bhavya Kailkhura, Pin- Yu Chen, Paishun Ting, Shiyu Chang, and Lisa Amini. Zeroth- order stochastic variance reduction for nonconvex optimization. In Advances in Neural Informa- fion Processing Systems (NeurIPS), pp. 3731-3741, 2018b.
12

---

Under review as a conference paper at ICLR 2022
Sijia Liu, Songtao Lu, Xiangyi Chen, Yao Feng, Kaidi Xu, Abdullah Al-Dujaili, Mingyi Hong, and Una-May O' O'Reilly. Min-max optimization without gradients: Convergence and applications to black-box evasion and poisoning attacks. In International Conference on Machine Learning (ICML), pp. 6282-6293, 2020b.
Jonathan Lorraine, Paul Vicol, and David Duvenaud. Optimizing millions of hyperparameters by implicit differentiation. International Conference Artificial Intelligence and Statistics (AIS- TATS), pp. 1540-1552, 2020.
Matthew Mackay, Paul Vicol, Jonathan Lorraine, David Duvenaud, and Roger Grosse, Self-tuning networks: Bilevel optimization of hyperparameters using structured best-response functions. In International Conference on Learning Representations (ICLR), 2018.
Dougal Maclaurin, David Duvenaud, and Ryan Adams. Gradient-barn based hyperparameter optimiza- tion through reversible learning. In International Conference on Machine Learning (ICML), pp- 2113-2122, 2015.
Gregory M Moore. Bilevel progranming algorithms for machine learning model selection. Rensse- laer Polytechnic Institute, 2010.
Yurii Nesterov and Vladimir Spokoiny. Random gradient-free minimization of convex functions. Foundations of Computational Mathematics, pp. 527-566, 2017.
Fabian Pedregosa. Hyperparameter optimization with approximate gradient. In International Con- ference on Machine Leurning (ICML), pp. 737-746, 2016
Joelle Pineau, Philippe Vincent-Lamarre, Koustuv Sinha, Vincent Lariviere, Alina Beygelzimer, Florence d' Alché Buc, Emily Fox, and Hugo Larochelle. Improving reproducibility in machine learning research. Journal of Machine Learning Research, 22:1-20: 2021.
Aniruddh Raghu, Maithra Raghu, Samy Bengio, and Oriol Vinyals. Rapid learning or feature reuse? towards understanding the effectiveness MAML. International Conference on Learning Rep- resentations (ICLR), 2019.
Aravind Rajeswaran, Chelsea Finn, Sham M Kakade, and Sergey Levine. Meta-learning with im- plicit gradients. In Advances in Neural Information Processing Systems (NeurIPS), pp. 113-124, 2019.
Olga Russakovsky, Jia Deng, Hao Su, Jonathan Krause, Sanjeev Satheesh, Sean Ma. Zhiheng Huang. Andrej Karpathy, Aditya Khosla, Michael Bernstein, Alexander C. Berg, and Fei- Fei. Imagenet large scale visual recognition challenge. International Journal of Computer Vision, 3(115):211-25) 2015.
Amirreza Shaban, Ching-An Cheng, Nathan Hatch, and Byron Boots. Truncated back proppprrpag for bilevel optimization. In International Conference Artificial Intelligence and Statistics (AIS- TATS), pp. 1723-1732, 2019.
Chenggen Shi, Jie Lu, and Guangquan Zhang. An extended kuhn-tucker approach for linear bilevel programming. Applied Mathematics and Computation, 162(1):51-63 2005.
Jake Snell, Kevin Swersky, and Richard Zemel. Prototypical networks for few-shot learning. In Advances in Neural Information Processing Systems (NIPS), 2017.
Xingyou Song, Wenbo Gao, Yuxiang Yang, Krzysztof Choromanski, Aldo Pacchiano, and Yunhao Tang. Es-maml: Simple hessian-free meta learning. In International Conference on Learning Representations (ICLR), 2019.
Konstantinos Varelas, Anne Auger, Dimo Brockhoff, Nikolaus Hansen, Ouassim Ait ElHara, Yann Semet, Rami Kassab, and Frédéric Barbaresco. A comparative study of large-scale variants of cma-es. In International Conference Parallel Problem Solving from Nature, pp. 3-15. Springer, 2018.
Oriol Vinyals, Charles Blundell, Timothy Lillicrap, and Daan Wierstra. Matching networks for one shot learning. In Advances in Neural Information Processing Systems (NIPS), 2016.
13

---

Under review as a conference paper at ICLR 2022
Tengyu Xu, Zhe Wang, Yingbin Liang, and H Vincent Poor. Gradient free minimax optimization: Variance reduction and faster convergence. arXiv preprint arXiv:2006.0936/ 2020.
Junjie Yang, Kaiyi Ji, and Yingbin Liang, Provably faster algorithms for bilevel optimization. arXiv preprint arXiv:2106.040000 2021.
Miao Zhang, Steven Su, Shirui Pan, Xiaojun Chang, Ehsan Abbasnejad, and Reza Haffari. idarts: Differentiable architecture search with stochastic implicit gradients. arXiv preprint arXiv:2106. 10784, 2021.
Daniel Zügner and Stephan Günnemann. Adversarial attacks on graph neural networks via meta learning. In International Conference Learning Representations (ICLR), 2019.
14

---

Under review as a conference paper at ICLR 2022
Supplementary Material
A STOCHASTIC BILEVEL OPTIMIZER BASED ON ES METHOD
In the following, we present the algorithm specification for our proposed stochastic bilevel optimizer based on ES method, which we call ESJ-S.
Algorithm Stochastic bilevel optimizer via ES Jacobians (ESJ-S)
1: Input: lower- and upper-level stepsizes a, 0, initializations 2g Rs and yo E R inner and outer iterations numbers K and N, and number of Gaussian vectors Q.
2: for k = 0. ,,.,, do
3: Set Y = yo: Y = yo. j = 1,...,Q Q
4: Generate Uj.j (O, I) IR,, 1=1,..
5: for = 1,2, N do
6: Draw a sample batch S
7:
Update Y - V G( V (1-1))
8:
for j 1,..., ....QQQG do
9:
V Update Yh. k.j
v [xk + pun.j 11-1; Sc-1
end for
11: end for
-Y, Yp
Compute dj-
j = 1,.,.
13: Draw a sample batch Dr
14: Compute ( = V YE Dr) + E9 (4, F(r. Y :r)) wn.j
Update I++1 = It (7
16:
end for
B ADDITIONAL RELATED WORK
In this section, we provide further related work on the applications of bilevel optimization and
methods.
Bilevel optimization applications. Bilevel optimization has been employed various applications such few-shot meta- meta-learning (Snell et al. 2017 Franceschi et al. 2018 Rajeswaran al.) 2019 Zügner& & Günnemann, 2019 et al. 2020aa b). hyperparameter optimization (Franceschi et al,
Mackay et al. et al. 2021). etc. This paper demonstrates the superior performance of the proposed ES-based bilevel optimizer in meta-learning and hyperparameter optimization.
2017
2018 8 Shaban et al. (2019), neural architecture search (Liu et al. 201 8a /hhan
ES applications. Evolution Strategies (ES) method has been studied for more than two decades (see review papers (Hansen 2006 Varelas al. 2018)). In particular, Nesterov Spokoiny (2017) pro- posed 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 metaaremffreemment learning (Song et al., (019). hy hyperparameter optimization (Gu al., (021). adversarial machine learning (Ji al. 2019 Liu et al. 2018b), min- imax optimization (Liu et al. 2020b, Xu et al.) (220), etc. This paper proposes 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- 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.
15

---

Under review as a conference paper at ICLR 2022
C FURTHER SPECIFICATIONS FOR EXPERIMENTS IN SECTION 4
We note that the smoothing parameter (in Algorithms and was easy to set and value of 0.1 or 0.01 yields a good starting point across all our experiments. The batch size Q (in Algorithms and 2 is fixed to 1 (i.e., we use one Jacobian oracle) in all our experiments.
C.I SPECIFICATIONS ON BASELINE BILEVEL APPROACHES IN SECTION 4.1
We compare algorithm ESJ with the following baseline methods:
HOZOG (Gu et al. 2021): a hyperparameter optimization algorithm that uses evolution strate- gies 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 al., 2019): approximate implicit differentiation comm
with conjugate gradient. We use its implementation provided at
proleazner/bbyeettoree.
AID-FP (Grazzi et al., 2020): approximate implicit differentiation with fixed-point. We exper- com
imented with its implementation at the repositery http:///ggtuuu....mcoooeeeeeaa 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
com wcw./prolearreer/yyeettooo.
C.2 HYPERPARAMETERS DETAILS FOR SHALLOW HR EXPERIMENTS IN SECTION 4.1
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 & Bal 2014) with learning rate of 0.05. The value of y in eq. (12) is set to be 0.1. For 0.1 for ESJ and HOZOG. For all methods, set N 10, = 0.001, 3 = 0.001, and use Adam with a learning rate of 0.01 as the outer optimizer.
the two-layer net case,
C.3 SPECIFICATIONS ON BASELINE STOCHASTIC ALGORITTMMS SECTION
We compare our stochastic algorithm ESJ-S with the following baseline stochastic bilevel algo- rithms.
stocBiO: an approximate implicit differentiation method that uses Neumann Series to esti- mate the Hessian inverse. We use its implementation available at httttts://gtthh/ com/ Jur:JleYang37/Stoccl
AID-CG-S and AID-FP-S: stochastic versions of AID-CG and AID-FP, respectively. We use their implementations in the repository httpe://ggttaabb.cc//rrrreaaaaae
hypertorch
C.4 SPECIHICATIONS FOR FEW-SHOT
SECTION 4.3
The minilmageNet dataset (Vinyals et al. 2016) is a large-scale benchmark for few-shot learning generated from ImageNet (Russakovsky et al.J 2015) Russakovsky. The dataset consists of 100 classes with each class containing 600 images of size 84 84. Following (Arnold et al. 2019), we classes for meta-validation, and classes for
split the classes into 64 classes for meta-training, meta-testing. More specfically, we use 20000 tasks for meta-training, 600 tasks for meta-vall voca--vaidation, 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 learning rate of a = 0.01 as inner optimizer and Adam with a learning rate of 3 = 0.01 as outer (meta) optimizer.
16

---

Under review as a conference paper at ICLR 2022
025
-880
ES]
ES]
665
Validattoo
HOZOG
HOZOG
AID-FP
AID-FP
REVERSE
REVERSE
660
3.8
1/2
12
0.25
0.50
1.00
1.25
1.50
2.00
Running time
Running time
Figure 4: Classification results on 20 Newsgroup dataset. Left: number of inner GD step N 5.
Right: number of inner GD steps N = 10. Running time is in seconds.
D EXPERIMENTS ON HYPERPARAMETER OPTIMIZATION
Hyperparameter optimization (HO) is the problem finding the set of the best hyperparamters (ei- ther representational 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 training loss (usually regularized) for the given hyperparameters and then the outer problem mini- mizes over the hyperparameters. Hence, HO can be mathematically expressed as follows
min (A) Paal ΣEEDIN (A)())
(15)
1 s.t. w² (A(A)= = arg min (w, Darl ECED,, (C(w, x: + R(w, A)),
10
where C is loss function (c.g., logistic loss), R(w, A) is regularizer, and Dtr and Deel are respec- tively training and validation data. Note that the loss function used 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 D and Drr) to compute validation and training losses, and by adding regularizer term on the training loss.
Following (Franceschi al. 2017 Grazzi et al. 2020), we perform classification the 20 News- group dataset, where the classifier modeled by an affine transformation, the cost function L the cross-entrpy loss, and R(w, X) is strongly-conved regularizer. We set one (g-regularization hyperparameter for each weight w, so that À and w have the same size.
For ESJ and HOZOG, use GD with learning rate of 100 and momentum of 0.9 perform the inner updates. The outer learning rate is set to be 0.02. We set the smoothing parameter (μ in Algo- rithm to be 0.01. For AID-FP, AID-CG, and REVERSE in their implementations accompanying the paper Grazzi et al. (2020).
the suggested hyperparameters
can be from Figure 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, we have seen in the other experiments (not for HO), ESJ is more robust and stable bilevel optimizer than HOZOG, and it achieves good performance across many bilevel problems.
method ESJ slightly outperforms HOZOG and converges faster
E SUPPORTING TECHNICAL LEMMAS
In this section, provide auxiliary lemmas that used for proving the convergence results for the algorithms ESJ and ESJ-S.
In the following proofs, let I. max (Tgg [. and D be such that [x)] < D.
17

---

Under review as a conference paper at ICLR 2022
First recall that for any two matrices A Rmxr and B E R"*", have the following upper- bound on the Frobenius norm of their product,
|| F || B|
(16)
The following lemma follows directly from the Lipschitz properties Assumptions 2 and 3
Lemma 1. Suppose that Assumptions 2 and and3hood. hold. Then, the stochastic derivatives V y; E),
V. V Grr, E), and Ga, y: {) have bounded variances, ie.. for any (z, y) and & we have:
3 V FF - ff(z, 1) - M
E V, - v V ( v) 01 - L2;
V3 v: - 1| - L?
Using Lemma 2.2 Ghadimi & Wang 2018). the following lemma characterizes the Lipschitz property of the gradient of the total objective + x)) [ex ((0,v"((r) (1)).
Lemma 2. Suppose that Assumplions 7 and 3 hold. Then, we have:
(22 V La | 22 æ1l
Vas R2 22 R"
where the constant L I. +
MM pLM+L'+TML
eL A
Pg
We next provide some essential properties of the ES- (i.c., zeroth-order) gradient oracle in eq. (4) due to Nesterov & Spokoiny (2017).
Lemma 3. Let h R² R be differentiable function with L-Lipsctitz gradient. Define its Gaussian smooth approximation by(x) E,, [M(x + pu]), where μ > and u E R is standard Gaussian random vector. Then, hy is differentiable and have:
The gradient of by takes the form
Vhy(x) = Ey, h(x + pux) l((x)) n
For any a E R
S 122
[h, - Vh(x>)]
L(n + 338/22
2
For any x E R
h(x + pu) h(x) u 2 4(n (x)]
3 + z²" n=²,, 1. (n
Ey
Note the first item in Lemma implies that the oracle eq. is indeed an unbiased estimator of the gradient of the smoothed function hy.
Lemma 4. Suppose that Assumptions and 2 hold. The Jacobian J. [g"(x)
dx has bounded norm:
L
J
S
(17)
fg
Proof of Lemma From the first order optimality condition of y" (z), we have Vug (z, (x)) = 0.
Hence, the Implicit Function Theorem implies:
J. =
-V (z, y" (z))] V V Vug (z, y (x))
(18)
Taking norms and applying eq.
16 together with Assumptions [nd yield the desired result
|
< V V (z, (x)) (x, (I)l [ (x))
(19)
Pg
18

---

Under review as a conference paper at ICLR 2022
Lemma 5, Suppose that Assumptions
| and 2 hold. The Jacobian IN
ByC Jat
has bounded norm:
L
| || S
(20)
Pg
Proof of Lemma
The inner loop gradient descent updates writes
uk = u
Vug (It, 1) t == 1, I.....N.. N.
Taking derivatives w.r.t 2/ yields
Tt = Ji- aV V v2 (IR: Uk 12 0. Vig (Ik, Uk 1)
= Jr-1 (1- V. (IK, 3/ 11) V Bn (22. 1)
Telescoping over from 1 to N yields
N
N-1 (I a V gg (Tk, vi)) - a Σ
N-1 II 112+1
IN = To
II
Vx v Vyg Tky 8k)
(I- a (0k,3/)
0=0
V 99 (Ik. (
II
(I- a V gg (Irk (x))
(21)
mmt++1
Hence, we have
N
N
|J F
Sa a
V Va be (F
|ell
II
(1 - a (Tx,
m 1++11
||
2 (IxK, vi*)
m=t+1
(8)
N
5 aL
(1 (1-app)) (1p2))
t=
N-1
L
aL
apt))*
Pg
t=0
where (i) follows from Assumption 2 and (it) applies the completes the proof.
of function g(x, :). This
Lemma 6. Suppose that Assumptions 7 and 2 hold. Then, the Jacobian in the stochastic algorithm
8Y.5 has bounded norm, as shown below.
ESJ-S JN
0sg
| || <
(22)
Proof of Lemma
The proof follows similarly Lemma
E.1 PROOF OF PROPOSITION
Proposition (Re-stated), Suppose that Assumptions
7 [add] 2
hold. Define the constant
T = 1
L
pL
(23)
Mg
Pg
22
Then, the Jacobian TN(x) = (x) is L -ipsschhit with respect to x under the Frobenious norm:
dx
[x1) IN(e2)ll S Tog ||10
-xall
Van
R.. 32 R.P
(24)
19

---

Under review as a conference paper at ICLR 2022
Proof of Proposition
L
Using eq.
21 we have
JN(x) - a N-1
N-1
V (z
II m=1+1
(1- ([=,vgg(,
x 12
Hence, for 21 E Rp and x2 E R2, we have
[IJn(1)) - (12))l
N-1
N-1
=0
V. V Ve
II (7- ( (x)))
m=++1
N-1
V.Vyg 9 (12.
(T - a ( (12))
m=
N-1
Vgg
II (7-a a
m=++1
N-1
V V2V 60 (12))
II m=+++
(T - a V ( (x2))
V V bn (m)) 11 - A1
0
[ (2I, V
V wg (2 1|
(25)
where we define A(xx)
m
N-1 (I- (z. yrhh (x)) and (i) follows from eq.
(16)
Next we upper bound the quantity [A((a)) (x )|. as shown below.
|A((z)) (ax))) و (V/g (x2,5't1 (x2)) - Vigg (rz, (12,, (12))) (z1)
+ (1 - Vig (z2, (x2)) (Ac+1
A (4+))
Sal A++ (z)
(-1, (21)) - Vig (x2, (
+ [ - a (22, (42)) | (x1) - (z21
<(1 aH4 ||| (141 - (z2)
- ap,) -111 | -1 x22 ||
+ + +
(26)
where the last inequality follows from Lemma Telescoping eq. (26) over yields
5 and Assumptions I and
|| Ax(11) - A( 5(1 N- 11N---21AN-2(220) - AN-2(122 ||
N-1-3 3
ap(1+ + 1 )1 attg) ) N-12-m -aaty) ) - x2l
Σ
Mg
m=u
(1 (nd
N 1-2 N
3ig (21, -2 (m1)) Big (22 y (x2))
N-1-3
L - (stg) N - x2ll
+ Σ ap(1 +
Pg
m=0
(1)
L
Sap(1+
+ _ (11 - capty) N || - z22l
[1
L (1 ) og1g)
+ (N - 1. - 22)pp(1 +
N -1- 21 - x22
20

---

Under review as a conference paper at ICLR 2022
=ap(1 + - [N 1)(1 ap1)) - 20 x2
(27)
Mg
where (i) follows from Lemma and Assumption
2
Replacing eq.
27 in eq.
25 and using
Assumption 2 we have
[x - JN(22)]1
- ) -t 1)(1 - T2
$ a
Lap(1
1g
N-1
L ) aa41)) (uug))~11 N- || - 22ll
+ a
Σ
r(1 +
Pg
1=0
N-1
I a Lp(1 La²Lp(1+ + _ 0||21 x22
Σ
m(1 - a44g))
+
5 - -1+ _
=)|[21 x2
Ug
Mg
Hg
m=0
1 (1 + - |||0 - -02l +
PL
J (1
+
- ) x x22
L
(28)
g
Hg
Fg
Vg
where we use Σ N -m= mx m
28), which can be obtained by taking derivatives for the expression N- x" with respect to æ. Hence, rearranging and using the definition of L in
in eq.
a²
eq. (23) finishes the proof.
Lemma 7. Suppose that Assumptions 7 Dand
Define the constant
Lg = 1
1.
T
pl.
Mg
Pg
Then, the Jacobian [x(x)= 88 (m:
is L -ipschhhtt with respect to x under the Frobenius norm:
dx
[[
Ir S L || - x2
Va E E IRP.
Proof of Lemmal
The proof follows similarly that for Proposition T
F PROOFS FOR DETERMINISTIC BILEVEL OPTIMIZATION
For notation convenience, define the following quantities:
yf" (Ir +mn,) - [ank)
1
Oyk
Oyk 0zk
IN.J = (x)
IN
J.
(29)
0z
Ma (oo +pa,) (aFA)
where uy R 1, Q are standard Gaussian vectors. Let (TR)) be the Gaussian smooth approximation of yi" (Ix). We collect v (xx) for 1,...,dd d together vector y/ (Ix)), which the Gaussian approximation of the vector y (IL) If μ > O, s/ (xx) is differentiable and we let Jn be the Jocobian given by
Oy (œx)
Jμ
(30)
Ozy
We approximate aee using the average ES estimator given by IN = 2 E9 In.. The hypergra- dient is then approximated as
on
( = v.f(ZK v JJ v f v
Q
=Vaf(ER,
) + Q j=1
JN v Vyf(E.
(31)
21

---

Under review as a conference paper at ICLR 2022
Let 6, n" (ax ((x+))))))) (xx) and let d,g be the i-th component of dy. Hence, we have
82
IN-=
IN Vof = و1 و2 = (ōj, Vyf uk )) uj-
dagug Vof(En,
(32)
Using eq. 31 and eq. 32 the estimator for the hypergradient thus be computed
Q (ő.. oo ve us.
(Φ(()) = V v )
Q
F.1 PROOF OF PROPOSITION 2
Proposition 2 (Re-stated). Suppose that Assumptions 7 and hoold Then, the expectedestimation estimation enror can be upper-bounded as follows:
L
D2 at1)) + +4444p= + 15) I2
E (xx))) (<<x))| 52I - -app) +4
14
Qu4
12
12 M2y + 24T.2
a1p )NN 13
M 3)$
Q
1/3
48M²(TH, + Lp)²
(1 - (11))))) N D?2 ND??
(33)
where the expectation [[|] is conditioned on and yk N
Proof of Proposition Based on the definitions of ( and vo(2x) and conditioning on and y we have
(xx) - (Ix())
V / vx vN Vz vi 2B| V
y (x) - J.T.VV Vol(
<2L² vi vill + J 1| V ffff( vt v fo vill
+ 4B JN V v )11
ff
(1)
221. (2 - 01 +
- vill + 1--22B|1 J2
(0)
L
$ 22L2D2 aH,))* + 4
D2 #
any 4M2 JN J.
(34)
where (1) follows from Lemma 4 and Assumption and (i) and (ii) also use the following result for full GD (when applied to a strongly-convex function).
xill < (0,,, D2
Next, we upper-bound the last term Bl| - T. at the last line of eq. (3p. First note that
Bl J.11
3E|| Jnlle- + IN J. + 3 Jn IN 1|
(35)
22

---

Under review as a conference paper at ICLR 2022
We then upper- bound each term of the right hand side of eq. 35 For the first term, we have
Q
JN
Jullz
-B1
Σ
IN.S - JJulls
1=1
Q
1
Σ
(JN. -Jx) 1|r
Q2
i=1
Q
1 Q
E
Σ
|| -2 (Ix, - -n..Ns-- (Ja)
+
1
Σ
BIJNN, -JullE
Q2
|l| II Julle jE QL.
(36)
We next upper-bound the term B| B| Jyllir in eq. (36)
EIJNNs - Jull =E JN. 1 - -1|J
(37)
4(p + 4) V
3 5
V yi ||
3 (4p 15) V 2 + 22
(38)
where (i) follows by applying Lemma to the components of vector y (ax)) which have Lipschitz
gradients by Proposition Then, noting that | 12 < V ||= + +22²² rz (p+3)3 and replacing
in eq. (38) we have
E| ElJJNs -onl-- S Σ (2(4p+ 15) V +
52(4p 15) JN 1| D3
(1) 22(4p + 15)
I2
m + dPL
(39)
where (i) follows from Lemma and P4 a polynomial of degree 4 in p. Combining eq. (36) and eq.
338
yields
L
Jl Jn $2(4p + 15)
+ o d
(40)
Quz
We next upper-bound the second term the right hand side of eq. using eq. (41) in Ji et al. (2021), as shown below.
35 which can be upper-bounded
| J.11² 2L )22
4(TH, + Lp)²
(1 app)) N-
(41)
p3
mg
We finally upper-bound the last term the right hand side of eq. (35) using Lemma 3
1r, - Ix 1
V M - 2
1=1
12
L3
12dd(p + 3)
(42)
2
23

---

Under review as a conference paper at ICLR 2022
Substituting eq.
40 eq. (41) and eq. (22) into eq.
35 yields
L
+ 61 )2N
Bl| J. 1. <6(4p + 15)
L'
-
Quy
Q
19
12 + 2
3u²
(1-a N- +
Lí d
(43)
3
2
Finally, the bound for the expected estimation error in eq.
34 becomes
La D2 attlg) )N + 24(4p= + 15) L2
E| ((())) [ $2I -) +
Q14
2² L2 MdP + SAL2M22( (1-44))
att ) L3 M d(p + 3)3
4 (THg 2 + =p)? (1 (211)) D2
(44)
This completes the proof.
F.2 PROOF OF THEOREM
Theorem (Re-stated). Suppose that Assumptions 7 2 and hold. Choose the inner- and outer-
loop stepsizes respectively as a I 1
L and B3 = aLL," where Ls = -.+
+M + +LM+1++TME +
Pg
Hg
PLPM Then, the iterates for k = 0, ..., K 1 of ESJ in Algorithm 1
satisfy:
K
El| (x)))))))| < (6 xb)
+DD
(45)
K
K
with x- = inf, ((x) and D is the upper-bound established in Proposition 2
Land is given by
L2M
D -12L2D2 an,) (1,))))) +
- D2 D1 axp)) + 24(4p + 15)
Qu2
12
1i M2M) + 24L M )2p + 6,2 1.2213 M + 3)3
Q
11
(TH, + Lp)2
(1 app)) N-2 D2
(46)
Proof of Theorem
Using the Lipschitzness of function (xxx), we have
(xx+1) (xx) + V (TA) T641 .x) + Iod || - xx 1
2
< - 3(V (Tx), (.((x)) + Lg
2
< - (IV ( v=((x)) /|| (=(())))
+ L2 (1 (T)) ((x)) - v+(()))l
< - { LBB3 [+))) + +g + ((0 v (47)
Let E+1.] = Eu n Hov, v be the expectation over the Gaussian vectors uk......, uxQQ condi- tioned on 2y and 21 Applying the expectation [+1 to eq.
(47) yields
Es + (xx+1) (xx) - - - L.882)
V
(1x) 11 + ( 2 + 1 (xxx)) - (+(x)))]
2
<Φ(Ix) - (-
+
+ LBB
(48)
2
2
24

---

Under review as a conference paper at ICLR 2022
where D represent the upper-bound established for the estimation error in Proposition 2 total expectation over 24 = {21,1.Q,.... 4k.1 ) we have
Now taking
Lo [ + (2 + L
E Ekk
(49)
2
where En Bayy (xx)). Summing up the inequalities eq. (4) for k - ...... yields
E Eo (
K-1
{ + KTD.
LOB3
Σ
Bu,
[<()))]
(50)
2
k=0
Setting 8
"2." denoting by 4- inf & (z), and rearranging eq. (50). we have
16L (2((0)
x-)
Ben
[
K + 3D.
(51)
k=0
Hence, the proof is finished.
G PROOFS FOR STOCHASTIC BILEVEL OPTIMIZATION
Define the following quantities
Y. (zry -2u,n)) (2)
OY
Oy 0zk
IN.- = IN (Ix: uz)
IN
J.
YJ (uu+uay)))) (2)
where uy [R.,j 1...., Q are standard Gaussian vectors and Y is the output of SGD obtained with the minibatches (Su.... Sx-1).
Conditioning on 2ąK and Y and taking expectation over uj yields
YIN to,, Y (#2:8)
Euy JNJ = Eus
V (1)) VN (KNI)) Uj
1
(
V. m d. (ax:)
J.(S)
where ( (xx; is the i-th component of vector Y.N (I S), which is the entry-wise Gaussian smooth approximation of vector y (xx;:) Let E&H E[||zg Y EDr.ur. be the expecta- tion over the Gaussian vectors and the sample minibatch DF conditioned on 2 and Y
G.1 PROOF OF PROPOSITION 3
Proposition (Re-stated). Suppose that Assumptions 7 and hold. Choose the inner-loop 2 stepsize ItØj
Define the constants
a al C, =(1 - -44p) (1-
L
Czy =a (a+ 7/1 -asug) +
Cy =
Cav
7
7Mg
Mg
pgg
o²
1
r' 2=2(72 + b2,
+ (Cay Cy), == 2(72Can +²²,))D2
(52)
Lg
25

---

Under review as a conference paper at ICLR 2022
where y is such that y 2
Lt42, M
Then, we have:
La
X(I M (1)
0444 02 C.
E IN I.1 ICN
+ 1 (g+ - apt1) - (L -2g)²
+
1- C
Proof of Proposition 3 Based on the SGD updates, we have
= Y v (Ik, y S.-1), 1= 1=1,..,, N
Taking the derivatives w.r.t. Ze yields
- V (I Y (S-1) T. J. aV, (,, (mg, y (S--)
-
[kk
aJ V3G (mh, Yt (,-11 J, (0k, 9k))
+ (V= 69 (Zk,
- - J. a (V G ( , ; V V vy (Tx. vh))
- a (Ji-1 J.) V G [k ; + +JJ. (V2g (Ix,, v G (Ik, y
Hence, using the triangle inequality, we have
|| J- 2 (J - V (S:-1) llr
a V2 (,,, S.1) V.Veg veg (Ixz vi) lr
-. VV. (,x,Yk y - V2 (ox vii) I|
where (i) follows from Assumption 0
We then further have
[ - J-I z
appp (1 [i-1 - +o
(8 - (2x, (x) l1
V.V
+ 2 L?
V (xx, 1e-
(rk, vi) 113
1/
+ 2a(1 - -app))
| - J. 11 r V. (ze, Yt ;1-1) - V. V.g (ox mh) 11g
4 + +2a(1 aap)) - |i--1 J. 830 (Ie, 1
Vig (Iz, (k), 11
Mg
P2
+2a² 2 L Hg
| ( S(-1
Vig (E,, vil ll V. (EL, ; V. (Ex, wi) l|
The terms P1, Pg and P3 in the above inequality can be transformed as follows using the Peter-Paul version of Young's inequality.
P1
J J. + - V. ( -
V, V (8) || 7>0
27
|T-11 J. 1. 1:= + 2 |2G (Ix. yt (8-1)) - 11 12:0 7>0
P₂ <
27
P3 5
13G (IK, Yt (--1) V2g (It, v 1/
ア V. V G (IR, ; V (xx, vko ||
Note that the trade-off constant y controls the contraction coefficient (factor in front of J=1[3 Hence, we have
Ji-1
| J..1/²
26

---

Under review as a conference paper at ICLR 2022
(1 apg)) + (1 (pp)) + aL
aل (1 - apw)) |J--1 J.IlE
7
7M4
a² + ta7(1 - -apg) + aa² _ || V (Ix, 11 V u9 (Ik: yi) 1|
pq
1.2
-07 (1- (p1) (1-a + (1-
L
V (Ix. y S V (Tx: vi.) l
+
a
+
Let E. - [|| = [k. _ 1. Conditioning on ky and y and taking expectations yield
Ex-1 Jr J.11 z
Fr J.IIF= C. (kk, ( (8-1) V. un (3 vi) lFr
+ CyEee11 V (IR, Y : - Vig (Ik: vk: (zk. 1|:-
(53)
where Cyp Cay and Cy are defined as follows
C- = (1-(g) (1- (1-011- + a = + oL
M/g ) Cry = (a+ -rsp)) + a Cy =
L
Cay-
1g
pbg
Conditioning on 24k and yx-1, have
B V.
yt Sc-1) V Vug (Tx, wke) 11.
V V 29 (IN, Yx - V V va (Ih, l|
<2E,
2F 1 Vx V (TE, y S. - V2 V.V (FK, 1) l|?
( L
+ + 1
(54)
5
where (i) follows from Lemma and Assumption 2 Similarly we can derive
L2
E, V (IR, ( Vig (zx, 1l 52 52 2p2
-
will
(55)
Combining eq.
53 eq.
and eq. (5)) we obtain
E, J. J. SC J1-1 J.IIF= + 22(2" "any + +p²Co) | 1 -2ill
L
+ 5
+ (Czy (Cay + Cy).
(56)
Unconditioning on Ix and y and taking total expectations of eq. (56) yield
+2 122
B Jr J. sC, J J.112 + 2( + P²C,)
vill"
+2
(Cay Cy)
(
SC,, E J1-1 - 1. + +2( Cry + rp2)))
(T.4
L Hg
2(t-1)
D2 + σ²
11 L.S
+ 2 12 s (Cay + Cyy,
where (i) follows from the analysis of SGD for strongly-conver function. Let I = 2(72 + p²yry 0² Ls
244 (Cay + Ccy) and X 2(72Czg D² Then, we have
L Mgg
2(t-1)
Bl J.I1 SC,BI|Jl-1 -
+
+r.
(57)
I 1gg
Telescoping eq. (57 over from N down to 1 yields
N-1
21
N-1
E BN - J. Bl| J-113- - A Σ
- Mg
C2N CN-1-t
C
(58)
I + Mg
t=
1=
27

---

Under review as a conference paper at ICLR 2022
which, in conjunction with
)
S
aup and 2
Ltdg such that Cy S 1 M
app yields
LtHg/
BN|JN J- 0o
L
+ +1Og
N-1
Σ
(L - pg)
+ 1-
r
pg
(L + Mg) 2²) a(pg))
1-C
1=
<ON m2
1
X(L + Mg) (1 00
r
(59)
[L Hg) aHtg) - (L - ngg
1
C
The proof is then complete.
Lemma 8. Suppose that Assumptions 1 Then, we have
2
and 4 hold. Set the inner-loop stepsize as a = L.P4"
Eg (=())) - V
C.N I
XLI + (1 agg))
r
122
I8M²
+
+
L3d 3)
(L + ) -atta) -(L - (bg)
1 1-C, Cy
+ 222
2
((ITHt)
2N
12 o²
(60)
Mg
LgLS
where the expectation Eg[] is conditioned on Fy and Y.N
Proof of Lemma g
Conditioning and Y.N, we have
Es [ = Vel(ER, Y - Jy (o Y
Recall ()(())) V. v J. y vt). Thus, we have
Eg ( ( (1x)) 2
2 V. Yx V. +11 V (TK, Y - J..'VV/VV. vi)(
<21i|| - vill + J v fI Y J. V off X2")
+ T. Vy](y( Y - J. v (T 11
La (3rk) YE Vys (i)
22r.2|122" - vil 4M2 [ J.1.:
V |
<21'I|YE - u SM Jr, 1 - JN 12 + &M2 JN - J.1l2 + +4 |
will
14g
where (i) applies Lemma 4 and Assumption
3
Taking expectation of the above inequality yields
1 μ3
E E. (.(x))
-z)))) <21²
1
will + SM'E|JJN J. E Eliy
+ 8M2 JN 1|
(61)
Using the fact that Y.N (xx; has Lipschitz gradient (Proposition and Lemma the last term the right hand side of eq. (61 can be directly upper-bounded as
12 < 13 + 3)3
1| - -
()); VY, S)
(62)
2
where Ly is the Lipschitz constant of the Jacobian IN (and also of its rows VY.N (xx: as defined in Proposition]] [ Combining eq. (61). eq. (62). Proposition 3 and SGD analysis (as in eq. (60) in Ji
al. (202T)) yields
E E2 (x(x)) (x+x))
28

---

Under review as a conference paper at ICLR 2022
L
ALL + (1 aau -N
r
11²
<8M²
C
d 3)9)
L
(L + 12) -444
(L-
Cy
1
(4- 1g
2N
o²
+ 212
1
1 Ugg
D32 +
(63)
/
Ng L LS
This finishes the proof.
G.2 PROOF OF PROPOSITION
Proposition 4 (Re-stated), Suppose that Assumptions 7 and AadAhooid Ahold. Set the inner-loop stepsize
2 Then, we have: as a LtHe"
E B -V0()))-V(rT))) - B
(64)
1p+ where A = 8M² ((++
D,)
+ D,, 23
2
22 (1+
4N22 Dj
Q pddL3Pa(p)+ 2222
H P3L Ps(p)
Q
Dj
and Bi respresents the upper bound established in Lemma g
Proof of Proposition We have, conditioning on 2y and Y
Ex [()()))) v (xx) = -x
([x)) = - E. =|
+
-
(65)
Our next step is to boper-bound the first term in eq. 665
Ex -E (xxx))
<2E V F( YE Dr) Vef(EE, Y )|
+ 2 JN (Zk, Y Dr) - J, V n (IA. Y )1
M
4E 1y (akk Y ; ) |1 -Juli²
2 Dj
+ 4EA || J. | V V,F(aa F(rk Y (D) V. Y/")
M
1 1
.M2
s + 4M2 1+ D,
E& 1x Jnll + *D, ||
(66)
D
where the last steps follows from Lemmal I
Next, we upper-bound the term E IN Tnll?
Eel| - J1
Ex JN.J - T2llz j E (1......
Q
ER
JNN.J 1 - 11T
Q
1
Y,N (xx + pu Y," (xxxx) Uj 1| - -IVVY (67)
Ex
Q
f
Recall that for function h with I.-Lipschitz gradient, we have
h(x + pm) h(x) ulf s 4(p 4)/l
3 22 "u² 5)3.
E,
(68)
Then, applying eq. (68) to function Y. (=;)) yields
Y.N (xx + p) Y.N (xx:
Eg,
Ugy
P
5 4(p + 4)l (T) 2
11 5)8.
3 マ
29

---

Under review as a conference paper at ICLR 2022
Hence, eq. 66 becomes
4p + 15 d
(1:; S)l2 + 3p²dL3 (p+
Ex IIIN - Jull -
2Q
2(4p + 15)
VY, (IZ; S| 2 1d Pa(D)
Q
Q
24pp + 15)
+ pd
JN
Pa(p),
(69)
Q
Q
where P4 is a polynomial of degree p. Combining eq.
(66
eq. 669
and Lemma yields
BxllV((<xx- E& (xx))]
22
M
4M
(1
1
2(4p + 15)
JJN 1| + nddLz
Pa(p)
Dj
Dj
Q
Q
4M2
(2]]Jx 1'= +122d
Dj
$ 8M2 ( (1+ 1
4p + 15
) +2 M
D
1
1.2
+
Dj
Q
D,
Dj
(1 +
1 )
4M2
pd(l3
4M2
p
Q
Dj
A
(70)
where = SM² ((++
1p+ +
1
12
2 py
+ (1+ b
AM2
q 1e dL)))(())
+
pg
D, po
4322
Taking total expectations of both eq. 665 and eq. 700 and combining, we have
E ((()))))) (*xx))
= 1 (Fx) E. (xx))l Vv(xx) - (x))))
<A + B
(71)
where B1 representing the upper-bound established eq.
63 The proof is then complete.
O
G.3 PROOF OF THBOREM 2
Theorem 2 (Re-stated). Suppose that Assumptions 2 and and Aoold. Set the inner- and outer-loop stepsizes respectivelly as a L+ug and 1 ILo' where = max ( Tog, I.g) and the constant La is defined as in Theorem Then, the iterates Jk. k = O, .... K - 1 of the ESJ-S algorithm satisfy (x()( - Lo
EK-J B<<=B\/V(<>>> where A and B, have the following
3B, + A,
(72)
K
A =8
((1+)),
1 ).),))
4p +
+
L
M
+(1+
11)
1M2 pddLiyP.())
Q
Dj
12
D₁
Dj
Q
+ 1M2
D, 1² (ddL)Pa(pp)
C 1?
X(L + 21 - (1-- ax1))
r
"² L3d(p 1bd(p+33) 3)
B, 2
+
(L+ (2 - apt)) - (L- (L-Ha}²
1-
1,2
2N
+ 2L
1+2
D2
pgg
LS
and the constants r, X, and C.y are defined in Proposition 3
Proof of Theorem
Using the Lipschitzness of function ((Ix),, we have
Le
@ (Ix+1) < + (Ix),, IK+1 - xx) + 2
| - 2%
30

---

Under review as a conference paper at ICLR 2022
v)(x)) Le 32 - +
Le
< (zv) (mx))) + Lo3 (|V0(z))]*+ v
Hence, taking expectation over the above inequality yields
[D((++1) <Ed(xa) (E (x(())) + L [
+
(73)
Also, we have
-B(VØ(xx), (zxxx) = E(V4(x2))) - -
EE
V E. (2xx)) v=((2)) - EE (++x+)]1
( - v (zz=)
Ex ( [ {
which, in conjunction with eq.
73 yields
Eg (TTK) - ( (5 I.0
(0((+)1) SBD (3r)
E V (xx)
LaB (zx) [=*+))]
(74)
Setting 3 TLe
and using the bounds established in Lemma 8 and Proposition we have
E$ EΦ( <E 2 B. _ [ |
8 (B1 + A).
4
Summing up the above inequality over k from k = 0 K 11-1yin yields
8
8 K +
E$(xK) (ED(20) + KB
E
Rearranging the above inequality yields
K
4 (E - [T((K))
Σ E V
+ +3B, +
K
3K
16 (4(20)
4-) La
+3B + A.
K
The proof is then complete.
31


Reviews:
[{"review_id": "jG4umOlPbhr", "reviewer": "Reviewer_Papm", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "This paper tackles a very relevant problem. to the machine learning community. It introduces two interesting new algorithms, executing on a simple but elegant idea, which is to limit the black box to its bare minimum, and leverage whatever problem structure of the hypergradient possible.\nThe writing is clear and enjoyable to read, and the problem is well motivated.\n\nA number of significant issues remain, however.\n1. The paper references deep neural networks many times, but crucially relies on the objective being strongly convex with respect to y, which if I understand correctly are the parameters of the model. Deep learning optimisation landscapes are notoriously not even convex, let alone strongly convex. If y is actually only the final layer of a neural network, this needs to be explained right from the start to put the exposition in context. Further, this doesn't seem compatible with the claim that all the 12M resnet12 parameters are learned in the experiment. Some clarity there would go a long way.\n\n2. There are some overclaims, notable regarding 'large' neural networks,  when only fairly small ones are used in the experiments (that is in itself not an issue; claiming that they are large scale when they are not is a problem though).\n\n3. Another overclaim is that the proposed algorithms are 'efficient'. Only requiring gradient computation does not make an algorithm efficient, particularly when it requires NQ *full* gradient steps per outer step. A comprehensive cost analysis with AID and ITD which use second-order information is necessary here, to see when N and Q become too big for this approach to yield gains.\n\n4. The theoretical results, while interesting, are not properly put into context. First, the upper bounds contain a massive additive constant, which isn't smaller than any epsilon. So the results are only valid when this constant is smaller than epsilon, which should be made clear. Second, the settings needed to make the constant small are drastic, meaning Q needs to be in p/epsilon, which if p is the amount of weights of a neural network means using incredibly big Qs, and in turn massive amounts of full gradient steps.\nNo context is given here on the limitations of the theoretical results, which appear to be all but unusable from a practical perspective.\n\nOn a related note, in proposition 3 we need at least some context on how big the constants are.\n\n5.  The experiments are all using small scale neural networks, so claiming that the approach succeeds on large neural nets is misleading.\nFinally, there seems to be a discrepancy between the analysis on strongly convex objectives, and the experiments using non-strongly convex objectives.\n\nQuestion: is it necessary to use the same batch sampling path for every mu? Has this been ablated?\n\nAll told, this paper introduces promising new ideas, but at this stage the execution is a bit lacking.  At this stage I do not recommend it for publication at ICLR, but I encourage the authors to take this feedback to heart and improve the paper so that it can have all the more impact on the community.", "summary_of_the_paper": "The paper introduce two new algorithms to tackle bi-level optimisation, which has wide-ranging application in machine learning.\nCrucially, these two algorithms forego second order gradient computation completely, favoring instead an evolution strategies approach. \nWhile 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.\nTheoretical 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.", "summary_of_the_review": "Interesting new methods to tackle a very relevant problem.\nUnfortunately a few too many issues remain in the presentation of the paper to make it ready for publication this time around, but it absolutely can be with a bit of additional work.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "2: Several of the paper’s claims are incorrect or not well-supported.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "X3ZUyKQIH5", "reviewer": "Reviewer_Ws7X", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "**Things I like**: The paper proposes an alternative approach to bypass the need of computing the costly second-order information in bilevel optimization and provides theoretical guarantees for the proposed algorithms. The problem is important, and the algorithms are well-motivated.\n\n**Things can be improved**:\n\n**Limited technical novelty**: While this paper uses ES in a more careful manner in the bilevel optimization, it should be noted that the idea of using ES mainly follows from [Gu et al 2021] and the convergence analysis of the bilevel part mainly follows from [Ji et al 2021]. \n\nSpecifically, compared with [Gu et al 2021], the proposed algorithm only uses ES on Jacobian in (3) not the entire bilevel gradient. On the theoretical side, it is known in optimization literature that applying an ES method will only introduce a small bias that is proportional to \\mu, and the remaining of the analysis follows directly from the bilevel optimization [Ji et al 2021]. Also see a couple of earlier and important works such as\n\n[1] Ghadimi, S. and Lan, G., 2013. Stochastic first-and zeroth-order methods for nonconvex stochastic programming. SIAM Journal on Optimization, 23(4), pp.2341-2368.\n\nBuilding upon previous results in bilevel optimization and ES, the theoretical contribution of this work amounts to establish the estimation error of the proposed gradient estimator, which looks not as impressive as claimed in the introduction. \n\n**High actual complexity**: While the goal is to reduce the complexity of bilevel optimization, the proposed algorithm still incurs large computation overhead. By carefully looking at Theorems 1 and 2, to guarantee convergence, the theory requires the number of queries Q and the batch size S to grow with $\\epsilon^-1$, which makes the theory less favorable to practitioners. Such high complexity seems to be hidden in the simulations by choosing small Q and S, which makes the comparison with baselines less convincing. \n\n**Limited empirical verification**: One should provide a thorough empirical studies on the proposed algorithms and the existing ones. However, most of experiments are carried out on simple problems with small networks. And the running time (in second, minute, or hour) is not clearly marked. I have several questions:\n\ni) The baselines on MNIST only get a test accuracy of 93, which is too small. Have the authors tuned the parameters well? To my knowledge, it is not hard to achieve 97 or 98 test accuracies on reasonable networks.\n\nii) It looks like that the authors do not scale up the learning rate with e.g., 0.1, after several epochs. This learning rate trick is useful in practice, and some algorithm benefits significantly from this trick. \n\niii) ResNet12 is not big enough to justify the claim on scalability. \n\nAs a matter of fact, only a subset of experiments of [Gu et al 2021] has been tested with ESJ using synthetic data in this paper. This seems implausible to claim its practical advantage to HOZOG [Gu et al 2021]. \n", "summary_of_the_paper": "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. ", "summary_of_the_review": "Overall, first-order bilevel optimization is a worthy direction to pursue, but efforts on both technical and empirical novelty are needed to make this paper competitive.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "UmNoHSdrg7S", "reviewer": "Reviewer_bgEq", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The paper is generally well written. The core part of this work is the ES Jacobian algorithm, which is specialized in solving bilevel optimization problems. The assumptions made for analyzing the bilevel algorithms are standard, where the main theoretical contributions are Prop.1 and 2 that serve as the pillars of showing the convergence of the bilevel algorithms. I found that method is novel and interesting, but still, have the following concerns regarding the contribution and technical correctness of this work:\n\n1) The proposed ES algorithm is quite similar to the zeroth-order methods. What are the unique differences here? Comparing the first-order and zeroth-order methods, the theoretically guaranteed convergence rate of the zeroth-order ones is slower than the first-order algorithm with a dependency of problem dimension in general. Here, the authors try to argue that the proposed ES is computationally efficient of using the black-box type of methods, which contradicts the common intuition that ES would be slower w.r.t. rate. Please justify.\n\n2) The proposed algorithm is basically a double loop one, since $N$ is a function of $\\epsilon$. This again contradicts the intuition that single loop algorithms would be more efficient. To my understanding, the current algorithm cannot be improved as a single loop one.\n\n3) Regarding the convergence rate results, can authors compare the obtained one with existing works, e.g., single- or two-timescale bilevel algorithms? my understanding is that the proposed one for both the deterministic and stochastic cases needs more iterations resulting in higher sample complexity.\n\n4) The required step sizes for the ESJ-S are a little wired. Why do $\\alpha,\\beta$ be constants? I think the step sizes are either decreasing sequences or depending on $\\epsilon$. Please justify.\n\n5) In section 4.3, why not consider MAML and ANIL as baselines?  Also, it seems that all the compared methods have not been converged yet during the training.  Also, it would be great that more statistical measures can be reported such as variance.\n\n\n", "summary_of_the_paper": "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.", "summary_of_the_review": "Overall, I believe that it is an interesting work, but the proposed is only applicable for solving a class of problems given a special structure of hypergradient. Even the authors justify the computational efficiency of the proposed algorithm numerically w.r.t. running time, the theoretical justification of the efficiency of the iteration and sample complexities should be further addressed.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "ThYzZ_qehS", "reviewer": "Reviewer_Tjy4", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Pros:\n* The idea of using ES based method to approximate the response Jacobian matrix is interesting.\n* Overall, the paper is easy to follow. In particular, the main contribution section makes it easier to highlight the important contribution made in the paper and the related work is written well and makes the algorithm well motivated. \n\nConcerns:\n* I believe that the paper misses several important related works. There have been earlier works on using ES-based methods for bilevel optimization (e.g. [1]) that share a similar idea of approximating the Jacobian of the response function. Moreover, STNs ([2], [3]) use the ES gradient approximation method to approximate the response function (or response Jacobian) in an amortized manner. I believe that these should be discussed in the related work section and contribution section of the paper, and should be compared against in the experiment section since the core idea is very similar. \n* The key concern about the paper is the lack of rigorous experimentation to show that the proposed method is superior to other bilevel optimization approaches. The paper majorly focuses on the hyper-representation experiments and compares the proposed method to only three baselines (e.g. HOZOG, AID-CG, AID-FP). I believe that the paper should experiment with the algorithm in more diverse settings (e.g. hyperparameter optimization) and compare it to other gradient-based methods (e.g. [2], [3], [5] and various ITD based methods). At the moment, I am not convinced that the proposed algorithm is superior to other bilevel optimizers (both in terms of accuracy and computation it required).\n* In the first experiment, the authors do not properly explain why there is an improvement in the outer data accuracy. Does this mean that the proposed method approximates the response Jacobian more accurately than CG? \n* What are additional hyperparameters introduced by ESJ-S? I believe that they are Q (number of Gaussian samples) and N. It would be good to have an ablation study of how Q affects the quality of the algorithm. \n* In the first experiment, the authors mention that they follow the same experimental procedure as [4], but I couldn’t find the paper. I believe that there are typos in the reference. \n* Some mathematical notations are not properly defined throughout the paper and make the paper hard to follow in terms of the notation. For example, in Algorithm 1, the authors did not define variables such as K, N, and Q.\n\nSome minor comments: \n* In section “related work”, the paragraph header “Bilevel optimization” is missing “.” at the end.\n* In algorithm 1, it would be helpful to identify variables such as K, Q\n* In section 3.2, “Differently from …” → “Different from …”\n* I believe that some information in the Appendix is outdated. For example, in Appendix B, the authors say “paper demonstrates the superior performance of the proposed ES-based bilevel optimizer in meta-learning and hyperparameter optimization”. However, there was no experiment on hyperparameter optimization.\n\n[1] Sinha, Ankur, Pekka Malo, and Kalyanmoy Deb. \"An improved bilevel evolutionary algorithm based on quadratic approximations.\" 2014 IEEE congress on evolutionary computation (CEC). IEEE, 2014.\n\n[2] MacKay, Matthew, et al. \"Self-tuning networks: Bilevel optimization of hyperparameters using structured best-response functions.\" arXiv preprint arXiv:1903.03088 (2019).\n\n[3] Bae, Juhan, and Roger Grosse. \"Delta-STN: Efficient Bilevel Optimization for Neural Networks using Structured Response Jacobians.\" arXiv preprint arXiv:2010.13514 (2020).\n\n[4] Riccardo Grazzi, Luca Franceschi, Massimiliano Pontil, and Saverio Salzo. Optimizing millions of hyperparameters by implicit differentiation. International Conference on Machine Learning (ICML)), 2020.\n\n[5] Pedregosa, Fabian. \"Hyperparameter optimization with approximate gradient.\" International conference on machine learning. PMLR, 2016.", "summary_of_the_paper": "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. ", "summary_of_the_review": "The paper is easy to follow and I believe that the contents of the paper are technically correct. However, the paper doesn't discuss several important related works and the empirical analysis is not rigorous: the authors conduct the empirical investigation in a limited setting and don't compare the proposed method with competitive baselines. It is still unclear if the proposed algorithm is superior to the previous approaches. I believe that there should more in-depth empirical analysis to make the proposed method more convincing. \n", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "1: The contributions are neither significant nor novel."}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper presents an algorithm for approximating the hypergradient for bilevel optimization using a trick based on evolution strategies. It seems like an interesting approach, somewhat reminiscent of STNs, so it's interesting to see experiments with it.

I have a big concern about the proposed justification of the method, namely that each iteration is more efficient than methods based on HVPs. The authors claim that because they only require gradient computations and not HVPs, their method is more efficient. However, as various reviewers point out, the proposed method requires numerous inner optimization runs. By comparison, a method based on unrolled backprop simply requires a single inner run, followed by backprop on the trajectory; hence it should be about as expensive as 2-3 inner optimizations (or less if it is truncated BPTT). Similarly, each HVP has a small multiple of the cost of an inner optimization step, so methods based on HVPs ought to be cheaper unless they're doing quite a lot of HVPs.

It's conceivable the proposed method could be more efficient than AID, etc. if each hypergradient estimate is more accurate. However, this isn't shown, and it would seem surprising for an ES-based approximation to be more accurate than the exact gradient.

The authors claim in the rebuttal that the efficiency claims aren't based on the theoretical analysis, but rather on the experiments (which use Q=1); however, Section 3.2 still finishes with the conclusion that ESJ is more efficient, which seems problematic.

I encourage the authors to formulate their theoretical claims more carefully and to consider the reviewers' other feedback, and I think this could make an interesting submission for the next cycle.
Decision: Reject