Provably Faster Algorithms for Bilevel Optimization via Without-Replacement Sampling

Authors: junyi  li
heng  huang
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 05:25:05

Abstract: Bilevel Optimization has experienced significant advancements recently with the introduction of new efficient algorithms. Mirroring the success in single-level optimization, stochastic gradient-based algorithms are widely used in bilevel optimization. However, a common limitation in these algorithms is the presumption of independent sampling, which can lead to increased computational costs due to the unique hyper-gradient structure in bilevel problems. To address this challenge, we study the example-selection strategy for bilevel optimization in this work. More specifically, we introduce a without-replacement sampling based algorithm which achieves a faster convergence rate compared to its counterparts that rely on independent sampling. Beyond the standard bilevel optimization formulation, we extend our discussion to conditional bilevel optimization and also two special cases: minimax and compositional optimization. Finally, we validate our algorithms over both synthetic and real-world applications. Numerical results clearly showcase the superiority of our algorithms.

Main paper: 
 introduction: 
 Bilevel optimization [51, 47] has received a lot of interest recently due to its wide-ranging applicability in machine learning tasks, including hyper-parameter optimization [39] , meta-learning [57] and personalized federated learning [46] . A bilevel optimization problem is a type of nested two-level problems as follows:EQUATIONwhich includes an outer problem f (x, y) and a inner problem g(x, y). The outer problem f (x, y) relies on the solution y x of the inner problem g(x, y). Eq. ( 1) can be solved through gradient descent:x t+1 = x t -η∇h(x t ), with η be the stepsize and ∇h(x t ) be gradient at the state x t . Specially, the gradient ∇h(x) [17] is a function of the inner problem minimizer y x , first order derivatives ∇ x f (x, y), ∇ y f (x, y), ∇ y g(x, y) and second order derivatives ∇ xy g(x, y), ∇ y 2 g(x, y). In particular, y x can be solved though gradient descent: y t+1 = y t -γ∇ y g(x, y t ) with γ be the stepsize and ∇ y g(x t , y t ) be gradient at the iterate (x t , y t ). In the standard setup, the outer and inner problems are defined over two datasets D u = {ξ i , i ∈ [m]} and D l = {ζ j , i ∈ [n]}, respectively:f (x, y) = 1 m m i=1f (x, y; ξ i ), g(x, y) = 1 n n j=1 g(x, y; ζ j ),Naturally, we sample from the datasets to estimate the first and second-order derivatives in the hyper-gradient formulation. And to guarantee convergence, previous approaches require the examples Table 1 : Comparisons of the Bilevel Opt. & Conditional Bilevel Opt. algorithms for finding an ϵ-stationary point. We include comparison with stochastic gradient descent type of methods and a more comprehensive comparison including other acceleration methods can be found in Table 2 of Appendix. An ϵ-stationary point is defined as ∥∇h(x)∥ ≤ ϵ. Gc(f, ϵ) and Gc(g, ϵ) denote the number of gradient evaluations w.r.t. f (x, y) and g(x, y); JV (g, ϵ) denotes the number of Jacobianvector products; HV (g, ϵ) is the number of Hessian-vector products. m and n are the number of data examples for the outer and inner problems, in particular, n is the maximum number of inner problems for conditional bilevel optimization. Our methods have a dependence over example numbers (max(m, n)) q , q is a value decided by without-replacement sampling strategy and can have value in [0, 1] (A herding-based permutation [33] can let q = 0).be mutually independent [25] even at the same state (x, y). For example, two different examples ζ 1 and ζ 2 are sampled to estimate ∇ y g(x, y) and ∇ y 2 g(x, y), respectively. This leads to extra computational cost in practice. Indeed, two backward passes are required to evaluate ∇ x f (x, y) and ∇ y f (x, y), while five backward passes are needed to evaluate ∇ y g(x, y), ∇ y 2 g(x, y) and ∇ xy g(x, y) for a given state of (x, y). In contrast, one backward pass are needed to evaluate the former and two backward passes for the latter if not sampling independently for each property. This leads to notable computational cost difference for the current large-scale machine learning models [4] . Therefore, it is beneficial to explore other example-selection strategies beyond independent sampling. More specifically, we aim to answer the following question in this work: Can we solve the bilevel optimization problem without using independent sampling?Although example-selection is under-explored for bilevel optimization, it has been well studied in the single level optimization literature [2, 42, 44] . For single level problems: min x∈R p f (x) =1 n n i=1 f (x; ξ i ), we iteratively perform the stochastic gradient step x t+1 = x t -η∇f (x t , ξ t ), where the sample gradient is used to estimate the full gradient and the example ξ t is sampled according to some order. More specifically, the sampling schemes are roughly divided into two categories: withreplacement sampling and without-replacement sampling. Traditional stochastic gradient descent typically employs with-replacement sampling to ensure an unbiased estimation of the full gradient in each iteration. In contrast, the without-replacement sampling has biased estimation at each iteration. Despite this, when averaging the gradients of consecutive examples, without-replacement sampling can approximate the full gradient more efficiently [34] . In fact, any permutation-based withoutreplacement sampling has that the gradient estimation error contracts at rate of O(K -2 ) [34] where K is the number of steps, while stochastic gradient descent only has a rate of O(K -1 ). Inspired by the properties of without-replacement sampling demonstrated in the single-level optimization, we design the algorithm named WiOR-BO (and its variants) which performs without-replacement sampling for bilevel optimization. Naturally, WiOR-BO does not require the independent sampling as in previous methods. In fact, our algorithm enjoys favorable convergence rate. As shown in Table 1 , WiOR-BO converges faster compared with their counterparts (e.g. StocBiO [25] ) using independent sampling. Actually, WiOR-BO gets comparable convergence rate with methods using more complicated variance reduction techniques (e.g. MRBO [52] ). In Table 2 , we perform a more comprehensive comparison with other variance reduction techniques and in Section C of Appendix, we further compare WiOR-BO with these methods.The contributions of our work can be summarized as:1. We propose WiOR-BO which performs without-sampling to solve the bilevel optimization problem and converges to an ϵ stationary point with a rate of O(ϵ -3 ). This rate improves over the O(ϵ -4 ) rate of its counterparts (e.g. stocBiO [25] ) using independent sampling. 2. We propose WiOR-CBO for the conditional bilevel optimization problems and show that our algorithm converges to an ϵ stationary point with a rate of O(ϵ -4 ). This rate improves over the O(ϵ -6 ) rate of its counterparts (e.g. DL-SGD [23] ) using independent-sampling.3. We customize our algorithms to the special cases of minimax and compositional optimization problems and demonstrate similar favorable convergence rates.Notations. ∇ (∇ x ) denotes full gradient (partial gradient, over variable x), higher order derivatives follow similar rules.[n] represents sequence of integers from 1 to n. O(•) is the big O notation, and we hide logarithmic terms. 
related works: 
 Bilevel Optimization has its roots tracing back to the 1960s, beginning with the regularization method proposed by [51] , and then followed by many research works [14, 47, 43, 55] . In the field of machine learning, similar implicit differentiation techniques were used to solve Hyper-parameter Optimization [28, 6, 12] . Exact solutions for Bilevel Optimization solve the inner problem for each outer variable but are inefficient. More efficient algorithms solve the inner problem with a fixed number of steps, and use the 'back-propagation through time' technique to compute the hypergradient [13, 35, 16, 40, 45] . Recently, single loop algorithms [17, 22, 25, 26, 52, 9, 30, 24] are introduced to perform the outer and inner updates alternatively. Besides, other aspects of the bilevel optimization are also investigated: [23] studied a type of conditional bilevel optimization where the inner problem depends on the example of the outer problem; [53, 56] studied the Hessian-free bilevel optimization and proposed an algorithm with optimal convergence rate; [32] proposed an efficient single loop algorithm for the general non-convex bilevel optimization; [54, 31] studies the bilevel optimization in the federated learning setting.Sample-selection in stochastic optimization has been extensively studied for the case of single level optimization problems [2, 3, 19, 38] . Beyond the with-replacement sampling adopted by SGD, various without-replacement strategies were also studied in the literature. The random-reshuffling strategy shuffles the data samples at the start of each training epoch, and its property was first studied in [42] via the non-commutative arithmetic-geometric mean conjecture and followed by [20, 19] . [44, 36] studied the shuffling-once strategy, where the data samples are only shuffled at the start of the training once; [8] investigated data echoing where one data sample is repeatedly used. Quasi-Monte Carlo sampling is also applied in stochastic gradient descent [5, 34] . Recently, [34, 33] proposed to use average gradient error to study the convergence of various example order strategies and then proposed a herding-based sampling strategy [50] to proactively minimize the average gradient error. 
