Mathieu Dagréou
email: mathieu.dagreou@inria.fr

Pierre Ablin
email: pierre.ablin@cnrs.fr

Samuel Vaiter
email: samuel.vaiter@cnrs.fr

Thomas Moreau
email: thomas.moreau@inria.fr

A framework for bilevel optimization that enables stochastic and global variance reduction algorithms

T ) convergence rate, and that it achieves linear convergence under Polyak-Łojasciewicz assumption. This is the first stochastic algorithm for bilevel optimization that verifies either of these properties. Numerical experiments validate the usefulness of our method.

We also note that the computation of these directions does not require to compute the matrices ∇ 2 11 G(z, x) and ∇ 2 21 G(z, x): we only need to compute their product with a vector, which can be computed at a cost similar to that of computing a gradient.

Introduction

Bilevel optimization is attracting more and more attention in the machine learning community thanks to its wide range of applications. Typical examples are hyperparameters selection [START_REF] Bengio | Gradient-Based Optimization of Hyperparameters[END_REF][START_REF] Pedregosa | Hyperparameter optimization with approximate gradient[END_REF][START_REF] Franceschi | Bilevel programming for hyperparameter optimization and meta-learning[END_REF][START_REF] Bertrand | Implicit differentiation of lasso-type models for hyperparameter optimization[END_REF], data augmentation [START_REF] Ekin | AutoAugment: Learning Augmentation Strategies From Data[END_REF][START_REF] Rommel | CADDA: Classwise Automatic Differentiable Data Augmentation for EEG Signals[END_REF], implicit deep learning [START_REF] Shaojie Bai | Deep Equilibrium Models[END_REF] or neural architecture search [START_REF] Liu | Darts: Differentiable architecture search[END_REF]. Bilevel optimization aims at minimizing a function whose value depends on the result of another optimization problem: min

x∈R d h(x) = F (z * (x), x), such that z * (x) ∈ arg min z∈R p G(z, x) , (1)
where F and G are two real valued functions defined on R p × R d . G is called the inner function, F is the outer function and h is the value function. Similarly, z is the inner variable and x is the outer variable. In most cases, the function z * can only be approximated by an optimization algorithm, which makes bilevel optimization problems challenging. Under appropriate hypotheses, the function h is differentiable, and the chain rule and implicit function theorem give for any

x ∈ R d ∇h(x) = ∇ 2 F (z * (x), x) + ∇ 2 21 G(z * (x), x)v * (x) , (2)
36th Conference on Neural Information Processing Systems (NeurIPS 2022).

where v * (x) ∈ R p is the solution of a linear system

v * (x) = -∇ 2 11 G(z * (x), x) -1 ∇ 1 F (z * (x), x) . (3)
In the light of ( 2) and (3), it turns out that the derivation of the gradient of h at each iteration is cumbersome because it involves two subproblems: the resolution of the inner problem to find an approximation of z * (x) and the resolution of a linear system to find an approximation of v * (x). It makes the practical implementation of first order methods like gradient descent for (1) challenging. As is the case in many machine learning problems, we suppose in this paper that F and G are empirical means:

F (z, x) = 1 m m j=1 F j (z, x), G(z, x) = 1 n n i=1 G i (z, x) .
This structure suggests the use of stochastic methods to solve [START_REF] Akhtar | Projection-Free Algorithm for Stochastic Bi-level Optimization[END_REF]. For single-level problems (that is, classical optimization problems where one function should be minimized), using Stochastic Gradient Descent (SGD; [START_REF] Robbins | A stochastic approximation method[END_REF][START_REF] Bottou | Large-Scale Machine Learning with Stochastic Gradient Descent[END_REF]) and variants is natural because individual gradients are straightforward unbiased estimators of the gradient. In the bilevel framework, we want to develop algorithms that make progress on problem (1) by using only a few functions F j and G i at a time. However, since ∇h involves the inverse of the Hessian of G, building such stochastic algorithms is quite challenging, one of the difficulties being that there is no straightforward unbiased estimator of ∇h. Still, in settings where m or n are large, where computing even a single evaluation of F or G is extremely expensive, stochastic methods are the only scalable algorithms.

Variance reduction [START_REF] Johnson | Accelerating stochastic gradient descent using predictive variance reduction[END_REF][START_REF] Defazio | SAGA: A Fast Incremental Gradient Method With Support for Non-Strongly Convex Composite Objectives[END_REF][START_REF] Schmidt | Minimizing finite sums with the stochastic average gradient[END_REF][START_REF] Fang | SPIDER: Near-Optimal Non-Convex Optimization via Stochastic Path Integrated Differential Estimator[END_REF][START_REF] Cutkosky | Momentum-based variance reduction in non-convex SGD[END_REF] is a popular technique to obtain fast stochastic algorithms. In a single-level setting, these methods build an approximation of the gradient of the objective function using only stochastic gradients. Contrary to SGD, the variance of the approximation goes to 0 as the algorithm progresses, allowing for faster convergence. For instance, the SAGA method [START_REF] Defazio | SAGA: A Fast Incremental Gradient Method With Support for Non-Strongly Convex Composite Objectives[END_REF] achieves linear convergence if the objective function satisfies a Polyak-Łojasciewicz inequality, and O( 1 T ) convergence rate on smooth non-convex functions [START_REF] Sashank | Fast Incremental Method for Nonconvex Optimization[END_REF]. The extension of these methods to bilevel optimization is a natural idea to develop faster algorithms. However, this idea is hard to implement because it is hard to derive unbiased estimators of ∇h, let alone variance reduction ones. Contributions. We introduce a novel framework for bilevel optimization in Section 2, where the inner variable, the solution of the linear system (3) and the outer variable evolve jointly. The evolution directions are written as sums of derivatives of F j and G i , which allows us to derive simple unbiased stochastic estimators. In this framework, we propose SOBA, an extension of SGD (Section 2.1), and SABA (Section 2.2), an extension of the variance reduction algorithm SAGA [START_REF] Defazio | SAGA: A Fast Incremental Gradient Method With Support for Non-Strongly Convex Composite Objectives[END_REF]. In Section 3 we analyse the convergence of our methods. SOBA is shown to achieve

inf t≤T E[∥∇h(x t )∥ 2 ] = O(log(T )T - 1
2 ) with decreasing step sizes. We prove that SABA with fixed step sizes achieves 1

T T t=1 E[∥∇h(x t )∥ 2 ] = O( 1 T )
. SABA is therefore, to the best of our knowledge, the first stochastic bilevel algorithm that matches the convergence rate of gradient descent on h. We also prove that SABA achieves linear convergence under the assumption that h satisfies a Polyak-Łojasciewicz inequality. To the best of our knowledge, SABA is also the first stochastic bilevel algorithm to feature such a property. Importantly, these rates match the rates of the single level counterparts of each algorithm in non-convex setting (SGD for SOBA and SAGA for SABA). Finally, in Section 4, we provide an extensive benchmark of many stochastic bilevel methods on hyperparameters selection and data hyper-cleaning, and illustrate the usefulness of our approach. Related work. The bilevel optimization problem has a strong history in the optimization community, taking root in game theory [START_REF] Heinrich Von | Theory of the market economy[END_REF]. Gradient-based algorithms to solve (1) can be mainly classified in two different categories depending on how ∇h is computed, by automatic or implicit differentiation. Since the solution of the inner problem z * (x) is approximated by the output of an iterative algorithm, it is possible to use automatic differentiation [START_REF] Wengert | A simple automatic derivative evaluation program[END_REF][START_REF] Linnainmaa | Taylor expansion of the accumulated rounding error[END_REF] to approximate ∇h(x). It consists in differentiating the different steps of the inner optimization algorithm -see [START_REF] Gunes Baydin | Automatic differentiation in Machine Learning: A survey[END_REF] for a review -and has been applied successfully to several bilevel problems arising in machine learning [START_REF] Domke | Generic methods for optimization-based modeling[END_REF][START_REF] Franceschi | Forward and reverse gradient-based hyperparameter optimization[END_REF]. One of the main drawbacks of this approach is that it requires to store in memory each iterate of the inner optimization algorithm, although this problem can sometimes be overcome using invertible optimization algorithms [START_REF] Maclaurin | Gradient-based hyperparameter optimization through reversible learning[END_REF] or truncated backpropagation [START_REF] Shaban | Truncated Backpropagation for Bilevel Optimization[END_REF]. The use of the implicit function theorem to obtain (2) and ( 3) is known as implicit differentiation [START_REF] Bengio | Gradient-Based Optimization of Hyperparameters[END_REF]. While the cost of computing exactly (2) can be prohibitive for large scale problems, Pedregosa [START_REF] Pedregosa | Hyperparameter optimization with approximate gradient[END_REF] showed that we can still converge to a stationary point of the problem by using approximate solutions of the inner problem and linear system (3), if the approximation error goes to 0 sufficiently quickly. The complexity of approximate implicit differentiation has been studied in [START_REF] Grazzi | On the iteration complexity of hypergradient computation[END_REF]. Ramzi et al. [START_REF] Ramzi | SHINE: SHaring the INverse Estimate from the forward pass for bi-level optimization and implicit models[END_REF] propose to reuse the computations done in the forward pass to approximate the solution of the linear system (3) when the inner problem is solved thanks to a quasi-Newton method. In the last few years, several works have proposed different strategies to solve [START_REF] Akhtar | Projection-Free Algorithm for Stochastic Bi-level Optimization[END_REF] in a stochastic fashion. A first set of methods relies on two nested loops: one inner loop to solve the inner problem with a stochastic method, and one outer loop to update the outer variable with an approximate gradient direction. In [START_REF] Ghadimi | Approximation Methods for Bilevel Programming[END_REF][START_REF] Ji | Bilevel optimization: Convergence analysis and enhanced design[END_REF][START_REF] Chen | Closing the Gap: Tighter Analysis of Alternating Stochastic Gradient Methods for Bilevel Problems[END_REF] the authors use several SGD iterations for the inner problem and then use stochastic Neumann approximations to get an estimate solution of the linear system, which provides them with an approximation of ∇h used to update x. The analysis of this kind of method was refined by Chen et al. [START_REF] Chen | Closing the Gap: Tighter Analysis of Alternating Stochastic Gradient Methods for Bilevel Problems[END_REF], allowing to achieve the same convergence rates as those of SGD. The convergence of the hypergradient when using stochastic solvers for the inner problem and the linear system has been studied in [START_REF] Grazzi | Convergence properties of stochastic hypergradients[END_REF]. Arbel and Mairal [START_REF] Arbel | Amortized Implicit Differentiation for Stochastic Bilevel Optimization[END_REF] replace the Neumann approximation by SGD steps to estimate [START_REF] Shaojie Bai | Deep Equilibrium Models[END_REF]. Other authors have proposed single loop algorithms, alternating steps in the inner and the outer problem. Hong et al. [START_REF] Hong | A Two-Timescale Framework for Bilevel Optimization: Complexity Analysis and Application to Actor-Critic[END_REF] propose to perform Neumann approximations of the inverse Hessian and use a single SGD step for the inner problem. It was refined in [START_REF] Guo | On Stochastic Moving-Average Estimators for Non-Convex Optimization[END_REF] and [START_REF] Yang | Provably Faster Algorithms for Bilevel Optimization[END_REF] where the optimization procedure uses a momentum acceleration. Other variations around this idea include [START_REF] Huang | BiAdam: Fast Adaptive Bilevel Optimization Methods[END_REF][START_REF] Khanduri | A Near-Optimal Algorithm for Stochastic Bilevel Optimization via Double-Momentum[END_REF][START_REF] Chen | A Single-Timescale Stochastic Bilevel Optimization Method[END_REF][START_REF] Guo | Randomized Stochastic Variance-Reduced Methods for Multi-Task Stochastic Bilevel Optimization[END_REF][START_REF] Li | A Fully Single Loop Algorithm for Bilevel Optimization without Hessian Inverse[END_REF]. We refer to Table 1 in appendix for a detailed comparison of these methods. Notation. The set of integers between 1 and n (included) is denoted [n]. For f : R p × R d → R, we denote ∇ i f (z, x) its gradient w.r.t. the i th variable. The Hessian of f with respect to the first variable is denoted ∇ 2 11 f (z, x) ∈ R p×p , and the cross-derivatives matrix is

∇ 2 21 f (z, x) ∈ R d×p . If v is a vector, ∥v∥ is its Euclidean norm. If M is a matrix, ∥M ∥ is its spectral norm.
A function is said to be L-smooth, for L > 0, if it is differentiable, and its gradient is L-Lipschitz.

Proposed framework Algorithm 1 General framework

Input:

initializations z 0 ∈ R p , x 0 ∈ R d , v 0 ∈ R p , number of iterations T , step size sequences (ρ t ) t<T and (γ t ) t<T . for t = 0, . . . , T -1 do Update z: z t+1 = z t -ρ t D t z , Update v: v t+1 = v t -ρ t D t v , Update x: x t+1 = x t -γ t D t x , where D t z , D t v and D t x are unbiased esti- mators of D z (z t , v t , x t ), D v (z t , v t , v t ) and D x (z t , v t , x t ).

end for

In this section, we introduce our framework in which the solution of the inner problem, the solution of the linear system (3) and the outer variable all evolve at the same time, following directions that are written as a sum of derivatives of F j and G i . We define

D z (z, v, x) = ∇ 1 G(z, x) , (4)
D v (z, v, x) = ∇ 2 11 G(z, x)v + ∇ 1 F (z, x) , (5) D x (z, v, x) = ∇ 2 21 G(z, x)v + ∇ 2 F (z, x) . ( 6
) These directions are motivated by the fact that we have ∇h

(x) = D x (z * (x), v * (x), x), with z * (x) the minimizer of G(•, x) and v * (x) the solution of ∇ 2 11 G(z * (x), x)v = -∇ 1 F (z * (x), x).
When x is fixed, we approximate z * by doing a gradient descent on G, following the direction -D z (z, v, x). Finally, when z and x are fixed, we find v * by following the direction -D v (z, v, x), which corresponds to a gradient descent on v → 1 2 ⟨∇ 2 11 G(z, x)v, v⟩ + ⟨∇ 1 F (z, x), v⟩. The rest of the paper is devoted to the study of the global dynamics where the three variables z, v and x evolve at the same time, following stochastic approximations of D z , D v and D x . The next proposition motivates the choice of these directions. Proposition 2.1. Assume that for all

x ∈ R d , G(•, x) is strongly convex. If (z, v, x) is a zero of (D z , D v , D x ), then z = z * (x), v = v * (x) and ∇h(x) = 0.
The framework we propose is summarized in Algorithm 1. It consists in following a joint update rule in (z, v, x) that follows directions D t z , D t v and D t x that are unbiased estimators of D z , D v , D x . The first and most important remark is that whereas ∇h cannot be written as a sum over samples, the directions D z , D v and D x involve only simple sums, since their expressions are "linear" in F and G:

D z (z, v, x)= 1 n n i=1 ∇ 1 G i (z, x) , (7)
D v (z, v, x)= 1 n n i=1 ∇ 2 11 G i (z, x)v + 1 m m j=1 ∇ 1 F j (z, x) , (8)
D x (z, v, x)= 1 n n i=1 ∇ 2 21 G i (z, x)v + 1 m m j=1 ∇ 2 F j (z, x) . (9)
It is therefore straightforward to derive unbiased estimators of these directions. In [START_REF] Li | A Fully Single Loop Algorithm for Bilevel Optimization without Hessian Inverse[END_REF], the authors considered one particular case of our framework, where each direction is estimated by using the STORM variance reduction technique (see [START_REF] Cutkosky | Momentum-based variance reduction in non-convex SGD[END_REF]). Taking a step back by proposing the framework summarized in Algorithm 1 opens the way to potential new algorithms that implement other techniques that exist in stochastic single level optimization. In what follows, we study two of them.

First example: the SOBA algorithm

The simplest unbiased estimator is obtained by replacing each mean by one of its terms chosen uniformly at random, akin to what is done in classical single-level SGD. We call the resulting algorithm SOBA (StOchastic Bilevel Algorithm). To do so, we choose two independent random indices i ∈ [n] and j ∈ [m] uniformly and estimate each term coming from G using G i and each term coming from F using F j . This gives the unbiased SOBA directions

D t z = ∇ 1 G i (z t , x t ) , D t v = ∇ 2 11 G i (z t , x t )v t + ∇ 1 F j (z t , x t ) , D t x = ∇ 2 21 G i (z t , x t )v t + ∇ 2 F j (z t , x t ) . (10a) (10b) (10c)
This provides us with a first algorithm, SOBA, where we plug Equations (10a) to (10c) in Algorithm 1.

We defer its analysis to the next section. Importantly, we use different step sizes for the update in (z, v) and for the update in x. We use the same step size in z and in v since the inner problem and the linear system have similar conditioning, which is that of ∇ 2 11 G(z t , x t ). The need for a different step size for the outer and inner problem is clear: both problems can have a different conditioning. An important remark for SOBA is that all the stochastic directions used are computed at the same point z t , v t and x t with the same indices (i, j). The update of z, v and x can thus be performed in parallel instead of sequentially, benefiting from hardware parallelism. Moreover, this enables to share the computations between the different directions. This is the case in hyperparameters selection where G i (z, x) = ℓ i (⟨z, d i ⟩) + x 2 ∥z∥ 2 , with d i a training sample, and ℓ i that measures how good is the prediction ⟨z, d i ⟩. In this setting, we have

