Projection-Free Methods for Stochastic Simple Bilevel Optimization with Convex Lower-level Problem
Authors: Jincheng Cao
Ruichen Jiang
Nazanin Abolfazli
Erfan Yazdandoost Hamedani
Aryan Mokhtari
Venue: NeurIPS 2023 poster (2023)
Keywords: Bilevel optimization
stochastic optimization

Abstract: In this paper, we study a class of stochastic bilevel optimization problems, also known as stochastic simple bilevel optimization, where we minimize a smooth stochastic objective function over the optimal solution set of another stochastic convex optimization problem. We introduce novel stochastic bilevel optimization methods that locally approximate the solution set of the lower-level problem via a stochastic cutting plane, and then run a conditional gradient update with variance reduction techniques to control the error induced by using stochastic gradients. For the case that the upper-level function is convex, our method requires $\mathcal{O}(\max\\{1/\epsilon_f^{2},1/\epsilon_g^{2}\\}) $ stochastic oracle queries to obtain a solution that is $\epsilon_f$-optimal for the upper-level and $\epsilon_g$-optimal for the lower-level. This guarantee improves the previous best-known complexity of $\mathcal{O}(\max\\{1/\epsilon_f^{4},1/\epsilon_g^{4}\\})$. Moreover, for the case that the upper-level function is non-convex, our method requires at most $\mathcal{O}(\max\\{1/\epsilon_f^{3},1/\epsilon_g^{3}\\}) $ stochastic oracle queries to find an $(\epsilon_f, \epsilon_g)$-stationary point. In the finite-sum setting, we show that the number of stochastic oracle calls required by our method are  $\mathcal{O}(\sqrt{n}/\epsilon)$ and $\mathcal{O}(\sqrt{n}/\epsilon^{2})$ for the convex and non-convex settings, respectively, where $\epsilon=\min \\{\epsilon_f,\epsilon_g\\}$.