without-replacement sampling in bilevel optimization: 
 In this work, we focus on the following finite-sum bilevel optimization problem:min x∈R p h(x) := f (x, y x ) = 1 m m i=1 f (x, y x ; ξ i ), s.t. y x = arg min y∈R d g(x, y) = 1 n n j=1 g(x, y; ζ j ) (2)where both the outer and the inner problems are defined over a finite dataset, and we denote them asD u = {ξ i , i ∈ [m]} and D l = {ζ j , j ∈ [n]}, respectively.Under mild assumptions, the hypergradient of h(x) is:EQUATIONOur objective is to minimize h(x) by exploiting Eq. ( 3) with examples from D u and D l at each step. More specifically, for a given example order π denoting the following example-selection sequence as {ξ π t ∈ D u } and {ζ π t ∈ D l } for t ∈ [T ], where T denotes the length of the sequence and can be arbitrarily large. We then perform the following update steps iteratively for t ∈ [T ]:EQUATIONwhere Eq. ( 4) adopts the single-loop [30] update. More specifically, Eq. (4a) performs gradient descent to estimate the minimizer y x ; Eq. (4b) performs gradient descent to estimate u x defined in Eq. ( 3); and Eq. (4c) perform the gradient descent step given the estimation of y x and u x . for i = 0 to I -1 do 6:y r i+1 = y r i -γ r i ∇ y g(x r i , y r i ; ζ π i ); u r i+1 = u r i -ρ r i (∇ y 2 g(x r i , y r i ; ζ π i )u r i -∇ y f (x r i , y r i ; ξ π i )); 7: x r i+1 = x r i -η r i (∇ x f (x r i , y r i ; ξ π i ) -∇ xy g(x r i , y r i ; ζ π i )u r i ); 8:end for For large-scale machine learning models where back-propagation is expensive, our method can lead to significant computation saving in practice.Equation ( 4) is independent of the specific order of samples. It can accommodate both withreplacement and without-replacement sampling strategies. In particular, we consider two most widely-used without-replacement sampling methods: random-reshuffling and shuffle-once. For random-reshuffling, we set T = R × lcm(m, n), where lcm(m, n) denotes the least common multiple of m and n and R is a constant. More specifically, we generate R × lcm(m, n)/m random permutations of examples in D u then concatenate them to get {ξ π t ∈ D u }, and we follow similar procedure to generate {ζ π t ∈ D u }. In algorithm 1, we show our WiOR-BO algorithm for solving Eq. ( 2) with random-reshuffling sampling. Note that in Line 4 of Algorithm 1, P ι (•) is the projection operation to an ι-size ball and we perform projection to make u t be bounded during updates. Shuffleonce is another widely used without-replacement sampling strategy, where a single permutation is generated and reused for both D u and D l , and we can modify Line 3 of Algorithm 1 to incorporate shuffle-once. Finally, compared to independent sampling-based algorithm such as stocBiO [25] , we require extra space to generate and store the orders of examples, but this cost is negligible compared to memory and computational cost of training large scale models.In Eq. ( 2), we assume the outer dataset D u and the inner dataset D l are independent with each other. However, it is also common in practice that the inner problem not only relies on the outer variable x, but also the current outer example. This type of problems is known as conditional (contextual) bilevel optimization problems [23] in the literature and we consider the finite setting as follows:EQUATIONNote that in Eq. ( 5), the outer problem is defined on a finite dataset D u = {ξ i , i ∈ [m]}, and for each example ξ i ∈ D u , we have a inner problem g (ξi) (x, y) defined on a dataset D l,ξi = {ζ ξi,j , j ∈ [n ξi ]}.The hyper-gradient of Eq. ( 5) is:EQUATIONThe key difference of Eq. ( 6) with Eq. ( 3) is that y ξi x and u ξi x are defined for each example ξ i ∈ D u . To minimize h(x) in Eq. ( 5), we assume a sample order π denoting the following example selection sequence as {ξ π t ∈ D u , t ∈ [T ]}, where T denotes the example sequence length and can be arbitrarily Randomly sample a permutation of outer dataset to have {ξ π i , i ∈ [m]} and set x r 0 = x r-1 m or x 0 if r = 1.for i = 0 to m -1 do 5:Initial states y 0 and u 0 , learning rates {γ s j , ρ s j }, j ∈ [n i ], s ∈ [S].for s = 1 to S do 7:Sample a permutation of inner dataset to have {ζ π ξi,j , j ∈ [n i ]} , set y s 0 = y s-1 ni and u s 0 = P ι (u s-1 ni ) or y s 0 = y 0 and u s 0 = u 0 if s = 1.for j = 0 to n i -1 do 9:y s j+1 = y s j -γ s j ∇ y g(x r i , y s j ; ζ π ξi,j ); 10:u s j+1 = u s j -ρ s j (∇ y 2 g(x r i , y s j ; ζ π ξi,j )u s j -∇ y f (x r i , y s j ; ξ π i )); 11:end for 12:end for 13:x r i+1 = x r i -η r i (∇ x f (x r i , y S ni ; ξ π i ) -∇ xy g ξ π i (x r i , y S ni )u S ni ); 14:end for 15: end for large. We then perform the following update steps iteratively:EQUATIONwhere ŷt and ût are estimations of yξ π t xt u ξ π t xt .We get them by performing the following rules T l steps over the sample order {ζ πEQUATIONwhere both the outer variable state x t and sample ξ π t are fixed. We then set ŷt = y T l and ût = u T l in Eq. (7) . Note that Eq. (7) and Eq. ( 8) perform a double-loop update rule [23] . Since the inner problem relies on the example of the outer problem, we compute y ξi x and u ξi x for each new state x and sample ξ i instead of reusing them as in the single loop update of Eq. (4).Note that an important feature of our algorithm is that we select samples based on an example order π at each step instead of sampling independently. Similar to the unconditional bilevel optimization, various example orders can be incorporated. For random-reshuffling, we set T = Rm with R be a constant. Then we generate R random permutations of examples in D u and concatenate them to get the example order {ξ π t ∈ D u , t ∈ [T ]}. Similarly, we concatenate S random permutations of the inner dataset D l,ξt to get {ζ π ξ π t ,t l ∈ D l,ξ π t , t l ∈ [T l ]}} with a sequence length T l = S × n ξi . We summarize this in Algorithm 2 and denote it as WiOR-CBO. We slightly abuse the notation to replace n ξi with n i in the inner loop update for better clarity.Our discussion of without-replacement sampling for bilevel optimization can be customized for two important nested optimization problems: compositional optimization and minimax optimization problems, as they can be viewed as a special type of bilevel optimization problem. Suppose we let the outer problem only rely on y, set the inner problem as g(x, y) = 1 2 ∥y -r(x)∥ 2 . Consider the finite case of f (y) and r(x), we have the following compositional optimization problem:EQUATIONThen we can derive the hyper-gradient as: ∇h(x) = ∇ x r(x)u x , where u x = ∇ y f (r(x)), meanwhile, we have ∇ y g(x, y) = y -r(x). Then we can instantiate Algorithm 1 to get an algorithm for compositional optimization problem using without-replacement sampling (Algorithm 3 of Appendix). Similar to the SCGD algorithm [49] , we track the exponential moving average of the inner state y, besides we also track the exponential moving average of inner gradient (the update rule of u in Algorithm 3).Next, similar to the unconditional case, the conditional compositional optimization problem can be specialized from the conditional bilevel optimization problem Eq. ( 5) to have:EQUATIONand we can instantiate Algorithm 2 in the special case of Eq. ( 10) to get Algorithm 5 (in the Appendix).The minimax optimization can also be viewed as a special type of bilevel optimization. More specifically, if we have g(x, y) = -f (x, y) in Eq. ( 2), then we get the following minimax problem:EQUATIONIn the special case of Eq. ( 11), the hyper-gradient ∇h(x) in Eq. ( 3) degenerates to ∇h(x) = ∇ x f (x, y x ) due to the fact that ∇ y f (x, y x ) = -∇ y g(x, y x ) = 0 by the optimality condition. Meanwhile, we have ∇ y g(x, y) = -∇ y f (x, y). As a result, we can simplify Algorithm 1 to get an alternative update algorithm for minimax optimization using without-replacement sampling (Algorithm 4 in the Appendix). Note that due to the linearity of minimax problem w.r.t. examples, the conditional case for minimax optimization degenerates to the unconditional case. [11, 7] also considers the without-replacement sampling for the minimax problem, however, they focus on the strongly-convex-strongly-concave and two-sided Polyak-Łojasiewicz settings, in contrast, our WiOR-MiniMax can be applied to nonconvex-strongly-concave setting. 
theoretical analysis: 
 In this section, we provide theoretical analysis of our algorithms. We first state some assumptions: ζ i ) - ∇ 2 g(x, y)∥ ≤ A.As stated in The assumption 4.1 and 4.2, we consider the non-convex-strongly-convex bilevel optimization problems, this class of problems is widely studied in bilevel optimization literature [17, 25] .Assumption 4.3 and 4.4 guarantee that the example gradients are good estimation of full gradient and is widely used in the analysis of single-level finite-sum optimization problems [36] . Note that as the hyper-gradient (Eq. 3) involves second order derivatives, Assumption 4.4 also requires that the second order example gradients have bounded bias.We next make the following assumptions about the example selection sequence. We assume{ξ π t ∈ D u , t ∈ [T ]} and {ζ π t ∈ D l , t ∈ [T ]} satisfy: Assumption 4.5. Given sequences of samples: {ξ π t ∈ D u , t ∈ [T ]} and {ζ π t ∈ D l , t ∈ [T ]}, we assume there exists some constants α, C such that for any step t and any interval k > 0, we have:1 k t+k-1 τ =t ∇f (x t , y t ; ξ π τ ) -∇f (x t , y t ) 2 ≤ k -α C 2and similar bounds hold for ∇g and ∇ 2 g (See Appendix B.1 for the full version of this assumption.) Assumption 4.5 measures how well the sample gradients approximate the full gradient over an interval on average. Similar assumptions are used to study the sample order in single level optimization problems [34] . If samples are selected independently at each step, we get α = 1 and C = A based on Assumption 4.3 and 4.4. For any permutation-based order (e.g. random-reshuffling and shuffle-once), we can show the assumption holds with α = 2 and C = max(m, n) × A. Note that the C depends on the number of samples m (n), this dependence can be reduced to (max(m, n)) 0.5 for sufficiently small learning rates [34, 37] . More recently, [33] proposed a herding-based algorithm which explicitly optimizes the gradient error in Assumption 4.5 and shows that we can reach α = 2 and C = O(A).Bilevel Optimization. We are ready to show the convergence of Algorithm 1. Firstly, we denote the following potential functionG t := G t = h(x t ) + ϕ y ∥y t -y xt ∥ 2 + ϕ u ∥u t -u xt ∥ 2, where ϕ y and ϕ u are constants that relates to the smoothness parameters of h(x). Then we have: Theorem 4.6. Suppose Assumptions 4.1-4.4 are satisfied, and Assumption 4.5 is satisfied with some α and C for the example order we use in Algorithm 1. We choose learning rates , we have:η t = η, γ t = c 1 η, ρ t = c1 E E-1 e=0 ∥∇h(x ke )∥ 2 ≤ 2∆ kEη + 2 L2 k -α C 2where c 1 , c 2 , L0 , L are constants related to the smoothness parameters of h(x), ∆ is the initial sub-optimality, R and I are defined in Algorithm 1.To reach an ϵ stationary point, we choose η = 1 k L0, and for (E, k), we choose:E = 4 L0∆ ϵ 2 and k = 4 L2 C 2 ϵ 2 1/α, which means we need to choose the number of epochsR = T I = 4 L0(4 L2 C 2 ) 1/α ∆ Iϵ 2+2/α. Specially, if we choose examples independently at each step, we have C = A and α = 1, then we need T = O(ϵ -4 ) steps to reach an ϵ-stationary point. This recovers the rate of stocBiO algorithm [25] and SOBA [9] . Next, if we use some permutation-based without-replacement sampling strategy witch has C = (max(m, n)) q × A (q ∈ [0, 1]) and α = 2, then we need T = O((max(m, n)) q ϵ -3 ). In particular, by adopting a herding-based permutation as in [33] , we get q = 0, and our WiOR-BO strictly outperforms the independent sampling-based stocBiO and SOBA algorithms.Conditional Bilevel Optimization. Next, we study the convergence rate of Algorithm 2. Besides Assumptions 4.3 and 4.4, we need the bias of sample hyper-gradient be bounded too: Assumption 4.7. For all examples ξ i ∈ D u and states x ∈ R p , there exists an upper bound A such that the sample gradient errors satisfy ∥∇h(x; ξ i ) -∇h(x)∥ ≤ A.We also augment Assumption 4.5 to include the case of ∇h(x) (Seem Appendix B.1 for more details.). Then we are ready to show the convergence rate, we use the potential function G t = h(x t ) and have: Theorem 4.8. Suppose Assumptions 4.1-4.4 and 4.7 are satisfied, and Assumption 4.5 is satisfied with some α and C for the example order we use in Algorithm 2. We choose learning ratesη t = η = 1 8k L , γ t = γ = 1 256kLκ , ρ t = ρ = 1 512kLκand denote T = R × m be the total number of outer steps and T l = S × max(n i ) be the maximum inner steps. Then for any pair of values (E, k) which has T = E × k and T l = E l × k, we have:1 E E-1 e=0 ∥∇h(x ke )∥ 2 ≤ 8∆ h kEη + C u (1 - kµρ 2 ) E l ∆ u + C y (1 - kµγ 2 ) E l ∆ y + (C ′ ) 2 C 2 k -αwhere C u , C y , C ′ , L are constants related to the smoothness parameters of h(x), ∆ h , ∆ u , ∆ y are the initial sub-optimality, R and S are defined in Algorithm 2.To reach an ϵ stationary point, we choose:E = 128 L∆ ϵ 2 , and , k = 2(C ′ ) 2 C 2 ϵ 2 1/α and E l = O(log(ϵ -1 )). Then the sample complexity of Algorithm 2 is EE l k 2 = O(C 4/α ϵ -(2+4/α) ), where we omit the logarithmic term E l . Specially, if we choose examples independently at each step, then we have sample complexity O(ϵ -6 ) steps for an ϵ-stationary point. This recovers the rate of DL-SGD algorithm [23] . Next, if we use some permutation-based without-replacement sampling strategy, then we have sample complexity O((max(m, n)) 2q ϵ -4 ), which can match the state-of-art algorithm RT-MLMC [23] if we choose an example order with q = 0 [33] .  