∇ 1 G i (z, x) = ℓ ′ i (⟨z, d i ⟩)d i + xz and ∇ 2 11 G i (z, x)v = ℓ ′′ i (⟨z, d i ⟩)⟨v, d i ⟩d i .
The prediction ⟨z, d i ⟩ can thus be computed only once to obtain both quantities. For more complicated models, where automatic differentiation is used to compute the different derivatives and Jacobian-vector products, we can store the computational graph only once to compute at the same time

∇ 1 G i (z, x), ∇ 2 11 G i (z, x)v and ∇ 2 21 G i (z, x)v
, requiring only one backward pass, thanks to the R technique [START_REF] Barak | Fast exact multiplication by the hessian[END_REF]. Finally, like all single loop bilevel algorithms, our method updates at the same time the inner and outer variable, avoiding unnecessary optimization of the inner problem when x is far from the optimum.

Global variance reduction with the SABA algorithm

In classical optimization, SGD fails to reach optimal rates because of the variance of the gradient estimator. Variance reduction algorithms aim at reducing this variance, in order to follow directions that are closer to the true gradient, and to achieve superior practical and theoretical convergence. In our framework, since the directions D z , D v and D x are all written as sums of derivatives of F j and G i , it is easy to adapt most classical variance reduction algorithms. We focus on the celebrated SAGA algorithm [START_REF] Defazio | SAGA: A Fast Incremental Gradient Method With Support for Non-Strongly Convex Composite Objectives[END_REF]. The extension we propose is called SABA (Stochastic Average Bilevel Algorithm). The general idea is to replace each sum in the directions D by a sum over a memory, updating only one term at each iteration. To help the exposition, we denote y = (z, x, v) the vector of joint variables. Since we have sums over i and over j, we have two memories for each variable: w t i for i ∈ [n] and wt j for j ∈ [m], which keep track of the previous values of the variable y.

At each iteration t, we draw two random independent indices i ∈ [n] and j ∈ [m] uniformly and update the memories. To do so, we put w t+1 i = y t and w t+1 i ′ = w t i ′ for i ′ ̸ = i, and wt+1 j = y t and wt+1 j ′ = wt j ′ for j ′ ̸ = j. Each sum in the directions D is then approximated using SAGA-like rules:

given n functions ϕ i ′ for i ′ ∈ [n], we define S[ϕ, w] t i = ϕ i (w t+1 i ) -ϕ i (w t i ) + 1 n n i ′ =1 ϕ i ′ (w t i ′ )
. This is an unbiased estimators of the average of the ϕ's since

E i S[ϕ, w] t i = 1 n n i=1 ϕ i (y t ).
With a slight abuse of notation, we call ∇ 2 11 Gv the sequence of functions

(y → ∇ 2 11 G i (z, x)v) i∈[n] and ∇ 2 21 Gv the sequence of functions (y → ∇ 2 21 G i (z, x)v) i∈[n]
. We define the SABA directions as

D t z = S[∇ 1 G, w] t i , D t v = S[∇ 2 11 Gv, w] t i + S[∇ 1 F, w] t j , D t x = S[∇ 2 21 Gv, w] t i + S[∇ 2 F, w] t j . (11a) (11b) (11c)
These estimators are unbiased estimators of the directions D z , D v and D x . The SABA algorithm corresponds to Algorithm 1 where we use Equations (11a) to (11c) as update directions. When taking a step size γ t = 0 in the outer problem, hereby stopping progress in x, we recover the iterations of the SAGA algorithm on the inner problem. In practice, the sum in S is computed by doing a rolling average (see Appendix B for precision), and the quantities ϕ i (w t i ) are stored rather than recomputed: the cost of computing the SABA directions is the same as that of SGD. It requires an additional memory for the five quantities, of total size n × p + (n + m) × (p + d) floats that can be reduced by using larger batch sizes. Indeed, if b in and b out are respectively the inner and the outer batch sizes, the memory load is reduced to

n b × p + (n b + m b ) × (p × d) with n b = ⌈ n binn ⌉ and m b = ⌈ m
bout ⌉ which are smaller than the number of samples. This memory load can also be reduced in specific cases, for instance when G and F correspond to linear models, where the individual gradients and Hessian-vector products are proportional to the samples. In this case, we only store the proportionality ratio, reducing the memory load to 3n + 2m floats. Like for SOBA, the computations of the new quantities ϕ i (w t+1 i ) are done in parallel, thus benefiting from hardware acceleration and shared computations. Despite this memory load, using SAGA-like variance reduction instead of STORM as done in [START_REF] Li | A Fully Single Loop Algorithm for Bilevel Optimization without Hessian Inverse[END_REF][START_REF] Yang | Provably Faster Algorithms for Bilevel Optimization[END_REF][START_REF] Khanduri | A Near-Optimal Algorithm for Stochastic Bilevel Optimization via Double-Momentum[END_REF] has the advantage to bring the variance of the estimate directions to zero, enabling faster O( 1 T ) convergence. In the next section, we show that SABA is fast. It essentially has the same properties as SAGA: despite being stochastic, it converges with fixed step sizes, and reaches the same rate of convergence as gradient descent on h.

Theoretical analysis

In this section, we provide convergence rates of SOBA and SABA under some classical assumptions. Note that, unlike most of the stochastic bilevel optimization papers, we work in finite sample setting rather than the more general expectation setting. Actually, SABA does not make any sense for functions that don't have a finite sum structure. However, we stress that SOBA could be studied in a more general setting to obtain the same bounds as here. Also, the finite sum setting is still interesting since doing empirical risk minimization is very common in practice in machine learning. The proofs and the constants in big-O are deferred in Appendix C.

Background and assumptions

We start by stating some regularity assumptions on the functions F and G. Assumption 3.1. The function F is twice differentiable. The derivatives ∇F and ∇ 2 F are Lipschitz continuous in (z, x) with respective Lipschitz constants L F 1 and L F 2 . Note that the above assumption is typically verified in the machine learning context, e.g., when F is the ordinary least squares (OLS) loss or the logistic loss. Assumption 3.2. The function G is three times continuously differentiable on

R p × R d . For any x ∈ R d , G( • , x) is µ G -strongly convex. The derivatives ∇G, ∇ 2 G and ∇ 3 G are Lipschitz continuous in (z, x) with respective Lipschitz constants L G 1 , L G 2 and L G 3 .
Strong convexity and smoothness with respect to z of G are verified when G is a regularized leastsquares/logistic regression with a full rank design matrix, when the data is not separable for the logistic regression. Moreover, the strong convexity ensures the existence and uniqueness of the inner optimization problem for any x ∈ R d . Assumption 3.3. There exists C F > 0 such that for any x we have

∥∇ 1 F (z * (x), x)∥ ≤ C F .
This assumption, combined with the strong convexity of G( • , x), shows boundedness of v * . This assumption holds, for instance, in the case of hyperparameters selection for a Ridge regression problem. Note that in Assumptions 3.1 and 3.2, we assume more regularity of F and G than in stochastic bilevel optimization literature (see for instance [START_REF] Ghadimi | Approximation Methods for Bilevel Programming[END_REF][START_REF] Hong | A Two-Timescale Framework for Bilevel Optimization: Complexity Analysis and Application to Actor-Critic[END_REF][START_REF] Ji | Bilevel optimization: Convergence analysis and enhanced design[END_REF][START_REF] Arbel | Amortized Implicit Differentiation for Stochastic Bilevel Optimization[END_REF]). It is necessary to get the smoothness of v * which will allow to adapt the proof of Chen et al. [START_REF] Chen | Closing the Gap: Tighter Analysis of Alternating Stochastic Gradient Methods for Bilevel Problems[END_REF] and get tight convergence rates. The following lemma gives us some smoothness properties of the considered directions that will be useful to derive convergence rates of our methods. Lemma 3.4. Under the Assumptions 3.1 to 3.3, there exist constants L z , L v and L x such that

∥D z (z, v, x)∥ 2 ≤ L 2 z ∥z -z * (x)∥ 2 , ∥D v (z, v, x)∥ 2 ≤ L 2 v (∥z -z * (x)∥ 2 + ∥v -v * (x)∥ 2 ) and ∥D x (z, v, x) -∇h(x)∥ 2 ≤ L 2 x (∥z -z * (x)∥ 2 + ∥v -v * (x)∥ 2
). In first order optimization, a fundamental assumption on the objective function is the smoothness assumption. In the case of vanilla gradient descent applied to a function f , it allows to get a convergence rate of ∥∇f (x t )∥ 2 in O(1/T ), i.e. convergence to a stationary point [START_REF] Iu | Introductory Lectures on Convex Optimization: A Basic Course[END_REF]. The following lemma proved by Ghadimi and Wang [18, Lemma 2.2] ensures the smoothness of h. Lemma 3.5. Under the Assumptions 3.1 to 3.3, the function h is L h -smooth for some L h > 0.

The constant L h is specified in Appendix C.3. As usual with the analysis of stochastic methods, we define the expected norms of the directions

V t z = E[∥D t z ∥ 2 ], V t v = E[∥D t v ∥ 2 ] and V t x = E[∥D t x ∥ 2 ]
, where the expectation is taken over the past. Thanks to variance-bias decomposition, they are the sum of the variance of the stochastic direction and the squared-norm of the unbiased direction. For SOBA, we use classical bounds on variances like those found for instance in [START_REF] Hong | A Two-Timescale Framework for Bilevel Optimization: Complexity Analysis and Application to Actor-Critic[END_REF]: Assumption 3.6. There exists B z , B v and B x such that for all t,

E t [∥D t z ∥ 2 ] ≤ B 2 z (1 + ∥D z (z t , v t , x t )∥ 2 ) and E t [∥D t v ∥ 2 ] ≤ B 2 v (1 + ∥D v (z t , v t , x t )∥ 2
) where E t denotes the expectation conditionally to (z t , v t , x t ).

For SOBA and SABA, we need to bound the expected norm of D t

x . For SABA, this assumption allows to get a the same sample complexity as SAGA for single level problems. Assumption 3.7. There exists B x such that for all t, E t [∥D t

x ∥ 2 ] ≤ B 2 x . Assumptions 3.6 and 3.7 are verified for instance, if all the G i and ∇ 1 G i have at most quadratic growth, and if F has bounded gradients. They are also verified if the iterates remain in a compact set. Note that we do not assume that G has bounded gradients, as this would contradict its strong-convexity. Finally, for the analysis of SABA, we need regularity on each G i and F j : Assumption 3.8. For all i ∈ [n] and j ∈ [m], the functions ∇G i , ∇F j , ∇ 2 11 G i and ∇ 2 21 G i are Lipschitz continuous in (z, x).

Fundamental descent lemmas

Our analysis for SOBA and SABA is based on the control of both

δ t z = E[∥z t -z * (x t )∥ 2 ] and δ t v = E[∥v t -v * (x t )∥ 2 ],
Strong convexity of G and smoothness of z * (x) and v * (x) allow to obtain the following lemma by adapting the proof of Chen et al. [START_REF] Chen | Closing the Gap: Tighter Analysis of Alternating Stochastic Gradient Methods for Bilevel Problems[END_REF]. In what follows, we drop the dependency of the step sizes ρ and γ in t for clarity. Lemma 3.9. Assume that γ 2 ≤ min

µ G L 2 * 4B 2 x L 2 zx , µ G L 2 * 8B 2 x L 2 vx ρ.
We have:

δ t+1 z ≤ 1 - ρµ G 4 δ t z + 2ρ 2 V t z + β zx γ 2 V t x + β zx γ 2 ρ E[∥D x (z t , v t , x t )∥ 2 ] δ t+1 v ≤ 1 - ρµ G 8 δ t v + β vz ρδ t z + 2ρ 2 V t v + β vx γ 2 V t x + β zx γ 2 ρ E[∥D x (z t , v t , x t )∥ 2 ]
where

β zx = β vx = 3L 2 * , β zx = 8L 2 * µ G , β vx = 16L 2 *
µ G , L * is the maximum between the Lipschitz constants of z * and v * (see Lemma C.1),

β vz = 1 µ 3 G (L F 1 µ G + L G 2 ) 2
, L zx and L vx are respectively the smoothness constants of z * and v * .

We insist that this result is obtained in general for Algorithm 1 with arbitrary unbiased directions. We can therefore invoke this lemma for the analysis of both SOBA and SABA. We use the smoothness of h to get the following lemma, which is similar to [9, Lemma 1]. Lemma 3.10. Let h t = E[h(x t )] and g t = E[∥∇h(x t )∥2 ]. We have

h t+1 ≤ h t - γ 2 g t - γ 2 E[∥D x (z t , v t , x t )∥ 2 ] + γ 2 L 2 x (δ t z + δ t v ) + L h 2 γ 2 V t x . If z t = z * (x t ), v t = v * (x t ), that is δ z , δ v both cancel and D x (z t , v t , x t ) = ∇h(x t
), we get an inequality reminiscent of the smoothness inequality for SGD on h.

Analysis of SOBA

The analysis of SOBA is based on Lemmas 3.5 and 3.9. We have the following theorem, with fixed step sizes depending on the number of iterations:

Theorem 1 (Convergence of SOBA, fixed step size). Fix an iteration T > 1 and assume that Assumptions 3.1 to 3.7 hold. We consider fixed steps ρ t = ρ √ T and γ t = ξρ t with ρ and ξ precised in the appendix. Let (x t ) t≥1 the sequence of outer iterates for SOBA. Then,

1 T T t=1 E[∥∇h(x t )∥ 2 ] = O(T -1 2 ) .
As opposed to [START_REF] Hong | A Two-Timescale Framework for Bilevel Optimization: Complexity Analysis and Application to Actor-Critic[END_REF], we do not need that the ratio γ ρ goes to 0, which allows to get a complexity (that is, the number of call to oracles to have an ϵ-stationary solution) in O(ϵ -2 ) better than the Õ(ϵ -5

2 ) they have. Also, note that this rate is the same as the one of SGD for non-convex and smooth objective [START_REF] Bottou | Optimization methods for large-scale Machine Learning[END_REF]. We obtain a similar rate using decreasing step sizes:

Theorem 2 (Convergence of SOBA, decreasing step size). Assume that Assumptions 3.1 to 3.7 hold. We consider steps ρ t = ρt -1 2 and γ t = ξρ. Let x t the sequence of outer iterates for SOBA. Then,

inf t≤T E[∥∇h(x t )∥ 2 ] = O(log(T )T -1 2 ) .
As for SGD, SOBA suffers from the need of decreasing step sizes to get actual convergence because of the variance of the estimation on each directions. On the other hand, the analysis of SABA leverages the dynamic of all three variables, resulting in fast convergence with fixed step sizes.

SABA: a stochastic method with optimal rates

In what follows, we denote N = n + m the total number of samples. The following theorem shows O(N

2 3 T -1
) convergence for the SABA algorithm in the general case where we only assume smoothness of h. Our analysis of SABA is inspired by the analysis of single-level SAGA by Reddi et al. [START_REF] Sashank | Fast Incremental Method for Nonconvex Optimization[END_REF].

Theorem 3 (Convergence of SABA, smooth case). Assume that Assumptions 3.1 to 3.3 and 3.7 to 3.8 hold. We suppose ρ = ρ ′ N -23 and γ = ξρ, where ρ ′ and ξ depend only on F and G and are specified in appendix. Let x t the iterates of SABA. Then,

1 T T t=1 E[∥∇h(x t )∥ 2 ] = O N 2 3 T -1 .
To prove the theorem, the idea is to control the distance from the memory to the current variables. We define