Projection-Free Methods for Stochastic Simple Bilevel Optimization with Convex Lower-level Problem
Jincheng Cao ECE Department UT Austin
Ruichen Jiang ECE Department UT Austin rjianguutexas edu
Nazanin Abolfazli
SIE Department The University of Arizona
j inchengcao@utexaa.com edu
mazaninabolfazz1@arrr
Erfan Yazdandoost Hamedani
Aryan Mokhtari
SIE Department
ECE Department
The University of Arizona
UT Austin
erfany@arizone wwa.yaarrizoa....m edu
mokhtari@austin utexas.eeo... edu
Abstract
In this paper, study class of stochastic bilevel optimization problems, also known as stochastic simple bilevel optimization, where we minimize smooth stochastic objective function the optimal solution set of another stochastic convex optimization problem. We introduce novel stochastic bilevel optimization methods that locally approximate the solution set of the lower-level problem via a stochastic cutting plane, and then run a conditional gradient update with variance reduction techniques control the error induced by using stochastic gradients. For the case that the upper-level function is convex, our method requires Ommax(1/e2 1/21) stochastic oracle queries to obtain a solution that is cy- optimal for the upper-level and Eg cp-optimal for the lower-level. This guarantee improves the previous best-known complexity O(max (armax(1//5) (1/4]- 1/( Moreover, for the case that the upper-level function is non-convex, our method requires most (mmax(1/e} }) stochastic oracle queries find an (e,, sttationny point. In the finite-sum setting, show that the number of stochastic oracle calls required by method are ((Vn/c) and ((vn/e2) for the convex and non-convex settings, respectively, where cig(
Introduction
An important class of bilevel optimization problems is simple bilevel optimization in which we aim to minimize upper-level objective function over the solution set of lower-level problem [ Recently this class problems has attracted great attention in machine learning society due to their applications in continual learning [5). syper-parameter optimization [6. meta-learning | 9. and reinforcement learning [10, [1]. Motivated by large-scale learning problems, in this paper, we are particularly interested in the stochastic variant of the simple bilevel optimization where the upper and lower-level objective functions are the expectations of some random functions with unknown distributions and accessible only through their samples. Hence, the computation of the objective function values or their gradients is not computationally tractable. Specifically, we focus on the stochastic simple bilevel problem defined
min
f(x) =xx)=xff(x.. Exf(x,
E arg g(z) E[ğ E)
(I)
ZEZ
37th Conference on Neural Information Processing Systems (NeurIPS 2023).

---

Table 1: Results on stochastic simple bilevel optimization. The abbreviations "SC", "C", and "NC" stand for "strongly convex", "convex" and "non-convex" respectively. Note that e = m
References
Type
Upper level
Lower level
Convergence
Sample Complexity
Objective f
Objective g Feusible
Upper level
Lower
aR-IP-SeG [6
Stochastic
C, Lipschitz
C. Lipschitz
Closed
O(max(1/c3.11c2)))
0(1/e")
Algorithml
Stnchsttic
smooth
smooth
Compact
((max(1/(33/133)
0(1/())
Algorithml
Stochastic
NC, smooth
smooth
Compact
0(max(1/e3:1/633)
0(1/e")
Algorithm2
Finite-sum
C, smooth
C, smooth
Compact
0(max(1/es.1/xa))
O(/6/e)
Algorithm2
Finite-sum
NC, smooth
smooth
Compact
0(max(1/3,.1/33)
Ö(yā/e")
where is compact and convex and f. R - R are continuously differentiable functions on an open set containing 2, and 0 and & are some independent random variables drawn from some possibly unknown probability distributions. As a result, the functions f. R - R are also continuously differentiable functions on Z. We assume that y is convex but necessarily strongly convex, and hence the solution set of the lower-level problem in ①T is in general not a singleton. We also study the finite sum version of the above problem where both functions can be written as the average of component functions, i.e., f(xx = (L/n) Ei f(x 0.) and g(z) (I/n)>. E(z
The main challenge in solving problem (I). which is inherited from its deterministic variant, is the absence of access to the feasible set, i.e., the lower-level solution set. This issue eliminates the or projection-frs free methods. There have been some efforts to overcome this issue in the deterministic setting (where access I and g and their gradients is possible), including [5) however, there is little done on the stochastic setting described above. In fact, the only work that addresses the stochastic problem in ① is [16). where the authors present stochastic extra gradient algorithm and show that it requires ((L/E)) and ())/e)) queries to the stochastic gradient of the upper-level and lower-level function, respectively, to obtain a solution that is é opttmmal for the upper-level and optomma for the lower-level. We improve these bounds and also extend our results to nonconvex settings.
possibility of using any projection-basee
an iterative regularization-bane
Contributions. In this paper, we present novel projection free stochastic bilevel optimization methods with tight mon-asymptotic guarantees for both upper and lower-level problems. At each iteration, the algorithms use small number of samples to build unbiased and low variance estimates and construct a cutting plane to locally approximate the solution set of the lower-level problem and then combine it with a Frank Soln-Wooffe-tpee update on the upper-level objective. Our methods require careful construction of the cutting plane so that with high probability it contains the solution set of the lower-level problem, which is obtained by selecting proper function and gradient estimators to achieve the obtained optimal convergence guarantees. Next, summarize our main theoretical results for the proposed Stochastic Bilevel Conditional Gradient methods for Infinite and Finite sample settings denoted by SBCGI and SBCGF, respectively.
(Stochastic setting) We show that SBCGI (Algorithm D. in the convex setting, finds a solution x that satisfies f(x) f < e and g(x) - I eg with probability 1- I-6 within Olog (gllog/ddg)) (d/ōc)/e²) stochastic oracle queries, where is the optimal value of the lower-level problem. Moreover, in the non-convex setting, it finds x satisfying (xx < ef and g()) g" S 6g with probability 1 5 within stochastic oracle queries, where (x) is the Frank-Woff (FW) gap.
minfer, f* is the optimal value of problem (1) and g
(Finite-sum setting) We show that SBCGF (Algorithm 2 in the convex setting, finds & that satis- fies f(x) f S C and g(x) g" I Cg with probability 1 - within O ((vT(log((/8E) 3/2 stochastic oracle queries, where n is the number of samples of finite-sum problem. Moreover, in the nonconvex setting, it finds & that satisfies (xx) S E and g(x)) 9" S Eg with probability -- å within O(ynn (o)()))))/e) stochastic oracle queries.

---

1.1 Related work
General stochastic bilevel. In a general format of stochastic bilevel problems, the upper-level function f also depends on an extra variable y E IR" which also affects the lower-level objective,
min XER YERP
yxx,y) =
S.t. arg min g(z. y) Elg(z,y, E)]-
(2)
zE2
There have been several works including 10 21 on solving the general stochastic bilevel problem (2) However, they only focus on the setting where the lower-level problem is strongly
convex, i.e., g(z y) strongly convex with respect to for any value of y. In fact. (2) with a convex lower-level problem is known to be NP-hard [22] Hence, the results of these works are not directly
comparable with our work as we focus on simpler setting, but our assumption on the lower-level objective function is weaker and it only requires the function to be convex.
Deterministic simple bilevel. There have been some recent results on con-asymptotic guarantees for the deterministic variant of problem (T) The BiG-SAM algorithm was presented in [1], and it was shown that its lower-level objective error converges to zero a rate of (1/tt), while the upper-level error asymptotically converges to zero. In [15], the authors achieved the first non-asymptotic rate for both upper- and lower-level problems by introducing an iterative regularization-batt which achieves an (g)r,())ppiimal solution after O (max (1/e, 1/c²)) iterations. In [1],, the authors proposed a projection-fiee free method for deterministic simple bilevel problems that has a complexity of (mmxx [1/e,, 1/2 for convex upper-level and complexity of (mmx[[1/e) 1/(E)e,))) for non-convex upper-level. Moreover, in 23 the authors presented a switching gradient method to solve simple bilevel problems with convex smooth functions for both upper- and lower-level problems with complexity O (1/e). However, all the above results limited the deterministic setting.
method
General bilevel without lower-level strong convexity. Recently, there are several recent works on general bilevel optimization problems without lower-level strong convexity including 224 23
However, they either have a weaker theoretical results like asymptotic convergence rate in 22 or have some additional assumptions. Specifically, in [25], the authors reformulated the problem as a constrained optimization problem and further assumes such problem to be a convex program. Moreover, in 226 (2), the authors in both papers assumed that the lower-level objective satisfies the PL inequality, while we assumed that the lower-level objective is convex. In (23). the authors used a looser convergence criterion that only guarantees convergence to Goldstein stationary point. Since these works consider a more general class of problems, we argue that their theoretical results when applied to our setting necessarily weaker.
Preliminaries
2.1 Motivating examples
Example 1: Over--aaannn pheeeerraeteri regression. A general form of problem (I) is when the lower-level problem represents training loss and the upper-level represents test loss. The goal is minimize the test loss by selecting one of the optimal solutions for the training loss [13. An instance of that is the constrained regression problem, where we intend to find an optimal parameter vector & Rd that minimizes the loss fer ((9) over the training dataset D To represent some prior knowledge, we usually constrain ,3 to be in some subsets C Rd e.g., Z = {33 || S }} for some to induce sparsity. To handle multiple global minima, we adopt the over- paraamtteeete approach, where the number of samples is less than the parameters. Although achieving one of these global minima is possible, not all optimal solutions perform equally other datasets. Hence, introduce an upper-level objective: the loss on a validation set D This helps select a training loss optimizer that performs well both training and validation sets. It leads the following bilevel problem:
min BERd
[(3) A (3)
s.t. , E argmin gzz) A (z)
zEZ
In this case, both the upper- and lower-level losses are smooth and convex if & is smooth and convex.
Example 2: Dictionary learning. Problem ① also appears in lifelong learning, where the learner takes a series of tasks sequentially and tries to accumulate knowledge from past tasks to improve performance in new tasks. Here we focus continual dictionary learning. The aim of dictionary
3

---

learning is to obtain compact representation of the input data. Let A = fa...... a, E R denote dataset points. We seek dictionary [d..... dp) Rmxp such that all data points a, can be represented by linear combination of basis vectors in D which can be cast as [1:
min min DEREXP XERPX
Σ EN
aaa - Dx,
s.t. ||dll22 S 1, = .: Ix,l S S0ii (4)
2n
Moreover, we denote X = X [x....., x,, R as the coefficient matrix. In practice, data points usually arrive sequentially and the representation evolves gradually. Hence, the dictionary must be updated sequentially as well. Assume that already have learned a dictionary D E TR= and the corresponding coefficient matrix X E |Rxx= for the dataset A. As new dataset A' : a,,) arrives, we intend to enrich our dictionary by learning D E Rm > p) and the coefficient matrix X E Roxn' for the new dataset while maintaining good performance D on the old dataset as well as the learned coefficient matrix X. This leads to the following stochastic bilevel problem:
min min X DeRmaxy XEReXH'
n': s.t. | S k I....,n' D grgming(D)) | S1
(5)
where f(D, X
Ek= |a
Dx &xxl represents the average reconstruction error on the new |a. D represents the error the old dataset A. Note
2nd
k==
dataset A', and g(D) that we denote &; as the prolonged vector Rº by appending zeros the end. In problem (5) the upper-level objective is non-convex, while the lower-level loss is convex with multiple minima.
Σ
2.2 Assumptions and definitions
Next, we formally state the assumptions required in this work. Assumption 2.1. 2 is convex and compact with diameter D, i.e., Vx,y E 2, we have ||x 1x-yl s D. Assumption 2.2. The upper-level stochastic function J satisfies the following conditions:
(i) V is Lipschitz with constant L,g, ie., Vx.y E 2, [V 0) Vf(y, 0) < T -yl yll-
(ii) The stochastic gradients noise is sub-Gaussian, Elexp[llV 0) fxx)))///3}1 exp(1).
Assumption 2.3. The lower-level stochastic function ğ satisfies the following conditions:
(i) is convex and Vg is Ly-Lipschitz, i.e., Vx,y E 2, VC, [V0(x,) Voly, E)ll - Lg [|||x - yll-
(ii) The stochastic gradients noise is nub-Gaussian, E[cxp(]]
12(x) S exp
(iii) The stochastic functions noise is sub-Gaussian, Ejexp(lo(x, - g(x)) /of - exp(1).
Remark 2.1. Assumptions 2.1 and and 2.3 imply that ğ is Lipschitz continuous on an open set containing Z with some constant L,, i.e., for all x,y E Z, we have [x(x) g(y)l L yl-
In the paper, we denote g" - mingez (z) and X A arg min gzz) as the optimal value and the optimal solution set the lower-level problem, respectively. Note that by Assumption [.3] the set X is nonempty, convex, and compact, but typically not singleton minima on Z. Furthermore, we denote f" and x* as the optimal value and an optimal solution of problem (T. which are assured to exist since f is continuous and XG is compact.
potentially has multiple
Definition 2.1. When f is convex, point E Z is (cr,sy)-optima if f(x) - -5 S Eg and g(x) -g" S Eg. When is non-convex, I E Z is (Ez, Eg)ooptimal if (Xx) Ef and g(x) g" S Eg.
where ((x) the FW gap 13 33) defined as ((x) A maxxex; (VVI(X) s)}.
Algorithms
Conditional gradient for simple bilevel optimization. A variant of the conditional gradient (CG) method for solving bilevel problems has been introduced in [2] which uses cutting plane idea 34

---

to approximate the solution set of the lower-level problem denoted by X" More precisely, if one has access to Ky is possible to run the FW update with stepsize "/. as
Xx+1 = (1 - 7x)x; + 7t5t:
where
S; = arg m ( x
(6)
However, the set X is not explicitly given and the above method is not implementable. In 12. the authors suggested the use of the following set: X, = {s E 2 (Vg(x,), xx.) S g ((xo) ) g (x)) instead of the set X in the FW update given in (6). Note that Xo is selected in a way that g(xo) 9 is smaller than cg/2, and such a point can be efficiently computed. A crucial property of the above set is that it always contains the solution set of the lower-level problem denoted by X This can be Z and
easily verified by the fact that for any x in X2 we have x E
(Vg(x)).x; xx) S g(xg)) g(xx) - g(xo) - g(x)).
(7)
where the second inequality holds as g(xo) 2 g(x²) As shown in [12], this condition is sufficient to show that if one follows the update in 66) with Xr instead of X the iterates will converge to the optimal solution. However, this framework is not applicable to the stochastic setting as we cannot
access the functions or their gradients. Next, we present our main idea to address this delicate issue.
Random set for the subproblem. A natural idea to address stochasticity is to replace all gradients and functions with their stochastic estimators for both the subproblem in i... i.c., V fx,,. as well as the construction of the cutting plane X,, i.e., g(x)), and Vg(x.). However, this simple idea fails since the set Xg may no longer contain the solution set X². More precisely, if ÿt and Vg, are unbiased estimators of g (xx)) and (xx),, respectively, for the following approximation set
X = {s 2 (V 8- - xx) S g(xo) gc3
(8)
we can not argue that it contains X as the second inequality in (7) does not hold, i.e. (V)(x.), X x²) S ((x,)) (xx) L 9(xo) - (xx). In the appendix, we numerically illustrate this point.
To address this issue, tune the cutting plane by only moving it but not rotating it, i.e., adding another term to tolerate the noise from stochastic estimates. We introduce the stochastic cutting plane
i = {s E Z (Vg.. s - xx) S (xx)) - 9e + K,},
where Vg, and Åe are gradient and function value estimators, respectively, that formally define later. In the above expression, the addition of the term K₁, which is a sequence of constants converging to zero as t + oo, allows us to ensure that with high probability the random set & contains all optimal solutions of the lover-level problem. Choosing suitable values for the sequence K, is crucial task. If we select a large value for K, then the probability of & containing X goes up the price allowing points with larger values in the set. As result, once we perform an update similar to the one (6). the lower-level function value could increase significantly. On the other hand, selecting small values for K, would allow us to show that the lower level objective function is not growing, while the probability of & containing X becomes smaller which could even lead to a case that the set becomes empty and the bilevel problem becomes infeasible.
Remark 3.1. How to compute g(xo)? In the finite sum setting, we can accurately compute g(xo). and the additional cost of function evaluations will be dominated by the overall complexity. In the stochastic setting, we could use a large batch of samples to compute g(xo) with high precision at the beginning of the process. This additional operation will not affect the overall sample complexity of the proposed method, as the additional cost is negligible compared to the overall sample complexity. Specifically, we need to take a batch size of b = O(E-2) to estimate (xx)). Using the Hoeffding inequality for subgaussian random variables, we have the following bound: (0(xo) g(xp) I (2a((T + 1)-2/2 / (o)() with a probability of at least 1 - 5, where T is the maximum number of iterations. Comparing this with Lemma 4.1.3, we can further derive: l0(x)) - g(x))l S (T & V2 + 3- (t)((t + 1)- (log(6/5). with a probability of at least 1- for all < t < T. Consequently, the introduced error term would be absorbed in Ko,, and will not affect any parts of the analysis.
1) 1log///
3"
Variance reduced estimators. As mentioned above, a key point in the design of our stochastic bilevel algorithms is to select K, properly such that Et contains X with high probability, for all t 2 0. To
5

---

Algorithm SBCGI
1: Input: Target accuracy: ff, > 0, probability > 0, step size: 0t,, St, St, Pt. 7 > 0
2: Initialization: Initialize Xo E 2 such that gxxo) 9" /g//
3: for (= 0,..., ,
4: if' = 0 then
5:
f. V fx,,, 01), Vg. = Vg(xe, Ec), 9t = 5(x², Et
6: else
7:
a1- Update the estimate of Vf. (1
8:
Update the estimate of Vg. Vg, = v ((
9:
Update the estimate of g. 91 (1 -1-18 1Pp))rr- g() Et) -1-
10:
end if
11:
Compute S. arg mi (V 8) where N₂ = E 3: (Vi,, 8-xx)
x1) S <g(x0) (0x0)-9)-00 K.}
12: Update the variable Xx+1 = (1 7t- )x + 7+188t
13:
end for
achieve such guarantee, we first need to characterize the error of our gradient and function value estimators. More preciscly, suppose that for our function estimator have that P( g S
Ko,{) 2 5' and for the gradient estimator we have P(IVV Vg(xx)) S K1,²) 2 1 8', for some 5' E (0, 1). Then, by setting K; = Ko, + 1 we can guarantee that the conditions required for the inequalities in 7 hold with probability at least (1 10--0 22")..
Using simple sample average estimators would not allow for the selection of diminishing K,, as the variance is not vanishing, but by using proper variance-reduced estimators the variance of the estimators vanishes over time and eventually, we can send K; to zero. In this section, we focus on two different variance reduction estimators. For the stochastic setting in (I we use the STOchastic Recursive Momentum estimator (STORM), proposed in [35], and for the finite-sum setting, we utilize the Stochastic Path-Inte hath-Integrated Differential EstimatoR (SPIDER) proposed in [3).. If v is the gradient estimator of STORM at time 1, the next estimator is computed as
a1)v, + Vi(x., 01)
a. Vix
(10)
where V fx(x, f) is the stochastic gradient evaluated at x with sample 8. The main advantage of the above estimator is that it can be implemented even with one sample per iteration. Unlike STORM, for the SPIDER estimator, we need larger batch of samples per update. More precisely, if we consider Vt-1 as the estimator of SPIDER for fxxx), it is updated according to
V. = fs(.).) - fi((x-1)) + Vt-1,
(11)
where fs()) = (1/5) iES V fx( 0,) is the average sub-sampled stochastic gradient computed using samples that are the set S. As we will discuss later, in the finite sum case that we use SPIDER, the size of batch S depends on which is the number of component functions. We delay establishing high probability error bound for these estimators section 4.1
3.1 Conditional gradient algorithms with random sets: stochastic and finite-sum
Next, we present our Stochastic Bilevel Conditional Gradient method for Infinite sample case abbreviated by SBCGI for solving ( and its finite sum variant denoted by SBCGF. In both cases, we first find point Xo that satisfies (xx) - g" S <g/2, for some accuracy 6g. The cost of finding such a point is negligible compared to the cost of the main algorithm as we discuss later. At each iteration f, we first update the gradient estimator of the upper-level and the function and gradient estimators of the lower-level problem. In SBCGI, we follow the STORM idea as described in steps 4-6 of Algorithm [ while in SBCGF, we use the SPIDER technique as presented in steps 7-10 of Algorithm 2 In the case of SBCGF, we need to compute the exact gradient and function values once every 9 iteration as presented in steps 4-6 of Algorithm 2 Once the estimators are updated, we can define the random set & as in ( and solve the following subproblem over the set 2;;
S; = arg min (V)
(12)
sEA;

---

Algorithm 2 SBCGF
1: Input: Target accuracy: ef, eg > O, probability accuracy: 50, 51 > 0, step size: 7t >
2: Initialization: Initialize Xo E 2 such that g(xo) 9" eg/
3: for (= -,...,, , do do
4: if' mod(t,,))) q) = 0 then
5: Set Vs.- V fxx)). Vg.- - 9:
6: else
7:
Draw S samples
8:
Update the estimate of f as Vf, =V + VVs(xx)) - (xs(x
9:
Update the estimate of g as Vg.- = -9g-1 + (gs((x)) (x
10:
Update the estimate of g as 91 -t-1 + (x( as (x(-1)
11: end
T
12:
Compute S, arg mi { f. 8) where X; {st E Z (Vg,, 8-,
xx) S (x(x)) -g
13:
Update the variable X++1 = (1 - 7--1) )x; + 7+18$t
14:
end for
where fe is the unbiased estimator of V (xx). Note that we implicitly assume that we have access to linear optimization oracle that returns a solution of the subproblem in (I2). which is standard for projection--ii creee methods [32 33.. In particular, if Z can be described by a system of linear inequalities, then problem (12) corresponds to linear program and can be solved by standard
solver as we show in our experiments. Once, S; is calculated we simply update the iterate
X++1 (1 -*++))x; + 7
(13)
with stepsize "--1 E [0, 1). The only missing part for the implementation of our methods is the choice of K2 in the random set and the stepsize parameters. We address these points in the next section.
Remark 3.2. SBCGI can be implemented with batch size as small as S = 1. However, this does not imply that the batch size "has be" 5 = In other words, the main advantage of SBCGI, compared to SBCGF, is its capability to be implemented with any mini-batch size, even small as S = 1. Therefore, for SBCGI, the batch size can be set arbitrarily, whereas for SBCGF, it must be vTi.
Remark 3.3. In the finite-sum setting, if the numbers of functions in the upper- and lower-level losses are different, we could simply modify SBCGF by choosing S, = 9a = /nu and S = qi = vnn, where ny and ny are the number of functions in the upper- and lower-level, respectively.
4 Convergence analysis
In this section, we characterize the sample complexity of our methods for stochastic and finite-sum settings. Before stating our results, we first characterize a high probability bound for the estimators of our algorithms, which are crucial in the selection of parameter K, and the overall sample complexity.
4.1 High probability bound for the error terms
To achieve a high probability bound, it is common to assume that the noise of gradient or function is uniformly bounded as in 33 37), but such assumptions may not be realistic for most machine learning applications. Hence, in our analysis, we consider a milder assumption and assume the noise of function and gradient are sub-Gaussian assumptions, we next establish a high probability error bound for the estimators in SBCGI.
in Assumptions 2.2 and 2.3 respectively. Given these
Lemma 4.1. Consider SBCGI in Algorithm 7 with parameters or 3, Pt Ttt where w E (0, 1]. If Assumptions 2.1 2.2 and 2.3 are satisfied, for any 2 1 and å E (0, 1), with probability least for some absolute constant c, (d the number of dimension), we have
3"
||V,, vx(x,)) S 2v2((2,, +
3 1- (t) 11- /oog/b//))))
(14)
3
VVg. - Vg(xi)ll S cV2(22,
+ ot + +1) -1)- 1 (log(6d/8)
(15)
3
3
l0. - -(xx( S cv2(2L,D+ +
3 - at((t + +1) -22 /To((6/8)
(16)
7

---

Lemma 4.1 shows that for any E (0, 11, if we set a, = = = 1)", then with high probability the gradient and function approximation errois converge zero sublinear rate of O(1/t=/2) Moreover, the above result characterizes the choice of K, Specifically, if we define K1, as the upper bound in (15) and Ko. as the upper bound in (16). by setting K2 = Ku, + D then with probability (1 5) the random set &g contains X. Later, show that w 1 leads the best complexity bound for the convex setting and 2/3 is the best choice for the nonconvex setting.
Next, we establish a similar result for the estimators in SBCGE.
Lemma 4.2. Consider SBCGF with stepsize n and 5 9 = vn. If Assumpptions 2.7 2.3 hold, for any 2 1 and 8 E (0, 1). with probability 1 8 we have || f [I) S 4L, Dy (22(1)/)), ||Vg Vg((x)] S ALg Dy /t///o/(////// (22/6), and 99 9(x,)l S 4L,D7v log (OPDDr/oog11//)) (12/0))
Similarly, for SBCGF, we set K1,t 4L, V llo((1288) and Ko.: / Lao.D//loo(1/// and choose K, = Ko,, + 1 then the random set XG contains Xÿ with probability Next, we formalize our claim about the random set with the above choice of Kt.
Lemma 4.3. If Xg is the solution set of the lower-level problem and Xe is the feasible set constructed by cutting plane at iteration t, then for any t 2 and E (0, 1), we have P(X; S &,) 2-6.
This lemma shows all X² is subset of the constructed feasible set &, with a high probability of 1 - 6. Indeed, using a union bound one can show that the above statement holds for all iterations up to time i with probability t5.
4.2 Convergence and complexity results for the stochastic setting
Next, we characterize the iteration and sample complexity of the proposed method in SBCGI for the stochastic setting. First, we present the result for the case that is convex.
Theorem 4.4 (Stochastic setting with convex upper-level), Consider SBCGI in Algorithm Mfor solving problem O Suppose Assumption 2.7 2.2 and 2.3 hold and J is convex. If the stepsizes of SBCGI are selected as ar = Bi = Pt = = (t I)- and the cutting plane parameter is 3
K, c(([L + ()) log (og + D(2L, + 2 20) (ogg(6d/8)){t + 1)-1/2, then after T iterations,
g(xT) 0" C1C VT
+
Lg Do
logT
Eg 2
S(XT) 8
C2f VT
+ f(xo) f + D2 log T
T
T
with probability 8 for some absolute constants C1 and C2 and Ç == Vlog (6Td/8).
Theorem 4.4 shows a convergence rate of O O(ylog(T//6)/VT
As a corollary, SBCGI re- after Ollo(d/66e){e) iterations, where mini(e, ) Since we use one sample per iteration, the overall sample complexity is also (log(d/fc)) Note that the iteration complexity and sample complexity of our method outper- form the ones in 16), as they require (11/c4) iterations and sample to achieve the same guarantee.
turns an -optimml solution with probability
Remark 4.1. The task of finding X, which is equivalent to a single-level stochastic optimization problem requires ((1/2)) iterations and samples. As a result, this additional cost does not affect the overall complexity of our method. The same argument also holds the non-convex case.
Theorem 4.5 (Stochastic setting with non-convex upper level). Consider SBCGI for solving problem [L. Suppose Assumptions 2.1 2.3 hold. f is nonconvex. and define f = minxez fxx). the stepsizes
of SBCGI are selected 04 = 3 Pt (t 1) 1) 2/3 2 (T 1)-2/3 and the cutting plane 3²/3
parameter is K2 = c((2,, 33/3_ a1) 22 (8) ((LL D+ 32/3_1 ;g) /2g)/22l//) log(oo//))))))))
32/3
1) 1/33, then after T iterations, there exists an iterate Xx- the set E {xx, X1, ....... x1 for which
g(xx-)) -9" Cas + L + tg
f(x0) - Cas D2
g(x.)
(T +11/3
(T + 1)1/3
with probability 1 1- for some absolute constants Cs and CA and ;= Vlog (6T'd/).
8

---

As corollary of Theorem 4.5 the number of iterations required to find an (g), -optimall solution can be upper bounded by (logg(d/8) (d//d)d)/)))) /c²), where min(er.Cg) We note that the dependence on the upper-level accuracy ef also matches that in the standard CG method for single-level non-
convex problem 33 38 Moreover, as we only need one stochastic oracle query per iteration, SBCGI only requires (log (d/ōe) (ea) stochastic oracle queries find an (ez, )o)opiiall.
4.3 Convergence and complexity results for the finite-sum setting
Similarly, we present iteration and sample complexity for algorithm 2 under the finite-sum setting.
Theorem 4.6 (Finite-sum sum setting with convex upper-level). Consider SBCGF presented Algorithm
D solving the finite-sum version of (I. Suppose Assumptions 2.7 2.2 and hold and f is convex. If we set the stepsizes of SBCGF as 7 = log T/T, S q - n, and the cutting plane parameter as K₁ = =AD(LI 12T 12T + Dg Vlo (2) log T /T, then after T itterations.
9(XT) - 9" (C35' + Dg log T
+ Eo 2' ((xT) f f(xo) f + CaL'log
T
T
with probability at least 5, for some absolute constant C₂ and Cs, and Ç' = /logg12T///
Theorem 4.6 implies the number of stochastic oracle queries is (log(1 /5e )3/2 ///2)///////). where minfe,, EgJ. which matches the optimal sample complexity of single-level problems [39].
Theorem 4.7 (Finite-sum setting with non-convex upper-level). Consider SBCGF presented in
Algorithm 2 for solving the finite-sum version of (I). Suppose Assumption 2.1 2.3 hold, and f is non-convex. Define f = minxez fxx). If the parameters of SBCGF are selected as 7 L/VT.
(og//122///// + L /0)/ q= =nn and the cutting plane parameter is K, = 4D(Li /
then after T iterations, there exists an iterate Xx. (Xx, (xo,X1,... , xT ,,,,.,,xr-11,,, for which,
g(xx-) - g- C75' C Lgg
Cg. g(x.) f(x0) - CaS
T1/2
2'
T1/2
with probability at least 1 5, for some absolute constants C7 and Cs. and Ç' /1og(12T/6).
As corollary of Theorem 4.7 the number of stochastic oracle queries is (log(1//e) /n/e²) where
= min }g which matches the state-of-the ahrt single-level result O( (2) in 40 SBCGF also improves the number of linear minimization oracle queries of SBCGI from 0(1 (11/²2) to O(1(4) for convex upper-level and from ((1/e) to O(1 (c2) for non-convex upper-level.
5 Numerical experiments
In this section, we test our methods on two different stochastic bilevel optimization problems with real and synthetic datasets and compare them with other existing stochastic methods in [16] and [3]
Over-parameterrzzd regression. We consider the bilevel problem corresponding to sparse linear regression introduced in (3). We apply the Wikipedia Math Essential dataset (30) which composes of a data matrix A E R with
1068 samples and d 730 features and an output vector b E R", To ensure the problem is over-parameterized, we assign 1/3 of the dataset as the training set (Atr, b 1/3 as the validation set (Aval byal ) and the remaining 1/3 as the test set (Arga Diest ). For both upper- and lower-level loss functions we use the least squared loss, and we set X = 10. We compare the performance of our methods with the aR-IP-SeG method by [16] and the stochastic version of DBGD introduced by [1). We employ CVX [41 42) to solve the lower-level problem and the reformulation of the bilevel problem to obtain g" and f", respectively. We also include the additional cost of finding Xo in SBCGI and SBCGF in our comparisons.
In Figure (b)),), we observe that SBCGF maintains a smaller lower-level gap than other methods and converges faster than the rest in terms of upper-level error. SBCGI has the second-best performance terms of lower- and upper-level gaps, while -R-IP-SeG performs poorly terms of both lower- and upper-level objectives. The performance of DBGD-sto for the upper-level objective is well, however, underperforms in terms of lower-level error. Figure 1(c), SBCGE, SBCGI, and DBGD--to achieve almost equally small test errors, while aR-IP SeG fails to achieve low test error. Note

---

10
10
10
((ok) 71
104
10
Iv/A]
Test
IP
10
10
10
samplex
x100
sxmpim
x10F
mumber semplrs
<10²
(a) Lower-level gap
(b) Upper-level gap
(c) Test error
Figure 1: Comparison of SBCGI, SBCGF, aR-IP-SeG, and DBGD-Sto for solving Problem
10
10
HELUP
104
rameber Gaeplis
exedee of cample
Lower-level gap
(b) Upper-level gap
(c) Recovery rate
Figure 2: Comparison of SBCGI, SBCGF, aR-IP-SeG, and DBGD-Sto for solving Problem (
that after the initial stage, SBCGI increases slightly in terms of all the performance criteria, because SBCGI ( only takes sample per iteration and uses a decreasing step-size while SBCGF takes vni samples per iteration and uses a small constant stepsize, demonstrating a more robust performance.
Dictionary learning. To test our methods on problems with non-convex upper-level we consider problem (5) on synthetic dataset with a similar setup to [1]. We first construct the true dictionary D' E R 25x50 comprising of 50 basis vectors in 25 All entries of these basis vectors are drawn from the standard Gaussian distribution and then normalized to have unit mn-nrm. We also generate two more dictionaries D* and D' consisting of 40 and 20 basis vectors in D", respectively (thus they share least 10 bases). These two datasets A {ai.... a250 } and A' (a'...... 2 are constructed as
D*x; + n,, for i = 1...., 250, and a D'x', + n for k ..... 250, where {x; Ti=1> {x², 1250
250
1250 are random Gaussian noises. J A==1 are coefficient vectors and {n, 250 [m6 J==l
As neither A nor A' includes all the elements of D it is important renew our dictionary by using the new dataset A' while maintaining the knowledge from the old dataset A.
In our experiment, we initially solve the standard dictionary learning problem employing dataset A, achieving the initial dictionary D and coefficient vectors {x]250- We define the lower-level objective the reconstruction error on A using [x],00, 250 and the upper-level objective as the error on new dataset A'. We compare our algorithms with aR- Se- and DBGD (stochastic version). measuring performance with the recovery rate of true basis vectors. Note that basis vector d; in D' is considered successfully recovered if there exists d in D such that (dd; d > 0.9 (for more details of the experiment setup see Appendix FF. In Figure 2(a), we observe SBCGF converges faster than any other method regarding the lower-level objective. While SBCGI has the second-best performance in terms of the lower-level gap. aR-IP-SeG and DBGD-sto perform poorly compared that SBCGI, SBCGF, and DBGD-sto
with SBCGI and SBCGF In Figures 2(b) and (c), achieve good results in terms of the upper-level objective and the recovery rate. However, aR-IP-SeG still performs poorly in terms of both criteria, which matches the theoretical results in Table [
Acknowledgements
The research of J. Cao, R. Jiang, and A. Mokhtari is supported in part by NSF Grant 2127697 and the NSF Institute for Foundations of Machine Learning (IFML.) at UT Austin. The research of N. Abolfazli and E. Yazdandoost Hamedani is supported by NSF Grant 2127696.
10

---

References
[1] Jerome Bracken and James McGill, Mathematical programs with optimization problems in the constraints. Operations research, 21(1):37-44, 1973.
[2] Stephen Dempe, Nguyen Dinh, and Joydeep Dutta. Optimality conditions for a simple con- vex bilevel programming problem. Variational Analysis and Generalized Differentiation in Optimization and Control: In Honor of Boris Mordukhovich, pages 149-161, 2010.. 2010.
[3] Joydeep Dutta and Tanushree Pandit. Algorithms for simple bilevel programming. Bilevel Optimization: Advances and Next Challenges, pages 253-291. 2020.
[4] Yekini Shehu, Phan Tu Vuong, and Alain Zemkoho. An inertial extrapolation method for convex simple bilevel optimization. Optimization Methods and Software. 36(1):1-19, 2021.
[5] Zalán Borsos, Mojmir Mutny, and Andreas Krause. Coresets via bilevel optimization for contin- ual learning and streaming. Advances in Neural Information Processing Systems, 33:14879- 14890, 2020.
[6] 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, pages 1568-1577. PMLR, 2018.
[7] Amirreza Shaban, Ching-An Cheng, Nathan Hatch, and Byron Boots. Truncated back- propagation for bilevel optimization. In The 22nd International Conference Artificial Intelligence and Statistics, pages 1723-1732. PMLR, 2019.
[8] Aravind Rajeswaran, Chelsea Finn, Sham M Kakade, and Sergey Levine. Meta-learning with implicit gradients. Advances in neural information processing systems, 32, 2019.
[9] Luca Bertinetto, Joao Henriques, Philip HS Torr, and Andrea Vedaldi. Meta-learning with differentiable closed-form solvers. arXiv preprint 180.:v:880.08/1111 2018.
[10] Mingyi Hong. Hoi-To Wai, Zhaoran Wang. and Zhuoran Yang. A two-timesee nootiimesclle framework for bilevel optimization: Complexity analysis and application to actor-critic. arXiv preprint 2000:00077.0700 2020.
[11] Haifeng Zhang. Weizhe Chen, Zeren Huang, Minne Li, Yaodong Yang, Weinan Zhang, and Jun Wang. Bi-level actor-critic for multi-agent coordination. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 34, pages 7325-7332, 2020.
[12] Ruichen Jiang, Nazanin Abolfazli, Aryan Mokhtari, and Erfan Yazdandoost Hamedani. A conditional gradient-basee method for simple bilevel optimization with convex lower-level problem. In International Conference on Artificial Intelligence and Statisties, pages 10305- 10323. PMLR, 2023.
[13] Chengyue Gong and Xingchao Liu. Bi-objective trade-off with dynamic barrier gradient descent. NeurlPS 2021, 2021.
[14] Shoham Sabach and Shimrit Shtern. A first order method for solving convex bilevel optimization problems. SIAM Journal on Optimization, 27(2):640-66 (00)))6006660.00. 2017.
[15] Harshal D Kaushik and Farzad Yousefian. A method with convergence rates for optimization problems with variational inequality constraints. SIAM Journal on Optimization, 31(3):2171- 2198, 2021.
[16] Afrooz Jalilzadeh and Farzad Yousefian. Stochastic approximation for estimating the price of stability in stochastic nash games. arXiv preprint arXiv:2003.01271, 2022.
[17] Saced Ghadimi and Mengdi Wang. Approximation methods for bilevel programming. arXiv preprint arXiv: arXiv:1802.089 0100.822002266 2018.
[18] Junjie Yang, Kaiyi Ji, and Yingbin Liang. Provably faster algorithms for bilevel optimization. Advances in Neural Information Processing Systems, 131103677-1622 2021.

---

[19] Prashant Khanduri, Siliang Zeng, Mingyi Hong, Hoi-To To Wai, Zhaoran Wang, and Zhuoran Yang. near-optimal algorithm for stochastic bilevel optimization via in neural information processing systems, 34:30271-30283 2021.
Advances
[20] Tianyi Chen, Yuejiao Sun, and Wotao Yin. Tighter analysis of alternating stochastic gradient method for stochastic nested problems. arXiv preprint 2wrX:v2206.... 1378/, 2021.
[21] Zeeshan Akhtar, Amrit Singh Bedi, Srujan Teja Thomdapu, and Ketan Rajawat. Projection-free stochastic bi-level optimization. arXiv preprint arXiv:2110.11221, 2021.
[22] Benoît Colson, Patrice Marcotte, and Gilles Savard. An overview of bilevel optimization. Annals of operations research, 153:235-2555 256, 2007.
[23] Lesi Chen, Jing Xu, and Jingzhao Zhang. On bilevel optimization without lower-level strong convexity. arXiv preprint arXtv:2301.00712, 2023.
[24] Risheng Liu, Yaohua Liu, Shangzhi Zeng, and Jin Zhang. Towards gradient -hased bilevel opti- mization with non-convex followers and beyond. Advances in Neural Information Processing Systems, 34:8662-6675 2021.
[25] Daouda Sow, Kaiyi Ji, Ziwei Guan, and Yingbin Liang. A primal-dual approach to bilevel optimization with multiple inner minima. arXiv preprint arXiv:2203.01123 2022.
[26] Han Shen and Tianyi Chen. On penalty- based bilevel gradient descent method. arXiv preprint wwXiv:2302.05185 2023.
[27] Feihu Huang. On momentum-based gradient methods for bilevel optimization with nonconvex lower-level. arXiv preprint arXiv:2303.039444 2023.
[28] Kenneth Kreutz-Delgado, Joseph F Murray, Bhaskar D Rao, Kjersti Engan, Te- Won Lee, and Terrence J Sejnowski. Dictionary learning algorithms for sparse representation. Neural computation, 15(2):349-33) -36,, 2003.
[29] Mehrdad Yaghoobi, Thomas Blumensath, and Mike E Davies. Dictionary learning for sparse approximations with the majorization method. IEEE Transactions on Signal Processing, 2766)2178- 2009. 2009.
[30] Benedek Rozemberczki, Paul Scherer, Yixuan He, George Panagopoulos, Alexander Riedel, Maria Astefanoaci, Oliver Kiss, Ferenc Beres, Guzmán López, Nicolas Collignon, et al. Pytorch geometric temporal: Spatiotemporal signal processing with neural machine learning models. In Proceedings of the 30th ACM International Conference Information & Knowledge Management, pages 4564 4573, 2021.
[31] Chenglong Bao, Hui Ji, Yuhui Quan, and Zuowei Shen. Dictionary learning for sparse coding: Algorithms and convergence analysis. IEEE transactions on pattern analysis and machine intelligence, 38(7): 2015.
[32] Martin Jaggi. Revisiting frank-wolfe:
Projection-free sparse convex optimization. In Interna- lional Conference on Machine Learning, pages 427-435. PMLR, 2013.
[33] Simon Lacoste-Julien. Convergence rate of frank-wolfe for non non-convex objectives. arXiv preprint arXiv: 1007.001607.00..00 0034,, 2016.
[34] Stephen Boyd and Lieven Vandenberghe. Localization and cutting-plan methods. From Stanford EE 364b lecture notes, 2007.
[35] Ashok Cutkosky and Francesco Orabona. Momentum bassed variance reduction in non-convex sgd. Advances in neural information processing systems, 32, 2019.
[36] Cong Fang, Chris Junchi Li, Zhouchen Lin, and Tong Zhang. Spider: Near-optimal non- convex optimization via stochastic path-integrated differential estimator. Advances Neural Information Processing Systems, 31, 2018.
12

---

[37] Jiahao Xie, Zebang Shen, Chao Zhang, Boyu Wang, and Hui Qian. Efficient projection-free free online methods with stochastic recursive gradient. Proceedings of the AAAI Conference on Artificial Intelligence, volume 34, pages 6446 6453, 2020.
[38] Aryan Mokhtari, Asuman Ozdaglar, and Ali Jadbabaie. Escaping saddle points in constrained optimization. Advances in Neural Information Processing Systems, 31, 2018.
[39] Aleksandr Beznosikov, David Dobre, and Gauthier Gidel. Sarah frank-wolfe: Methods for constrained optimization with best rates and practical features. arXiv preprint arXiv:2304.11737, 2023.
[40] Alp Yurtsever, Suvrit Sra, and Volkan Cevher. Conditional gradient methods via stochastic path-integrated differential estimator. In International Conference on Machine Learning, pages 7282-72 1291-729..11 PMLR, 2019.
[41] Michael C Grant and Stephen P Boyd. Graph implementations for nonsmooth convex programs. In Recent advances learning and control, pages 95-110. Springer, 2008.
[42] Michael Grant and Stephen Boyd. Cvx: Matlab software for disciplined convex programming, version 2.1, 2014.
[43] Roman Vershynin. High-diimnnsional probability: An introduction with applications in data science, volume 47. Cambridge university press, 2018.
[44] Iosif Pinelis. Optimum bounds for the distributions of martingales in banach spaces. The Annals of Probability, pages 1679-1706. 1994.
[45] Chi Jin, Praneeth Netrapalli, Rong Ge, Sham M Kakade, and Michael Jordan. short note concentration inequalities for random vectors with subgaussian norm. arXiv preprint arXiv:1902.033366 2019.
[46] Tong Zhang. Learning bounds for kernel regression using effective data dimensionality. Neural Computation, 17(9):2077-2098, 2005.
13

---

general form of the examples covered in the main Supplementary material
A
Additional Motivating Examples
The bilevel optimization problem in ( provides a versatile framework that covers a broad class of optimization problems. In addition to the motivating examples provided in the main body of the paper, here we also provide a generic example of stochastic convex constrained optimization that can be formulated as TT. We further present body.
Stochastic convex optimization with many conic constraints: Consider the following convex optimiza- tion problem
min B[f(x, 0)]
s.t. h(x,() E -K. Ag E Q
XER*
where K. C Rd is closed convex cone. This problem be formulated special case of ①I) E)) where d_x(-)) - -P-x(() denotes the distance funetion and P-x(-) denotes the projection map. Our proposed framework provides efficient method for solving this class of problems when the projections onto can be computed efficiently, while the projection onto the preimage h (-K, ) is not practical, c.g., when I. is the positive semidefinite cone, computing projection onto the preimage set requires solving a nonlinear SDP.
by letting (x. x
A.1 Lexicographic optimization
Example I (over-paranne parametrrmde regression) can be generalized as broader class of problem, which known as lexicographic optimization [13] and uses the secondary loss to improve generalization. The problem can be formulated as the following stochastic simple bilevel optimization problem,
min ((3)
s.t.
B E arg arg min fir (0) Ep,, [[ (x()
(17)
dEZ
In general, the lower-level problem could have multiple optimal solutions and be very sensitive to small perturbations. To tackle the issue, we use a secondary criterion L(-) to select some of the optimal solutions with desired properties. For instance, minimal f2-norm by letting L(3) = || which is also known as Lexicographic &a Regularization.
find the optimal solutions with
A.2 Lifelong learning
Example (dictionary learning) is an instance of a popular framework known as lifelong learning, which can be formulated follows,
min Σ ((x4,)))
Σ (((x,,3).0k) S
Σ ((txx, (18)
s.t.
(X.
(NWM))M
In this problem, the objective is the training loss the current tasks D, {x(, vt) v While the constraint enforces that the model parameterized by 3 performs no worse than the previous one on the episodic memory M (i.e., data samples from all the past tasks).
In the paper, we discuss variant of the problem above, where we slightly change the constraint and ensure that the current model also minimizes the error the past tasks. It can be formulated the following cum/sscchsm/cctticcccccc simple bilevel optimization problem (12).
min B
((x) .n
s.t.
B argmin
Σ
((xi, v3))
(19)
(X.M))EM
B Supporting lemmas
B.1 Proof of Lemma 4.1
Before we proceed the proof for Lemma 4.1 we present the following technical lemma, which gives an upper bound for complex term appearing in the following analysis.
14

---

Lemma B.1. Define P: 1/(t 1) where w E (0, and t 2 1. For all 2 2, let {si} be
sequence of real numbers given by
2
pk)
Pr
K=T
Then holds that
8t (1 1)
(20)
Proof. We prove the result by induction. For 2, we can verify that
3
2
1
$2
3
3
- 3
Now we suppose that the inequality in (20) holds when & T for some T 2 2, i.e.
ST Σ Pr
II pk))
11
First note that the sequence {s²} satisfies the following recurrence relation:
T+1
T+1
T+1
T
Σ
Pr II (1- mi) = -1- p1+1)+
II (1 ph))
ST+1
k=
7=2
k=T
2 pk)
= (1 (1-rr+1) pr+1)²
+ ph+
Pr
(1 pT+1| )2 (ST + pē+1).
Moreover, since w E (0, 1}. we have (T 23 - S (1 + 1) Therefore, we obtain
(T' 2) - 2
1
1
8T+1
(T + 2)
(1+
(T 2)22
((T + 2) - 1)(( + 1)t
1
1
(1+1)= +
(T + 2)2w
(T+
(T' + 2) - (T+ 1) +1
(T+ 2)2w
(T 1)
(T 2) (t
1)tH (T + 2) - (t + 1)
2) 2w (t + 1)x
(T 2) 11
1
[T 2 1)
(T+2) 2)
By induction, the inequality in 20 holds for all 2 2.
Now we proceed prove Lemma 4.1
We show the proof of part (i) here. The proof of part (ii) is very similar to part (i). The first step is to reformulate
Proof of Lemma 4.1
Vg. Vy(x,) as the sum of martingale difference sequence. For / 1, by unrolling the reucurrence we have
e: :0= -11 8.) (Vğ(x., {.) Vg(()))
+ +(1 (1-3)) (1133)(V0((....
Vg(x. (Vg(() -Vx())
II (1 3x e1 + ΣII (L) - Vo(x,- -(V(xx)xx))))))) - Vg(x--1))
2=2
k=
Σ B, II (1- B. -V)(x)))
7=2
k=T+1
(21)
15

---

Thus, we can write e, as the sum e, = Ei= Sr. where we define (1 = ITk=2 - Bx)e1 and
II (V - Vg((+-1))
(22)
B II (1 (1--Ba (V) Vg(x+))
(23)
km7+1
for r > 1. Recall that Vo(x,, (1) Vg(x1). We observe that B[s=| = 0 where Fr-1 the or-fiel generated by (I......... ,.al,..,,,,,,,,, 11 1--11 Therefore, {(r]+= martingale difference sequence.
Next, derive upper bounds of ||s/||- To begin with, observe that for any 1,2,...,f,
1 II
II
(k 11
ka
7u (1 (1+11 11
B.) -1-Bx)
(ks + 13 II
II
(24)
II
k=y
where we used the fact that (k + 1)- - s knd the last inequality, By using the above inequality, can bound ||s/]| as follows:
2
2 (1 1)
[Vg(x,, vg(xi)ll
|kll| II (1-Bx)lleil) k=2
S (1+ 1) EI)
(7+ then by Assumption 2.3 ii) we have Elexp (I<:112/c]) I exp (1). Moreover, for r > 1, by triangle inequality, ||<>|| can be bounded by
Define 61
2
4 II 3.) & - vo(x,-1)1)
1|+ll -
k=T
+ 3,
II (1 Bb) Vx),) - vo(x,)1
km7+1
$ $2bgl Xrr - X-- II (1 3k) + (+) Vg(x,) II (1 - 3x)
km
kmT+1
DD II (1 Bx) |V &+) - Vg(x,)13,
II (1- 3k)
(25)
R=T
k=
3
S D3,
II (1 - Bx) + 3" - Vg(x+, Vg(x,)ll ||3,
(1 3xx
k=
k=r
= 2hg D + 3
Vg( X,, {+) - vo(xx)1) Br II - 3k)
k=
-(2L, 22g D + 3d
Vg((x, ) Vg(x,) |||
3, II (1 (A)
--
Og
k=T
Define e, = (2L, + - )r, 11 Bz). Note that if we have Elexp(X3/<2) < 1 and
Elexp(X3/cc) (X3//3)) I then we have E(exp((X; X2)2/(a + c2)3)) < 43 Thus, we have Elexp (IlG, // < exp (1) for all < < t. Hence by proposition with probability 8
2d 8
||o/|| Sc.
Σ 7=1
c
log
(26)
16

---

where c is an absolute constant, d is the number of dimension, and Σx= of can be bounded by as follows,
Lemma B.I I
2
3
Σ = Σ
2
+ +2L +
0g)
( II 3x))
(T + 1)20
3
km
2 og
(2L, D + 3--1" )2
(T + 1)20
(t + 1)
(27)
((v2)"00)²
(2Lg D + 3 3- 0 2
3
(++
(1+1))
2 D + 3 3
11
where the last inequality follows from the fact that (v2) S 3 33/(3" - 1) for any w E (0,
Combining (26) and (27). we have with probability at least 1 - 6,
[Vg(x,) - gell S 2V2 D+ + 3 5l) + 1) w0/2 w/2 log(2d/8)) del K1,t (28)
3
Similarly with probability at least 1
6',
[g(xx) gel S c22((2,)
3
ot w/2 / /log(2d/8)) def Ko.
(29)
3 1
and with probability least 1 - 5',
[V - -..I S eV2 3 - (t+ /11g) (10g(22/8)) (2d/6))
3
def K2,
(30)
where is absolute constant and d is the dimension of vectors. We can union bound to obtain that these three inequalities hold for at least probability 1 - 35 = 1 - 5. For simplicity, we define constant Ay and Ağ such that,
AT + 1) w/2 / (log(6d/0) = Ki,t and Ab Vlo((l/o5)) Ko.t
(31)
and similarly Ag (t + w/2 Tog(6d/8) K2,t-
B.2 Proof of Lemma
[t/g] for any E {O...., T 1}, then whenever tog according the Algorithm a full batch of sample gradients are selected, hence, Vg. = Vg( x,): otherwise, the error of computing sample gradient can be expressed as follows
Proof. Let us define to A
E,ei 5 (V9s() (x:) Vgs(s
Vg(x,) + V)(x,-1)),
(32)
where i is the index with S(i) denoting the i-th random component function selected at iteration l. Furthermore, from the update rule of Ie we have ||x: - x|- 74 /=|1s:11 - xi-1l1 S Dy for any t > 0, therefore.
||c,,|| S 1 (1))) - Vg((x--1)l + VVxx), Vg(xx-1)I)
1
(33)
5 24g ||x x,- 5 Dy
for all / E (1o + +1...., lo + 9} and E (l,..., S}. On the other hand, from the update of Vg, and
332 we have that for any * tog. Vg. Vg((x)) = Vg.- (9(x+-1) + E Cr.i- Therefore, by
continuing the recursive relation and taking the norm from both sides of the equality obtain
|[Vg. Vg(xx)1 (x)().)) (g((ee))
Σ
gall
j=to1
(34)
j=to 1=1
17

---

where the last equality follows from Vg(x,,) = Vg((x,)).
Then by Proposition E.1
we have
A2
A2
P( Vg((x)] 2 X) v 4exp( 45(t - to) 413
S expp
(35)
16
g
where the last inequality follows from the fact S vn and to S q Vn.
By setting
(+LgD))/lo(((8") for some 8' (0, we have with probability at least 1 6',
IVg. Vg(xx))l S 4Lg ALgg4//00(4/8) log(4/8))
(36)
Similarly, with probability at least 1 1-
19: - g(x+)l S 4L,Dry log(4/8"),
(37)
and with probability 8',
VV. vf(x))] DLy 1og(4/8')
(38)
Then by union bound and å = 35", we show these three equalities hold with probability 1 - 6.
B.3 Proof of Lemma
Proof. Let x be any point in X7 i.e., any optimal solution of the lower-level problem. By definition, have (xx) =g" Since g is convex and g" 9 (Xo), we have
9 (xo) 9 (xx) g" g (xx) (x,)) g (xx) (V -xi)
(39)
Add and subtract terms in inequality above, we have,
(Vg, x x.) + Gr g(xo) N kV Vx((;)) X xx) + 10e - g(xx)l
(40)
Considering the random hyperplane we used in we want to prove the following inequality holds with high probability,
Vg.. xx) + g(xo) K.
(41)
Recall K, Ko + DK1,t And Ko and K1, were set as the high probability bounds of
Vg:
Vg(x))] and 9r -xxx)] in Lemma 4.1 for Algorithm Oor Lemma 4.2 for Algorithm 2 Then
compare the two inequalities above and use Jensen's inequality, (Vge, X + g(xo) Ke
holds with high probability 1 - 6 for all 2 0. Hence, Lemma 4.3 holds with probability 1 - - for all 0.
B.4 Improvement in one step
The following lemma characterizes the improvement of both the upper-level and lower-level objective values after one step of the algorithms.
Lemma B.2. If Assumptions
2.7
2.2 2.3
are satisfied,
(i) For all 2 0, assume that X C Nr. Then we have
12
7x19(xx) & f(x:) f(x++)) + 7+1 D(Vf()) - F/ +
(42)
2
As corollary, if f is convex, we further have
D2
f(x+11) f" s (1- (1-2411(((x))) -f + D(V -f.11 +
(43)
2
18

---

(ii) We have
g(xx+1) - (9x0) S (1 - -1+)))()((())) g(x0)) + Dr 1 (x))) - - Vghll 1 + Ki,t)
+ (1+1()))(x,) Bell + ko,r) + Lg DgD
(44)
2
Proof. (i) Based on the Ly-smoothness of the expected function f we show that f fxx++1) is bounded by
[(x++1) S [(x,) + V (xx) (Xe+1
xx)
X - xx
(45)
2
Replace the terms X - X by M+11 x1) and add and subtract the term 7+t1 V (s: - xx) to the right hand side obtain,
1x++1) S fxx) + 7 (e+i(V((())) (s. 7+1 VT, (s; x,) + L1
2 Xx++ xzI1? (46)
By Lemma 4.3 X C Xr with high probability 1 - 6, for all / = 1,....
T. Note that if we define
8/ = arg m [xV)(x), X - s)}. Recall that FW gap is G(X) = max (VV (x), & - -)).
We can replace the inner product V f St) by its upper bound (V st)). Applying this substitution leads to
1(x+11) S f((x)) + 7 (V/(x)) - VTt) (s: - xx) + Y++1 Vh. (5) - xx) + L1
x -xel?
2
= -((x) + + (VS(x.) - VT.) (s: - xx) + M+t1 (VI. - vx(xx)) (5/ 166--x xx)
-741 -74+V/(())) (x: ())) L1 ||x+11 - -xell²
f(xx)) + 7+ + (Vf(x,) - vf)) (s: -s') -n+19(x;)) + I,f 2
X xx,12
S f(x)) + R+1 - V (+))(()) + I. 1,
(47)
Rearrange the terms for the inequality above, we can obtain,
L D
(x+) S f(x)) - f(x++1) + 7x4 VftIl As simple corollary, since (xx;) f(xx) when f is convex, have,
(48)
2
f(x++1) - < (1 - + 7 DI Vf.ll Lg D2
(49)
2
(ii) Based on the L, smoothhhssss of the expected function g we show that g(x(+1) bounded by
Lg
g(x++1) S g(xx) + Vg((x)) (xx+1 - xx) + | x;
(50)
2
T
Replace the terms X X2 by (s: xx) and add and subtract the term 7 Vgi (s: xx) to the right-hand side to obtain,
g(x++1)= S g(x,) + (Vgx)) Vg.) T (s, x,) + 7 Vgt (s, x,) + Lg 2
xll? (51)
Now by definition of the set X,, using (Vit: 8; xx) g(xo) 5r + Ko,t + DK1, In addition, we could use Cauchy-Schward inequality to upper bound the second term. Then add and subtract (x0) on the right hand side to obtain,
g(x++1) S (xx:) + (g(1(((x)) g(x))) + 74 7+1 D(Vgg))) Fgell
7 (g(x:) 9.) + (Ko. + +KK,)+ + Lg ||x++1 =x,
(52)
Then subtract g(xo) on both sides,
g(xx-1)= g(xo) (1 (g(x) 9(x0))
+ 7 1 - F + (1g(x:) - gell + Ko,: + DK1,)) + Lg
X - 12
2
(53)
and the claim in the lemma follows.
0
19

---

C Proof of Theorem for Algorithm 1
C.1 Proof of Theorem 4.4
Proof. For lower-level, by Lemma B.2 have
s(x++1) ((xo) S (1 x+1) - g(xo)) + Dr++1 ([x))) Villl + Ki,t)
7++1((9((x))) dcll Ko,e) Lgg Di
(54)
2
By Lemma 4.1
have Vx(x,) g - 1,, and [19(x;) Dell S Ko with probability 1 8'.
Plug them in the inequality above to obtain,
L
g(xx-1) g(xo) 5(1 (1 (gx)) g(xu)) + 1 Ko,))
2
$ (1 (o(x)) g((o))
(55)
2 ( /oo((///)) + 1AbVo+g (log((//)))
+ D
+ 1)3/2
2(t + 1)2
with probability 1 for all t. Let C 4 A + Ab) and = T6'. Then, by applying the inequality recursively for .,...." we obtain that
g(XT) ((x0) < (1 (1 7) (o(xT-1) g(xo)) + C1/21 log (6d/8")
Lg D²
T3/2
272
1
C C1/2 log (6d/5")
1
(0(x0) g(xo)) +
II
(1
1+7
(1. + 1)3/2
i+1
(=++1
Lg
T-1
(1
2(1 + 1)2
i+1
i=++1
C1/ 02/2% /log (6d/8')
D2 T-1
< 0
Vt+1
2T
t t+1
Ci log (6d/8")
D²
(1
+ log 7)
VT
2T
C1 / log (6td/0)
D1 log
T
(56)
with probability 1-8.
For upper-level, by Lemmal B.2 we have
T D.pD2
(1- f*)
D ([V Fl +
(57)
By Lemma 4.1 we have ||v)(x;) - Ff.l < inequality above to obtain,
A1 /log(6d//8) (141)122
with probability 1-8'. Plug it in the
DA / log(6d/8))
D2
f(x(1)) - 5 8+1
)(((x)) - -f*)
(58)
1)3/2
2(t 1)2
20

---

with probability 1- 8' for all t. Let C2 = 2DA2 and ė = T8'. Then, by applying the inequality recursively fort .,...,T T 1, we obtain that
DA2 D llg(6d/88)
D
f(XT) - 1 T
(f(xT-1)
T3/2
2T
T-1
-+11 (f(x))- - -f)
1
T-1 Σ
DA v (6g(///
(t + 1)3/2
i+1
D2D
T-1
(1
(1- -+1
II
2(t + 1)²
i=t+1
fxx) 8 DAby log (6//6
T-1
D T-1
T
Vt+1
2T
t+1
1=1
((x0) -
2 log(6d/5))
D22
(1
loglog T)
T
VT
2T
[(x0) - S
2DA1 log
(log//dd/// 8)
D
log Dogg
T
VT
T
(59)
with probability 1 The theorem is obtained.
C.2 Proof of Theorem 4.5
Proof. For lower-level, by Lemma B.2 we have
g(x(+1) - g(xo) s (1 (1(x(((x))) - g(Xo)) + Dit+a((VV(xx) Fo.1l + K2,))
+ (+1((1 - dell + Ko.)) + DgD
(60)
2
By Lemma 4.1
we have [x,))
- < K1 and [|o(xx) - gell S Ko. with probability 1 1 6'.
Plug them the inequality above to obtain,
g(x++1) g(xo) S (1 - -++1) - g(x))) + (DT++(DK,, - Ko,e) + Lg DgD
2
(T (T + 1)2/3 (xxx)) - g(Xo)
(61)
(A (log(bi/88) 22 V vlog(6d/88) 1)+3)1/3(1 + 1)2/3
LgD
2(T + 1)443
with probability 8' for all t. Let C3 = 2(DA 2/3 + A 2/3 J.Then we can sum all the inequality up for all to obtain,
Ca loggob////8")
I.g12
g(x++1) - g(xo) (1
(T 1)2/3 1)2/2) (o(xi) g(xo)) +
(t 13T + 1)2/3 2(T + 1)4/3
1
Ca (og(6Td//) Lg D2/2
<1 (T + 1)2/8 )o(x;) g(xo))
(I + (TT + 1)2/3
(62)
By induction, we have for all 2
Cs (og(666)))))) + Lg D2
g(x(+1) g(x0)= -
(63)
(I + 11133
with probability 1- where 8 = T8'.
21

---

For upper-level, by Lemma B.2 we have
7x+19(x)) - fxx)) f(x++1) + 7411 (++1)/VV(x))) ffell L 1y D2
(64)
By Lemma 4.1 we have [V Vs,1 S
1/ vi'ving Vlog(ed/6") (141)173
with probability 1 1- 5
Plug it and
7++1 17 1)2/3 inequality above to obtain,
T-1
T-1
D2
T-1
(x+) - f(x)) - f(xx) + D
7441 V.f(())(() - Vs.I+ +
2
7++1
1=
=
1=0
2/3 (log(6d/8")
D2
s f(xx) f(XT) DD _
(65)
1) 1/3 (T 1)2/3
2
(T 1)4/3
t=0
3 2/3
D
S f(x)) f(XT) + DA V (og((6d/8))
( 1)
Let X₁- = arg minigisT (xx),, then
T-1
g(xx-)) S
Σ 7
2imo Yo41
$ (T 1)2/3 (f(x)) - f(xx) + = = Vlog(6Td/8) +
1
D
(66)
(T + 1)1/3)
2
(T (113)) (f(x)) - + 2/ D
1/3 (log(6Td/8) + D2 1
2
(+11)188
with probability 8, where 5 = Tō'. By letting C4 DA2/3 2/3 tte theorem is obtained.
L
D Proof of Theorem for Algorithm 2
D.1 Proof of Theorem 4.6
Proof. For lower-level By Lemma B.2 we have
g(x++1) - g(xo) S (1 744 (g(xx) - g(xo)) + D ( Vg.1l + K1,))
+ 7+41 ((9(()) dell + +Ko.)) + I.g T.gDPn?211
(67)
2
By Lemma 4.2 we have [V))(;)) - Fg,1l S 41g Dr Tog 12/88) and [ gell S
4L, AL. (DDy/ooo(/ (2/ with probability 1 - 5'. Let Cs 8D(DL, L)) and T6', Plug
them in inequality above and let 7 7 = log T/T to obtain,
9(xr+1) g(xa)= S (g(x)) g(xx)) + (Cs (12/8))))) + Lg /22²
(68)
with probability T/T. Sum up the inequalities for all 1 < t $ to get,
T
g(xT+1) - g(xo) (1 7) (g(xo) g(xo)) + (C's V 1og((2/88)) Ly D2 /22 > (1 7)t
k=1
(Cs (Csylog log 12T 12T + D2 /2) log T
< 0 (Co 1C (2)0))((1)))+)))) 12
<
T
(69)
with probability 1 1-6 6. For upper-level, by Lemma B.2 have,
f(xT) - fr s (1 -- ff -f" DV+(-(((7-I) Vl. D2
(70)
22

---

Now we proceed by replacing the terms (xx) Vlrll || by its upper bounds from Lemma 4.2 i.e.
[V+ -Vlll < , Drylog (12/8),
f(x)) - f s (1 (1-7) (1x)) + " (12/6') + 1/2)
(71)
with probability (1 (1-8 8') And choose 8 = 3T8 Then by telescope, with 7 log T obtain.
T
we can
f(xT) (1 - f*) + (4y 1/2)) D2 7)
$ - (11-)) (f(xx) - f*) + (4 log (12/8)+ 1/2)
(72)
S (exp(~>T)(f(x)) - + (4 log +
S (f(x)) f"/T (4y/og (12T/0) + 1/2)L, log T/T
with probability 1 - 8. Note that without loss of generality, less than O, we can bound it by O. By letting Ce - 5Ly D2, the theorem is obtained.
assume f(x)) fr 2 0. is
D.2 Proof of Theorem 4.7
Proof. For lower-level, by Lemma B.2 have
g(xe+1) - g(xo) s (1 )+1)(((x))) - g(xo)) + D4++ ( Vg(x)) - Vg.l + Kie))
+ 7+11(19((x) - dell + -Koo)) + L, D2g2
(73)
2
By Lemma 4.2 we have [[])(x;) lzlll S Dy
Dy Do/DD//oo////0000 and [[g(x;) gell S
4L.Dr//ogg(( (12//1/88))) with probability 1 - 6'. Let C7 = 8D(DLg - Li) and = T6 Plug them in inequality above and let 7 1/ /T to obtain,
T1/2 (g(x)) 9((x)) + Cr Toylog( Tog/oo/////// (22/8))
I.g
g(x++1) - gxx)) S S1
+
T
2T
(74)
$ 2 T1/2 (((x)) - g(xo)) + CT Crylog 12/8)) + L.g. D2/22/2
T
with probability 1 1 5/T. Sum up the inequalities for all 1 get,
T1/2 ([o((x)) -[xx)]] + (Crv (C7//og(12/66) Lg D2/2)
1
1
g(xx-1) g((x)) (1
T1/2
k=1
C7 CT CoyTog( (1)7))))) + 12/2
T1/2
(75)
with probability 1 - 5. For upper-level, by Lemma B.2 have
L L D2
7+19((x)) S (xx) f(x+11) + ((1 ff +
(76)
By Lemma 4.2 we have [x,)) &l.l S Dry A (12()) with probability 1 - 5'. Plug it
and 7+41
1/vT in inequality above to obtain.
D2 T-1
Σ ((x:) s f(xo) - f(XT) D
7041 [ Fl
2
2
7++1
(77)
$ f((o)) - f(xT) + Ly D2 (4 log( (28')))))))))))) + 1/2)
23

---

Divide both sides by vT, we can get, Let X,- arg mig SeST (xx). then
1 T-1
v f(x0) - f L, D2 log log( 1/2)
((x-)) -
Σ ((x,)
(78)
T
T1/2
1=0
with probability By letting Cs 5L,, the theorem is obtained.
E
Axuma-Hoefiding-yppe inequalities
In this section, we present two useful vector versions of Azuma-Hooffding-yyp concentration inequalities with uniform bound assumption or sub-gaussian assumption. They are crucial in our high probability analysis.
Proposition E.1. (Pinelis and other 1994 (44), Theorem 3.5) Lel CI...., Se IR be a vector-
valued martingale difference sequence w.r.f. filtration {Fi), i.e. for each T E 1,...,h, , we have B|sr = O. Suppose that ||<= S c, almost surely. Then 2 1,
X
P(I Σ Sil 4exp(
(79)
4 + c?
Proposition E.2. (Jin el al. 145), Corollary 7) Let Si,.... St E Bdd be vector-valued martingale difference sequence w.r.t. a filtration {Fi), i.e. for each 1,...,l, have [[ Suppose that E E[exp[]]<>]I /12 /c,] < exp(1). Then there exists a absolute constant e such that, any 8 > with probability at least 1 8,
2d c log
fll Se.
(80)
8
This proposition was also used in previous literature including [46] and [37]. is common to use such martingale inequality to obtain some eigh-probability results recently.
F Experiment details
In this section, we include more details about the numerical experiments in Section 5 For complete- ness, we briefly introduce the update rules of aR-IP-SeG in [16] and DBGD [13 In the following. we use the notation Il z-)(-) (-) to denote the Euclidean projection onto the set Z. The aR-IP-SeG algorithm given by,
y++1 Il 7o (V Vg
X 1Iz((: Vel (Vf(y, 01)) + P: vg(y..())
T441 П. + (rup))"
(81)
Yi+1 T. + ('e
P(+1
where 7 is the stepsize, Pt is the regularization parameter, and Yr is the output of the algorithm. In this experiment, we choose 7 7// (L 1/3/1 and pt po + 1)1/1 for some constants "Yo and Po- The DBGD-sto is stochastic version DBGD, which simply replaces stochastic gradients. Although the stochastic version of DBGD does not have a theoretical guarantee, it has been used to solve stochastic simple bilevel optimization problems [13), which worked pretty well empirically. Hence, we use it as baseline for solving stochastic simple bilevel problems and compare it with our proposed algorithms. The DBGD algorithm is given by
gradients in DBGD with
X = X Yk (V. (xx) + d Vg (x+))
where Y is the stepsize and we set λ
24

---

Φ (xx)
(V/ (xx() Vg (xx)) VV (xx)
.0} and 6(x) = mmnn {a(g(x) - 9). 311 Vg())
A = max
where a and 3 are hyperparameters and 9 is lower bound of g". In this experiment, we choose 9= We also note that [13 only considered unconstrained simple bilevel optimization, i.e. 2 R
We further project X2 onto Z for each iteration to ensure the constraints are satisfied.
F11 Over-parameterized regression
Dataset generation. The original Wikipedia Math Essential dataset 330 composes of data matrix of size 1068 x 731. We randomly select one of the columns as the outcome vector b E R 1000 and the rest to be
matrix A E IR x068x We set constraint parameter À 10 in this experiment. Initialization. We run the algorithm, SPIDER-FW [40], with stepsize chosen as 7: 0.1/(l 1) on the lower-level problem in ( We terminate the process to get Xo as the initial point for both SBCGI and SBCGF 2 after 10 stochastic oracle queries.
Implementation details. We query stochastic oracle 9 10y times with stepsize 7t- 0.01/ and 7 10-5 for SBCGI and SBCGF with K, 10-4 ////// /vt + respectively. In each iteration, we need to solve the following subproblem induced by the methods,
min f (Bx).s s)
s.t. ||s||2 S X, (Vg (Ba) s 3x) و (Bx)
(82)
Introduce s st.s 2 0 such that st -s Then we can reformulate the problem above as follows,
min VV (Bx) st -
(83)
-Bx) s.t. s.t.s s 2 0, (st 1) + (s 1) I X, Vg (Be), - S - -Bx) 9 ( - (Bx)
where 1 E Rd is the all-one vector.
For aR-IP-SeG, we choose To 10-7 and pa 10 For DBGD, we
and "/ 10-6
F22 Dictionary learning
Dataset generation. We generate 500 sparse coefficient vectors {x:) and {x() nonzero entries, whose absolute values are drawn uniformly from [0.2, 1]. The entries of the random
250 with 5 random
noise vectors {n,},-) 250 and (n') 250 [n() are drawn from i..d. Gaussian distribution with mean 0 and
standard deviation 0.01.
Initialization. We use a similar initialization procedure as (1), which consists of two phases. In the first phase, run the standard Toalk algorithm both the variables D R25x40 and X E 250 for 10a iterations with the stepsize fix the variable X and only update D using the Frank Wolfe algorithm with exact line search for additional 10 iterations to obtain D and X the initial point for the full bilevel problem.
1/vt=I1 Next, in the second phase, we
Implementation Details. We choose 5 = 3 in both problems (5). To be fair, all four algorithms start from the same initial point. We slightly modify the initial point by letting D E R225550 be the concatenation of D E 25x40 and 10 columns of all zeros vectors. Furthermore, we initialize another variable X randomly by choosing its entries from a standard Gaussian distribution and then normalizing each column have a & -nrmm of 5. We choose the stepsize as 7 = 0.1/(t + 1)2/3 and 7 10-3 for our SBCGI and SBCG12 with K; 0.01/(t 1)1/3, respectively. Empirically, we observe that taking one sample per iteration leads to a very unstable process in this problem. In this case, we choose mini-batch of size 8 for SBCGI, aR-IP-SeG, and the stochastic version of DBGD. For each iteration, we will solve the following subproblem,
min V
D X3) .))
s.t.
|ā.1|_ < 1, Vg (Dx) DD -DD s (D.) - (D(b)
(84)
The above problem can be reformulated by using the KKT condition, which equivalent to get root of the following one-dimensional nonlinear equation involving À 2 0:
D = 1l (V/b D Xx) AVg (Dx)): (Vg (Dx) D D2) = (Do) -9 (D:) (85)
25

---

10
10
10
ROOKG 3U9DD-90
Iv/A] V7
UtA) 11
104
10
Tes
SBCEP
10
10
(4)
(a) Lower-level gap
(b) Upper-level gap
(c) Test error
Figure 3: Comparison of SBCGI, SBCGF, aR-IP-SeG, and DBGD-Sto for the over-parameterized
regression problem
109
3100 -.
ALSEOPP Biirb-ik
10
20.00
12
1130:
(a) Lower-level gap
(b) Upper-level gap
(c) Recovery rate
Figure 4: Comparison of SBCGI, SBCGF, aR-IP-SeG, and DBGD-Sto for solving the dictionary
learning problem.
1, I,..., 50} is equivalent to project each column the Euclidean ball. In practice, the reformulated problem can be solved efficiently by MATLAB's s root-finding solver.
where the projection onto Z D
IR
25x50
For -R-IP--eed we choose Yu 10 and Po 1. For the stochastic version of DBGD, set a = 3 = 100 and "/ = 5 x 10-3
Additional plots illustrating the comparison of the studied methods in terms of runtime rather than the number of sample used are provided in Figure and Figure
F.3 Experiments with different random seeds
We further repeat the experiment times with different random seeds to see more realizations the stochastic algorithms. The results are reported in Figure and Figure The solid lines denote the average statistics over 10 trials of the algorithms. While the shaded regions surrounding each line reflect the span of all the random instances involved. Figure 5 and Figure Figure and Figure which eliminates the possibility of choosing a particularly good instance.
6 present similar results as
150
107
SHIMENSU
10
-10
Euappio
x107
musaplo
(a) Lower-level gap
(b) Upper-level
(c) Test error
Figure 5: Comparison of SBCGI, SBCGE, aR-IP-SeG, and DBGD-Sto for solving Problem 3 with 10 different random seeds
26

---

10
10"
104
auddec of Dezgli
craaler Suaplecs
Lower-level gap
(b) Upper-level gap
Recovery
Figure 6: Comparison of SBCGI, SBCGF, aR-IP-SeG, and DBGD-Sto for solving Problem (5) with 10 different random seeds
GEEAO SROSE SAIIEE PM
231077
10
1015
10
S110ES PM
12
cauglie
X00
sander docogio
txseler meaple
v10
(a) Lower-level gap
(b) Upper-level gap
(c) Recovery rate
Figure 7: Comparison of SBCGI, SBCGE, SBCGI-M, SBCGE-M, STORM-FW, and SPIDER-FW for solving Problem (5)
F44 Importance of the right cutting plane
In this section, we numerically illustrate the importance of choosing the right cutting plane on Example 2 (dictionary learning). Specifically, we compare our proposed methods with the ones without cutting plane and with an unregularized cutting plane (without additional term Kj). If we replace the stochastic cutting plane
p with the unregularized cutting plane (8) in SBCGI and then the algorithm usually fail at some point
SBCGF 2 parameters chosen, based on our experimental observations. More specifically, algorithms' failure means that the subproblem of dictionary learning (85) infeasible. we slightly modify by adding a checkpoint before solving the subproblem. If the subproblem is infeasible at the current
process, depending on the datasets and
iteration, then we choose the update direction s, Vge. This adjustment prevents unnecessary interruptions during the process and enforce the algorithms to focus only on the lower-level problem when the subproblem infeasible. We denote the modified algorithms SBCGI-M and SBCGF-M. Moreover, we also take SBCGI and SBCGF without cutting planes into consideration, denoted as STORM-FW and SPIDER-FV In fact, this case, the bilevel algorithms degenerate single-level projection-free algorithms similar to algorithms in 37 and 40
Figure 7 (a) indicates that SBCGI-M and SBCGF-M focus more on the lower-level problem due to the design of the algorithms and extremely unstable as we can see in Figure (b)) While STORM-FW and SPIDER-FTE only focus on the upper-level problem, which leads to terrible results on the lower-level gap and recovery rate.
27


Reviews:
[{"review_id": "yC5z35KYVr", "reviewer": "Reviewer_FAuV", "summary": "This paper proposes methods for solving nonconvex/convex-convex simple bilevel optimization problems. The proposed methods are proven to have better non-asymptotic convergence guarantee than existing works. Numerical experiments are done on two applications to verify the effectiveness of the proposed methods.", "questions": "Please see the weaknesses section.", "limitations": "The fixed batch size limits the proposed methods from applying on large-scale datasets.", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "The main contribution of this work is the improved convergence rates for solving nonconvex/convex-convex simple bilevel optimization problems. The proposed methods use a cutting plane to approximate the solution set of the lower-level problem, of which the solution may not be unique. In order to have a more accurate approximation, the work proposes to use variance reduced estimators (STORM for stochastic setting and SPIDER for finite-sum setting) for gradient and function value estimations.", "weaknesses": "1. In the related work section, line 69~71 says that existing works on solving general bilevel optimization only consider strongly convex lower-level problems. This statement is not true. There have several works done on studying general bilevel optimization without strong convexity condition of lower-level problems. Just to name a few: (Liu et al., 2021), (Sow et al., 2022), (Shen et al., 2023), etc. Since these works are solving general bilevel optimization, they are potentially applicable to the simple bilevel optimization. I believe it is reasonable to provide a more detailed discussion on how these existing works are/aren't applicable to the target formulation of this draft.\n2. I notice that the algorithm SBCGI uses only one sample data per iteration, and the sample data batch size in algorithm SBCGF is fixed to $\\sqrt{n}$. In general, fixed batch size is often not very practical, especially for large dataset. I wonder if the authors have any insights on why the batch size cannot be set arbitrarily and be a dependency parameter in the convergence rate.\n3. Typo: in line 73, one->on.\n\n\n\nReferences:\n\nRisheng Liu, Yaohua Liu, Shangzhi Zeng, Jin Zhang, Towards Gradient-based Bilevel Optimization with Non-convex Followers and Beyond, 2021.\n\nDaouda Sow, Kaiyi Ji, Ziwei Guan, Yingbin Liang, A Primal-Dual Approach to Bilevel Optimization with Multiple Inner Minima, 2022.\n\nHan Shen, Quan Xiao, Tianyi Chen, On Penalty-based Bilevel Gradient Descent Method, 2023.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank reviewer FAuV for the comments and feedback!\n\n**Q1 There are several works on general bilevel optimization problems without lower-level strong convexity. The authors should provide more detailed discussions on how these existing works are/aren't applicable to the target formulation of this draft.** \n\n**A1** We thank the reviewer for pointing out these references. Indeed, there are several recent works on general bilevel optimization problems without lower-level strong convexity, and we will add them to the discussions in related work. \nSince simple bilevel problems can be viewed as a special case of general bilevel problems (where the lower-level objective is fixed), in principle these existing methods could be applied to our setting as well. However, there will be some technical issues, which we detail below.  \n\n- Some of the works [Liu et al., 2021; Sow et al., 2022; Chen et al., 2023; Lu and Mei, 2023] only consider deterministic bilevel problems, where one has access to the exact gradient of the upper-level and lower-level objectives. In comparison, in this paper we focus on stochastic bilevel problems and one of the main challenges is to properly control the stochastic noises present in our oracle.  \n\n- To the best of our knowledge, the only existing works considering stochastic bilevel problems without lower-level strong convexity are [Shen et al., 2023] and [Huang, 2023]. Note that the authors in both papers assume that the lower-level objective satisfies the PL inequality, while we assume that the lower-level objective is convex. Due to the different assumptions, their results are not directly comparable to ours. \n\n- Finally, since these works consider a more general and challenging class of problems, we argue that their theoretical results when applied to our setting are necessarily weaker. \n\n\nWe will add the above discussion to the related work section of our revised paper. Thanks for your comment.\n\nReferences:\n\n- Risheng Liu, Yaohua Liu, Shangzhi Zeng, Jin Zhang, Towards Gradient-based Bilevel Optimization with Non-convex Followers and Beyond, 2021.\n- Daouda Sow, Kaiyi Ji, Ziwei Guan, Yingbin Liang, A Primal-Dual Approach to Bilevel Optimization with Multiple Inner Minima, 2022.\n- Lesi Chen, Jing Xu, and Jingzhao Zhang, On Bilevel Optimization without Lower-level Strong Convexity, 2023.\n- Zhaosong Lu and Sanyou Mei, First-Order Penalty Methods for Bilevel Optimization, 2023.\n- Han Shen, Quan Xiao, Tianyi Chen, On Penalty-based Bilevel Gradient Descent Method, 2023.\n- Feihu Huang, On Momentum-based Gradient Methods for Bilevel Optimization with Nonconvex Lower-level, 2023.\n\n--- \n\n**Q2 SBCGI uses a fixed batch size of 1 and SBCGF uses a fixed batch size of $\\sqrt{n}$. Why the batch sizes cannot be set arbitrarily and be a parameter in the convergence rate?**\n\n**A2** \nNote that SBCGI can be implemented with a batch size as small as $b=1$. However, this does not imply that the stepsize \"has to be\" $b=1$. In other words, the main advantage of SBCGI, compared to SBCGF, is its capability to be implemented with any mini-batch size, even as small as $b=1$. Therefore, for SBCGI, the batch size can be set arbitrarily, whereas for SBCGF, it must be $\\sqrt{n}$. We will emphasize this point in the revised paper\n\n---\n\n**Typos.**\n\nThank you for pointing out this typo. We will fix all typos mentioned by all reviewers in the revised paper."}, {"review_id": "07GdxbrloL", "reviewer": "Reviewer_iBKW", "summary": "This paper introduces novel projection-free optimization methods for solving a class of *simple* stochastic bilevel optimization problems. The proposed methods are based on the conditional gradient algorithm and use variance reduction techniques to approximate the solution set of the lower-level problem. The paper provides theoretical guarantees for the convergence of the proposed methods and analyzes their sample complexity. The methods are shown to be efficient for both convex (global convergence) and non-convex upper-level (local convergence) functions.", "questions": "Questions: Can you describe how we can see this optimization problem from a game theory perspective?", "limitations": "NA", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper introduces a new framework for solving a class of stochastic bilevel optimization problems, which is an important and challenging problem in optimization. The proposed methods are shown to be efficient in experiments and have theoretical guarantees.", "weaknesses": "I am concerned that the paper makes no reference to the game theory literature on the topic (general-sum Stackelberg games).  Related work on the topic would make the work accessible to a broader audience.\n\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank reviewer iBKW for the comments and feedback!\n\n**Q1 Can you describe how we can see this optimization problem from a game theory perspective?**\n\n**A1** As the reviewer has correctly pointed out, a Stackelberg game (SCG) can be viewed as a general bilevel program in which a leader aims to maximize their own gain by anticipating and influencing the equilibrium state at which followers settle by engaging in a game [Liu et al., 2021; Li et al., 2022].\n\nHowever, in our paper, we are addressing a (stochastic) \"simple\" bilevel optimization problem that cannot be viewed as a two-player game. To be more precise, in a general bilevel optimization problem, the action of the leader player can influence the problem faced by the follower, since the lower-level objective function depends on a secondary variable. In contrast to this scenario, in simple bilevel optimization, the lower-level problem remains fixed and does not vary based on the actions of another player. As a result, it is not clear how one can view this optimization problem as a two-player game.\n\n\nWe will add the above discussion to the revised paper to better highlight the connection between our setting and game theory literature. \n\nReferences:\n\n- Risheng Liu, Jiaxin Gao, Jin Zhang, Deyu Meng, Zhouchen Lin, Investigating Bi-Level Optimization for Learning and Vision from a Unified Perspective: A Survey and Beyond, 2021.\n- Jiayang Li, Jing Yu, Qianni Wang, Boyi Liu, Zhaoran Wang, Yu (Marco) Nie, Differentiable Bilevel Programming for Stackelberg Congestion Games, 2022"}, {"review_id": "8Na8GxyI40", "reviewer": "Reviewer_Xgt8", "summary": "The paper presents a novel approach for solving stochastic simple bilevel optimization problems. The authors develop new methods that locally approximate the solution set of the lower-level problem using a stochastic cutting plane and then employ a conditional gradient update with variance reduction techniques to manage error from stochastic gradients.\n\nThe specific contributions are:\n- For the convex upper-level function case, the proposed method requires ${\\tilde{\\mathcal{O}}}(\\max ( 1/\\epsilon_f^2, 1/\\epsilon_g^2 ) )$ stochastic oracle queries to find a solution that is $\\epsilon_f$-optimal for the upper-level and $\\epsilon_g$-optimal for the lower-level. This is a substantial improvement over the previous complexity of $\\mathcal{O}(\\max(1/\\epsilon_f^4, 1/\\epsilon_g^4))$.\n- For the non-convex upper-level function case, the method requires at most ${\\tilde{\\mathcal{O}}}(\\max(1/\\epsilon_f^3, 1/\\epsilon_g^3))$ stochastic oracle queries to find an $(\\epsilon_f, \\epsilon_g)$-stationary point.\n- In the finite-sum setting, the method requires ${\\tilde{\\mathcal{O}}}(\\sqrt{n}/\\epsilon)$ and ${\\tilde{\\mathcal{O}}}(\\sqrt{n}/\\epsilon^2)$ stochastic oracle calls for the convex and non-convex settings, respectively, where $\\epsilon= \\min(\\epsilon_f,\\epsilon_g)$.\n", "questions": "\n- Line 15 and elsewhere: The authors discard the logarithmic terms in the complexity results. I suggest the authors use $\\tilde{\\mathcal{O}}$ to hide the logarithmic terms.\n- Line 46: The authors state that their approach is ``tight\". Why?\n- Example 1: The paper doesn't clarify where the element of randomness lies, nor does it explain why the problem specifically requires a stochastic simple bilevel optimization model.\n- Line 156, the formula (8): How to compute the exact $g(x_0)$?\n- Line 158: Why is there $\\hat{g}(x_0)$ instead of $g(x_0)$?\n- Algorithm 1, Line 7; Line 160, the formula (9);  Line 163 and elsewhere: The notations $\\mathcal{X}_t$, $\\hat{\\mathcal{X}}_t$, $\\mathcal{X}_t^{S}$ should be unified.\n- Line 177: The citation (7) should be (9).\n- Line 226, the formula (16): The term $\\sqrt{\\log(6/\\delta)}$ should be changed to $\\sqrt{\\log(6d/\\delta)}$.\n- The authors state that $\\omega=1$ and $\\omega=2/3$ are the ``best'' choice for the convex and nonconvex settings, respectively. How to ensure that the choices are the best?\n- Lemma 4.2 and Lemma 4.3: The authors could swap the places of these two lemmas and rewrite Lemma 4.2 to be more precise.\n- Line 238, the last sentence: add \"we\" before \"establish\".\n- In numerical experiments: Why do the authors use the number of samples be the abscissa axis? The paper doesn't include a comparative analysis of the time efficiencies across different algorithms, an omission that needs clarification.\n- Line 316: The matrix $D^{'*}$ is not clearly defined.\n- In Supplementary material: Many formulas are numbered but not cited, please consider removing these numbers. Sometimes the paper doesn't provide an explanation for the validity of two adjacent formulas. An appropriate clarification would be helpful for understanding.\n- Line 479: The $S$ in the formula should be changed to $s_2$.\n- Line 480 and elsewhere: The use of both \"t\" and \"T\" in the paper could lead to confusion. To avoid this, the authors could consider replacing one of them, such as changing \"T\" to \"M\".\n- Line 482 and elsewhere: The \"$\\alpha$\" should be replaced by \"$\\omega$\".\n- Line 526: the citation (47) is not correct.\n- Line 563, the formula (55) and elsewhere: The term $g(x_t)-g(x_0)$ should be $(g(x_t)-g(x_0))$.\n- Line 565, the formulas and elsewhere: The subscript $t+1$ in $g(x_{t+1})$ in the formula does not match the subscript $t$ in $g(x_t)$ in the theorem. This situation occurs many times in the proofs of theorems.\n- Line 567, the formula (57) and elsewhere: The Lipchitz constant $L_g$ should be replaced by $L_f$.\n- Line 677: The number $40$ should be replaced by $50$. \t\n", "limitations": "NA", "rating": 4, "confidence": 4, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "- Originality: In this paper, the idea of using cutting planes to approximate the optimal solution set for solving deterministic simple bilevel optimization is extended to the stochastic version. Due to the stochastic nature of the problem, this generalization is not trivial.\n- Quality: From a theoretical standpoint, the quality of the research is notable, with the authors demonstrating a substantial improvement over previous complexity results.\n- Clarity: The overall structure of the paper is clear.\n- Significance: This research provide lower complexity guarantees for stochastic simple bilevel optimization.", "weaknesses": "- The authors make an assumption that the constraint set $\\mathcal Z$ is compact and convex, which may not hold for many real applications. This strong assumption limits the method's applicability. \n\n- In their experiments, the authors employ the number of samples as the measure on the abscissa axis without providing a clear rationale for this choice. Furthermore, a comparative analysis of the time efficiency across various methods is missing, which should ideally be included to give a more comprehensive evaluation of their approach. \n\n- The paper contains numerous typographical errors which, although not affecting the technical content, detract from the overall readability and professionalism of the paper.\n\n- In Example 1, the authors do not clearly delineate how randomness is incorporated and why it necessitates the use of stochastic simple bilevel optimization as a modeling approach for this particular problem.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank reviewer Xgt8 for the comments and feedback!\n\n**Q1 The authors assume that the constraint set $\\mathcal{Z}$ is compact and convex, which may not hold for many real applications.** \n\n**A1** While this assumption may not apply to certain real-world scenarios, it remains a common assumption in many constrained simple bilevel optimization papers. These include works like [Jalilzadeh et al. (2023)] and [Jiang et al. (2023)], also in constrained general bilevel optimization papers, exemplified by [Liu et al. (2021)] and [Shen et al. (2023)]. Generally, assuming a compact constraint set is a prevalent practice when dealing with projection-free (Frank-Wolfe) methods, as demonstrated in works such as [Zhang et al. (2020)] and [Yurtsever et al. (2019)].  Since our work also employs a projection-free framework, it necessitates the assumption.\n\nReferences: \n\n- Jalilzadeh et al., Stochastic Approximation for Estimating the Price of Stability in Stochastic Nash Games, 2023\n- Jiang et al., A conditional gradient-based method for simple bilevel optimization with convex lower-level problem, 2023\n- Liu et al., Towards Gradient-based Bilevel Optimization with Non-convex Followers and Beyond, 2021\n- Shen et al., On Penalty-based Bilevel Gradient Descent Method, 2023\n- Zhang et al., One Sample Stochastic Frank-Wolfe, 2020\n- Yurtsever et al., Conditional Gradient Methods via Stochastic Path-Integrated Differential Estimator, 2019\n---\n**Q2 Why do the authors use the number of samples as the abscissa axis? The paper doesn't include a comparative analysis of the time efficiencies across different algorithms.**\n\n**A2** A widely used metric for comparing iterative methods in stochastic optimization is the number of required sample queries to achieve a specific accuracy. However, as the reviewer suggested, it is also possible to compare these methods based on their runtime if sample complexity is not the primary point of comparison. To address the reviewer's concern, we have included new plots in the attached PDF illustrating the comparison of the studied methods in terms of runtime. \n\n--- \n**Q3 In Example 1, the paper doesn't clarify where the element of randomness lies.**\n\n**A3** Example 1 presents a bilevel optimization problem, where the lower-level problem pertains to training loss, while the upper-level problem corresponds to testing loss. In both instances, we are confronted with an Empirical Risk Minimization loss, which can be regarded as stochastic optimization problems. This involves a uniform distribution across individual elements of the loss associated with each sample point. Therefore, this scenario can be viewed as a specific instance of the stochastic simple bilevel problem described in equation (1). We will clarify it in the revised paper. \n\n---\n**Q4 In Line 46, the authors state that their approach is \"tight\". Why?**\n\n**A4** It is essential to note that the sample complexity bounds of our proposed methods in both convex and nonconvex settings match the optimal sample complexity of \"stochastic single-level\" optimization problems. Consequently, these bounds are also tight for the more generalized setting of \"stochastic bilevel\" optimization problems considered in this paper. We will add the point to the revised paper. \n\n--- \n**Q5 How to compute the exact $g(x_0)$?**\n\n**A5** In the finite sum setting, we can exactly compute $g(x_{0})$. In the stochastic setting, we could use a batch of samples to estimate $g(x_{0})$ with high precision at the beginning. This additional operation will not affect the overall sample complexity of the proposed method, as the additional cost is negligible compared to the overall sample complexity. Specifically, we need to take a batch size of $b = \\tilde{\\mathcal{O}}(\\epsilon^{-2})$ to form the estimate $\\hat{g}(x_{0})$. By Hoeffding's inequality for subgaussian random variables, we have the following bound: $|\\hat{g}(x_0) - g(x_{0})|\\leq\\sqrt{2}\\sigma_{l}(T+1)^{-\\omega/2}\\sqrt{\\log(2/\\delta)}$, with a probability of at least $1-\\delta$, where $T$ is the maximum number of iterations. Comparing this with Lemma 4.1.3, we can further derive: $|\\hat{g}(x_0) - g(x_{0})|\\leq \\sqrt{2}\\sigma_{l}(T+1)^{ \\omega/2}\\sqrt{\\log(2/\\delta)} \\leq \\sqrt{2}(2L_{l}D+\\frac{3^{\\omega}}{3^{\\omega}-1}\\sigma_{l})(t+1)^{-\\omega/2}\\sqrt{\\log(6/\\delta})$, with a probability of at least $1-\\delta$ for all $0\\leq t \\leq T$. Consequently, the introduced error term would be absorbed in $K_{0, t}$ (Line 175) and will not affect any parts of the analysis. We will add the point to the revised paper.\n\n---\n**Q6 Line 158: Why is there $\\hat{g}(x_0)$ instead of $g(x_0)$?** \n\n**A6** This is a typo and it should be $g(x_0)$. We will fix it in the revised paper. \n\n---\n**Q7 Line 226, the formula (16): The term $\\sqrt{\\log(6/\\delta)}$ should be changed to $\\sqrt{\\log(6d/\\delta)}$.**\n\n**A7** Since function values are real values, the dimension is $d=1$ and hence the bound $\\sqrt{\\log(6d/\\delta)}$ can be simplified as $\\sqrt{\\log(6/\\delta)}$.\n\n---\n**Q8 The authors state that $\\omega = 1$ and $\\omega = 2/3$ are the ``best'' choices for the convex and nonconvex settings. How to ensure that the choices are the best?**\n\n**A8** What we meant to convey was that the best overall iteration and sample complexity, based on our analysis, for the convex setting is achieved when $\\omega = 1$. Moreover, again based on our analysis, the best upper bound for the nonconvex setting is obtained for $\\omega = 2/3$. We will clarify it in the revised paper. \n\n---\n**Q9 Line 316: The matrix $D'^\\*$ is not clearly defined.**\n\n**A9** As mentioned in Line 315-317, the matrix $\\tilde{D}^* \\in \\mathbb{R}^{25 \\times 50}$ contains all 50 true basis vectors and we generate matrices $D^* \\in \\mathbb{R}^{25 \\times 40}$ and $D'^\\* \\in \\mathbb{R}^{25 \\times 20}$ by randomly selecting $40$ and $20$ columns of  $\\tilde{D}^*$, respectively. \n\n---\n**Other typos.** Thanks for pointing out these typos. We will make sure to fix them in our revision."}, {"review_id": "HkDHFwTtcC", "reviewer": "Reviewer_oeu1", "summary": "A nice and clean result for stochastic simple bilevel optimization. The manuscript extends the algorithm design from the deterministic simple bilevel optimization with convex lower-level problem to the stochastic setting and demonstrate the sample complexity, which seems to be near-optimal (ignoring the log factor) when restricting the problem to classical stochastic optimization setting. ", "questions": "1. Throughout your discussion of related literature and contributions, please accurately denote if the reference or contribution pertains to simple stochastic bilevel optimization or the broader general stochastic bilevel optimization. This applies to numerous instances, such as Line 139 and Line 45.\n\n2. Regarding Line 103, is it $A\\subset \\mathbb{R}^m$ or $a_i\\in\\mathbb{R}^m$? There appears to be a typographical error.\n\n3. How are $p$ and $q$ determined in continual dictionary learning? Do the old dataset $A$ and the new dataset $A^\\prime$ have any correlation? If uncorrelated, why merge them? If correlated, would the results still hold?\n\n4. In the 'Random set for the subproblem' subsection, is there a missing $\\hat $ on $g(x_0)$? How is $g(x_0)$ known? This query also applies to equations (8)(9) and Algorithms 1 and 2.\n\n5. Lemma 4.1 seems to miss a space in \"c, (d is \". Additionally, why does the estimation error in Lemma 1 depend on the dimension? In other variance reduction papers, like STORM and SPIDER, the estimation error is generally dimension-independent.\n\n6. In Theorem 1 and 2, $\\mathcal{G}$ is not defined anywhere in the main context.\n\n7. My understanding of general bilevel optimization implies variance reduction aims for $O(\\epsilon^{-4})$ complexity bound in the nonconvex regime ($O(\\epsilon^{-3})$ with variance reduction). However, Chen et al., 2021 [NeurIPS, closing the gap, tighter analysis] indicated that variance reduction isn't essential for achieving optimal complexity in strong and convex settings. Moreover, for nonconvex settings, an $O(\\epsilon^{-4})$ rate can still be obtained without variance reduction. Therefore, I question whether variance reduction is genuinely necessary when $\\hat{\\nabla g}$ can be adaptively estimated using a moving average estimator, possibly achieving the same rate. This could be an avenue for future consideration.", "limitations": "The authors have adequately addressed the limitations and there is no societal impact involved as it is a pure theoretical work. ", "rating": 6, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "New results with good quality that are demonstrated in a nice and clean way. The motivation is very clear. It may push forward the study of general stochastic bilevel optimization when the lower-level is non-singleton. ", "weaknesses": "The manuscript has no major flaws. The idea to add a $K_t$ and to estimate g using samples are rather straightforward though. The motivation for using variance reduction are not well-explained. Using variance reduction for lower level is to ensure that decaying $K_t$ makes sense. This step I doubt if it is possible to only use moving average estimator [Chen et al., 2021 NeurIPS, Closing the Gap: Tighter Analysis of Alternating Stochastic Gradient Methods for Bilevel Problems] on the lower-level problem to achieve the same effect. Using variance reduction for upper level is to obtain the current complexity bound. Otherwise, using vanilla Frank-wolfe on the upper level will result in larger complexity. But nevertheless, the manuscript already provides a full story. So that the algorithm would requires less hyperparameter to tune. See my questions below.\n\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank reviewer oeu1 for the comments and feedback!\n\n**Q1 Why do we use STORM and SPIDER rather than a simpler moving average estimator?**\n\n**A1**  This is an excellent point. Since our framework utilizes a projection-free approach for stochastic bilevel optimization, we have endeavored to incorporate the best variance reduction techniques studied in the literature when applying the projection-free method to stochastic single-level optimization problems.\n\nIn fact, [Mokhtari et al. 2020] developed a stochastic Frank-Wolfe method using a moving average for stochastic optimization. However, their method's complexity was worse than the approaches achieved by unbiased variance reduction techniques like STORM [Zhang et al. 2020] or SPIDER [Yurtsever et al. 2019]. Given that, for stochastic single-level projection-free methods, the use of a biased moving average estimator led to suboptimal guarantees, we decided to employ unbiased variance reduction techniques such as STORM and SPIDER.\n\nThat being said, we concur with the reviewer's observation that utilizing a simpler moving average estimator with fewer parameters to tune would be advantageous. Nevertheless, we are not aware of any analysis, even in the context of single-level stochastic settings, that matches the bounds established by [Zhang et al. 2020] and [Yurtsever et al. 2019] for projection-free methods.  Perhaps the fact that  [Chen et al. 2021] considered an unconstrained problem was crucial in allowing them to establish a strong complexity bound for the moving average estimator.  Exploring the feasibility of attaining similar results for constrained problems using a projection-free method remains an intriguing avenue of research that warrants further investigation.\n\nReferences:\n\n- Mokhtari et al., Stochastic Conditional Gradient Methods: From Convex Minimization to Submodular Maximization, 2020.\n- Zhang et al., One Sample Stochastic Frank-Wolfe, 2020\n- Yurtsever et al., Conditional Gradient Methods via Stochastic Path-Integrated Differential Estimator, 2019\n\n---\n\n**Q2 How are $p$ and $q$ determined in continual dictionary learning? Do the old dataset $A$ and the new dataset $A'$ have any correlation? If uncorrelated, why merge them? If correlated, would the results still hold?**\n\n**A2** We used synthetic data for this problem, with $p$ and $q$ able to take any values, provided that $p<q$. The old dataset $A$ and the new dataset $A'$ are both generated from the same source, which in our example is a uniform distribution. The primary concept behind lifelong learning is that we only have access to a portion of the data at the outset, with new data gradually becoming available over time. The correlation between the datasets is not a critical factor.\n\n---\n\n**Q3 How is $g(x_{0})$ known?**\n\n**A3** In the finite sum setting, we can exactly compute $g(x_{0})$, and the additional cost of $n$ function evaluations will be dominated by the overall complexity.\n\nIn the stochastic setting, we could use a large batch of samples to estimate $g(x_{0})$ with high precision at the beginning of the process. This additional operation will not affect the overall sample complexity of the proposed method, as the additional cost is negligible compared to the overall sample complexity.\nSpecifically, we need to take a batch size of $b = \\tilde{\\mathcal{O}}(\\epsilon^{-2})$ to form the estimate $\\hat{g}(x_{0})$. Using Hoeffding's inequality for subgaussian random variables, we have the following bound:\n$|\\hat{g}(x_0) - g(x_{0})| \\leq \\sqrt{2}\\sigma_{l}(T+1)^{-\\omega/2}\\sqrt{\\log(2/\\delta)}$, with a probability of at least $1-\\delta$, where $T$ is the maximum number of iterations. Comparing this with Lemma 4.1.3, we can further derive:\n$   \n|\\hat{g}(x_0) - g(x_{0})| \\leq \\sqrt{2}\\sigma_{l}(T+1)^{-\\omega/2}\\sqrt{\\log(2/\\delta)}\n\\leq \\sqrt{2}(2L_{l}D+\\frac{3^{\\omega}}{3^{\\omega}-1}\\sigma_{l})(t+1)^{-\\omega/2}\\sqrt{\\log(6/\\delta}) \n$, with a probability of at least $1-\\delta$ for all $0\\leq t \\leq T$. Consequently, the introduced error term would be absorbed in $K_{0, t}$ (Line 175) and will not affect any parts of the analysis. We will add this point to the revised paper. Thanks for your comment.\n\n---\n\n**Q4 Why does the estimation error in Lemma 4.1 depend on the dimension? In other variance reduction papers, like STORM and SPIDER, the estimation error is generally dimension-independent.**\n\n**A4** The main reason our bounds depend on dimension is due to the fact that we provide high-probability bounds. This contrasts with most works that establish guarantees in expectation, which are dimension independent. More precisely, our results demonstrate convergence in high probability under the assumption of sub-gaussian noise (Assumptions 2.2.3 and 2.3.3 in our paper).\n\nIt is worth noting that there are other works that provide guarantees with high probability and are dimension independent; however, they require stronger assumptions, such as the boundedness of random variables with probability 1, as seen in [Fang et al. (2018)] and [Xie et al. (2019)]. In Appendix E of our paper, Proposition E.1 is utilized to prove results under the assumption of bounded random variables, while Proposition E.2 is used to prove results under the sub-gaussian noise assumption.\n\nReferences:\n\n- Fang et al., SPIDER: Near-Optimal Non-Convex Optimization via Stochastic Path Integrated Differential Estimator, 2018.\n- Xie et al., Efficient Projection-Free Online Methods with Stochastic Recursive Gradient, 2019.\n---\n\n**Q5 In Theorem 1 and 2, $\\mathcal{G}$ is not defined anywhere in the main context.**\n\n**A5** Thanks for your comment. The definition of $\\mathcal{G}(\\hat{x})$ is provided in Definition 2.1 which states that \n  $  \\mathcal{G}(\\hat{x}) \\triangleq \\max _{s \\in \\mathcal{X}_g^*}\\{\\langle\\nabla f(\\hat{x}), \\hat{x}-s\\rangle\\}$. We will add a pointer to this definition in the statements of both theorems.\n\n---\n\n**Other typos.**\nThank the reviewer for pointing out the typos. We will fix them."}]

Author rebuttal:
In response to Reviewer Xgt8, we included some additional plots in the attached PDF file comparing the studied methods in terms of runtime.

Meta-review: There was a bit of disagreement between the reviewers on this submission, but overall, the interest of this submission was underlined by all, and the reviewer in favor of the rejection did not participate in the discussion.
Decision: Accept (poster)