applications and numerical experiments: 
 In this section, we verify the effectiveness of the proposed algorithm through a synthetic invariant risk minimization task and two real world bilevel tasks: Hyper-Data Cleaning and Hyper-representation Learning. The code is written in Pytorch. Our experiments were conducted on servers equipped with 8 NVIDIA A5000 GPUs. Experiments are averaged over five independent runs, and the standard deviation is indicated by the shaded area in the plots. In this section, we consider a synthetic invariant risk minimization task [21] . More specifically, we perform logistic regression over a set of examples{(c i , b i ), i ∈ [M ]} with cbe the input and b be the target. However, we does not have access to c but a set of noisy observations {c j,i , j ∈ n i } for each c i . To learn the coefficient x between c and b, we optimize the following objective:min x∈R p 1 m m i=1 log(1 + exp(-b i y (i) x )), s.t. y (i) x = arg min y∈R d 1 n i ni j=1 1 2 ∥y -c j,i x∥ 2This is a conditional bilevel optimization problem and we can solve it using our WiOR-CBO. We generate synthetic data to test the effects of different sample order (sample generation process is described in Appendix A). More specifically, we compare independent-sampling (WiR-CBO), shuffleonce (WiOR-CBO (SO)) and random-reshuffling (WiOR-CBO (RR)). As shown in Figure 1 , the two without-replacement sampling methods outperforms the independent sampling one.In the Hyper-Data Cleaning task, we are given noisy training dataset whose labels are corrupted by noise, meanwhile, we have a validation set whose labels are not corrupted. Our target is then using the clean validation samples to identify corrupted samples in the training set. More specifically, we learn optimal weights for training samples such that a model learned over the weighted training set performs well on the validation set. We can formulate this task as a bilevel optimization problem (Eq. ( 2)). A mathematical formulation of the task is included in Appendix A.Dataset and Baselines. We construct datasets based on MNIST [29] . For the training set, we randomly sample 40000 images from the original training dataset and then randomly perturb a fraction of labels of samples. For the validation set, we randomly select 5000 clean images from the original training dataset. In our experiments, we test our WiOR-BO algorithm (Algorithm 1), with the shuffle-once (WiOR-BO-SO) and random-reshuffling (WiOR-BO-RR) sampling; additionally, we also consider the following non-adaptive bilevel algorithms as baselines: reverse [16] , stocBiO [25] , BSA [17] , AID-CG [18] , MRBO [52] , VRBO [52] and WiR-BO (the variant of our WiOR-BO using independent sampling, which is similar to the single loop algorithms such as SOBA [9] , AmIGO [1] and FLSA [30] ). We perform grid search for hyper-parameters and report the best results.We summarize the results in Figure 2 . In the figure, we compare the validation loss and F1 score w.r.t. both Hyper-iterations and running time. As shown in the figure, the two variants of our algorithm WiOR-BO-SO and WiOR-BO-RR get similar performance and they both outperform other baselines.In particular, note that WiR-BO differs with WiOR-BO-SO (RR) only over the example order, but the In this section, we consider Hyper-representation learning task. In this task, we learn a hyperrepresentation of the data such that a linear classifier can be learned quickly with a small number of samples. We consider the Omniglot [27] and MiniImageNet [41] data sets. This task can be viewed as a conditional bilevel optimization problem (Eq. ( 2)) and a mathematical formulation of the task is included in Appendix A.Dataset and Baselines. The details of the datasets are included in Appendix A and we consider N-way-K-shot classification task following [48] . In our experiments, we test our WiOR-CBO (Algorithm 2), using the shuffle-once (WiOR-BO-SO) and random-reshuffling (WiOR-BO-RR) sampling. Besides, we compare with the following baselines: DL-SGD [23] and RT-MLMC [23] . Note that our WiOR-CBO can also use independent sampling, and it has similar performance as DL-SGD. We perform grid search of hyper-parameters for each method and report the best results.We summarize the experimental results for the Omniglot dataset in Figure 3 and we defer the results for MiniImageNet to the Appendix A. As shown in the figure, our WiOR-CBO SO/RR outperforms DL-SGD and is comparable to the state-of-the-art algorithm RT-MLMC. 
conclusion: 
 In this work, we investigated example-selection of bilevel optimization. Beyond the classical independent sampling, we assumed an example-order based on without-replacement sampling, such as random-reshuffling and shuffle-once. We proposed the WiOR-BO algorithm for bilevel optimization and show the algorithm converges to an ϵ-stationary point with rate O(ϵ -3 ). After that, we also discussed the conditional bilevel optimization problems and introduced an algorithm with convergence rate of O(ϵ -4 ). As special cases of bilevel optimization, we studied the minimax and compositional optimization problems. Finally, we validated the efficacy of our algorithms through one synthetic and two real-world tasks; the numerical results show the superiority of our algorithms.Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?Answer: [Yes] Justification: Most of the experimental details can be found in Section 5, while the rest is included in the Appendix A.Guidelines:• The answer NA means that the paper does not include experiments.• The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them. • The full details can be provided either with the code, in appendix, or as supplemental material. 
experiment statistical significance: 
 Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?Answer: [Yes] Justification: The numerical experiments in the paper are averaged over 5 independent runs and we report the mean and standard deviation (shadowed area).Guidelines:• The answer NA means that the paper does not include experiments.• The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper.• The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions). • The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors). • It should be clear whether the error bar is the standard deviation or the standard error of the mean. Justification: We state the compute resource to run our experiments at the start of Section 5.Guidelines:• The answer NA means that the paper does not include experiments.• The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage. • The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute. • The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper). 
code of ethics: 
 Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https://neurips.cc/public/EthicsGuidelines? 
broader impacts: 
 Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?Answer: [NA] Justification: This paper presents work related to the theoretical aspect of Machine Learning. The topic of the paper, bilevel optimization, has many applications in machine learning. So there are some potential societal consequences of our work when our algorithms are applied to various machine learning tasks, but none which we feel is directly related to our algorithm.Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: The paper poses no such risks.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?Answer: [Yes]Justification: All the datasets used in our experiments are publicly available Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. • If this information is not available online, the authors are encouraged to reach out to the asset's creators. 
new assets: 
 Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?Answer: [NA]Justification: The paper does not release new assets.Guidelines:• The answer NA means that the paper does not release new assets.• Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc. • The paper should discuss whether and how consent was obtained from people whose asset is used. • At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file. 
crowdsourcing and research with human subjects: 
 Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?Answer: [NA]Justification: The paper does not involve crowdsourcing nor research with human subjects Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should be included in the main paper. • According to the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector. 