S t = 1 n n i=1 ∥y t -w t i ∥ 2 + 1 m m j=1 ∥y t -wt j ∥ 2 .
In appendix, we show that we can find scalars ϕ s , ϕ z , ϕ v > 0 such that the quantity L t = h t +ϕ s S t +ϕ z δ t z +ϕ v δ t v satisfies L t+1 ≤ L tγ 2 g t . Summing these inequalities for t = 1 . . . T and using the fact that L t is lower bounded demonstrates the theorem. Note that the step sizes are constant with respect to the time, but they scale with N -2 3 . As a consequence, the sample complexity is O(N non-convex single level problems [START_REF] Sashank | Fast Incremental Method for Nonconvex Optimization[END_REF]. This is better than the sample complexity of Algorithm 1 with full batch directions, which is O(N ϵ -1 ). Hence, with SABA, we get the best of both worlds: the stochasticity makes the scaling in N of the sample complexity goes from N in full batch mode to N 2 3 for SABA, and the variance reduction makes the scaling in ϵ goes from ϵ -2 for SOBA to ϵ -1 for SABA. Our experiments in Section 4 confirm this gain. Furthermore, if we assume that h satisfies a Polyak-Łojasiewicz (PL) inequality, we recover linear convergence. Recall that h has the PL property if there exists µ h > 0 such that for all x ∈ R d ,

1 2 ∥∇h(x)∥ 2 ≥ µ h (h(x) -h * ) with h * the minimum of h.
Theorem 4 (Convergence of SABA, PL case). Assume that h satisfies the PL inequality and that Assumptions 3.1 to 3.3 and 3.7 to 3.8 hold. We suppose ρ = ρ ′ N -2 3 and γ = ξρ ′ N -1 , where ρ ′ and ξ depend only on F and G and are specified in appendix. Let x t the iterates of SABA and c ′ ≜ min µ h , 1 16P ′ with P ′ specified in the appendix. Then,

E[h T ] -h * = (1 -c ′ γ) T (h 0 -h * + C 0 )
where C 0 is a constant specified in appendix that depends on the initialization of z, v, x and memory.

The proof is similar to that of the previous theorem: we find coefficients ϕ s , ϕ z , ϕ v such that

L t = h t + ϕ s S t + ϕ z δ t z + ϕ v δ t v satisfies the inequality L t+1 ≤ (1 -c ′ γ)L t
, which is then unrolled. Note that in the case where we initialize z and v with z 0 = z * (x 0 ), v 0 = v * (x 0 ), and the memories w 0 i = w 0 , w0 j = w 0 for all i, j, the constant C 0 cancels and the bound simplifies to

E[h(x T )] -h * ≤ (1 -c ′ γ) T (h(x 0 ) -h * ).
Just like classical variance reduction methods in single-level optimization, this theorem shows that our method achieves linear convergence under PL assumption on the value function. To the best of our knowledge, our method is the first stochastic bilevel optimization method that enjoys such property. We note that the PL hypothesis is more general than µ h -strong convexity of h -it is a necessary condition for strong convexity. We see here the importance of global variance reduction. Indeed, using variance reduction only on z and SGD on x would lead to sub-linear convergence in x. This would be the case even with a perfect estimation of z * (x). Similarly, using variance reduction only on x and SGD on z would lead to sub-linear convergence in z, and hence in x. Using global variance reduction with respect to each variable as we propose here is the only way to achieve linear convergence. We now turn to experiments, where we find that our method is also promising from a practical point of view.

Experiments

Here we compare the performances of SOBA and SABA with competitor methods on different tasks. The different methods being compared are stocBiO [START_REF] Ji | Bilevel optimization: Convergence analysis and enhanced design[END_REF], AmiGO [START_REF] Arbel | Amortized Implicit Differentiation for Stochastic Bilevel Optimization[END_REF], FSLA [START_REF] Li | A Fully Single Loop Algorithm for Bilevel Optimization without Hessian Inverse[END_REF], MRBO [START_REF] Yang | Provably Faster Algorithms for Bilevel Optimization[END_REF], TTSA [START_REF] Hong | A Two-Timescale Framework for Bilevel Optimization: Complexity Analysis and Application to Actor-Critic[END_REF], BSA [START_REF] Ghadimi | Approximation Methods for Bilevel Programming[END_REF] and SUSTAIN [START_REF] Khanduri | A Near-Optimal Algorithm for Stochastic Bilevel Optimization via Double-Momentum[END_REF]. A detailed account of the experiments is provided in Appendix B.

Hyperparameters selection

The first task we perform is hyperparameters selection to choose regularization parameters on ℓ2 logistic regression. Let us denote ((d train i , y train i )) 1≤i≤n and ((d val i , y val i )) 1≤i≤m the training and the validation sets. In this case, the inner variable θ corresponds to the parameters of the model, and the outer variable λ to the regularization. The functions F and G of the problem (1) are the logistic loss, with ℓ 2 penalty for G, that is to say

F (θ, λ) = 1 m m i=1 φ(y val i ⟨d val i , θ⟩) and G(θ, λ) = 1 n n i=1 φ(y train i ⟨d train i , θ⟩) + 1 2 p k=1 e λ k θ 2
k where φ(u) = log(1 + e -u ). We fit a binary classification model on the IJCNN11 dataset. Here, n = 49 990, m = 91 701 and p = 22. The suboptimality gap is plotted in Figure 2a for each method. The lowest values are reached by AmIGO and SABA, but much quicker for SABA. Moreover, SABA is the only single-loop method that reaches a suboptimality below 10 -3 . Among all methods, SOBA is the first to reach it best value, but this value is also the highest. The gap between SOBA and SABA highlights the benefits of variance reduction: it gives us a lower plateau and the fixed step sizes enable faster convergence.

Data hyper-cleaning

The second task we perform is data hyper-cleaning introduced in [START_REF] Franceschi | Forward and reverse gradient-based hyperparameter optimization[END_REF] on the MNIST 2 dataset. The data is patitioned into a training set (d train i , y train i ), a validation set (d val i , y val i ), and a test set. The training set contains 20000 samples, the validation set 5000 samples and the test set 10000 samples. The targets y take values in {0, . . . , 9} and the samples x are in dimension 784. Each sample in the training set is corrupted with probability p: a sample is corrupted when we replace its label y i by a random label in {0, . . . , 9}. Samples in the validation and test sets are not corrupted. The goal of datacleaning is to train a multinomial logistic regression on the train set and learn a weight per training sample, that should go to 0 for corrupted samples. This is formalized by the bilevel optimization problem [START_REF] Akhtar | Projection-Free Algorithm for Stochastic Bi-level Optimization[END_REF] with

F (θ, λ) = 1 m m i=1 ℓ(θd val i , y val i ) and G(θ, λ) = 1 n n i=1 σ(λ i )ℓ(θd train i , y train i ) + C r ∥θ∥ 2
where ℓ is the cross entropy loss and σ is the sigmoid function. The inner variable θ is a matrix of size 10 × 784, and the outer variable λ is a vector in dimension n train = 20000. For the estimated parameters θ during optimization, we report in Figure 2b the test error, i.e., the percent of wrong predictions on the testing data. We use for this experiment a corruption probability p = 0.5. In general, the error decreases quickly until it reaches a final value. We observe that our method SABA outperforms all the other methods by reaching faster its smallest error, which is smaller than the ones of the other methods. For SOBA, it reaches a lower final error than stocBiO and BSA. In appendix, we provide other convergence curves, and find that for higher values of p, SABA is still the fastest algorithm to reach its final accuracy, but the algorithm MRBO can reach a lower test accuracy. Overall, we find that among all methods, even those that implement variance reduction (that is FSLA, MRBO, SUSTAIN, SABA), SABA is the one that demonstrates the best empirical performance.

Conclusion

In this paper, we have presented a framework for bilevel optimization that enables the straightforward development of stochastic algorithms. The gist of our framework is that the directions in Equations ( 4) to ( 6) are all written as simple sums of samples derivatives. We leveraged this fact to propose SOBA, an extension of SGD to our framework, and SABA, an extension of SAGA to our framework, which both achieve similar convergence rates as their single level counterparts. Finally, we think that our framework opens a large panel of potential methods for stochastic bilevel optimization involving techniques of extrapolation, variance reduction, momentum and so on.

A Extensive comparison between stochastic methods for bilevel optimization

We provide here tables summarizing other methods in stochastic bilevel optimization. They are grouped between methods that are based on two nested loops and methods that use only one loop.

In the following tables, the inner iterations are referred with the variable k and the outer iterations are referred with the variable t (or T for the total number of iterations).

In the literature, there are three main ways to perform Hessian inversion. The HIA, first proposed in [START_REF] Ghadimi | Approximation Methods for Bilevel Programming[END_REF], and SHIA, proposed in [START_REF] Ji | Bilevel optimization: Convergence analysis and enhanced design[END_REF], procedures used for Hessian inversion are precised in Algorithm 2 and 3. These methods are based on Neumann approximation of the inverse of a matrix. SGD for Hessian inversion refers to Stochastic Gradient Descent on v → 1 2 ⟨∇ 2 11 G(z, x)v, v⟩-⟨∇ 1 F (z, x), v⟩. The complexity refers to the number of call to the oracles to get an ϵ-stationary solution. In these complexities, the notation Õ hide polynomial factors in log ϵ -1 .

Algorithm 2 Hessian Inverse Approximation (HIA)

Input: variables z ∈ R p , x ∈ R d , gradient ∇ 1 F (z, x) ∈ R p , maximum number of iterations b, a parameter η. Set v 0 = ∇ 1 F (z, x) Choose p ∈ {0, . . . , b -1} randomly. for k = 1, . . . , p do Sample i ∈ [n] Update v : v k+1 = (I -η∇ 2 11 G(z, x))v k end for Return: bηv p+1 Algorithm 3 Summed Hessian Inverse Approximation (SHIA) Input: variables z ∈ R p , x ∈ R d , gradient ∇ 1 F (z, x) ∈ R p , maximum number of iterations b, a parameter η. Set v 0 = ∇ 1 F (z, x) Set s 0 = v 0 for k = 0 . . . , b -1 do Sample i ∈ [n] Update v: v k+1 = (I -η∇ 2 11 G(z, x))v k Update s: s k+1 = s k + v k+1 end for Return: ηs b
The momentum column refers to the use of STORM [START_REF] Cutkosky | Momentum-based variance reduction in non-convex SGD[END_REF] momentum in the inner loop or the outer loop. This momentum can be applied to either the inner or the implicit gradient estimate. If we consider the current estimate y t = (z t , v t , x t ) and the previous estimate y t-1 = (z t-1 , v t-1 , x t-1 ), and we apply STORM to the quantity ϕ(y t ) with the memory φt , the momentum update rule reads

φ(t+1) = ηϕ(y t ) + (1 -η)( φt + ϕ(y t ) -ϕ(y t-1 )) .
Note that this update requires to evaluate the quantity ϕ twice per iteration, once in y t and once in y t-1 . The memory is need to store the previous estimates y t-1 as well as the running estimate of the gradient φ. )

MRBO [46] SHIA SGD with STORM Yes (STORM) O(t -1/3 ) O(t -1/3 ) Õ(ϵ -3/2 ) STABLE [10] Direct SGD No O(T -1/2 ) O(T -1/2 ) O(ϵ -2 ) SUSTAIN [27] HIA SGD with STORM Yes (STORM) O(t -1/3 ) O(t -1/3 ) O(ϵ -3/2
) SVRB [START_REF] Guo | Randomized Stochastic Variance-Reduced Methods for Multi-Task Stochastic Bilevel Optimization[END_REF] Direct + momentum

SGD with momentum

Yes

O(t -1/3 ) O(t -1/3 ) Õ(ϵ -3 ) SBFW [1] HIA SGD No O(t -1/2 ) O(T -3/4 ) Õ(ϵ -4
) FSLA [START_REF] Li | A Fully Single Loop Algorithm for Bilevel Optimization without Hessian Inverse[END_REF] SGD with STORM SGD with STORM

Yes (STORM) O(t -1/2 ) O(T -1/2 ) O(ϵ -2 ) SOBA SGD step SGD No O(t -1/2 ) O(t -1/2
)

O(ϵ -1/2 ) SABA SAGA step SAGA No Constant Constant O((n + m) 2/3 ϵ -1
)

Table 1: Comparison of the stochastic bilevel optimization solvers in the literature.

The complexity represents the number of oracle calls necessary to attain an ϵ accurate stationary point.

B Details on experiments

We provide here additional informations on the experiments.

B.1 Generalities

All the experiments are performed with Python, using the package Benchopt [START_REF] Moreau | Benchopt: Reproducible, efficient and collaborative optimization benchmarks[END_REF]. For each problem, we use oracles for a function given function f that (f (z, x),

∇ 1 f (z, x), ∇ 2 11 f (z, x)v, ∇ 2 21 f (z, x)v
) avoiding duplicate computation of intermediate results for these quantities. We find that using mini-batches instead of individual samples to compute the stochastic estimates allowed for much faster computations, thanks to hardware acceleration and vectorization of the computations. We use continuous batches to avoid random memory access that slow down the computations. Concretely, if i b is the index of the current batch and B is the batch-size, the indices of the corresponding samples are those in the set {i b × B, . . . , (i b + 1) × B -1}. By doing so, the samples in a same batch are contiguous in memory, which facilitates the access. We use a batch-size of 64 in all experiments. For the methods involving an inner loop (stocBiO, BSA, AmIGO), we perform 10 inner steps at each outer iteration as proposed in the papers which introduced these methods. For the approximate Hessian vector product, we perform 10 steps per outer iteration for each methods using HIA (BSA, TTSA, SUSTAIN), SHIA (MRBO, stocBiO) or SGD (AmIGO) for the inversion of the linear system. For the step sizes, they all have the form ρ t = α/t a and γ t = β/t b . For the pair of exponents (a, b), we choose the theoretical one from the original papers, that is (1/2, 1/2) for BSA and FSLA, (1/3, 1/3) for MRBO and SUSTAIN, (0, 0) for SABA, AmIGO and stocBiO, (2/5, 3/5) for TTSA and SOBA. For (α, β), we perform a grid search (the grid is precised in the subsection dedicated to each experiment) and we keep for each method, the pair (α, β) that gives the lowest value of h (for the hyperparameters) or the lowest test accuracy (for the data cleaning task) in median over 10 runs for each possible pair. When we use HIA or SHIA for the Hessian inversion, we set η = α since the Hessian inversion problem has the same conditioning as the inner optimization problem.

For the STORM's momentum parameter in MRBO and SUSTAIN, we take 0.5/t 2/3 .

For SABA, we have to maintain the estimate S[ϕ, w] i t = ϕ i (w t+1 i

) -ϕ i (w t i ) + 1 n n i ′ =1 ϕ i ′ (w t i ′ ) of 1
n n i=1 ϕ i (y t ) (see Section 2.2 for the notations). The sum inside S is maintain by performing a rolling mean on the past gradients computed. More precisely,

A t = 1 n n i ′ =1 ϕ i ′ (w t i ′ ).
To get A t+1 , instead of computing the summing all the gradients stored, which has O(n) computational complexity, we do

A t+1 = A t + 1 n (ϕ i (w t+1 i ) -ϕ i (w t i ))
, which is equivalent mathematically but has O(1) computational complexity.

B.2 Hyperparameter selection on a toy problem

The Figure 1 corresponds to the methods SABA et SOBA applied to an hyperparameter selection problem for a Ridge regression. We generate 1000 samples x 1 , . . . , x 1000 ∈ R 1 0 for N (0, I 1 0). We generate a parameter β ∼ N (0, I 10 ) and do y = (X ⊙ W )β + ϵ where ϵ ∼ N (0, 0.01I 10 ) and the entries of W have the form

W i,j = 1 + u j v i,j with v i,j ∼ U([0, 1]) and u j ∼ U([0, 1]) if 1 ≤ j ≤ 5 or u j ∼ U([0, 10]) if 6 ≤ j ≤
F (θ, λ) = 1 2n val n val i=1 ((x val i ) ⊤ θ -y val i ) 2 and G(θ, λ) = 1 2n train ntrain i=1 ((x train i ) ⊤ θ -y train i ) 2 + λ 2 ∥θ∥ 2
with n train = 750 and n val = 250.

B.3 Hyperparameters selection on IJCNN1

In this experiment, we select the parameters regularization for a multiregularized logistic regression model precised in Equations ( 12) and ( 13) where we have one hyperparameter per feature

F (θ, λ) = 1 m m i=1
φ(y val i ⟨d val i , θ⟩) and ( 12)

G(θ, λ) = 1 n n i=1 φ(y train i ⟨d train i , θ⟩) + 1 2 θ ⊤ diag(e λ1 , . . . , e λp )θ . ( 13
)
Note that the parametrization in e λ of the penalty instead of λ can be surprising at first glance, but it is classical in the bilevel optimization literature [START_REF] Pedregosa | Hyperparameter optimization with approximate gradient[END_REF][START_REF] Ji | Bilevel optimization: Convergence analysis and enhanced design[END_REF][START_REF] Grazzi | Convergence properties of stochastic hypergradients[END_REF] because it avoids positivity constraints on λ. In order to choose the select proper parameters (α, β) for each algorithm, we perform a grid search. We search α in a set of 9 values between 2 -5 and 23 spaced on a log scale. For β, we choose r in a set of 7 values between 10 -2 and 10 spaced on a logarithmic scale and we set β = α r . For this experiments, we use Just-In-Time (JIT) compilation thanks to the package Numba [START_REF] Siu Kwan Lam | Numba: A llvm-based python jit compiler[END_REF], to decrease the python overhead in the iteration loop. To evaluate the value function h, we use L-BFGS [START_REF] Liu | On the limited memory BFGS method for large scale optimization[END_REF] to solve compute z * (x t ) and then evaluate the function h(x t ) = F (z * (x t ), x t ). We use Python code with Numba [START_REF] Siu Kwan Lam | Numba: A llvm-based python jit compiler[END_REF] for fast implementation of stochastic methods, and the Python package Benchopt [START_REF] Moreau | Benchopt: Reproducible, efficient and collaborative optimization benchmarks[END_REF] to perform the benchmark 3 . For each algorithm, we plot the median performance over 10 runs. In both plots, SABA achieves the best performance. The dashed lines are for one loop competitor methods, the dotted lines are for two loops methods and the solid lines are the proposed methods. Left: performance in running time, Right: performance in number of gradient/Hessian-vector products sampled.

B.4 Data hyper-cleaning

For the regularization parameter C r , we choose C r = 0.2 after a manual search in order to get the best final test accuracy. In this experiment, the selection of the good pair (α, β) is also performed by grid search. The parameter α is picked in a set of 11 numbers between 10 -3 and 100 spaced on a logarithmic scale.