institutional review board (irb) approvals or equivalent for research with human subjects: 
 Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?Answer: [NA]Justification: The paper does not involve crowdsourcing nor research with human subjects.Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.• Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 
 
 Appendix: In this section, we introduce more details of the experiments.The data generation process is as follows: we first randomly sample a ground truth x * , then we randomly generate m = 1000 input variables c, and calculate the ground truth label b = cx, then for each c i , we generate a set of n = 100 noisy observations by adding Guassian noise of scale 0.1. We also add a L 2 regularization term of strength 0.1 for the outer problem. During trainging, we use both the inner and outer learing rates of 0.001.The formulation of the problem is as follows:min x∈R p h(x) := f (x, y x ) = 1 N (val) N (val) n=1 Θ(y x ; ξ val n ) s.t. y x = arg min y∈R d g(x, y) = N (tr) n=1x n Θ(y; ξ tr n )In the above formulation, we have a pair of (noisy) training set {ξ tr n } N (tr) n=1 and validation set {ξ val n } N (val) n=1 , and x n , n ∈ [N (tr) ] are weights for training samples, y is the parameter of a model, and we denote the model by Θ. Note that y x is the model learned over the weighted training set. We fit a model with 3 fully connected layers for the MNIST dataset. We also use L 2 regularization with coefficient 10 -3 to satisfy the strong convexity condition. In the Experiments, we choose inner learning rate (γ, ρ) as 0.1 and outer learning rate η 1000. . Θ defines the model, x is the parameter of the backbone model and y is the parameter of the linear classifier. In summary, the lower level problem is to learn the optimal linear classifier y given the backbone x, and the upper level problem is to learn the optimal backbone parameter x.The Omniglot dataset includes 1623 characters from 50 different alphabets and each character consists of 20 samples. We follow the experimental protocols of [48] to divide the alphabets to train/validation/test with 33/5/12, respectively. We perform N -way-K-shot classification, more specifically, for each task, we randomly sample N characters and for each character, we sample K samples for training and 15 samples for validation. We augment the characters by performing rotation operations (multipliers of 90 degrees). We use a 4-layer convolutional neural network where each convolutional layer has 64 filters of 3×3 [15] . For the MiniImageNet, it has 64 training classes and 16 validation classes. Similar to Omniglot, we also perform the N -way-K-shot classification. We use a 4-layer convolutional neural network where each convolutional layer has 64 filters of 3×3 [15] for experiments. For the experiments, we use inner learning rates 0.4 and outer learning rates 0.1 for Omniglot related experiments and inner learning rates 0.01 and outer learning rates 0.05 for MiniImageNet-related experiments. We perform 4 inner gradient descent steps and set K m ax = 6 for the RT-MLMC method. The experimental results for the MiniImageNet dataset is shown in Figure 4 .In this section, we provide proof for the convergence results in the main text. First, we restate all assumptions needed in our proof below:Assumption B.1 (Assumption 4.1). The function f (x, y) is possibly non-convex and uniformly L-smooth, i.e. for anyx 1 , x 2 ∈ X , y 1 , y 2 ∈ R d , any ξ ∈ D u . Denote z 1 = (x 1 , y 1 ), z 2 = (x 2 , y 2 ),then we have:f (z 1 ; ξ) ≤ f (z 2 ; ξ) + ⟨∇f (z 2 , ξ), z 1 -z 2 ⟩ + L 2 ||z 1 -z 2 || 2 .or equivalently: ||∇f (z 1 ; ξ) -∇f (z 2 ; ξ)|| ≤ L||z 1 -z 2 ||. We also assume for any x ∈ X and any y ∈ R d , and we denote z = (x, y), then we have||∇ y f (z)|| ≤ C f .Assumption B.2 (Assumption 4.2). Function g(x, y) is uniformly µ-strongly convex w.r.t y for any given x and uniformly L-smooth, i.e. for any y 1 , y 2 ∈ R d and ζ ∈ D l , we have:g(x, y 1 ; ζ) ≥ g(x, y 2 ; ζ) + ⟨∇ y g(x, y 2 ; ζ), y 2 -y 1 ⟩ + µ 2 ||y 2 -y 1 || 2 .and we have:g(z 1 ; ζ) ≤ g(z 2 ; ζ) + ⟨∇g(z 2 ; ζ), z 1 -z 2 ⟩ + L 2 ||z 1 -z 2 || 2 . equivalently: ||∇g(z 1 ; ζ) -∇g(z 2 ; ζ)|| ≤ L||z 1 -z 2 ||. Furthermore, we have ∇ xy g(x, y) and ∇ y 2 g(x, y) are Lipschitz continuous with constant L xy and L y 2 respectively, i.e. we have: [T ]}, we assume there exists some constants α, C such that for any step t and any k > 0, we have:||∇ xy g(z 1 ; ζ) -∇ xy g(z 2 ; ζ)|| ≤ L xy ||z 1 -z 2 || and ||∇ y 2 g(z 1 ; ζ) -∇ y 2 g(z 2 ; ζ)|| ≤ L y 2 ||z 1 -z 2 ||.1 k t+k-1 τ =t ∇f (x t , y t ; ξ π τ ) - ∇f (x t , y t ) 2 ≤ k -α C 2 , 1 k t+k-1 τ =t ∇g(x t , y t ; ζ π τ ) -∇g(x t , y t ) 2 ≤ k -α C 2 , 1 k t+k-1 τ =t ∇ 2 g(x t , y t ; ζ π τ ) -∇ 2 g(x t , y t ) 2 ≤ k -α C 2 , 1 k t+k-1 τ =t ∇h(x t ; ξ π τ ) -∇h(x t ) 2 ≤k -α C 2 (Only for the conditional bilevel optimization case). for i = 0 to I -1 do 6:y r i+1 = (1 -γ r i )y r i + γ r i r(x r i , ζ π i ); u r i+1 = (1 -ρ r i )u r i + ρ r i ∇ y f (y r i ; ξ π i ); 7:x r i+1 = x r i -η r i ∇r(x r i ; ζ π i )u r i ; for i = 0 to m -1 do 6:y r i+1 = y r i + γ r i ∇ y f (x r i , y r i ; ξ π i ); 7: x r i+1 = x r i -η r i ∇ x f (x r i , y r i ; ξ π i ); 8:end for 9: end for B.1 Proof for Bilevel Optimization Problems Suppose we define:Φ(x, y) =∇ x f (x, y) -∇ xy g(x, y) × [∇ y 2 g(x, y)] -1 ∇ y f (x, y)Then we have ∇h(x) = Φ(x, y x ). Furthermore, we have the following proposition: Proposition B.7. Suppose Assumptions 4.1 and 4.2 hold, the following statements hold: a) y x is Lipschitz continuous in x with constant ρ = κ, where κ = L µ is the condition number of g(x, y).b) ∥Φ(x 1 ; y 1 ) -Φ(x 2 ; y 2 )∥ 2 ≤ L2 (∥x 1 -x 2 ∥ 2 + ∥y 1 -y 2 ∥ 2 ), where L = O(κ 2 ). c) h(x) is smooth in x with constant L i.e., for any given x 1 , x 2 ∈ X, we have ∥∇h(x 2 ) - ∇h(x 1 )∥ ≤ L∥x 2 -x 1 ∥ where L = O(κ 3 ).This is a standard results in bilevel optimization and we omit the proof here.The main technique we use to analysis the convergence of Algorithm 1 is aggregated analysis, in other words, we perform analysis over a block of length k instead of the more common per step analysis in stochastic bilevel optimization. Note that this technique is essential in analyzing without-replacement sampling for single-level optimization problems, and we extend them to the bilevel setting. ρ τ q τ , where ν t = ∇ x f (x t , y t ; ξ t ) -∇ xy g(x t , y t ; ; ζ t )u t , w t = ∇ y g(x t , y t ; ζ t ) and q t = ∇ y 2 g(x t , y t ; ζ t )u t -∇ y f (x t , y t ; ξ t ). In particular, we denote r x (u) as: Randomly sample a permutation of outer dataset to have {ξ π i , i ∈ [m]} and set x r 0 = x r-1 m or x 0 if r = 1.r x (u) = 1 2 u T ∇ y 2 g(x, y x )u -∇ y f (x, y x ) T ufor i = 0 to m -1 doInitial states y 0 and u 0 , learning rates {γ s j , ρ s j }, j ∈ [n i ], s ∈ [S].for s = 1 to S do 7:Sample a permutation of inner dataset to have {ζ π ξi,j , i ∈ [n i ]}, set y s 0 = y s-1 ni and u s 0 = P ϵ (u s-1 ni ) or y s 0 = y 0 and u s 0 = u 0 if s = 1.for j = 0 to n i -1 do 9:y s j+1 = (1 -γ s j )y s j + γ s j r(x r i ; ζ π ξi,j ); 10:u s j+1 = (1 -ρ s j )u s j + ρ s j ∇f (y s j ; ξ π i ); 11:end for 12:end for 13:x r i+1 = x r i -η r i ∇r ξ π i (x r i )u S ni ;14:end for 15: end for {ζ t,1 }, {ζ t,2 } and {ζ t,3 } for ∇ y f (x, y), ∇ x f (x, y),∇ y g(x, y), ∇ y 2 g(x, y), ∇ xy g(x, y) respectively. The two orders used by Algorithm 1 is a special case of the proof and is practically appealing as they can reuse computation as we discussed in the introduction.Lemma B.8. For τ > t ≥ 0, the bias of ∥ τ -1 τ =t η τ ν τ -∇h(x t ) ∥ 2 , τ -1 τ =t ρ τ ∇r xt (u t )-q τ 2 and ∥ τ -1 τ =t γ τ ∇ y g(x τ , y τ , ζ τ ,1) -∇ y g(x t , y t ) ∥ 2 are bounded by Eq. (15), Eq. ( 16) and Eq. (17) .EQUATIONFor the first and the third terms above, we have:EQUATIONandEQUATIONCombine Eq ( 13) and ( 14) with Eq. ( 12) and denoteL2 1 = L 2 + 2L 2 xy C 2 f µ 2to have:EQUATIONNext, for the term ∥ τ -1 τ =t ρ τ (∇r xt (u t ) -q τ )∥ 2 , we have:τ -1 τ =t ρ τ ∇r xt (u t ) -q τ 2 ≤ τ -1 τ =t ρ τ ∇ y 2 g(x τ , y τ , ζ τ ,2 )u τ -∇ y f (x τ , y τ , ξ τ ,1 ) -∇ y 2 g(x t , y xt )u t -∇ y f (x t , y xt ) 2 ≤ 2∥ τ -1 τ =t ρ τ ∇ y f (x t , y xt ) -∇ y f (x τ , y τ ; ξ τ ,1 ) ∥ 2 + 2∥ τ -1 τ =t ρ τ ∇ y 2 g(x t , y xt )u t -∇ y 2 g(x τ , y τ ; ζ τ ,2 )u τ ∥ 2 ≤ 2∥ τ -1 τ =t ρ τ ∇ y f (x t , y xt ) -∇ y f (x τ , y τ ; ξ τ ,1 ) ∥ 2 + 4∥ τ -1 τ =t ρ τ ∇ y 2 g(x τ , y τ ; ζ τ ,2 ) u t -u τ ∥ 2 + 4C 2 f µ 2 ∥ τ -1 τ =t ρ τ ∇ y 2 g(x t , y xt ) -∇ y 2 g(x τ , y τ ; ζ τ,2 ) ∥ 2 ≤ 2∥ τ -1 τ =t ρ τ ∇ y f (x t , y xt ) -∇ y f (x τ , y τ ; ξ τ ,1 ) ∥ 2 + 4L 2 ρτ t τ -1 τ =t ρ τ ∥(u t -u τ )∥ 2 + 4C 2 f µ 2 ∥ τ -1 τ =t ρ τ ∇ y 2 g(x t , y xt ) -∇ y 2 g(x τ , y τ ; ζ τ,2 ) ∥ 2The first and the third terms can be bounded similarly as in Eq. ( 13) and Eq. ( 14), and denoteL2 2 = L 2 + 2L 2 y 2 C 2 f µ 2, then we have:EQUATIONFinally, for ∥ τ -1 τ =t γ τ (∇ y g(x τ , y τ , ζ τ ,1 ) -∇ y g(x t , y t ))∥ 2 , we have:EQUATIONThis completes the proof.Lemma B.9. For all t ≥ 0 and any constant m > 0, suppose ηt < 1 2 L , the iterates generated satisfy:h(x t+k ) ≤ h(x t ) - ηt 2 ∥∇h(x t )∥ 2 - ηt 4 ∥ν t ∥ 2 + 3 L2 1 ηt ∥y xt -y t ∥ 2 + 3 L2 1 t+k-1 τ =t η τ ∥z t -z τ ∥ 2 + 4L 2 ηt ∥u xt -u t ∥ 2 + 4L 2 t+k-1 τ =t η τ ∥u t -u τ ∥ 2 + 3 ηt ∥ t+k-1 τ =t η τ ∇ x f (x t , y t ) -∇ x f (x t , y t ; ξ τ,2 ) ∥ 2 + 6C 2 f ηt µ 2 ∥ t+k-1 τ =t η τ ∇ xy g(x t , y t ) -∇ xy g(x t , y t ; ζ τ,3 ) ∥ 2Proof. By the smoothness of f , we have:EQUATIONwhere equality (a) uses ⟨a, b⟩ = 1 2 [∥a∥ 2 + ∥b∥ 2 -∥a -b∥ 2 ]; (b) follows the assumption that ηt < 1/2 L. Next, for the third term, by using Eq. ( 15), we have:ηt 2 ∥∇h(x t ) -νt ∥ 2 ≤ 3 L2 1 ηt ∥y xt -y t ∥ 2 + 3 L2 1 t+k-1 τ =t η τ ∥z t -z τ ∥ 2 + 4L 2 ηt ∥u xt -u t ∥ 2 + 4L 2 t+k-1 τ =t η τ ∥u t -u τ ∥ 2 + 3 ηt ∥ t+k-1 τ =t η τ ∇ x f (x t , y t ) -∇ x f (x t , y t ; ξ τ,2 ) ∥ 2 + 6C 2 f ηt µ 2 ∥ t+k-1 τ =t η τ ∇ xy g(x t , y t ) -∇ xy g(x t , y t ; ζ τ,3 ) ∥ 2whereL2 1 = L 2 + 2L 2 xy C 2 f µ 2. This completes the proof.Lemma B.10. When γt < 1 4L , the error of the inner variable y t are bounded through Eq. (19) .Proof. then by proposition B.19 and choose γt < 1 4L , we have:∥y t+k -y xt ∥ 2 ≤ (1 - µγ t 2 )∥y t -y xt ∥ 2 - γ2 t 2 ∥∇ y g(x t , y t )∥ 2 + 4γ t µ ∥∇ y g(x t , y t ) -wt ∥ 2Furthermore, by the generalized triangle inequality, we have:EQUATIONwhere the second inequality is due to γt < 1 4L < 1 µ and uses Eq. ( 17). This completes the proof.Lemma B.11. When ρt < 1 4L , the error of variable u t are bounded through Eq. (20) .Proof. suppose we choose ρt < 1 4L , then we have:∥u t+k -u xt ∥ 2 ≤ (1 - µρ t 2 )∥u t -u xt ∥ 2 - (ρ t ) 2 2 ∥∇r xt (u t )∥ 2 + 4ρ t µ ∥∇r xt (u t ) -qt ∥ 2Furthermore, by the generalized triangle inequality, we have:∥u t+k -u x t+k ∥ 2 ≤ (1 + µρ t 4 )∥u t+k -u xt ∥ 2 + (1 + 4 µρ t )∥u x t+k -u xt ∥ 2 ≤ (1 + µρ t 4 )∥u t+k -u xt ∥ 2 + 5 L2 µρ t ∥x t -x t+k ∥ 2 ≤ (1 - µρ t 4 )∥u t -u xt ∥ 2 - ρ2 t 2 ∥∇r xt (u t )∥ 2 + 5 L2 µρ t ∥x t -x t+k ∥ 2 + 30 L2 2 ρt µ ∥y xt -y t ∥ 2 + 30 L2 2 µ t+k-1 τ =t ρ τ ∥z t -z τ ∥ 2 + 20L 2 µ t+k-1 τ =t ρ τ ∥u t -u τ ∥ 2 + 30 µρ t ∥ t+k-1 τ =t ρ τ ∇ y f (x t , y t ) -∇ y f (x t , y t ; ξ τ,1 ) ∥ 2 + 60C 2 f µ 3 ρt ∥ t+k-1 τ =t ρ τ ∇ y 2 g(x t , y t ) -∇ y 2 g(x t , y t ; ζ τ,2 ) ∥ 2 (20)where in the last inequality, we use the condition that ρt < 1 4L and Eq. ( 16). This completes the proof.Lemma B.12. Suppose ηt ≤ min1 96 L2 1 , 1 32L 2 c 2 1 , 1 128L 2 , 1 64L 2 c 2 2 , 1 96 L2 2 c 2 2, then the drift of z τ = (x τ , y τ ) and u τ can be bounded as in Eq. ( 25) and Eq. (24) .Proof. We first have:EQUATIONFor the term ∥z t -z τ ∥ 2 , we have ∥z t -z τ ∥ 2 = ∥x t -x τ ∥ 2 + ∥y t -y τ ∥ 2 . Combine Eq. ( 21) with Eq. ( 15) and ( 17), sum τ over [t, t + k -1] and use Proposition B.18 to have:t+k-1 τ =t η τ ∥z t -z τ ∥ 2 ≤ t+k-1 τ =t η τ τ -1 τ =t 12 L2 1 ητ t η τ + 4L 2 γτ t γ τ ∥z t -z τ ∥ 2 + 2 ᾱt η2 t ∥∇h(x t )∥ 2 + 2 ᾱt γ2 t ∥∇ y g(x t , y t )∥ 2 + 12 L2 1 ᾱt (η t ) 2 ∥y xt -y t ∥ 2 + 16L 2 ᾱt (η t ) 2 ∥u xt -u t ∥ 2 + 16L 2 ηt t+k-1 τ =t η τ τ -1 τ =t η τ ∥u t -u τ ∥ 2 + 12η t η 2 t k 2-α C 2 + 24C 2 f µ 2 ηt η 2 t k 2-α C 2 + 4η t γ 2 t k 2-α C 2By the fact that τ ≤ t + k, we have:EQUATIONFor ∥u t -u τ ∥ 2 , we combine Eq. ( 21) with Eq. ( 16) and sum τ over [t, t + k -1] and use the fact that τ ≤ t + k to have:EQUATIONSuppose we haveγ τ = c 1 η τ = 20κ L3 µ η τ , ρ τ = c 2 η τ = 40κLη τ and set:(η t ) 2 < min 1 96 L2 1 , 1 32L 2 c 2 1 , 1 128L 2 , 1 64L 2 c 2 2 , 1 96 L2 2 c 2 2Then we combine Eq. ( 22) and Eq. ( 23) to have:EQUATIONandEQUATIONThis completes the proof.Suppose we define the potential function G(t):G t = h(x t ) + ϕ y ∥y t -y xt ∥ 2 + ϕ u ∥u t -u xt ∥ 2where ϕ u = 32L 2 ηt µ ρt , ϕ y = , we have:1 E E-1 e=0 ∥∇h(x ke )∥ 2 ≤ 2∆ kEη + 2 L2 k -α C 2where c 1 , c 2 , L0 , L are constants related to the smoothness parameters of h(x), ∆ is the initial sub-optimality, R and I are defined in Algorithm 1.Proof. First we combine Lemma B.9-Lemma B.11, and use the condition that γ t = 20κ L3 µ η t , ρ t = 40κ Lη t , and by the assumption about gradient error:EQUATION)Suppose we denoteL2 4 = 960κ 2 L2 2 + 3 L2 1 + 80κ 2 L2 3 + 640κ 2 L 2 + 4L 2, and set:EQUATIONand we bound the terms ∥u xt -u t ∥ 2 and ∥z t -z τ ∥ 2 through Lemma B.12 to have:G t+k -G t ≤ - ηt 2 ∥∇h(x t )∥ 2 + 24 + 48C 2 f µ 2 + 8c 2 1 + 24c 2 2 + 48c 2 2 C 2 f µ 2 L2 4 ηt η 2 t k 2-α C 2 + 3 + 6C 2 f µ 2 + 80 L2 3 µ 2 + 960κ 2 + 1920κ 2 C 2 f µ 2 η 2 t k 2-α C 2 ηtNext suppose we set L4 ηt < 1, then we have:G t+k -G t ≤ - ηt 2 ∥∇h(x t )∥ 2 + 24 + 48C 2 f µ 2 + 8c 2 1 + 24c 2 2 + 48c 2 2 C 2 f µ 2 + 3 + 6C 2 f µ 2 + 80 L2 3 µ 2 + 960κ 2 + 1920κ 2 C 2 f µ 2 η 2 t k 2-α C 2 ηtFor ease of notation we denote L2 = 24 +48C 2 f µ 2 + 8c 2 1 + 24c 2 2 + 48c 2 2 C 2 f µ 2 + 3 + 6C 2 f µ 2 + 80 L2 3 µ 2 + 960κ 2 + 1920κ 2 C 2 f µ 2, then we have:G t+k -G t ≤ - ηt 2 ∥∇h(x t )∥ 2 + L2 η 2 t k 2-α C 2 ηtSum the above inequality over E phases to have:E-1 e=0 ηke 2 ∥∇h(x ke )∥ 2 ≤ G 0 -G kE + E-1 e=0 L2 η 2 ke k 2-α C 2 ηke ≤ ∆ h + ϕ y ∆ y + ϕ u ∆ u + E-1 e=0 L2 η 2 ke k 2-α C 2 ηkewe define ∆ h = h(x 0 ) -h * as the initial sub-optimality of the function, ∆ y = ∥y 0 -y x0 ∥ 2 as the initial sub-optimality of the inner variable estimation, ∆ u = ∥u 0 -u x0 ∥ 2 as the initial sub-optimality of the hyper-gradient estimation. Then suppose we choose η t = η be some constant, then we have:kη 2 E-1 e=0 ∥∇h(x ke )∥ 2 ≤ ∆ h + ϕ y ∆ y + ϕ u ∆ u + L2 Ek 1-α C 2 ηDivide by kEη/2 on both sides to have:1 E E-1 e=0 ∥∇h(x ke )∥ 2 ≤ 2 kEη ∆ h + ϕ y ∆ y + ϕ u ∆ u + 2 L2 k -α C 2Suppose we choose the largest η = 1 k L0such that the conditions of Lemma B.9-Lemma B.12 are satisfied, then we have:1 E E-1 e=0 ∥∇h(x ke )∥ 2 ≤ 2 L0 E ∆ + ϕ y ∆ y + ϕ u ∆ u + 2 L2 k -α C 2Then to reach an ϵ-stationary point, we need to have:E ≥ 4 L0 ∆ ϵ 2 , and , k ≥ 4 L2 C 2 ϵ 2 1/α in other words T = Ek ≥ 4 L0(4 L2 C 2 ) 1/α ∆ ϵ 2+2/α. This completes the proof.In this section we study the convergence rate of Algorithm 2. Lemma B.14. For all t ≥ 0 and any constant k > 0, suppose ηt < 1 4 L , the iterates generated satisfy:h(x t+k ) ≤ h(x t ) - kη 8 ∥∇h(x t )∥ 2 - kη 4 ∥ν t ∥ 2 + 3 1 + 4κ 2 (192C 2 f + 96 * 24 L2 2 ) µ 2 + 48 L2 1 µ 2 ηk 1-α C 2 + 12ηkL 2 (1 - kµρ 2 ) E l ∆ u + 3ηk(2 L2 1 + 4 * 96κ 2 L2 2 )(1 - kµγ 2 ) E l ∆ yProof. By the smoothness of f , follow similar derivation as in Eq. ( 18), we have:h(x t+k ) ≤ h(x t ) - ηt 2 ∥∇h(x t )∥ 2 + ηt 2 ∥∇h(x t ) -νt ∥ 2 - ηt 4 ∥ν t ∥ 2In particular, for the term ∥ τ -1 τ =t η τ (ν τ -∇h(x t ))∥ 2 , we have:EQUATIONFor the third term, we follow similar derivation as in Eq. ( 12)-Eq. ( 15), we have:EQUATIONwhereL2 1 = L 2 + 2L 2 xy C 2 f µ 2. Combine Eq. ( 27) and Eq. ( 28), we have:EQUATIONFor the term t+k-1 τ =t η τ x t -x τ 2 , by Eq. ( 27) and Eq. ( 28), we have:x t -x τ 2 ≤ 2∥ τ -1 τ =t η τ ν τ -∇h(x t ) ∥ 2 + 2(η τ t ) 2 ∥∇h(x t )∥ 2 ≤ 6 τ -1 τ =t η τ ∇h(x t ) -∇h(x t ; ξ τ ) 2 + 6 L2 ητ t τ -1 τ =t η τ x t -x τ 2 + 6η τ t τ -1 τ =t η τ 2 L2 1 ∥y ξτ xτ -y τ ,T l ∥ 2 + 4L 2 ∥u ξτ xτ -u τ ,T l ∥ 2 + 2(η τ t ) 2 ∥∇h(x t )∥ 2 ≤ 6η 2 t k 2-α C 2 + 6 L2 ηt t+k-1 τ =t η τ x t -x τ 2 + 6η t t+k-1 τ =t η τ 2 L2 1 ∥y ξτ xτ -y τ ,T l ∥ 2 + 4L 2 ∥u ξτ xτ -u τ ,T l ∥ 2 + 2(η t ) 2 ∥∇h(x t )∥ 2Multiply η τ on both sides and sum over [t, t + k -1], we have:(1 -6 L2 (η t ) 2 ) t+k-1 τ =t η τ x t -x τ 2 ≤ 6(η t ) 2 t+k-1 τ =t η τ 2 L2 1 ∥y ξτ xτ -y τ ,T l ∥ 2 + 4L 2 ∥u ξτ xτ -u τ ,T l ∥ 2 + 2(η t ) 3 ∥∇h(x t )∥ 2 + 6η t η 2 t k 2-α C 2By the condition that ηt < 1 4 L , we have:EQUATIONCombine Eq. ( 30) with Eq. ( 29), then we have:EQUATIONwhere the second inequality uses the conditions that ηt < 1 4 L . By Lemma B.16,we have:EQUATIONwhere ∆ u and ∆ y denotes the upper bounds of the initial estimation errors of y x and u x . This completes the proof.Lemma B.15. Suppose we have γl < 1 2L and ρl < 1 2L , thenl+k-1 l ′ =l γ l ′ ∥y l ′ -y l ∥ 2 and l+k-1 l ′ =l ρ l ′ ∥u l -u l ′ ∥ 2 canbe bounded as in Eq. (33) and Eq. (Proof. We first bound l+k-1 l ′ =l γ l ′ ∥y l ′ -y l ∥ 2 , in fact, we have:∥y l -y l ′ ∥ 2 = ∥ l ′ -1 l′ =l γl′ wl′ ∥ 2 ≤ 2∥ l ′ -1 l′ =l γl′ ∇ y g ζ l′ (x, yl′ ) -∇ y g(x, y l ) ∥ 2 + 2(γ l ′ l ) 2 ∥∇ y g(x, y l )∥ 2 ≤ 4L 2 γl ′ l l ′ -1 l′ =l γl′ ∥yl′ -y l ∥ 2 + 2(γ l ′ l ) 2 ∥∇ y g(x, y l )∥ 2 + 4∥ l ′ -1 l′ =l γl′ ∇ y g ζ l′ (x, y l ) -∇ y g(x, y l ) ∥ 2 ≤ 4L 2 γl l+k-1 l′ =l γl′ ∥yl′ -y l ∥ 2 + (γ l ) 2 ∥∇ y g(x, y l )∥ 2 + 4γ 2 l k 2-α C 2Multiple γ l ′ on both sides, and then sum l′ in [l, l + k -1],we have:(1 -4L 2 γ2 l ) l+k-1 l ′ =l γ l ′ ∥y l ′ -y l ∥ 2 ≤ 2(γ l ) 3 ∥∇ y g(x, y l )∥ 2 + 4γ l γ 2 l k 2-α C 2Since we have γl < 1 2L , we have:EQUATIONNext for l+k-1 l ′ =l ρ l ′ ∥u l -u l ′ ∥ 2 , we have:∥u l -u l ′ ∥ 2 = ∥ l ′ -1 l′ =t ρl′ ql′ ∥ 2 ≤ 2∥ l ′ -1 l′ =l ρl′ ql′ -∇r x (u l ) ∥ 2 + 2(ρ l ′ l ) 2 ∥∇r x (u l )∥ 2 ≤ 12 L2 2 (ρ l ′ l ) 2 ∥y x -y T l ∥ 2 + 8L 2 ρl ′ l l ′ -1 l′ =l ρl′ ∥u l -ul′ ∥ 2 + 2(ρ l ′ l ) 2 ∥∇r x (u l )∥ 2 + 24C 2 f µ 2 ∥ l ′ -1 l′ =l ρl′ ∇ y 2 g(x, y T l ) -∇ y 2 g ζ l′ (x, y T l ) ∥ 2 ≤ 12 L2 2 (ρ l ) 2 ∥y x -y T l ∥ 2 + 8L 2 ρl l+k-1 l′ =l ρl′ ∥u l -ul′ ∥ 2 + 2(ρ l ) 2 ∥∇r x (u l )∥ 2 + 24C 2 f µ 2 ρ 2 l k 2-α C 2Multiple ρ l ′ on both sides, and then sum l ′ in [l, l + k -1], we have:(1 -8L 2 ρ2 l ) l+k-1 l ′ =l ρ l ′ ∥u l -u l ′ ∥ 2 ≤ 12 L2 2 (ρ l ) 3 ∥y x -y T l ∥ 2 + 2(ρ l ) 3 ∥∇r x (u l )∥ 2 + 24C 2 f µ 2 ρl ρ 2 l k 2-α C 2Since we have ρl < 1 2L , we have:l+k-1 l ′ =l ρ l ′ ∥u l -u l ′ ∥ 2 ≤ 24 L2 2 (ρ l ) 3 ∥y x -y T l ∥ 2 + 4(ρ l ) 3 ∥∇r x (u l )∥ 2 + 48C 2 f µ 2 ρl ρ 2 l k 2-α C 2 (34)Lemma B.16. When γt < 1 128Lκ and ρt < 1 256Lκ , the error of the inner variable y l and the variable u l are bounded through Eq. (35) and Eq. (36) .Proof. Since the outer iteration t is fixed for each inner loop, we omit t and ξ t in the notation for clarity. Furthermore, we assume performing the updates to y and u separately in Algorithm 2. Although the alternative updates in Algorithm 2 is appealing in practice as we can reuse ∇ y g for ∇ y 2 g, update y first and u next can avoid treating complicated higher order terms. Follow similar derivation of Lemma B.10, if γt < 1 4L , we have:∥y l+k -y x ∥ 2 ≤ (1 - µγ l 2 )∥y l -y x ∥ 2 - γ2 l 2 ∥∇ y g(x, y l )∥ 2 + 8Lκ l+k-1 l ′ =l γ l ′ ∥y l ′ -y l ∥ 2 + 8 µγ l ∥ l+k-1 l ′ =l γ l ′ ∇ y g(x, y l ) -∇ y g(x, y l ; ζ l ′ ) ∥ 2Combine with Eq. ( 33) and use the condition that γl < 1 128Lκ < 1 4L , we have:∥y l+k -y x ∥ 2 ≤ (1 - µγ l 2 )∥y l -y x ∥ 2 - γ2 l 4 ∥∇ y g(x, y l )∥ 2 + 64Lκγ l γ 2 l k 2-α C 2 + 8 µγ l γ 2 l k 2-α C 2Suppose we perform E l rounds, then by telescoping, we have:EQUATIONNext, follow similar derivations as in Lemma B.11 and choose ρt < 1 4L , then we have:∥u l+k -u x ∥ 2 ≤ (1 - µρ l 2 )∥u l -u x ∥ 2 - (ρ l ) 2 2 ∥∇r x (u l )∥ 2 + 24 L2 2 (ρ t ) µ ∥y x -y T l ∥ 2 + 16Lκ l+k-1 l ′ =l ρ l ′ ∥u l -u l ′ ∥ 2 + 48C 2 f µ 3 ρl ∥ l+k-1 l ′ =l ρ l ′ ∇ y 2 g(x, y T l ) -∇ y 2 g(x, y T l ; ζ l ′ ) ∥ 2Combine with Eq. ( 34) and use the condition that ρl < 1 256Lκ < 1 4L , we have:∥u l+k -u x ∥ 2 ≤ (1 - µρ l 2 )∥u l -u x ∥ 2 - (ρ l ) 2 4 ∥∇r x (u l )∥ 2 + 48 L2 2 (ρ l ) µ ∥y x -y T l ∥ 2 + 96C 2 f µ 3 ρl ρ 2 l k 2-α C 2we perform E l rounds, then by telescoping, we have:EQUATIONThis completes the proof.Theorem B.17. Suppose Assumptions 4.1-4.4 and 4.7 are satisfied, and Assumption 4.5 is satisfied with some α and C for the example order we use in Algorithm 2. We choose learning ratesη t = η = 1 8k L , γ t = γ = 1 256kLκ , ρ t = ρ = 1 512kLκand denote T = R × m be the total number of outer steps and T l = S × max(n i ) be the maximum inner steps. Then for any pair of values (E, k) which has T = E × k and T l = E l × k, we have:1 E E-1 e=0 ∥∇h(x ke )∥ 2 ≤ 8∆ h kEη + C u (1 - kµρ 2 ) E l ∆ u + C y (1 - kµγ 2 ) E l ∆ y + (C ′ ) 2 C 2 k -αwhere C u , C y , C ′ , L are constants related to the smoothness parameters of h(x), ∆ h , ∆ u , ∆ y are the initial sub-optimality, R and S are defined in Algorithm 2.Proof. First, by Lemma B.14, we have:h(x t+k ) ≤ h(x t ) - kη 8 ∥∇h(x t )∥ 2 - kη 4 ∥ν t ∥ 2 + 3 1 + 4κ 2 (192C 2 f + 96 * 24 L2 2 ) µ 2 + 48 L2 1 µ 2 ηk 1-α C 2 + 12ηkL 2 (1 - kµρ 2 ) E l ∆ u + 3ηk(2 L2 1 + 4 * 96κ 2 L2 2 )(1 - kµγ 2 ) E l ∆ ySum the above inequality over E phases to have:kη 8 E-1 e=0 ∥∇h(x ke )∥ 2 ≤ ∆ h + 12ηkEL 2 (1 - kµρ 2 ) E l ∆ u + 3ηkE(2 L2 1 + 4 * 96κ 2 L2 2 )(1 - kµγ 2 ) E l ∆ y + 3E 1 + 4κ 2 (192C 2 f + 96 * 24 L2 2 ) µ 2 + 48 L2 1 µ 2 ηk 1-α C 2Divide by kEη/8 on both sides to have:1 E E-1 e=0 ∥∇h(x ke )∥ 2 ≤ 8∆ h kEη + 96L 2 (1 - kµρ 2 ) E l ∆ u + 24(2 L2 1 + 4 * 96κ 2 L2 2 )(1 - kµγ 2 ) E l ∆ y + 24 1 + 4κ 2 (192C 2 f + 96 * 24 L2 2 ) µ 2 + 48 L2 1 µ 2 k -α C 2By the condition η < 1 4k L , γ < 1 E E-1 e=0 ∥∇h(x ke )∥ 2 ≤ 64 L∆ h E + 96L 2 (1 - kµρ 2 ) E l ∆ u + 24(2 L2 1 + 4 * 96κ 2 L2 2 )(1 - kµγ 2 ) E l ∆ y + 24 1 + 4κ 2 (192C 2 f + 96 * 24 L2 2 ) µ 2 + 48 L2 1 µ 2 k -α C 2Then to reach an ϵ-stationary point, we need to have: τ =t η τ ∇ y g(x t , y xt ) -∇ y g(x τ , y τ ; ζ τ ,1 ) ∥ 2 ≤ η 2 k 2-α C 2 for any τ ≤ t + k -1. The same upper bound is achieved for other properties: ∇h(x), ∇ x f (x, y),∇ y g(x, y), ∇ y 2 g(x, y), ∇ xy g(x, y) Proof. This proposition can be adapted from Lemma 2 of [34] .E ≥ 128 L∆ ϵ 2 , and , k α ≥ 48C 2 ϵ 2 1 + 4κ 2 (192C 2 f + 96 * 24 L2 2 ) µ 2 + 48 L2Proposition B.19. Suppose we have function g(y), which is L-smooth and µ-strongly-convex, then suppose γ < 1 2L , the progress made by one step of gradient descent is:∥y t+1 -y * ∥ 2 ≤ 1 - µγ 2 ∥y t -y * ∥ 2 - γ 2 2 ∥∇ y g(y t )∥ 2 + 4γ µ ∥∇ y g(y t ) -w t ∥ 2where y * is the minimum of g(y) and we have update rule y t+1 = y t -γω t .Proof. First, by the strong convexity of of function g(y), we have:g(y * ) ≥ g(y t ) + ⟨∇ y g(y t ), y * -y t ⟩ + µ 2 ∥y * -y t ∥ 2= g(y t ) + ⟨∇ y g(y t ), y * -y t+1 ⟩ + ⟨∇ y g(y t ), y t+1 -EQUATIONThen by L-smoothness, we have:EQUATIONTable 2 : Comparisons of the Bilevel Opt. & Conditional Bilevel Opt. algorithms for finding an ϵ-stationary point. An ϵ-stationary point is defined as ∥∇h(x)∥ ≤ ϵ. Gc(f, ϵ) and Gc(g, ϵ) denote the number of gradient evaluations w.r.t. f (x, y) and g(x, y); JV (g, ϵ) denotes the number of Jacobian-vector products; HV (g, ϵ) is the number of Hessian-vector products. m and n are the number of data examples for the outer and inner problems, in particular, n is the maximum number of inner problems for conditional bilevel optimization. Our methods have a dependence over example numbers (max(m, n)) q , q is a value decided by without-replacement sampling strategy and can have value in [0, 1] (A herding-based permutation [33] can let q = 0).Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental material?Answer: [Yes] Justification: the datasets used in experiments are publicly available and we include the code implementation in the supplementary material.Guidelines:• The answer NA means that paper does not include experiments requiring code. • Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted.Justification: We follow all of NeurIPS Code of Ethics Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 


Decision: The paper introduces WiOR-BO, a new algorithm that leverages without-replacement sampling to enhance the efficiency of bilevel optimization. The method is designed to address the limitations of traditional SGD-type bilevel algorithms, which rely on independent sampling. WiOR-BO is applicable to standard bilevel optimization, conditional bilevel optimization, and specific cases such as minimax and compositional optimization problems. The authors provide theoretical convergence guarantees and demonstrate the algorithm's effectiveness through empirical validation on multiple tasks.

Pros:

Novel Application: Introduces without-replacement sampling to bilevel optimization, a relatively unexplored area that offers potential benefits in terms of computational efficiency.
Theoretical Contribution: Provides a comprehensive theoretical analysis, demonstrating improved convergence rates compared to existing methods.
Empirical Validation: Strong empirical results on various tasks, showing the effectiveness of the proposed algorithms.
Practical Relevance: The proposed method offers a practical alternative to existing bilevel optimization algorithms, particularly in large-scale machine learning scenarios.

Cons:

Limited Novelty in Algorithmic Design: The method is perceived as a direct combination of existing techniques (e.g., SOBA) with without-replacement sampling, limiting the novelty.
Unclear Theoretical Advantage: Despite the theoretical contributions, the actual advantage of the proposed method over deterministic and stochastic SOBA remains unclear, especially when considering both gradient and Hessian terms.
Complexity Considerations: The improvement in convergence rates is not as significant as initially claimed, leading to questions about the overall theoretical superiority of the method.
Presentation Issues: The clarity of some figures and explanations in the paper could be improved, particularly in terms of technical details and parameter settings.

Thins that should be addressed when preparing the final version of the paper.

-  There are lingering concerns about the theoretical advantage of the proposed method over existing methods like deterministic and stochastic SOBA, particularly in terms of the trade-off between gradient and Hessian complexities.

-  While the empirical results are strong, the theoretical justification for why the proposed method is preferable in all scenarios remains insufficiently addressed.

The paper offers a meaningful contribution to bilevel optimization, particularly through its introduction of without-replacement sampling. However, the novelty and theoretical justification of the method require further clarification. The empirical results are compelling, but the paper would benefit from additional theoretical insights and clearer presentation of its advantages over existing methods. These issues can be addressed when preparing a camera ready version.

Reviews: review: paper_summary: paper_summary: The paper introduces a new algorithm called WiOR-BO, which leverages without-replacement sampling to achieve faster convergence rates compared to traditional independent sampling methods. The algorithm is applied to standard bilevel optimization, conditional bilevel optimization, and specific cases such as minimax and compositional optimization problems. The authors provide a theoretical analysis demonstrating that WiOR-BO achieves a significantly improved convergence rate over existing methods. Empirical validation on synthetic and real-world tasks, including hyper-data cleaning and hyper-representation learning, shows the superior performance and computational efficiency of the proposed algorithms.
strength_weakness: strength_weakness: The authors have selected learning rates and other tuning parameters to optimize performance, but they do not report any sensitivity analysis regarding these choices. This lack of analysis leaves open questions about the robustness of the algorithm's performance under different parameter settings. Additionally, there is no provided methodology or algorithm for efficiently selecting these parameters without the risk of overfitting to the validation set, which could undermine the generalizability of the results.

Furthermore, the proposed algorithms are designed under the assumption of a static dataset, which restricts their applicability in scenarios where data arrives sequentially, such as in dynamic or online learning environments. This limitation could be mitigated by extending the algorithms to handle streaming data or datasets that evolve over time, making the methods more versatile and applicable to real-world situations where data distribution may change. Addressing these issues would enhance the robustness and broader applicability of the proposed solutions.
questions: questions: 1. Comparative Benchmarks: The paper provides empirical validation, but could the authors include more comprehensive comparisons with state-of-the-art methods, especially in dynamic environments? How does the proposed method perform in cases where data distributions shift over time?
2. Computational Efficiency: While the method reduces the number of backward passes, how does it scale with very large datasets or high-dimensional data?
3. Hyperparameter Tuning: Could the authors elaborate on the process of selecting hyperparameters, such as learning rates and other tuning parameters? Specifically, are there any strategies employed to avoid overfitting during hyperparameter tuning? Additionally, how sensitive is the performance of the algorithm to these hyperparameter choices, and have the authors conducted any sensitivity analyses to assess this? Addressing these points could provide more clarity on the robustness and practical implementation of the proposed algorithms.
limitations: limitations: Are there any specific limitations that the authors anticipate in applying this method to different domains or tasks? What future directions do the authors envision for improving the algorithm, particularly concerning hyperparameter tunning, adaptivity and scalability?
paper_summary: paper_summary: This paper investigates stochastic bilevel optimization. One common practice herein is that the data samples used to calculate all the derivatives are mutually independent, which may however introduce more computational cost compared to the case where the samples are reused. Motivated by this, the authors explore stochastic bilevel optimization without using independent sampling, particularly under the without-replacement sampling scheme. Based on previous single-loop algorithm, two with-replacement bilevel optimization algorithms are proposed for standard bilevel optimization and conditional bilevel optimization, respectively, which can also be applied to several special cases. The authors show that faster convergence rates can be achieved by the proposed algorithms compared to the counterparts with independent sampling. Experiments on multiple tasks are conducted to further justify the performance of the proposed algorithms.
strength_weakness: strength_weakness: 1. It is not clear what the key technical contributions in the theoretical analysis are.

2. As the authors claim the extension to the special cases of minimax and compositional optimization problems as a contribution, it is better to compare with the baselines in these problems.

3. The presentation can be further improved. For examples, the figures in the experiments are not clear because of the small font size.
questions: questions: 1. How was the bias due to the without-replacement sampling handled in the analysis? Using window averaging? In this case, did you need to store the gradients of previous samples?

2. When you compared the convergence rate of algorithm 1 with baselines, e.g., stocBiO, how did you justify that the faster convergence is because of the sampling but not the single-loop update? Given that stocBiO is not a single-loop algorithm.
limitations: limitations: Yes
paper_summary: paper_summary: This paper improves the computational inefficiencies in bilevel optimization algorithms that rely on independent sampling. It proposes a novel without-replacement sampling-based algorithm, which achieves a faster convergence rate than existing methods.
strength_weakness: strength_weakness: - The paper claimed that "Compared to independent sampling-based algorithm such as stocBiO, we require extra space to generate and store the orders of examples, but this cost is negligible compared to memory and computational cost of training large scale models." Is there any empirical result to show the memory costs?

- Projection is applied to $u_i^r$, but I found other bilevel papers (such as StocBio, BSA) don't use that. Is there any intuition or explanation about that?

- In Algorithm 2, the generation of $\\{\zeta_{\xi_t^{\pi}, t_l}^{\pi} \in \mathcal{D}_{l, \xi_t^{\pi}}. t_l\in [T_l]\\}$ is not very clear. It seems that there are $S\times n_i$ samples for inner loops given a outer sample $\xi_i$. How do you execute the permutation for the inner dataset.

- In experiments, what is the best hyper-parameters for each baselines, such inner and outer learning rates?

- What are the details of data sampling for other baselines? For other baselines, they usually need the independent training and validation set for upper and lower-level update, how did you construct data for them?

- Typos: Line 107, $w_x$ should be $u_x$.
questions: questions: Please refer to weaknesses.
limitations: limitations: Yes.
paper_summary: paper_summary: This paper introduces new algorithms that leverage without-replacement sampling to enhance bilevel optimization. The main contributions are:
1. Introducing more practical bilevel optimization algorithms using without-replacement sampling.
2. Providing comprehensive theoretical analysis demonstrating improved convergence rates and efficiency.
Novel Algorithms: Introduces simpler, more practical bilevel optimization algorithms using without-replacement sampling.
3. Validating the effectiveness through applications in hyper-data cleaning and hyper-representation learning.
strength_weakness: strength_weakness: 1. Why are the upper bounds the same in Assumptions 4.3, 4.4, and 4.7?
questions: questions: Line 71: Please use $O(\cdot)$ to denote the big O notation and $\tilde{O}(\cdot)$ to indicate that logarithmic terms are hidden.
LInes 131, 536, etc.: Some formulas exceed the line width.
limitations: limitations: NA
paper_summary: paper_summary: This paper proposes WiOR-(C)BO for without-replacement sampling algorithms for (conditional) bilevel optimization. The authors prove the convergence rate for both algorithms under regular assumptions, which are improved upon existing works on the complexity of $\epsilon$.
strength_weakness: strength_weakness: 1. The improvement from order $\mathcal{O}(\epsilon^{-4})$ to $\mathcal{O}(\max\\{m,n\\}^p\epsilon^{-3})$ is not as great as stated to be $\mathcal{O}(\epsilon^{-3})$. If we can easily ignore the existence of $m$ and $n$, why not calculate the true gradients and Hessian matrices (instead of stochastic estimates) all time and reach a $\mathcal{O}(\max\\{m,n\\}\epsilon^{-2})$ rate of SOBA? Please clarify the influence of term $\\{m,n\\}^p$ so that the degree of improvement is clearer. If $p$ can be forced to 0 without additional assumptions, the result will be better.
2. The $p$ values in both $\mathcal{O}(\max\\{m,n\\}^p\epsilon^{-3})$ for bilevel optimization and $\mathcal{O}(\max\\{m,n\\}^{2p}\epsilon^{-4})$ for conditional bilevel optimization are not clear. Please clarify how $p$ is determined. By the way, the authors should have used another character since $p$ has already been used as the dimension of $x$.
3. The motivation of using without-replacement sampling is not solid. In the introduction section, the author stated that using without-replacement sampling achieves more accurate gradient estimations than with-replacement sampling. However, SOBA implemented by with-replacement sampling can achieve $\mathcal{O}(\epsilon^{-2})$ complexity which is better than $\mathcal{O}(\epsilon^{-3})$. Could you give some explanations on this?
4. The novelty in algorithmic design is limited. It looks quite like a direct combination of SOBA and natural without-replacement sampling strategies.
questions: questions: 1. The paper regards the special cases as one of the main contributions. Can the special cases be fully covered by the general cases? Are convergence results or algorithms for these cases different from what we could achieve by applying the general case?
2. How is $p$ (or $q$) in the convergence rate determined? If it can be forced to 0 by using certain sampling strategy, why not directly use that strategy? Will it contradict with other assumptions?
3. Why should we use without-replacement sampling rather than with-replacement sampling? I'm asking because I believe SOBA with moving average and with-replacement sampling has a better convergence rate of $\mathcal{O}(\epsilon^{-2})$ under regular assumptions. Please correct me if I'm wrong.
4. Can you compare the emprical performance with SOBA?
limitations: limitations: The limitations are stated as assumptions.

Comments: title: Any further concerns?
Additional Experiments
comment: Thank you for the fair comment. Here's another question I'm concerned about. Lets assume we have $p<1$ for some sampling methods, then the proposed method has a better rate on $\max\\{m,n\\}$ but a worse rate on $\epsilon$ than deterministic SOBA, while the stochastic SOBA has a even better rate on $\max\\{m,n\\}$ but worse rate on $\epsilon$ than the proposed method. Thus, it seems that no matter which part dominates the convergence rate, the proposed method is not the best one in theory. Taking both $\max\\{m,n\\}$ and $\epsilon$ into considerations, it still remains unclear why the proposed method are theoretically better, and it just appears to be a trade-off choice. It would be better if the authors can theoretically justify the advantage on the complexities with any of the sampling strategies mentioned.
Sorry for the confusion! Since in your original review, it reads "SOBA with moving average and with-replacement sampling", we interpret this as meaning the standard stochastic case (which involves sampling).

Firstly, we agree with the reviewer that the value of $p$ is important and we indeed include it in Table 1 of the manuscript and also discuss them in the main text. As stated in the rebuttal, the dependence over $max(m,n)^p$ comes from the sample gradient error $||\frac{1}{k} \sum_{\tau = t}^{t+k-1} \nabla f (x_t,y_t; \xi_{\tau}^{\pi}) - \nabla f (x_t,y_t)||^2$, while the value of $p$ depends on specific sampling orders. In particular, random-reshuffling and the GraB method can force $p < 1$, which has a better dependence over $max(m,n)$ compared to deterministic algorithms. Additionally, we want to add that similar dependence over $max(m,n)^p$  exists for  single level algorithms that adopt without-replacement sampling.

Next, we believe our work makes a valuable contribution to the study of bilevel optimization. Specifically, we demonstrate that the independent sampling assumption, commonly used in SGD-type bilevel algorithms, can be replaced by without-replacement sampling, leading to a reduced per-iteration cost (less number of backward passes) and an improved convergence rate (with the improvement depending on specific example orders). Furthermore, the new algorithms (WiOR-BO/WiOR-CBO) still perform SGD-type updates, making them practical for modern machine learning, in contrast to other variance-reduction-based or deterministic algorithms.
Thanks for your rebuttal. 

The overall response looks good yet I'm still not satisfied with the part related to SOBA. In the rebuttal the authors say the convergence rate of SOBA is $\mathcal{O}(\epsilon^{-4})$, which is true for the stochastic case. However, in my original review, I'm talking about SOBA with full-batch true gradient, which is a deterministic case. With a single modification on SOBA (projecting variable $z^k$ to a ball of radius no greater than $C_f/\mu_g$), the convergence rate of deterministic SOBA can achieve $\mathcal{O}(\epsilon^{-2})$, thus $\mathcal{O}(\max\\{m,n\\}\epsilon^{-2})$ computation complexity in total. That's why I believe the dependence on $\max\\{m,n\\}^p$ is also of great importance and we cannot only compare the order of $\epsilon$. In view of this, when using shuffle-once where $p=1$, I do not see any advantage in the complexity bound compared with deterministic SOBA.
Thanks for your reply. Could you please elaborate any concerns that may still be unclear? We would be very happy to discuss them further.
Dear reviewer,

Thank you once again for spending time reviewing our manuscript. With the discussion period ending soon, we wanted to kindly check if your concerns have been fully resolved. If you have any further questions, please don't hesitate to reach out and we are very happy to discuss them.

Best regards,

Authors
Below, we show experimental results for one more task:  Risk-averse portfolio management. In this task, we are given a set of assets with their returns at different time slots and then we need to design a portfolio with the objective of maximum mean return and minimum variance across time slots. This task can be formulated as a *compositional optimization problem* (See [1] for a more detailed description). We compare our WiOR-Comp-RR and WiOR-Comp-SO with the SCGD [2], for our methods, we set $\gamma=0.99$, $\rho=0.99$ and $\eta=0.001$, and for SCGD, we set $\alpha=0.001$ and $\beta=0.99$. Table below shows the objective loss over number of iterations (we set batch-size as 100):

|  Iters     | 1     | 500   | 1000  | 1500  | 2000  |
|-------|-------|-------|-------|-------|-------|
| SCGD  | 0.985 | 0.379 | 0.137 | 0.057 | 0.029 |
| WiOR-Comp-SO | 0.982 | 0.335 | 0.122 | 0.041 | 0.015 |
| WiOR-Comp-RR | 0.913 | 0.331 | 0.118 | 0.037 | 0.014 |

As shown by the experimental results, our algorithms outperforms the baseline SCGD.

Finally, we would like to once again thank the reviewer for their time and effort in reviewing our manuscript. Please do not hesitate to reach out if you have any further questions.

Reference

[1]. Chen, Tianyi, Yuejiao Sun, and Wotao Yin. "Solving stochastic compositional optimization is nearly as easy as solving stochastic optimization." IEEE Transactions on Signal Processing 69 (2021): 4937-4948.

[2]. Wang, Mengdi, Ethan X. Fang, and Han Liu. "Stochastic compositional gradient descent: algorithms for minimizing compositions of expected-value functions." Mathematical Programming 161 (2017): 419-449.
Thank you for your replies. This addresses some of the points I raised and I decide to maintain my score.
Thanks for your reply!
Thanks for your Reply!
Thanks to the authors' responses, I have no further questions and will keep my score.
Thanks for your rebuttal, I would like to keep my rating.