For β, we choose r in a set of 11 values between 10 -5 and 1 spaced on a logarithmic scale and we set β = α r . Note that in this case, we could not use JIT from Numba since at the moment of the experiment, the softmax function coming from Scipy was not compatible with Numba. We report in Figure B.2 some additional convergence curves with different corruption probabilities p ∈ {0.5, 0.7, 0.9} (the figure in the main text corresponds to p = 0.5). SABA is always the fastest algorithm to reach its final accuracy. For p = 0.9 and p = 0.7, the algorithm MRBO reaches an accuracy roughly 1% smaller than that of SABA, but takes roughly 10 times longer to get there. where θ ∈ R p×C and λ ∈ R C . As for the other experiments, we performed and grid search over 63 pairs (α, β) to set the step sizes. The parameter α is chosen among values between 2 -5 and 2 3 spaced in log scale. For β, we choose it in a set of values between 10 -2 and 10 spaced in log scale. We used a batch size of 64. The experiment took 525 CPU hours. We show in Figure B.3 the error on the test samples with respect to the running time and the number of gradients/Hessian-vector products sampled. We observe that SABA and SOBA achieve the best performances. The initial gap between the first and the second plot for SABA is due to the overhead of the initialization of the memory. This gap can be reduced by increasing the batch size.

Proof. Let (z, v, x) a zero of (D z , D v , D x ). For D z , this means that ∇ 1 G(z, x) = 0. Since G( • , x) is strongly convex, z is the minimizer of G( • , x), i.e. z = z * (x). The fact that (z, v, x) is a zero of D v implies that ∇ 2 11 G(z, x)v = -∇ 1 F (z, x). Replacing z by its value, we get v = -∇ 2 11 G(z * (x), x) -1 ∇ 1 F (z * (x), x) which is v * (x) by definition.
Putting all together and using the expression of ∇h(x) given by ( 2), we get

D x (z, v, x) = ∇ 2 F (z * (x), x) + ∇ 21 G(z * (x), x)v * (x) = ∇h(x) .
On the other hand, D x (z, v, x) = 0 so ∇h(x) = 0.

C.2 Proof of Lemma 3.4

Proof. Let (z, v, x) ∈ R p ×R p ×R d . Using the fact that ∇ 1 G(z * (x), x) = 0 and the L G 1 -smoothness of G( • , x), we have

∥D z (z, v, x)∥ 2 = ∥∇ 1 G(z, x) -∇ 1 G(z * (x), x)∥ 2 ≤ L 2 G ∥z -z * (x)∥ 2 . For D v , since ∇ 2 11 G(z * (x), x)v * (x) = -∇ 1 F (z * (x), x), we write ∥D v ∥ = ∥(∇ 2 11 G(z, x)v + ∇ 1 F (z, x)) -(∇ 2 11 G(z * (x), x)v * (x) + ∇ 1 F (z * (x), x))∥ (14) ≤ ∥[∇ 2 11 G(z, x) -∇ 2 11 G(z * (x), x)]v * (x)∥ + ∥∇ 2 11 G(z, x)[v -v * (x)]∥ (15) + ∥∇ 1 F (z, x) -∇ 1 F (z * (x), x)∥ .
For the first term, we use the Lipschitz continuity of ∇ 2 11 G:

∥[∇ 2 11 G(z, x) -∇ 2 11 G(z * (x), x)]v * (x)∥ ≤ L G 2 ∥z -z * (x)∥∥v * (x)∥ . Then, since G in µ G -strongly convex w.r.t. z, ∇ 1 F (z * (•), •) is bounded and v * (x) = -[∇ 2 11 G(z * (x), x)] -1 ∇ 1 F (z * (x), x), we have ∥[∇ 2 11 G(z, x) -∇ 2 11 G(z * (x), x)]v * (x)∥ ≤ L G 2 C F µ G ∥z -z * (x)∥ . ( 16
)
For the second term, we use the L G 1 -smoothness of G( • , x) and for the third term, we use the L F 1 -smoothness of F and we finally get

∥D v ∥ ≤ L G 2 C F µ G + L F 1 ∥z -z * (x)∥ + L G 1 ∥v -v * (x)∥ . (17)
Then, taking

L v = √ 2 max L G 2 C F µ G + L F , L G 1 , we get ∥D v (z, v, x)∥ 2 ≤ L 2 v (∥z -z * (x)∥ 2 + ∥v -v * (x)∥ 2 ) . (18)
For D x (z, v, x) -∇h(x) we start by writing

∥D x (z, v, x) -∇h(x)∥ ≤ ∥∇ 2 F (z, x) -∇ 2 F (z * (x), x)∥ + ∥∇ 2 21 G(z, x)v -∇ 2 21 G(z * (x), x)v * (x)∥ (19)
≤ ∥∇ 2 F (z, x) -∇ 2 F (z * (x), x)∥ + ∥∇ 2 21 G(z, x)∥∥v -v * (x)∥ (20)
+ ∥v * (x)∥∥∇ 2 21 G(z, x) -∇ 2 21 G(z * (x), x)∥ .
We bound the first term using the fact that ∇ 2 F is L F 1 -Lipschitz continuous. For the second term, the fact that ∇ 2 21 G is bounded thanks to the Lipschitz continuity of ∇ 1 G(z, • ). For the third term, we use that ∇ 2 21 G( • , x) is L G 2 -Lipschitz continuous and the same derivation as Equation ( 16). We finally get

∥D x -∇h(x)∥ ≤ L F 1 + C F L G 2 µ G ∥z -z * (x)∥ + L G 1 ∥v -v * (x)∥ . (21)
Taking

L x = √ 2 max L F 1 + C F L G 2 µ G , L G 1 yields ∥D x (z, v, x) -∇h(x)∥ 2 ≤ L 2 x (∥z -z * (x)∥ 2 + ∥v -v * (x)∥ 2 ) . (22)

C.3 Smoothness constant of h

From Ghadimi and Wang [18, Lemma 2.2], we get the Lemma 3.10 which states the L h -smoothness of h with Lemma C.1. There exists a constant L * > 0 such that for any x 1 , x 2 ∈ R d we have

L h = L F 1 + 2L F 1 L G 2 + C 2 F L G 2 µ G + L G 11 L G 1 C F + L G 1 L G 2 C F + (L G 1 ) 2 L F 1 µ 2 G + (L G 1 ) 2 L G 2 C F µ 3 G . C.4
∥z * (x 1 ) -z * (x 2 )∥ ≤ L * ∥x 1 -x 2 ∥, ∥v * (x 1 ) -v * (x 2 )∥ ≤ L * ∥x 1 -x 2 ∥ . Proof. Let x ∈ R d . The Jacobian of z * is given by dz * (x) = -[∇ 2 11 G(z * (x), x)] -1 ∇ 2 1,2 G(z * (x), x). Thanks to the µ G -strong convexity of G and the fact that ∇ 2 21 G is bounded, we have ∥dz * (x)∥ ≤ L G 1 µ G . Thus, z * is Lipschitz continuous. For ∥v * (x 1 ) -v * (x 2 )∥, we start from the definition v * : ∥v * (x 1 ) -v * (x 2 )∥ = ∥[∇ 2 11 G(z * (x 1 ), x 1 )] -1 ∇ 1 F (z * (x 1 ), x 1 ) -[∇ 2 11 G(z * (x 2 ), x 2 )] -1 ∇ 1 F (z * (x 2 ), x 2 )∥ (23) ≤ ∥([∇ 2 11 G(z * (x 1 ), x 1 )] -1 -[∇ 2 11 G(z * (x 2 ), x 2 )] -1 ∇ 1 F (z * (x 1 ), x 1 )∥ (24) + ∥[∇ 2 11 G(z * (x 2 ), x 2 )] -1 (∇ 1 F (z * (x 2 ), x 2 ) -∇ 1 F (z * (x 1 ), x 1 ))∥ .
For the first term, we use that for any invertible matrix A and B we have

A -1 -b -1 = A -1 (B - A)B -1 to get ∥[∇ 2 11 G(z * (x 1 ), x 1 )] -1 -∇ 2 11 G(z * (x 2 ), x 2 )] -1 ∥ = ∥[∇ 2 11 G(z * (x 1 ), x 1 )] -1 (∇ 2 11 G(z * (x 2 ), x 2 )]- ∇ 2 11 G(z * (x 1 ), x 1 )])[∇ 2 11 G(z * (x 2 ), x 2 )] -1 ∥ ≤ 1 µ 2 G ∥∇ 2 11 G(z * (x 1 ), x 1 ) -∇ 2 11 G(z * (x 2 ), x 2 )∥ ≤ L G 2 µ 2 G ∥(z * (x 1 ), x 1 ) -(z * (x 2 ), x 2 )∥ ≤ L G 2 µ 2 G [∥z * (x 1 ) -z * (x 2 )∥ + ∥x 1 -x 2 ∥] ≤ L G 2 µ 2 G 1 + L G 1 µ G ∥x 1 -x 2 ∥ .
And then, since

∇ 1 F (z * ( • ), • ) is bounded: ([∇ 2 11 G(z * (x 1 ), x 1 )] -1 -[∇ 2 11 G(z * (x 2 ), x 2 )] -1 ∇ 1 F (z * (x 1 ), x 1 ) ≤ C F L G 2 µ 2 G 1 + L G 1 µ G ∥x 1 -x 2 ∥ .
For the second term, the strong convexity of G( • , x) and the fact that

∇ 1 F is Lipschitz continuous lead to ∥[∇ 2 11 G(z * (x 2 ), x 2 )] -1 (∇ 1 F (z * (x 2 ), x 2 ) -∇ 1 F (z * (x 1 ), x 1 ))∥ ≤ 1 µ G ∥∇ 1 F (z * (x 2 ), x 2 ) -∇ 1 F (z * (x 1 ), x 1 )∥ (25)
≤ L F 1 µ F ∥(z * (x 1 ), x 1 ) -(z * (x 2 ), x 2 )∥ (26)
≤ L F 1 µ G [∥z * (x 1 ) -z * (x 2 )∥ + ∥x 1 -x 2 ∥] (27)
≤ L F 1 µ G 1 + L G 1 µ G ∥x 1 -x 2 ∥ . (28)
Then we get

∥v * (x 1 ) -v * (x 2 )∥ ≤ C F L G 2 µ 2 G 1 + L G 1 µ G + L F 1 µ G 1 + L G 1 µ G ∥x 1 -x 2 ∥ . ( 29
)
We conclude by setting

L * = max L G 1 µ G , C F L G 2 µ 2 G 1 + L G 1 µ G + L F 1 µ G 1 + L G 1 µ G .
In what follows, we denote by E t [ • ] the expectation conditionally on z t , v t and x t . We have the smoothness property of z * provided in [9, Lemma 2].

Lemma C.2. Under the Assumptions 3.1, 3.2 and 3.3, the function z * : R d → R p is L zx -smooth with

L zx = L G 2 (1 + L * ) µ G + L G 1 L G 11 (1 + L * ) µ 2 G . ( 30
)
We establish the same result for v * . To this, we need more regularity on G and F .

We can now bound each term using the regularity assumptions on G and F :

∥[∇ 2 11 G(z * (x), x)] -1 (B(x, ϵ) -B(y, ϵ))∥ ≤ 1 µ G (∥∇ 2 11 F (z * (x), x)dz * (x) -∇ 2 11 F (z * (y), y)dz * (y)∥ (34) + ∥∇ 2 12 F (z * (x), x) -∇ 2 12 F (z * (y), y)∥) ≤ 1 µ G (∥∇ 2 11 F (z * (x), x) -∇ 2 11 F (z * (y), y)∥∥dz * (x)∥ (35)
+ ∥dz * (x) -dz * (y)∥∥∇ 2 11 F (z * (y), y)∥ + L F 2 (∥z * (x) -z * (y)∥ + ∥x -y∥) ≤ 1 µ G (L F 2 L * (1 + L * ) + L zx L F 1 + L F 2 (1 + L * ))∥x -y∥ (36)
(37) For the second term:

∥([∇ 2 11 G(z * (x), x)] -1 -[∇ 2 11 G(z * (y), y)] -1 )B(y, ϵ)∥ ≤ 1 µ 2 G ∥∇ 2 11 G(z * (x), x) -∇ 2 11 G(z * (y), y)∥∥B(y, ϵ)∥ (38) ≤ 1 µ 2 G ∥∇ 2 11 G(z * (x), x) -∇ 2 11 G(z * (y), y)∥ (39)
× (∥∇ 2 11 F (z * (x), x)∥∥dz * (x)∥ + ∥∇ 2 12 F (z * (x), x)∥) ≤ (L G 2 + L F 1 )(L * + 1) µ 2 G ∥x -y∥ (40)
For the third term, we have:

∥A(x, ϵ)(∇ 1 F (z * (x), x) -∇ 1 F (z * (y), y))∥ ≤ L F 1 (1 + L * ) µ 2 G ∥(∇ 3 111 G(z * (x), x)|dz * (x).ϵ) (41) + (∇ 3 112 G(z * (x), x)|ϵ)∥∥x -y∥ ≤ (L F 1 + L G 2 )(1 + L * ) µ 2 G ∥x -y∥ (42)
And finally, for the forth term:

∥(A(x, ϵ) -A(y, ϵ))∇ 1 F (z * (y), y)∥ ≤ C F {∥[∇ 2 11 G(z * (x), x)] -1 ∥ (43) × ∥(∇ 3 111 G(z * (x), x)|dz * (x).ϵ) + (∇ 3 112 G(z * (x), x)|ϵ)∥ × ∥[∇ 2 11 G(z * (x), x)] -1 -[∇ 2 11 G(z * (y), y)] -1 ∥ + ∥[∇ 2 11 G(z * (x), x)] -1 -[∇ 2 11 G(z * (y), y)] -1 ∥ × ∥(∇ 3 111 G(z * (x), x)|dz * (x).ϵ) + (∇ 3 112 G(z * (x), x)|ϵ)∥ × ∥[∇ 2 11 G(z * (y), y)] -1 ∥ + ∥[∇ 2 11 G(z * (y), y)] -1 ∥ 2 × (∥(∇ 3 111 G(z * (x), x)|dz * (x).ϵ) -(∇ 3 111 G(z * (y), y)|dz * (y).ϵ)∥ ∥(∇ 3 112 G(z * (x), x)|ϵ) -(∇ 3 112 G(z * (y), y)|ϵ)∥)} (44) ≤ C F 2 2L G 2 (1 + L * ) µ 3 G + L G 3 (1 + L * ) µ 2 G ∥x -y∥ Thus v * is L vx -smooth with L vx = L F 2 L * (1 + L * ) + L zx L F 1 + L F 2 (1 + L * ) µ G +2 (L G 2 + L F 1 )(L * + 1) µ 2 G + C F L G 3 (1 + L * ) µ 2 G +4 C F L G 2 (1 + L * ) µ 3 G .
C.5 Proof of Lemma 3.9

We now provide the proof of Lemma 3.9.

Proof. Inequality for δ z . We start by expanding the square:

∥z t+1 -z * (x t+1 )∥ 2 = ∥z t+1 -z * (x t )∥ 2 + ∥z * (x t+1 ) -z * (x t )∥ 2 (45) -2⟨z t+1 -z * (x t ), z * (x t+1 ) -z * (x t )⟩
We study each member, using the unbiasedness of D t z and the µ G -strong convexity of G( • , x t ):

E t [∥z t+1 -z * (x t )∥ 2 ] = E t [∥z t -z * (x t )∥ 2 ] -2ρE t [⟨D t z , z t -z * (x t )⟩] + ρ 2 E t [∥D t z ∥ 2 ] ( 46
)
= ∥z t -z * (x t )∥ 2 -2ρ⟨∇ 1 G(z t , x t ), z t -z * (x t )⟩ + ρ 2 E t [∥D t z ∥ 2 ] (47) ≤ (1 -ρµ G )∥z t -z * (x t )∥ 2 + ρ 2 E t [∥D t z ∥ 2 ] . ( 48
)
Taking the total expectation yields

E[∥z t+1 -z * (x t )∥ 2 ] ≤ (1 -ρµ G )δ t z + ρ 2 V t z . ( 49
)
The second member is bounded using Lipschitz continuity of z * :

E[∥z * (x t+1 ) -z * (x t )∥ 2 ] ≤ L 2 * E[∥x t+1 -x t ∥ 2 ] = L 2 * γ 2 V t x .
For the remaining scalar product, we have

-2⟨z t+1 -z * (x t ), z * (x t+1 ) -z * (x t )⟩ = -2[⟨z t -z * (x t ), z * (x t+1 ) -z * (x t )⟩ -ρ⟨D t z , z * (x t+1 ) -z * (x t )⟩] . (50)
The second term can be bounded using Cauchy-Schwarz inequality, the Lipschitz-continuity of z * and Young inequality:

E[ρ⟨D t z , z * (x t+1 ) -z * (x t )⟩] ≤ E[ρ∥D t z ∥∥z * (x t+1 ) -z * (x t )∥] (51) ≤ ρL * E[∥D t z ∥∥x t+1 -x t ∥] (52) ≤ ρ 2 2 V t z + L 2 * 2 ∥x t+1 -x t ∥ 2 (53) ≤ ρ 2 2 V t z + L 2 * γ 2 2 V t x . (54)
For -2⟨z t -z * (x t ), z * (x t+1 ) -z * (x t )⟩, we follow the proof of [START_REF] Chen | Closing the Gap: Tighter Analysis of Alternating Stochastic Gradient Methods for Bilevel Problems[END_REF] which consists in making appear the "unbiased part of z * (x t+1 -z * (x t ) by a linear approximation. More precisely, we have

⟨z t -z * (x t ), z * (x t+1 ) -z * (x t )⟩ = ⟨z t -z * (x t ), dz * (x t )(x t+1 -x t )⟩ A (55) ⟨z t -z * (x t ), z * (x t+1 ) -z * (x t ) -dz * (x t )(x t+1 -x t )⟩ B .
For A, we use the unbiasedness of D t x , Cauchy-Schwarz inequality, the Lipschitz continuity of z * (Lemma C.1) and the identity ab ≤ ηa 2 + b 2 η for any η > 0:

-2E[A] = -2γE[⟨z t -z * (x t ), dz * (x t )D t x ⟩] (56) = -2γE[⟨z t -z * (x t ), dz * (x t )E t [D t x ]⟩] (57) = -2γE[⟨z t -z * (x t ), dz * (x t )D x (z t , v t , x t )⟩] (58) ≤ 2γE[∥z t -z * (x t )∥∥dz * (x t )D x (z t , v t , x t )∥] (59) ≤ 2L * γE[∥z t -z * (x t )∥∥D x (z t , v t , x t )∥] (60)
≤ 2ηδ t z + 2L 2 * η γ 2 E[∥D x (z t , v t , x t )∥ 2 ] . (61)
We take η = ρµ G 4 and we get

-2E[A] ≤ ρµ G 2 δ t z + 8L 2 * µ G γ 2 ρ E[∥D x (z t , v t , x t )∥ 2 ] . (62)
For B, we use Cauchy-Schwarz inequality, the smoothness of z * (Lemma C.2), Young inequality and the boundedness of

E t [∥D t x ∥ 2 ] to get -2E[B] ≤ 2E[∥z t -z * (x t )∥∥z * (x t+1 ) -z * (x t ) -dz * (x t )(x t+1 -x t )∥] (63) ≤ L zx E[∥z t -z * (x t )∥∥x t+1 -x t ∥ 2 ] (64) ≤ L zx νE[∥z t -z * (x t )∥ 2 ∥x t+1 -x t ∥ 2 ] + L zx ν E[∥x t+1 -x t ∥ 2 ] (65) ≤ L zx νγ 2 E[∥z t -z * (x t )∥ 2 E t [∥D t x ∥ 2 ]] + L zx γ 2 ν V t x (66) ≤ L zx B 2 x νγ 2 δ t z + L zx γ 2 ν V t x . ( 67
)
We take ν = Lzx L 2 * and we get

-2E[B] ≤ L 2 zx B 2 x γ 2 L 2 * δ t z + L 2 * γ 2 V t x (68) Now, using γ 2 ≤ ρµ G L 2 * B 2 x L 2 zx
, we end up with

δ t+1 z ≤ (1 - ρµ G 4 )δ t z + 2ρ 2 V t z + β zx γ 2 V t x + β zx γ 2 ρ E[∥D x (z t , v t , x t )∥ 2 ] , (69)
with

β zx = 3L 2 * and β zx = 8L 2 *
µ G . Inequality for δ v . We proceed in a similar way for v:

δ t+1 v ≤ E[∥v t+1 -v * (x t )∥ 2 ]+E[∥v * (x t+1 )-v * (x t )∥ 2 ]-2E[⟨v t+1 -v * (x t ), v * (x t+1 )-v * (x t )⟩] .
(70) For the first term, we have

E t [∥v t+1 -v * (x t )∥ 2 ] = ∥v t -v * (x t )∥ 2 -2ρ⟨D v (z t , v t , x t ), v t -v * (x t )⟩ + ρ 2 E t [∥D t v ∥ 2 ] (71) Now, using that D v (z * (x t ), v * (x t ), x t ) = 0: ⟨D v (z t , v t , x t ), v t -v * (x t )⟩ = ⟨D v (z t , v t , x t ) -D v (z * (x t ), v * (x t ), x t ), v t -v * (x t )⟩ (72) = ⟨∇ 2 11 G(z t , x t )(v t -v * (x t )), v t -v * (x t )⟩ (73) + ⟨(∇ 2 11 G(z t , x t ) -∇ 2 11 G(z * (x t ), x t ))v * (x t ), v t -v * (x t )⟩ + ⟨(∇ 1 F (z t , x t ) -∇ 1 F (z * (x t ), x t )), v t -v * (x t )⟩ ≥ µ G ∥v t -v * (x t )∥ 2 - L G 2 C F µ G ∥z t -z * (x t )∥∥v t -v * (x t )∥ (74) -L F 1 ∥z t -z * (x t )∥∥v t -v * (x t )∥ ≥ µ G ∥v t -v * (x t )∥ 2 -ω∥z t -z * (x t )∥∥v t -v * (x t )∥ ( 75
)
where

ω = L F 1 + L G 2 C F µ G . We then use ω∥z t -z * (x t )∥∥v t -v * (x t )∥ ≤ 1 2 c∥v t -v * (x t )∥ 2 + ω 2 2c ∥z t - z * (x t )∥ 2 with c = µ G to get -⟨D v (z t , v t , x t ), v t -v * (x t )⟩ ≤ - 1 2 µ G δ t v + ω 2 2µ G δ t z .
We get the overall inequality by taking the total expectation

E[∥v t+1 -v * (x t )∥ 2 ] ≤ 1 - ρµ G 2 δ t v + ρ ω 2 2µ G δ t z + ρ 2 V t v .
We also use Lipschitz on v * to bound the other term

E[∥v * (x t+1 ) -v * (x t )∥ 2 ] ≤ L 2 * γ 2 V t x .
As previously, the scalar product is bounded by: -

E[⟨v t+1 -v * (x t ), v * (x t+1 ) -v * (x t )⟩] = -E[⟨v t -v * (x t ), v * (x t+1 ) -v * (x t )⟩] -ρE[⟨D t v , v * (x t+1 ) -v * (x t )⟩] (76) ≤ E[⟨z t -z * (x t ), v * (x t+1 ) -v * (x t )⟩] + ρ 2 2 V t v + L 2 * γ 2 2 V t x (77)
We do similar manipulations pour v * , thanks to Lemma C.3. We have as for z from Lemma C.1 for any η > 0:

-

E[⟨v t -v * (x t ), dv * (x t )(x t+1 -x t )⟩] ≤ ηδ t v + L 2 * γ 2 η E[∥D x (z t , v t , x t )∥ 2 ] . (78)
We take η = ρµ G 8 and we get

-E[⟨v t -v * (x t ), dv * (x t )(x t+1 -x t )⟩] ≤ ρµ G 8 δ t v + 8L 2 * γ 2 µ G ρ E[∥D x (z t , v t , x t )∥ 2 ] (79)
(80) Then smoothness of v * for any η > 0 gives us

-E[⟨v t -v * (x t ), v * (x t+1 ) -v * (x t ) -dv * (x t )(x t+1 -x t )⟩] ≤ L vx B 2 x ν 2 γ 2 δ t v + L vx 2ν γ 2 V t x . (81)
With

ν = Lvx L 2 * we get -E[⟨v t -v * (x t ), v * (x t+1 ) -v * (x t ) -dv * (x t )(x t+1 -x t )⟩] ≤ L 2 vx B 2 x 2L 2 * γ 2 δ t v + L 2 * 2 γ 2 V t x . (82)
With the assumption

γ 2 ⩽ ρµ G L 2 * 8L 2 vx B 2 x
, we get

δ t+1 v ≤ 1 - ρµ G 2 + ρµ G 4 + L 2 vx B x L 2 * δ t v + ρβ vz δ t z + 2ρ 2 V t z + 3L 2 * γ 2 V t x + 16L 2 * γ 2 µ G ρ E[∥D x (z t , v t , x t )∥ 2 ] (83) ≤ 1 - ρµ G 8 δ t v + ρβ vz δ t z + 2ρ 2 V t z + 3L 2 * γ 2 V t x + 16L 2 * γ 2 µ G ρ E[∥D x (z t , v t , x t )∥ 2 ] . (84)
And finally we have

δ t+1 v ≤ 1 - ρµ G 8 δ t v + ρβ vz δ t z + 2ρ 2 V t z + β vx γ 2 V t x + β vx γ 2 ρ E[∥D x (z t , v t , x t )∥ 2 ] (85)
with

β vz = ω 2 2µ G , β vx = 3L 2 * and β vx = 16L 2 * γ 2 µ G .
C.6 Proof of Lemma 3.10

Proof. We use smoothness of h to get

E t [h(x t+1 )] ≤ h(x t ) -γ⟨D x (z t , v t , x t ), ∇h(x t )⟩ + L h 2 γ 2 E t [∥D t x ∥ 2 ] (86) ≤ h(x t ) - γ 2 (∥∇h(x t )∥ 2 + ∥D x (z t , v t , x t )∥ 2 -∥∇h(x t ) -D x (z t , v t , x t )∥ 2 ) + L h 2 γ 2 E t [∥D t x ∥ 2 ] (87)
where the last inequality comes from the identity ⟨a, b⟩ = 1 2 (∥a∥ 2 + ∥b∥ 2 -∥a -b∥) 2 . We take the total expectation and use the previous Lemma 3.4 to get

h t+1 ≤ h t - γ 2 g t - γ 2 E[∥D x (z t , v t , x t )∥ 2 ] + γL 2 x 2 (δ t z + δ t v ) + L h 2 γ 2 V t x (88)

C.7 Proof of Theorem 1

This section is devoted to the proof of Theorem 1 that we recall here.

Theorem 1 (Convergence of SOBA, fixed step size). Fix an iteration T > 1 and assume that Assumptions 3.1 to 3.7 hold. We consider fixed steps ρ t = ρ √ T and γ t = ξρ t with ρ and ξ precised in the appendix. Let (x t ) t≥1 the sequence of outer iterates for SOBA. Then,

1 T T t=1 E[∥∇h(x t )∥ 2 ] = O(T -1 2 ) .
The values of the differents constants are

ϕ ′ z = 1 8β zx , ϕ ′ v = min 1 8β vx , µ G ϕ ′ z 32β vz , ρ = min 16 µ G , µ G 16L 2 z B 2 z , µ G 32L 2 v B 2 v , β vz L 2 v B 2 v ,
and

ξ 2 = µ G 4 min min 1 L 2 zx , 1 L 2 vx L 2 * B 2 x ρ , min (ϕ ′ v , ϕ ′ z ) 1 2L 2 x .
Before, one has to adapt our descent lemmas to the case of SOBA.

Lemma C.4. Assume that the step sizes ρ and γ verify ρ ≤ min

µ G 16L 2 z B 2 z , µ G 32L 2 v B 2 v , βvz L 2 v B 2 v and γ 2 ≤ min ρµ G L 2 * 4B 2 x L 2 zx , ρµ G L 2 * 8B 2 x L 2 vx . Then it holds δ t+1 z ≤ 1 - ρµ G 8 δ t z + 2ρ 2 B 2 z + β zx γ 2 B 2 x + β zx γ 2 ρ E[∥D x (z t , v t , x t )∥ 2 ] (89)
δ t+1 v ≤ 1 - ρµ G 16 δ t v + 2β vz ρδ t z + 2ρ 2 B 2 v + β vx γ 2 B 2 x + β vx E[∥D x (z t , v t , x t )∥ 2 ] . (90)
Proof. From Assumption 3.6 and Lemma 3.4, we have

V t z ≤ B 2 z (1 + D z (z t , v t , x t )) ≤ B 2 z (1 + L 2 z δ t z )
. Plugging this into Equation (69) and using V t

x ≤ B 2 x yields

δ t+1 z ≤ 1 - ρµ G 4 + 2L 2 z B 2 z ρ 2 δ t z + 2ρ 2 B 2 z + β zx γ 2 B 2 x + β zx γ 2 ρ E[∥D x (z t , v t , x t )∥ 2 ] . (91)
Since by assumption ρ

≤ µ G 16L 2 z B 2 z
, we have

δ t+1 z ≤ 1 - ρµ G 8 δ t z + 2ρ 2 B 2 z + β zx γ 2 B 2 x + β zx γ 2 ρ E[∥D x (z t , v t , x t )∥ 2 ] . (92)
For δ t v , Assumption 3.3 and Lemma 3.4 provide us

V t v ≤ B 2 v (1 + L 2 v (δ t z + δ t v ))
. Since the assumptions of Lemma 3.9 are verified, we can plug the previous inequality into Equation (85) to get

δ t+1 v ≤ 1 - ρµ G 8 + 2L 2 v B 2 v ρ 2 δ t v + (β vz ρ + 2L 2 v ρ 2 B 2 v )δ t z + 2ρ 2 B 2 v + β vx γ 2 B 2 x + β vx E[∥D x (z t , v t , x t )∥ 2 ] (93)
which can be simplified using ρ ≤ min

µ G 32L 2 v B 2 v , βvz L 2 v B 2 v to get finally δ t+1 v ≤ 1 - ρµ G 16 δ t v + 2β vz ρδ t z + 2ρ 2 B 2 v + β vx γ 2 B 2 x + β vx E[∥D x (z t , v t , x t )∥ 2 ] . (94)
We can now prove Theorem 1.

Proof. Consider the Lyapunov function L

t = h t + ϕ z δ t z + ϕ v δ t v .
Using the Equations ( 88), ( 69) and (85), we can bound L t+1 -L t :

L t+1 -L t ≤ - γ 2 g t - γ 2 -ϕ z β zx γ 2 ρ -ϕ v β vx γ 2 ρ E[∥D x (z t , v t , x t )∥ 2 ] (95) -ϕ z µ G 8 ρ - L 2 x 2 γ -2ϕ v β vz ρ δ t z -ϕ v µ G 16 ρ - L 2 x 2 γ δ t v + L h 2 + ϕ z β zx + ϕ v β vx B 2 x γ 2 + 2(ϕ z B 2 z + ϕ v B 2 v )ρ 2 . Let ϕ ′ z = ϕ z γ ρ and ϕ ′ v = ϕ v γ ρ
, so that:

L t+1 -L t ≤ - γ 2 g t - γ 2 -ϕ ′ z β zx γ -ϕ ′ v β vx γ E[∥D x (z t , v t , x t )∥ 2 ] (96)
-ϕ ′ z µ G 8 ρ 2 γ - L 2 x 2 γ -2ϕ ′ v β vz ρ 2 γ δ t z -ϕ ′ v µ G 16 ρ 2 γ - L 2 x 2 γ δ t v + L h 2 + ϕ ′ z β zx ρ γ + ϕ v β vx ρ γ B 2 x γ 2 + 2 ϕ ′ z B 2 z ρ γ + ϕ ′ v B 2 v ρ γ ρ 2 .
In order to get a decrease, ϕ ′ z , ϕ ′ v , ρ and γ must verify

     ϕ ′ z β zx + ϕ ′ v β vx ≤ 1 2 L 2 x 2 γ + 2ϕ ′ v β vz ρ 2 γ ≤ ϕ ′ z µ 8 ρ 2 γ L 2 x 2 γ ≤ ϕ ′ v µ G 16 ρ 2 γ (97) Let us take ϕ ′ z = 1 8β zx and ϕ ′ v = min 1 8β vx , µ G ϕ ′ z 32βvz . We have ϕ ′ z β zx + ϕ ′ v β vx ≤ 1 4 < 1 2 and L 2 x 2 γ + 2ϕ ′ v β vz ρ 2 γ ≤ L 2 x 2 γ + ϕ ′ z µ G 16
ρ 2 γ .
If we impose

L 2 x 2 γ + ϕ ′ z µ G 16 ρ 2 γ ≤ ϕ ′ z µ G 8 ρ 2
γ , this combined with the third condition in Equation (97) gives the condition

L 2 x 2 γ 2 ≤ min (ϕ ′ v , ϕ ′ z ) µ G 16 ρ 2 .
We also have the conditions coming from the assumptions of C.4, that is

ρ ≤ ρ = min 16 µ G , µ G 16L 2 z B 2 z , µ G 32L 2 v B 2 v , β vz L 2 v ( 98
)
and

γ 2 ≤ min 1 L 2 zx , 1 L 2 vx µ G L 2 * 4B 2 x ρ.
Let us take ρ = ρ √ T with γ = ξρ where ξ is defined as

ξ 2 ≜ µ G 4 min min 1 L 2 zx , 1 L 2 vx L 2 * B 2 x ρ , min (ϕ ′ v , ϕ ′ z ) 1 2L 2 x . (99)
From now, we have

L t+1 -L t ≤ - γ 2 g t + L h 2 B 2 x γ 2 + (ϕ ′ z β zx + ϕ ′ v β vx ) B 2 x ργ + 2 ϕ ′ z B 2 z + ϕ ′ v B 2 v ρ 3 γ . (100)
Summing and telescoping yields

1 T T t=1 g t ≤ 2L 1 T γ + L h B 2 x γ + 2 (ϕ ′ z β zx + 2ϕ ′ v β vx ) B 2 x ρ + 4 ϕ ′ z B 2 z + ϕ ′ v B 2 v ρ 3 γ 2 (101) ≤ 2L 1 √ T ξρ + L h B 2 x ξα √ T + (ϕ ′ z β zx + 2ϕ ′ v β vx ) B 2 x α √ T + 4 ϕ ′ z B 2 z + ϕ ′ v B 2 v α ξ 2 √ T (102) (103)
and so

1 T T t=1 g t = O 1 √ T .
C.8 Proof of Theorem 2

Proof. In the decreasing step size case, we take ρ t = ρ √ t and γ t = ξρ t where ρ is defined in Equation ( 98) and ξ is defined in Equation (99). We recall the integral majorization:

T t=1 t -1 ≤ 1 + T 1 t -1 dt = 1 + log(T ) .
With such definition of ρ t and γ t , Equation (100) is still valid for any t ⩾ 1. The only difference is that the step sizes decrease with t. Hence, by summing and rearranging in Equation (100), we get

T t=1 γ t g t ≤ 2L 1 + L h + 2 (ϕ ′ z β zx + 2ϕ ′ v β vx ) B 2 x 1 ξ + 4 ϕ ′ z B 2 z + ϕ ′ v B 2 v 1 ξ 3 T t=1 (γ t ) 2 (104)
The left-hand-side in Equation ( 104) can be lower bounded by

T t=1 γ t g t ≥ inf t∈[T ] g t ξρ T t=1 t -1 2 ≥ inf t∈[T ] g t ξρT 1 2 . ( 105
)
Also we have

T t=1 (γ t ) 2 = ξ 2 ρ 2 T t=1 t -1 ≤ ξ 2 ρ 2 (1 + log(T )) . (106)
Plugging Equations ( 105) and (106) into Equation (104) and rearranging give

inf t∈[T ] g t ≤ 2L 1 ξρ √ T + ξρ L h + 2 (ϕ ′ z β zx + 2ϕ ′ v β vx ) B 2 x 1 ξ + 4 ϕ ′ z B 2 z + ϕ ′ v B 2 v 1 ξ 3 1 + log(T ) √ T (107) that is to say inf t∈[T ] g t = O 1 √ T + log(T ) √ T . (108)
C.9 Proof of Theorem 3

In this section, we prove Theorem 3 that we recall here Theorem 3 (Convergence of SABA, smooth case). Assume that Assumptions 3.1 to 3.3 and 3.7 to 3.8 hold. We suppose ρ = ρ ′ N -2 3 and γ = ξρ, where ρ ′ and ξ depend only on F and G and are specified in appendix. Let x t the iterates of SABA. Then,

1 T T t=1 E[∥∇h(x t )∥ 2 ] = O N 2 3 T -1 .

29

The constants ρ ′ and ξ are given by

ρ ′ = min K 1 K 5 , K 2 K 5 2 5 , K 3 K 5 5 7 , K 4 K 5 1 3 , µ G 64L 2 z , β zx 2β zx , µ G 128(L 2 v + L ′′ v ) , β vz 8(L 2 v + L ′′ v ) , β vx 2β vx (109) and ξ = min(K 1 , K 2 (ρ ′ ) -1 2 , K 3 (ρ ′ ) -3 2 , K 4 (ρ ′ ) -1 ) (110)
where

ϕ ′′ z = 1 32β zx , ϕ ′′ v = min 1 32β vx , ϕ ′′ z µ G 128β vz , K 1 = min ϕ ′′ z µ G 32L 2 x , ϕ ′′ v µ G 48L 2 x , L ′ z 2L ′ x β zx , L ′ v 2L ′ x β vx , K 2 = min µ G 64β zx L ′′ x , µ G 128β vx L ′′ x , β vz 4L ′′ x β vx , K 3 = ϕ ′′ v µ G 384ϕ ′′ z L ′′ x , K 4 = min 1 4L h , L 2 x 2L h L ′′ x , Γ ′ 6L h L ′ x , 1 8P ′ , ϕ ′′ z µ G 32β ′ sz , ϕ ′′ v µ G 48β ′ sv and K 5 = 15(ϕ ′′ z L ′ z + ϕ ′′ v L ′ v ) Γ ′ .

C.9.1 Control of distance from memory to iterates

We can view our method has having two "parallel" memories for each variable (z t i , v t i , x t i ) for i ∈ 1[n] corresponding to calls in G and (z ′t j , v ′t j , x ′t j ) for j ∈ [m] corresponding to calls to F . At each iteration, we sample i at random uniformly and do (z

t+1 i , v t+1 i , x t+1 i ) = (z t , v t , x t ) and (z t+1 i ′ , v t+1 i ′ , x t+1 i ′ ) = (z t i ′ , v t i ′ , x t i ′
) for i ′ ̸ = i, and similarly for the other memory. In what follows, we focus on controlling the error between the iterates and the memories. We define to make things simpler

E t z = 1 n n i=1 E[∥z t -z t i ∥ 2 ] , E t v = 1 n n i=1 E[∥v t -v t i ∥ 2 ] , E t x = 1 n n i=1 E[∥x t -x t i ∥ 2 ] ,
and similarly E ′t x , E ′t v and E ′t x . Lemma C.5. We have the following inequalities:

E t+1 z ≤ 1 - 1 2n E t z + ρ 2 E∥D t z ∥ 2 + 2nρ 2 E[∥D z (z t , v t , x t )∥ 2 ] , E t+1 v ≤ 1 - 1 2n E t v + ρ 2 E∥D t v ∥ 2 + 2nρ 2 E[∥D v (z t , v t , x t )∥ 2 ] , E t+1 x ≤ 1 - 1 2n E t x + γ 2 E∥D t x ∥ 2 + 2nγ 2 E[∥D x (z t , v t , x t )∥ 2 ] , E ′t+1 z ≤ 1 - 1 2m E t z + ρ 2 E∥D t z ∥ 2 + 2mρ 2 E[∥D z (z t , v t , x t )∥ 2 ] , E ′t+1 v ≤ 1 - 1 2m E t v + ρ 2 E∥D t v ∥ 2 + 2mρ 2 E[∥D v (z t , v t , x t )∥ 2 ] ,
and

E ′t+1 x ≤ 1 - 1 2m E t x + γ 2 E∥D t x ∥ 2 + 2mγ 2 E[∥D x (z t , v t , x t )∥ 2 ] .
Proof. We provide the detailed proof for E t z . The approach for the five others is similar.

Let i ∈ [n]. Taking the expectation of ∥z t+1 -z t+1 i ∥ 2 conditionaly to z t , v t , x t yields E t [∥z t+1 -z t+1 i ∥ 2 ] = 1 n E t [∥z t+1 -z t ∥ 2 ] + n -1 n E t [∥z t+1 -z t i ∥ 2 ] .
Then, using the fact that E t [D t z (z t , v t , x t )] = D z (z t , v t , x t ), we have

E t [∥z t+1 -z t i ∥ 2 ] = E t [∥z t+1 -z t ∥ 2 ] + ∥z t -z t i ∥ 2 -2ρ⟨D z (z t , v t , x t ), z t -z t i ⟩ . (111)
We then upper-bound crudely the scalar product by Cauchy-Schwarz and Young inequalities with parameter β:

E t [∥z t+1 -z t i ∥ 2 ] ≤ E t [∥z t+1 -z t ∥ 2 ] + ρβ -1 ∥D z (z t , v t , x t )∥ 2 + (1 + ρβ)∥z t -z t i ∥ 2
As a consequence, by taking the total expectation and summing for all i ∈ [n], we find

E t+1 z ≤ ρ 2 E[∥D t z ∥ 2 ] + ρβ -1 1 - 1 n E[∥D z (z t , v t , x t )∥ 2 ] + (1 + ρβ) 1 - 1 n E t z .
Finally, we take β = 1 2nρ to obtain

E t+1 z ≤ 1 - 1 2n E t z + ρ 2 E∥D t z (z t , v t , x t )∥ 2 + 2nρ 2 E[∥D z (z t , v t , x t )∥ 2 ] . (112)

C.9.2 Bounds on the variances

The following lemma gives us upper-bounds for

E[∥D t z (z t , v t , x t )∥ 2 ], E[∥D t v (z t , v t , x t )∥ 2 ]
, and E[∥D t

x (z t , v t , x t )∥ 2 ]. Lemma C.6. For SABA, there are constants

L ′ z , L ′ v , L ′ x > 0 such that E[∥D t z (z t , v t , x t )∥] 2 ≤ 2E[∥D z (z t , v t , x t )∥ 2 ] + 2L ′ z (E t z + E t x ) , E[∥D t v (z t , v t , x t )∥ 2 ] ≤ 2E[∥D v (z t , v t , x t )∥ 2 ] + 2L ′ v (E t z + E t x + E t v + E ′t z + E ′t x ) + 2L ′′ v (δ t z + δ t v ) and E[∥D t x (z t , v t , x t )∥ 2 ] ≤ 2E[∥D x (z t , v t , x t )∥ 2 ] + 2L ′ x (E t z + E t x + E t v + E ′t z + E ′t x ) + 2L ′′ x (δ t z + δ t v ) .
Proof. For SABA, if we consider i sampled from [n] at iteration t, we have

D t z = ∇ 1 G i (z t , x t ) -∇ 1 G i (z t i , x t i ) + 1 n n i ′ =1 ∇ 1 G i ′ (z t i ′ , x t i ′ ) .
Hence we get

E t [∥D t z (z t , v t , x t )∥ 2 ] = E t [∥∇ 1 G i (z t , x t ) -∇ 1 G i (z t i , x t i ) + 1 n N i ′ =1 ∇ 1 G i ′ (z t i ′ , x t i ′ ) -∇ 1 G(z t , x t ) + ∇ 1 G(z t , x t )∥ 2 ] ≤ 2∥∇ 1 G(z t , x t )∥ 2 + 2E t [∥∇ 1 G i (z t , x t ) -∇ 1 G i (z t i , x t i ) (113) + 1 n N i ′ =1 ∇ 1 G i ′ (z t i ′ , x t i ′ ) -∇ 1 G(z t , x t )∥ 2 ] .
The second term is the variance of ∇ 1 G i (z t , x t ) -∇ 1 G i (z t i , x t i ), which is therefore upper-bounded by

E t [∥[∇ 1 G i (z t , x t ) -∇ 1 G i (z t i , x t i )∥ 2 ] = 1 n n i=1 ∥[∇ 1 G i (z t , x t ) -∇ 1 G i (z t i , x t i )∥ 2 ≤ L ′ z n n i=1 (∥z t -z t i ∥ 2 + ∥x t -x t i ∥ 2 ) ( 114
)
where the inequality comes from the Lipschitz continuity of each

∇ 1 G i with L ′ z = max i∈[n] L Gi 1 .
Then, by plugging (114) into (113) and taking the total expectation, we get

E[∥D t z (z t , v t , x t )∥] 2 ≤ 2E[∥D z (z t , v t , x t )∥ 2 ] + 2L ′ z (E t z + E t x ) . (115)
Things are quite similar for the other variables, albeit a bit more difficult.

In v, it holds

E t [∥D t v (z t , v t , x t )∥ 2 ] =E t [∥∇ 1 F j (z t , x t ) -∇ 1 F j (z ′t j , x ′t j ) + 1 m m j ′ =1 ∇ 1 F j ′ (z ′t j ′ , x ′t j ′ ) (116)
+ ∇ 2 11 G i (z t , x t )v t -∇ 2 11 G i (z t i , x t i )v t i + 1 n n i ′ =1 ∇ 2 1 G i ′ (z t i ′ , x t i ′ )v i ′t -D v (z t , v t , x t ) + D v (z t , v t , x t )∥ 2 ] ≤2[∥D v (z t , v t , x t )∥ 2 (117)
+ 2E t [∥∇ 1 F j (z t , x t ) -∇ 1 F j (z ′t j , x ′t j ) + 1 m m j ′ =1 ∇ 1 F j ′ (z ′t j ′ , x ′t j ′ ) + ∇ 2 11 G i (z t , x t )v t -∇ 2 11 G i (z t i , x t i )v t i + 1 n n i ′ =1 ∇ 1 G 2 i ′ (z t i ′ , x t i ′ )v t i ′ -D v (z t , v t , x t )∥ 2 ] (118)
Here, we see that we need to control the variance of

∇ 1 F j (z t , x t ) -∇ 1 F j (z ′t j , x ′t j ) + ∇ 2 11 G i (z t , x t )v t -∇ 2 11 G i (z t i , x t i )v t i .
Since i and j are independent, this is a sum of two independent random variables, hence its variance is the sum of the variances, which is upper-bounded by

E t [∥∇ 1 F j (z t , x t ) -∇ 1 F j (z ′t j , x ′t j )∥ 2 ] + E t [∥∇ 2 11 G i (z t , x t )v t -∇ 2 11 G i (z t i , x t i )v t i ∥ 2 ]
.

For E t [∥∇ 1 F j (z t , x t ) -∇ 1 F j (z ′t j , x ′t j )∥ 2 ]
we use the lipschitz continuity of the ∇ 1 F j :

E t [∥∇ 1 F j (z t , x t ) -∇ 1 F j (z ′t j , x ′t j )∥ 2 ] ≤ max j∈[m] L Fj 1 E t [∥z t -z t j ∥ 2 + ∥x t -x t j ∥ 2 ] (119) ≤ max j∈[m] L Fj 1 1 m m j=1 (∥z t -z t j ∥ 2 + ∥x t -x t j ∥ 2 ) . ( 120
)
The control of

E t [∥∇ 2 11 G i (z t , x t )v t -∇ 2 11 G i (z t i , x t i )v t i ∥ 2
] is a bit harder without assuming the boundness of v beforehand. But, we can bypass the difficulty by introducing ∇ 2 11 G i (z * (x t ), x t )v * (x t ):

E t [∥∇ 2 11 G i (z t , x t )v t -∇ 2 11 G i (z t i , x t i )v t i ∥ 2 ] ≤ 4{E t [∥∇ 2 11 G i (z t , x t )(v t -v * (x t ))∥ 2 ] (121) + E t [∥(∇ 2 11 G i (z t , x t ) -∇ 2 11 G i (z * (x t ), x t ))v * (x t )∥ 2 ] + E t [∥(∇ 2 11 G i (z * (x t ), x t ) -∇ 2 11 G i (z t i , x t i ))v * (x t )∥ 2 ] + E t [∥∇ 2 11 G i (z t i , x t i )(v * (x t ) -v t i )∥ 2 ]} ≤ 4((max i∈[n] L Gi 1 )∥v t -v * (x t )∥ 2 + (max i∈[n] L Gi 2 ) C F µ G ∥z t -z * (x t )∥ 2 (122)
+ (max i∈[n] L Gi 2 ) C F µ G (∥x t -x t i ∥ 2 + 2(∥z t -z * (x t )∥ 2 + ∥z t -z t i ∥ 2 )) + (max i∈[n] L Gi 1 )(∥x t -x t i ∥ 2 + 2(∥v t -v * (x t )∥ 2 + ∥v t -v t i ∥ 2 )) Let L ′ v = 4 max 2 max i∈[n] L Gi 1 , 2 max i∈[n] L Gi 2 C F µ G , max j∈[m] L Fj 1 and L ′′ v = 4 max 3 max i∈[n] L Gi 1 ), 3 max i∈[n] L Gi 2 ) C F µ G .
Taking the total expectation and putting all together yields

E[∥D t v (z t , v t , x t )∥ 2 ] ≤ 2E[∥D v (z t , v t , x t )∥ 2 ] + 2L ′ v (E t z + E t x + E t v + E ′t z + E ′t x ) + 2L ′′ v (δ t z + δ t v ) . (123) In x we have similarly E[∥D t x (z t , v t , x t )∥ 2 ] ≤ 2E[∥D x (z t , v t , x t )∥ 2 ] + 2L ′ x (E t z + E t x + E t v + E ′t z + E ′t x ) + 2L ′′ x (δ t z + δ t v ) . (124)
E[∥D t z (z t , v t , x t )∥ 2 ] ≤ 2L 2 z δ t z + 2L ′ z S t , (125)
E[∥D t v (z t , v t , x t )∥ 2 ] ≤ 2(L 2 v + L ′′ v )(δ t z + δ t v ) + 2L ′ v S t (126)
and

E[∥D t x (z t , v t , x t )∥ 2 ] ≤ 2E[∥D x ∥ 2 ] + 2L ′ x S t + 2L ′′ x (δ t z + δ t v ) . (127)
We have the following lemma Lemma C.7.

If 4ρ 2 (L ′ z + L ′ v ) + 4γ 2 L ′ x ≤ Γ 2 and 4L ′′ x γ 2 ≤ ρ 2 (L 2 v + 4L ′′ v ), it holds S t+1 ≤ 1 - Γ 2 S t + β sz ρ 2 δ t z + β sv ρ 2 δ t v + P γ 2 E[∥D x ∥ 2 ]
for some L s , β sz , P > 0.

Proof. It holds following eq. (112) (and omitting the dependencies in (z t , v t , x t ) in the direction for simplicity)

S t+1 ≤ (1 -Γ) S t + E 2ρ 2 (∥D t z ∥ 2 + ∥D t v ∥ 2 ) + 2γ 2 ∥D t x ∥ 2 +2(m + n)[ρ 2 (∥D z ∥ 2 + ∥D v ∥ 2 ) + γ 2 ∥D x ∥ 2 ] .
Using the previous bounds (115), ( 123) and (124), we get

S t+1 ≤ 1 -Γ + 4ρ 2 (L ′ z + L ′ v ) + 4γ 2 L ′ x S t + (2(m + n) + 4)E[ρ 2 (∥D z ∥ 2 + ∥D v ∥ 2 ) + γ 2 ∥D x ∥ 2 ] + 4L ′′ v ρ 2 (δ t z + δ t v ) + 4L ′′ x γ 2 (δ t z + δ t v ) . Next, using 4ρ 2 (L ′ z + L ′ v ) + 4γ 2 L ′ x ≤ Γ 2 and letting P = (2(m + n) + 4) we get S t+1 ≤ 1 - Γ 2 S t +P E[ρ 2 (∥D z ∥ 2 +∥D v ∥ 2 )+γ 2 ∥D x ∥ 2 ]++4L ′′ v ρ 2 (δ t z +δ t v )+4L ′′ x γ 2 (δ t z +δ t v ) .
To finish, we use Lemma 3.4 to get

S t+1 ≤ 1 - Γ 2 S t + P [ρ 2 ((L 2 z + L 2 v )δ t z + L 2 v δ t v ) + (4L ′′ v ρ 2 + 4L ′′ x γ 2 )(δ t z + δ t v ) + γ 2 E[∥D x ∥ 2 ]] .
Then, using that

4L ′′ x γ 2 ≤ ρ 2 (L 2 v + 4L ′′ v ), we get the bound, letting L sz = L 2 z + L 2 v + 4L ′′ v and L sv = L 2 v + 4L ′′ v : S t+1 ≤ 1 - Γ 2 S t + β sz ρ 2 δ t z + β sv ρ 2 δ t v + P γ 2 E[∥D x ∥ 2 ] with β sz = 2P L sz , β sv = 2P L sv C.9.

Putting it all together

Recall that we denote g t = E[∥∇h(x t )∥ 2 ] and h t = E[h(x t )]. In the following lemma, we adapt Lemma 3.9 and Lemma 3.10 to the SABA algorithm.

Lemma C.8. If ρ ≤ min µ G 64L 2 z , β zx 2β zx , µ G 128(L 2 v + L ′′ v ) , β vz 8(L 2 v + L ′′ v ) , β vx 2β vx and γ ≤ min ρµ G 64β zx L ′′ x , L ′ z 2L ′ x β zx ρ, ρµ G 128β vx L ′′ x , ρβ vz 4L ′′ x β vx , L ′ v 2L ′ x β vx ρ, 1 4L h , L 2 x 2L h L ′′ x then it holds δ t+1 z ≤ 1 - ρµ G 8 δ t z + 2L ′′ x β zx γ 2 δ t v + 5L ′ z ρ 2 S t + 2β zx γ 2 ρ E[∥D x (z t , v t , x t )∥ 2 ] , (128)
δ t+1 v ≤ 1 - ρµ G 16 δ t v + 3β vz ρδ t z + 5L ′ v ρ 2 S t + 2β vx γ 2 ρ E[∥D x (z t , v t , x t )∥ 2 ] (129) and h t+1 ≤ h t - γ 2 g t - γ 4 E[∥D x (z t , v t , x t )∥ 2 ] + L 2 x γ(δ t z + δ t v ) + L h L ′ x γ 2 S t . (130)
Proof. We start from Lemma 3.9 and plug the bounds of Equations ( 125) and (126).

δ t+1 z ≤ 1 - ρµ G 4 + 4L 2 z ρ 2 + 4β zx L ′′ x γ 2 δ t z + 2L ′′ x β zx γ 2 δ t v (131) + (4L ′ z ρ 2 + 2L ′ x β zx γ 2 )S t + 2β zx γ 2 + β zx γ 2 ρ E[∥D x (z t , v t , x t )∥ 2 ] Since ρ ≤ µ G 64L 2 z and γ 2 ≤ ρµ G 64βzxL ′′ x , we have - ρµ G 4 + 4L 2 z ρ 2 + 4β zx L ′′ x γ 2 ≤ - ρµ G 8 . (132)
The condition

γ 2 ≤ L ′ z 2L ′ x βzx ρ 2 gives us 4L ′ z ρ 2 + 2L ′ x β zx γ 2 ≤ 5L ′ z ρ 2 . ( 133
)
With ρ ≤ β zx 2βzx , we get

2β zx γ 2 + β zx γ 2 ρ ≤ 2β zx γ 2 ρ . (134)
We can plug Equations ( 132), ( 133) and (134) into Equation (131) and we end up with

δ t+1 z ≤ 1 - ρµ G 8 δ t z + 2L ′′ x β zx γ 2 δ t v + 5L ′ z ρ 2 S t + 2β zx γ 2 ρ E[∥D x (z t , v t , x t )∥ 2 ] .
The proof for δ t v is quite similar. From Lemma 3.9, Equations ( 126) and (127).

δ t+1 v ≤ 1 - ρµ G 8 δ t v + β vz ρδ t z + 2ρ 2 V t v + β vx γ 2 V t x + β vx γ 2 ρ E[∥D x (z t , v t , x t )∥ 2 ] (135) ≤ 1 - ρµ G 8 + 4(L 2 v + L ′′ v )ρ 2 + 4L ′′ x β vx γ 2 δ t v + (4(L 2 v + L ′′ v )ρ 2 + 2L ′′ x β vx γ 2 + β vz ρ)δ t z + (136) + 4L ′ v ρ 2 + 2L ′ x β vx γ 2 S t + 2β vx γ 2 + β vx γ 2 ρ E[∥D x (z t , v t , x t )∥ 2 ] . Using ρ ≤ µ G 128(L 2 v +L ′′ v ) and γ 2 ≤ ρµ G 128L ′′
x βvx , we get

- ρµ G 8 + 4(L 2 v + L ′′ v )ρ 2 + 4L ′′ x β vx γ 2 ≤ - ρµ G 16 . ( 137
) With γ 2 ≤ ρβvz 4L ′′ x βvx and ρ ≤ βvz 8(L 2 v +L ′′ v ) , we have 4(L 2 v + L ′′ v )ρ 2 + 2L ′′ x β vx γ 2 + β vz ρ ≤ 3β vz ρ . (138)
The condition

γ 2 ≤ L ′ v 2L ′ x βvx ρ 2 yields 4L ′ v ρ 2 + 2L ′ x β vx γ 2 ≤ 5L ′ v ρ 2 . (139)
With ρ ≤ β vx 2βvx we get

2β vx γ 2 + β zx γ 2 ρ ≤ 2β vx γ 2 ρ . (140)
As a consequence of Equations ( 135), (137), (138), ( 139) and (140), we have

δ t+1 v ≤ 1 - ρµ G 16 δ t v + 3β vz ρδ t z + 5L ′ v ρ 2 S t + 2β vx γ 2 ρ E[∥D x (z t , v t , x t )∥ 2 ] .
For the inequality on h t , we start from Equations ( 88) and (127)

h t+1 ≤ h t - γ 2 g t - γ 2 -L h γ 2 E[∥D x (z t , v t , x t )∥ 2 ] (141) + L 2 x 2 γ + L h L ′′ x γ 2 (δ t z + δ t v ) + L h L ′ x γ 2 S t . Assuming γ ≤ min 1 4L h , L 2 x 2L h L ′′ x leads h t+1 ≤ h t - γ 2 g t - γ 4 E[∥D x (z t , v t , x t )∥ 2 ] + L 2 x γ(δ t z + δ t v ) + L h L ′ x γ 2 S t . (142)
We are now ready to prove Theorem 3.

Proof. We consider the Lyapunov function

L t = h t + ϕ s S t + ϕ z δ t z + ϕ v δ t v (143)
for some constants ϕ s , ϕ z and ϕ v .

We have

L t+1 -L t ≤ - γ 2 g t - γ 4 -2ϕ z β zx γ 2 ρ -2ϕ v β vx γ 2 ρ -ϕ s P γ 2 E[∥D x (z t , v t , x t )∥ 2 ] -ϕ z µ G 8 ρ -L 2 x γ -8ϕ v β vz ρ -ϕ s β sz ρ 2 δ t z -ϕ v µ G 16 ρ -L 2 x γ -2ϕ z L ′′ x γ 2 -ϕ s β sv ρ 2 δ t v -ϕ s Γ 2 -5ϕ z L ′ z ρ 2 -5ϕ v L ′ v ρ 2 -L h L ′ x γ 2 S t .
To get a decrease, ϕ z , ϕ v and ϕ s , ρ and γ must be such that:

2ϕ z β zx γ 2 ρ + 2ϕ v β vx γ 2 ρ + ϕ s P γ 2 ≤ γ 4 L 2 x γ + 8ϕ v β vz ρ + ϕ s β sz ρ 2 ≤ ϕ z µ G 8 ρ L 2 x γ + 8ϕ z L ′′ x γ 2 + ϕ s β sv ρ 2 ≤ ϕ v µ G 16 ρ 5ϕ z L ′ z ρ 2 + 5ϕ v L ′ v ρ 2 + L h L ′ x γ 2 ≤ ϕ s Γ 2 .
In order to take into account the scaling of the quantities with respect to N = n + m, we take

ρ = ρ ′ N nρ , γ = γ ′ N nγ , ϕ z = ϕ ′ z N nz , ϕ v = ϕ ′ v N nv and ϕ s = ϕ ′ s N ns . Since Γ = O(N -1 ), P = O(N ), β sz = O(N ) and β sv = O(N ), we also define Γ ′ = ΓN , P ′ = P N -1 , β ′ sz = β sz N -1 and β ′ sv N -1
. Now, the previous Equations read (after slight simplifications):

(2ϕ ′ z β zx + 2ϕ ′ v β vx ) γ ′ ρ ′ N nz+nγ -nρ + ϕ ′ s P ′ γ ′ N ns+nγ +1 ≤ 1 4 L 2 x γ ′ N nγ + 8ϕ ′ v β vz ρ ′ N nv+nρ + ϕ ′ s β ′ sz (ρ ′ ) 2 N 2nρ+ns+1 ≤ ϕ ′ z µ G 8 ρ ′ N nz+nρ L 2 x γ ′ N nγ + 8ϕ ′ z L ′′ x (γ ′ ) 2 N 2nγ +nz + ϕ ′ s β ′ sv (ρ ′ ) 2 N ns+2nρ+1 ≤ ϕ ′ v µ G 16 ρ ′ N nv+nρ 5ϕ ′ z L ′ z (ρ ′ ) 2 N nz+2nρ + 5ϕ ′ v L ′ v (ρ ′ ) 2 N 2nρ+nv + L h L ′ x (γ ′ ) 2 N nγ ≤ ϕ s Γ ′ 2 N ns-1 .
In order to ensure that the exponents on N are lower in the left-hand-side than those on the righthand-side, we take n z = n v = 0, n ρ = n γ = -2 3 and n s = -1 3 . The Equations become

(2ϕ ′ z β zx + 2ϕ ′ v β vx ) γ ′ ρ ′ + ϕ ′ s P ′ γ ′ ≤ 1 4 L 2 x γ ′ N -2/3 + 8ϕ ′ v β vz ρ ′ N -2/3 + ϕ ′ s β ′ sz (ρ ′ ) 2 N -2/3 ≤ ϕ ′ z µ G 8 ρ ′ N -2/3 L 2 x γ ′ N -2/3 + 8ϕ ′ z L ′′ x (γ ′ ) 2 N -4/3 + ϕ ′ s β ′ sv (ρ ′ ) 2 N -2/3 ≤ ϕ ′ v µ G 16 ρ ′ N -2/3 5ϕ ′ z L ′ z (ρ ′ ) 2 N -4/3 + 5ϕ ′ v L ′ v (ρ ′ ) 2 N -4/3 + L h L ′ x (γ ′ ) 2 N -4/3 ≤ ϕ ′ s Γ ′ 2 N -4/3 .
We can replace the penultimate equation by the stronger

L 2 x γ ′ N -2/3 + 8ϕ ′ z L ′′ x (γ ′ ) 2 N -2/3 + ϕ ′ s β ′ sv (ρ ′ ) 2 N -2/3 ≤ ϕ ′ v µ G 16 ρ ′ N -2/3
so that we can simplify all the equations by dropping the dependencies in N :

(2ϕ ′ z β zx + 2ϕ ′ v β vx ) γ ′ ρ ′ + ϕ ′ s P ′ γ ′ ≤ 1 4 L 2 x γ ′ + 8ϕ ′ v β vz ρ ′ + ϕ ′ s β ′ sz (ρ ′ ) 2 ≤ ϕ ′ z µ G 8 ρ ′ L 2 x γ ′ + 8ϕ ′ z L ′′ x (γ ′ ) 2 + ϕ ′ s β ′ sv (ρ ′ ) 2 ≤ ϕ ′ v µ G 16 ρ ′ 5ϕ ′ z L ′ z (ρ ′ ) 2 + 5ϕ ′ v L ′ v (ρ ′ ) 2 + L h L ′ x (γ ′ ) 2 ≤ ϕ ′ s Γ ′ 2 .
Let us take ϕ

′ s = 1, ϕ ′ z = ϕ ′′ z ρ ′ γ ′ and ϕ ′ v = ϕ ′′ v ρ ′ γ ′ with ϕ ′′ z = 1 32β zx and ϕ ′′ v = min 1 32β vx , ϕ ′′ z µ G 128βvz .
The equations become

P ′ γ ′ ≤ 1 8 L 2 x γ ′ + β ′ sz (ρ ′ ) 2 ≤ ϕ ′′ z µ G 16 (ρ ′ ) 2 γ ′ L 2 x γ ′ + 8ϕ ′′ z L ′′ x γ ′ ρ ′ + β ′ sv (ρ ′ ) 2 ≤ ϕ ′′ v µ G 16 (ρ ′ ) 2 γ ′ 5ϕ ′′ z L ′ z (ρ ′ ) 3 γ ′ + 5ϕ ′′ v L ′ v (ρ ′ ) 3 γ ′ + L h L ′ x (γ ′ ) 2 ≤ Γ ′ 2 .
The condition γ ′ ≤ 1 8P ′ ensures that the first equation is verified. With γ ′ ≤ min

ϕ ′′ z µ G 32L 2 x ρ ′ , ϕ ′′ z µ G 32β ′ sz ,
the second equations is verified. With γ ′ ≤ min

ϕ ′′ v µ G 48L 2 x ρ ′ , ϕ ′′ v µ G 48β ′ sv , ϕ ′′ v µ G 384ϕ ′′ z L ′′ x ρ ′ , the third is verified. With γ ′ ≤ Γ ′ 6L h L ′ x
, the last can be simplified:

(5ϕ ′′ z L ′ z + 5ϕ ′′ v L ′ v )(ρ ′ ) 3 ≤ Γ ′ 3 γ ′ .
Let us write γ ′ = ξρ ′ . If we want that equation does no contradict the previous upper bound on γ ′ involving ρ ′ and the conditions of Lemma C.8, that is

γ ′ ≤ min ϕ ′′ z µ G 32L 2 x , ϕ ′′ v µ G 48L 2 x , L ′ z 2L ′ x β zx , L ′ v 2L ′ x β vx K1 ρ ′ γ ′ ≤ min µ G 64β zx L ′′ x , µ G 128β vx L ′′ x , β vz 4L ′′ x β vx K2 ρ ′ γ ′ ≤ ϕ ′′ v µ G 384ϕ ′′ z L ′′ x K3 1 √ ρ ′ γ ′ ≤ min 1 4L h , L 2 x 2L h L ′′ x , Γ ′ 6L h L ′ x , 1 8P ′ , ϕ ′′ z µ G 32β ′ sz , ϕ ′′ v µ G 48β ′ sv K4 γ ′ ≥ 15(ϕ ′′ z L ′ z + ϕ ′′ v L ′ v ) Γ ′ K5 ρ 3 ξ must verify ξ ≤ K 1 ξ ≤ K 2 (ρ ′ ) -1 2 ξ ≤ K 3 (ρ ′ ) -3 2 ξ ≤ K 4 (ρ ′ ) -1 ξ ≥ K 5 (ρ ′ ) 2 which is possible if ρ ′ satisfies ρ ′ ≤ min K 1 K 5 , K 2 K 5 -3 2 , K 3 K 5 -5 2 , K 4 K 5 -2
.

Proof. For simplicity, we assume that h * = 0 and so for any x ∈ R d the PL inequality reads:

1 2 ∥∇h(x)∥ 2 ≥ µ h h(x) . (146)
Then, eq. ( 130) gives

h t+1 ≤ 1 - γµ h 2 h t - γ 4 E[∥D x (z t , v t , x t )∥ 2 ] + γL 2 x (δ t z + δ t v ) + L h L ′ x γ 2 S t .
We take L t the Lyapunov function given in Equation (143). We find

L t+1 -L t ≤ -γµ h h t - γ 4 -2ϕ z β zx γ 2 ρ -2ϕ v β vx γ 2 ρ -ϕ s P γ 2 E[∥D x (z t , v t , x t )∥ 2 ] -ϕ z µ G 8 ρ -L 2 x γ -8ϕ v β vz ρ -ϕ s β sz ρ 2 δ t z -ϕ v µ G 16 ρ -L 2 x γ -2ϕ z L ′′ x γ 2 -ϕ s β sv ρ 2 δ t v -ϕ s Γ 2 -5ϕ z L ′ z ρ 2 -5ϕ v L ′ v ρ 2 -L h L ′ x γ 2 S t .
We now try to find linear convergence, hence we add to this cL t to get

L t+1 -(1 -c)L t ≤ -(γµ h -c)h t - γ 4 -2ϕ z β zx γ 2 ρ -2ϕ v β vx γ 2 ρ -ϕ s P γ 2 -c E[∥D x (z t , v t , x t )∥ 2 ] -ϕ z µ G 8 ρ -L 2 x γ -8ϕ v β vz ρ -ϕ s β sz ρ 2 -cϕ z δ t z -ϕ v µ G 16 ρ -L 2 x γ -2ϕ z L ′′ x γ 2 -ϕ s β sv ρ 2 -cϕ v δ t v -ϕ s Γ 2 -5ϕ z L ′ z ρ 2 -5ϕ v L ′ v ρ 2 -L h L ′ x γ 2 -cϕ S S t .
Hence, the set of inequations for decrease becomes

c ≤ γµ h 2ϕ z β zx γ 2 ρ + 2ϕ v β vx γ 2 ρ + ϕ s P γ 2 + c ≤ γ 4 L 2 x γ + 8ϕ v β vz ρ + ϕ s β sz ρ 2 + ϕ z c ≤ ϕ z µ G 8 ρ L 2 x γ + 8ϕ z L ′′ x γ 2 + ϕ s β sv ρ 2 + ϕ v c ≤ ϕ v µ G 16 ρ 5ϕ z L ′ z ρ 2 + 5ϕ v L ′ v ρ 2 + L h L ′ x γ 2 + ϕ s c ≤ ϕ s Γ 2 .
We see that it is more convenient to write c = γc ′ . As previously, we write γ = γ ′ N nγ , ρ = ρ ′ N nρ ,

ϕ z = ϕ ′ z N nz , ϕ v = ϕ ′ v N nv , ϕ s = ϕ ′ s N ns , P = P ′ N , Γ = Γ ′ N -1 , β sx = β ′
sx N and β sv = β ′ sv N . The equations read:

c ′ ≤ µ h 2ϕ ′ z β zx γ ′ ρ ′ N nz+nγ -nρ + 2ϕ ′ v β vx γ ′ ρ ′ N nv+nγ -nρ + ϕ ′ s P ′ γ ′ N ns+1+nγ + c ′ ≤ 1 4 L 2 x γ ′ N nγ + 8ϕ ′ v β vz ρ ′ N nv+nρ + ϕ ′ s β ′ sz (ρ ′ ) 2 N ns+2nρ+1 + ϕ ′ z c ′ γ ′ N nz+nγ ≤ ϕ ′ z µ G 8 ρ ′ N nρ+nz L 2 x γ ′ N nγ + 8ϕ ′ z L ′′ x (γ ′ ) 2 N nz+2nγ + ϕ ′ s β ′ sv (ρ ′ ) 2 N ns+1+2nρ + ϕ ′ v c ′ γ ′ N nv+nγ ≤ ϕ ′ v µ G 16 ρ ′ N nv+nρ 5ϕ ′ z L ′ z (ρ ′ ) 2 N nz+2nρ + 5ϕ ′ v L ′ v (ρ ′ ) 2 N nv+2nρ + L h L ′ x (γ ′ ) 2 N 2nγ + ϕ ′ s c ′ γ ′ N ns+nγ ≤ ϕ ′ s Γ ′ 2 N ns-1 .
In order to ensure that the exponents on N are lower in the left-hand-side than those on the righthand-side, we take n z = n v = 0, n ρ = -2 3 , n γ = -1 and n s = -1 3 . The Equations become

c ′ ≤ µ h 2ϕ ′ z β zx γ ′ ρ ′ N -1 3 + 2ϕ ′ v β vx γ ′ ρ ′ N -1 3 + ϕ ′ s P ′ γ ′ N -1 3 + c ′ ≤ 1 4 L 2 x γ ′ N -1 + 8ϕ ′ v β vz ρ ′ N -2 3 + ϕ ′ s β ′ sz (ρ ′ ) 2 N -2 3 + ϕ ′ z c ′ γ ′ N -1 ≤ ϕ ′ z µ G 8 ρ ′ N -2 3 L 2 x γ ′ N -1 + 8ϕ ′ z L ′′ x (γ ′ ) 2 N -2 + ϕ ′ s β ′ sv (ρ ′ ) 2 N -2 3 + ϕ ′ v c ′ γ ′ N -1 ≤ ϕ ′ v µ G 16 ρ ′ N -2 3 5ϕ ′ z L ′ z (ρ ′ ) 2 N -4 3 + 5ϕ ′ v L ′ v (ρ ′ ) 2 N -2 + L h L ′ x (γ ′ ) 2 N -2 + ϕ ′ s c ′ γ ′ N -4 3 ≤ ϕ ′ s Γ ′ 2 N -4 3 .
Now we have to find ρ ′ , γ ′ , ϕ ′ z , ϕ ′ v and ϕ ′ s that verifies the following conditions (which are a bit stronger than thoose in the previous Equations):

c ′ ≤ µ h 2ϕ ′ z β zx γ ′ ρ ′ + 2ϕ ′ v β vx γ ′ ρ ′ + ϕ ′ s P ′ γ ′ + c ′ ≤ 1 4 L 2 x γ ′ + 8ϕ ′ v β vz ρ ′ + ϕ ′ s β ′ sz (ρ ′ ) 2 + ϕ ′ z c ′ γ ′ ≤ ϕ ′ z µ G 8 ρ ′ L 2 x γ ′ + 8ϕ ′ z L ′′ x (γ ′ ) 2 + ϕ ′ s β ′ sv (ρ ′ ) 2 + ϕ ′ v c ′ γ ′ ≤ ϕ ′ v µ G 16 ρ ′ 5ϕ ′ z L ′ z (ρ ′ ) 2 + 5ϕ ′ v L ′ v (ρ ′ ) 2 + L h L ′ x (γ ′ ) 2 + ϕ ′ s c ′ γ ′ ≤ ϕ ′ s Γ ′ 2 .
As previously, we take ϕ ′ s = 1 and we denote ϕ ′ z = ϕ ′′

c ′ ≤ µ h P ′ γ ′ + c ′ ≤ 1 8 L 2 x (γ ′ ) 2 + β ′ sz (ρ ′ ) 2 γ ′ + ϕ ′′ z c ′ ρ ′ γ ′ ≤ ϕ ′′ z µ G 16 (ρ ′ ) 2 L 2 x (γ ′ ) 2 + 8ϕ ′′ z L ′′ x ρ ′ (γ ′ ) 2 + β ′ sv (ρ ′ ) 2 γ ′ + ϕ ′′ v c ′ ρ ′ γ ′ ≤ ϕ ′′ v µ G 16 (ρ ′ ) 2 5ϕ ′′ z L ′ z (ρ ′ ) 3 + 5ϕ ′′ v L ′ v (ρ ′ ) 3 + L h L ′ x (γ ′ ) 3 + c ′ (γ ′ ) 2 ≤ Γ ′ 2 γ ′ .
Since c ′ ≤ 1 16 and γ ′ ≤ 1 16P ′ , the second equation is verified. With γ ′ ≤ min

ϕ ′′ z µ G 48L 2 x ρ ′ , ϕ ′′ z µ G 48βsv
and c ′ ≤ µ G ρ ′ 48γ ′ the third is verified. The conditions γ ′ ≤ min

ϕ ′′ v µ G 64L 2 x ρ ′ , ϕ ′′ v µ G 512ϕ ′′ z L ′′ x ρ ′ , ϕ ′′ v µ G 64β ′ sv and c ′ ≤ µ G ρ ′ 64γ ′ ensure that the forth is verified. With γ ′ ≤ Γ ′ 8L h L ′
x and c ′ ≤ Γ ′ 8γ ′ , the fifth is simplified in

5ϕ ′′ z L ′ z (ρ ′ ) 3 + 5ϕ ′′ v L ′ v (ρ ′ ) 3 ≤ Γ ′ 4 γ ′ .
The main difference with Lemma 3.9 is that we have O( γ 2 ρ ) in factor of V t x instead of O(γ 2 ). As a consequence, we need that the ratio γ ρ goes to zero to get convergence, as in [START_REF] Hong | A Two-Timescale Framework for Bilevel Optimization: Complexity Analysis and Application to Actor-Critic[END_REF]. This prevent us in getting rates that match rates of single level algorithms.

Hence, for SOBA, we have to choose γ = O(T - 3 5 ) and ρ = O(T -2 5 ) and we end up with a convergence rate in O(T - 2 5 ). For SABA, we get a O((n + m)ϵ -1 ) sample complexity, which is actually the sample complexity of SOBA used with full batch estimated directions.

Figure 1 :

 1 Figure 1: Convergence curves of the two proposed methods on a toy problem. SABA is a stochastic method that achieves fast convergence on the value function.

Figure 2 :

 2 Figure 2: Comparison of SOBA and SABA with other stochastic bilevel optimization methods. For each algorithm, we plot the median performance over 10 runs. In both experiments, SABA achieves the best performance. The dashed lines are for one loop competitor methods, the dotted lines are for two loops methods and the solid lines are the proposed methods. Left: hyperparameter selection for ℓ 2 penalized logistic regression on IJCNN1 dataset , Right: data hyper-cleaning on MNIST with p = 0.5 corruption rate.

Checklist 1 .

 1 For all authors... (a) Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope? [Yes] See Section 3 and Section 4 (b) Did you describe the limitations of your work? [Yes] (c) Did you discuss any potential negative societal impacts of your work? [N/A] (d) Have you read the ethics review guidelines and ensured that your paper conforms to them? [Yes] 2. If you are including theoretical results... (a) Did you state the full set of assumptions of all theoretical results? [Yes] See Section 3.1 (b) Did you include complete proofs of all theoretical results? [Yes] See Appendix C 3. If you ran experiments... (a) Did you include the code, data, and instructions needed to reproduce the main experimental results (either in the supplemental material or as a URL)? [Yes] See supplementary material (b) Did you specify all the training details (e.g., data splits, hyperparameters, how they were chosen)? [Yes] See Appendix B (c) Did you report error bars (e.g., with respect to the random seed after running experiments multiple times)? [Yes] (d) Did you include the total amount of compute and the type of resources used (e.g., type of GPUs, internal cluster, or cloud provider)? [Yes] See Appendix B 4. If you are using existing assets (e.g., code, data, models) or curating/releasing new assets... (a) If your work uses existing assets, did you cite the creators? [Yes] (b) Did you mention the license of the assets? [Yes] (c) Did you include any new assets either in the supplemental material or as a URL? [Yes] (d) Did you discuss whether and how consent was obtained from people whose data you're using/curating? [N/A] (e) Did you discuss whether the data you are using/curating contains personally identifiable information or offensive content? [N/A] 5. If you used crowdsourcing or conducted research with human subjects... (a) Did you include the full text of instructions given to participants and screenshots, if applicable? [N/A] (b) Did you describe any potential participant risks, with links to Institutional Review Board (IRB) approvals, if applicable? [N/A] (c) Did you include the estimated hourly wage paid to participants and the total amount spent on participant compensation? [N/A]

10 .

 10 Then we use 750 pairs (x train i , y train i ) 1≤i≤750 as training samples and the remaining pairs (x val i , y val i ) 1≤i≤250 as validation samples. Finally, we solve (1) with

Figure B. 1 :

 1 Figure B.1: Comparison of SOBA and SABA with other stochastic bilevel optimization methods in a problem of hyperparameter selection for ℓ 2 penalized logistic regression on IJCNN1 dataset.For each algorithm, we plot the median performance over 10 runs. In both plots, SABA achieves the best performance. The dashed lines are for one loop competitor methods, the dotted lines are for two loops methods and the solid lines are the proposed methods. Left: performance in running time, Right: performance in number of gradient/Hessian-vector products sampled.

Figure B. 2 :

 2 Figure B.2: Datacleaning experiment, with different corruption probability (higher means that more data are contamined). Top: Performance with respect to the number of gradient/Hessian-vector product sampled, Bottom: Performance with respect to running time

Figure B. 3 :

 3 Figure B.3: Comparison of SOBA and SABA with other stochastic bilevel optimization methods in a problem of hyperparameter selection for ℓ 2 penalized multical logistic regression on covtype dataset. For each algorithm, we plot the median performance over 10 runs. The dashed lines are for one loop competitor methods, the dotted lines are for two loops methods and the solid lines are the proposed methods. Left: performance in running time, Right: performance in number of gradient/Hessian-vector products sampled.

  Lemmas on the regularity of z * and v * We start by showing the Lipschitz continuity of z * and v * .

ϵ -1 ) which is analogous of the one of SAGA for

https://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/binary.html

http://yann.lecun.com/exdb/mnist/

All the code is provided in supplementary materials

https://scikit-learn.org/stable/modules/generated/sklearn.datasets.fetch_ covtype.html

Acknowledgments and Disclosure of Funding

We thank Othmane Sebbouh, Zaccharie Ramzi and Benoît Malézieux for their precious comments. The authors acknowledge the support of the ANER RAGA BFC. SV acknowledges the support of the ANR GraVa ANR-18-CE40-0005. This research was supported by DATAIA convergence institute as part of the " Programme d'Investissement d'Avenir ", (ANR-17-CONV-0003) operated by Inria.

Lemma C.3. The function v * : R d → R p is differentiable and its differential is defined for any x, ϵ ∈ R d by: dv * (x).ϵ = [∇ 2 1 G(z * (x), x)] -1 [∇ 2 11 F (z * (x), x)dz * (x).ϵ + ∇ 2 12 F (z * (x), x).ϵ] (31)

where for any z, α ∈ R p and x ∈ R d , (∇ 3 111 G(z, x)|α) ∈ R p×p is defined by

and for any β ∈ R d , (∇ 3 112 G(z, x)|β) ∈ R p×p is defined by

Using the differentiability of ∇ 2 11 G, ∇ 1 F and of the matrix inversion, we have

that proves [START_REF] Liu | On the limited memory BFGS method for large scale optimization[END_REF]. Now, let x, y, ϵ ∈ R d with ∥ϵ∥ = 1. Let us denote

We have

and therefore, summing and telescoping yields

T .

Since with respect to N we have

T .

C.10 Proof of Theorem 4

We are now going to prove Theorem 4 that we recall here:

Theorem 4 (Convergence of SABA, PL case). Assume that h satisfies the PL inequality and that Assumptions 3.1 to 3.3 and 3.7 to 3.8 hold. We suppose ρ = ρ ′ N -2 3 and γ = ξρ ′ N -1 , where ρ ′ and ξ depend only on F and G and are specified in appendix. Let x t the iterates of SABA and c ′ ≜ min µ h , 1 16P ′ with P ′ specified in the appendix. Then,

where C 0 is a constant specified in appendix that depends on the initialization of z, v, x and memory.

Here, we have

) . where

As in the proof of Theorem 3, let us denote γ ′ = ξρ ′ . To verify this equation and the previous bounds on γ ′ and c ′ , we need

So, ξ, ρ ′ and c ′ must verify

D Convergence rates with weaker regularity assumptions

To get our rates, we need stronger assumptions than in the stochastic bilevel optimization literature [START_REF] Ghadimi | Approximation Methods for Bilevel Programming[END_REF][START_REF] Hong | A Two-Timescale Framework for Bilevel Optimization: Complexity Analysis and Application to Actor-Critic[END_REF][START_REF] Ji | Bilevel optimization: Convergence analysis and enhanced design[END_REF][START_REF] Arbel | Amortized Implicit Differentiation for Stochastic Bilevel Optimization[END_REF]. In this section, we shortly present the convergence rates we can expect if we replace Assumptions 3. . With these assumptions, we are not ensured that v * is smooth, and so the descent lemmas take the form of Lemma D.3. Lemma D.3. Assume that ρ ≤ 2 µ G . We have:

where L * is the maximum between the Lipschitz constants of z * and v * (see Lemma C.1) and

Proof. Inequality for δ z .

Instead of expanding the square as done in the proof of Lemma 3.9 in Equation ( 45), we use Young's inequality for some a > 0

Treating

In order to keep a decrease in δ z , we might want to use a = 1 2 ρµ G , which gives the bound

with β zx = 4

L 2 * µ G . Indeed, this gives (1 + 1 2 ρµ G )(1 -ρµ G ) ≤ 1 -1 2 ρµ G . We have a ≤ 1 since ρ ≤ 2 µ G , so (1 + a)ρ 2 ≤ 2ρ 2 . Finally, we also have 1 + a -1 ≤ 2a -1 = 4 ρµ G . Inequality for δ v . As for δ z , the difference with the proof of Lemma 3.9 is that we use we use Young's inequality for some b > 0 to get

(150) The remaining part of the proof is similar to the proof of Lemma 3.9.