domain: math info

text:
Mathieu Dagréou

Thomas Moreau

Samuel Vaiter

Pierre Ablin

A Lower Bound and a Near-Optimal Algorithm for Bilevel Empirical Risk Minimization

Bilevel optimization problems, which are problems where two optimization problems are nested, have more and more applications in machine learning. In many practical cases, the upper and the lower objectives correspond to empirical risk minimization problems and therefore have a sum structure. In this context, we propose a bilevel extension of the celebrated SARAH algorithm. We demonstrate that the algorithm requires O((n + m) 1 2 ε -1 ) oracle calls to achieve ε-stationarity with n + m the total number of samples, which improves over all previous bilevel algorithms. Moreover, we provide a lower bound on the number of oracle calls required to get an approximate stationary point of the objective function of the bilevel problem. This lower bound is attained by our algorithm, making it optimal in terms of sample complexity.

Introduction

In the last few years, bilevel optimization has become an essential tool for the machine learning community thanks to its numerous applications. Among them, we can cite hyperparameter selection [START_REF] Bengio | Gradient-Based Optimization of Hyperparameters[END_REF]Pedregosa, 2016;Franceschi et al., 2017;Lorraine et al., 2020), implicit deep learning [START_REF] Shaojie Bai | Deep Equilibrium Models[END_REF], neural architecture search (Liu et al., 2019;Zhang et al., 2021), data augmentation (Li et al., 2020;Rommel et al., 2022) or meta-learning (Franceschi et al., 2018;Rajeswaran et al., 2019). Bilevel optimization consists in minimizing a function under the constraint that (1)

The function F is called the outer function and the function G is the inner function. Likewise, we refer to x as the outer variable and z as the inner variable. The function h is the value function and it can be minimized using gradient descent. To compute its gradient, we use implicit differentiation which yields

∇h(x) = ∇ 2 F (z * (x), x) + ∇ 2 21 G(z * (x), x)v * (x) (2)
where v * (x) is the solution of a linear system v * (x) = -∇ 2 11 G(z * (x), x) -1 ∇ 1 F (z * (x), x) .

(3)

When we have exact access to z * (x), solving (1) boils down to a smooth nonconvex optimization problem which can be solved using solvers for single-level problems. However, computing exactly z * (x) and v * (x) is often too costly, and implicit differentiation-based algorithms rely on approximations of z * (x) and v * (x) rather than their exact value. Depending on the precision of the different approximations, we are not ensured that the approximate gradient used is a descent direction. Results by Pedregosa (2016) characterized the approximation quality for z * (x) and v * (x) required to ensure convergence, opening the door to various algorithms to solve bilevel optimization problems (Lorraine et al., 2020;Ramzi et al., 2022).

In many applications of interest, the functions F and G correspond to Empirical Risk Minimization (ERM), and as a consequence have a finite sum structure

F (z, x) = 1 m m j=1 F j (z, x), G(z, x) = 1 n n i=1 G i (z, x) .
For instance, in hyperparameter selection, F is the validation loss which is an average on the validation set and G is the training loss which is an average on the training set. In single-level optimization, the finite sum structure has been widely leveraged to produce fast first-order algorithms that provably converge faster than gradient descent. Among them, we can cite stochastic methods such as stochastic gradient descent (SGD) (Robbins and Monro, 1951;[START_REF] Bottou | Large-Scale Machine Learning with Stochastic Gradient Descent[END_REF] and its variance-reduced variants such as SAGA (Defazio et al., 2014), STORM (Cutkosky and Orabona, 2019) or SPIDER/SARAH (Fang et al., 2018;Nguyen et al., 2017) that use only a handful of samples at a time to make progress. To get faster methods than full-batch approaches, it is natural to extend these methods to the bilevel setting. The main obstacle comes from the difficulty of obtaining stochastic approximations of ∇h(x) because of its structure (2) which involves a Hessian inversion. Several strategies have been proposed to overcome this obstacle, and some works demonstrate that stochastic implicit differentiation-based algorithms for solving (1) have the same complexity as single-level analogous algorithms. For instance, ALSET from Chen et al. ( 2021) and SOBA from Dagréou et al. (2022) have the same convergence rate as SGD for nonconvex singlelevel problems (Ghadimi and Lan, 2013;[START_REF] Bottou | Optimization methods for large-scale machine learning[END_REF]. Also, Dagréou et al. (2022) show that SABA, an adaptation of SAGA (Defazio et al., 2014), has an analogous sample complexity to its single-level counterparts for nonconvex problems (Reddi et al., 2016).

Yet, in classical single-level optimization, it is known that neither of these algorithms is optimal: the SARAH algorithm (Nguyen et al., 2017) achieves a better sample complexity of O(m 1 2 ε -1 ) with m the number of samples. Furthermore, this algorithm is near-optimal (i.e. optimal up to constant factors) because the lower bound for single-level nonconvex optimization is Ω(m 1 2 ε -1 ) as proved by Zhou and Gu (2019). It is natural to ask if we can extend these results to bilevel optimization.

Contributions

In Section 2, we introduce SRBA, an adaptation of the SARAH algorithm to the bilevel setting. We then demonstrate in Section 3 that, similarly to the single-level setting, O (n + m) 1 2 ε -1 ∨ (n + m) oracle calls are sufficient to reach an ε-stationary point. As shown in Table 1, it achieves the best-known complexity in the regime n + m ≲ O(ε -2 ). In Section 4, we analyze the lower bounds for such problems. We show that we need at least Ω(m 1 2 ε -1 ) oracle calls to reach an ε-stationary point (see Definition 3.1), hereby matching the previous upper-bound in the case where n ≍ m and ε ≤ m -1 2 . SRBA is thus near-optimal in that regime. Even though our main contribution is theoretical, we illustrate the numerical performances of the algorithm in Section 5. Related work There are several strategies to solve (1) with a stochastic method. The first one is the value-function-based method which consists in recasting Problem 1 as a single-level constrained optimization problem as done with F 2 SA (Kwon et al., 2023) or BOME (Ye et al., 2022). The second way is to use first-order methods on h with approximate gradients. The approximate gradient of h can be estimated using two approaches: iterative differentiation (ITD) and approximate implicit differentiation (AID). On the one hand, in ITD algorithms, the Jacobian of z * is estimated by differentiating the different steps used to compute an approximation of z * . On the other hand, AID algorithms leverage the implicit gradient given by (2) replacing z * and v * by some approximations z and v. In the class of ITD algorithms, Maclaurin et al. (2015) propose to approximate the Jacobian of the solution of the inner problem by differentiating through the iterations of SGD with momentum. The complexity of the hypergradient computation in ITD solvers is studied in Franceschi et al. (2017); Grazzi et al. (2020); [START_REF] Ablin | Super-efficiency of automatic differentiation for functions defined as a minimum[END_REF]. For AID algorithms, Ghadimi and Wang (2018); Chen et al. (2021);Ji et al. (2021) propose to perform several SGD steps in the inner problem and then use Neumann approximations to approximate v * (x) defined in (3). A method consisting of alternating steps in the inner and outer variables was proposed in Hong et al. (2023). These methods can be improved by using a warm start strategy for the inner problem (Ji et al., 2021;Chen et al., 2021) and for the linear system [START_REF] Arbel | Amortized Implicit Differentiation for Stochastic Bilevel Optimization[END_REF]. Some works adapt variance reduction methods to like STORM (Cutkosky and Orabona, 2019;Khanduri et al., 2021;Yang et al., 2021) or SAGA (Defazio et al., 2014;Dagréou et al., 2022). We take a similar approach and extend the SARAH variance reduction method to the bilevel setting. Recent works propose to approximate the Jacobian of z * by stochastic finite difference (Sow et al., 2022) or to use Bregman divergence-based methods (Huang et al., 2022).

In single-level optimization, the problem of finding complexity lower bound has been widely studied since the seminal work of Nemirovsky and Yudin (1983). On the one hand, [START_REF] Agarwal | A Lower Bound for the Optimization of Finite Sums[END_REF] provided a lower bound to minimize strongly convex and smooth finite sum with deterministic algorithms that have access to individual gradients. These results were extended to randomized algorithms for (strongly) convex finite sum objective by Woodworth and Srebro (2016). On the other hand, [START_REF] Carmon | Lower bounds for finding stationary points I[END_REF] provided a lower bound for minimizing nonconvex functions with deterministic and randomized algorithms. The nonconvex finite sum case is treated in Fang et al. (2018); Zhou and Gu (2019). In the bilevel case, Ji and Liang (2023) showed a lower bound for deterministic algorithms. However, this result is restricted to the case where the value function h is convex or strongly convex, which andC 3,3 L Table 1: Comparison between the sample complexities and the Assumptions of some stochastic bilevel solvers. It corresponds to the number of calls to gradient, Hessian-vector products, and Jacobian-vector product sufficient to get an ε-stationary point. The tilde on the Õ hide a factor log(ε -1 ). "SC" means "strongly convex". C p,k L means p-times differentiable with Lipschitz kth order derivatives for k ≤ p.

F 2 SA (Kwon et al., 2023) O(ε -7 2 ) General expectation C 2,2 L SC and C 2,2 L SRBA O((n + m) 1 2 ε -1 ) Finite sum C 2,2 L SC
is not the case with most ML-related bilevel problems. Our results are instead in a nonconvex setting.

Notation The quantity A • refers to A z , A v , or A x , depending on the context. If f : R p ×R d → R is a twice differentiable function, we denote ∇ i f (z, x) its gradient w.r.t. its i th variable. Its Hessian with respect to z is denoted ∇ 2 11 f (z, x) ∈ R p×p and its cross derivative matrix

∂ 2 f ∂zi∂xj i∈[p] j∈[d] is denoted ∇ 2 12 f (z, x) ∈ R p×d .
We denote Π C the projection on a closed convex set C.

SRBA: a Near-Optimal Algorithm for Bilevel Empirical Risk Minimization

In this section, we introduce SRBA (Stochastic Recursive Bilevel Algorithm), a novel algorithm for bilevel empirical risk minimization which is provably nearoptimal for this problem. This algorithm is inspired by the algorithms SPIDER (Fang et al., 2018) and SARAH (Nguyen et al., 2017[START_REF] Arbel | Amortized Implicit Differentiation for Stochastic Bilevel Optimization[END_REF] which are known for being near-optimal algorithms for nonconvex finite sum minimization problems. It relies on a recursive estimation of directions of interest, which is restarted periodically. Proofs are deferred to the appendix.

Assumptions

Before presenting our algorithm, we formulate regularity Assumptions on the functions F and G.

Assumption 2.1. For all j ∈ [m], F j is twice differentiable and L F 0 -Lipschitz continuous. Its gradient is L F 1 -Lipschitz continuous and its Hessian is L F 2 -Lipschitz continuous.

Assumption 2.2. For all i ∈ [n], G i is three times differentiable. Its first, second, and third order derivatives are respectively

L G 1 -Lipschitz continuous, L G 2 -Lipschitz continuous, and L G 3 -Lipschitz continuous. For x ∈ R d , the function G i ( . , x) is µ G -strongly convex.
The strong convexity and the smoothness with respect to z hold for instance when we consider an ℓ 2 -regularized logistic regression problem with nonseparable data. These regularity assumptions up to first-order for F and second-order for G are standard in the stochastic bilevel literature [START_REF] Arbel | Amortized Implicit Differentiation for Stochastic Bilevel Optimization[END_REF]Ji et al., 2021;Yang et al., 2021). The secondorder regularity for F and third-order regularity for G are necessary for the analysis of the dynamics of v, as is the case in Dagréou et al. (2022). As shown in Ghadimi and Wang (2018, Lemma 2.2), these assumptions imply the smoothness of h, which is a fundamental property to get a descent.

Proposition 2.3. Under Assumptions 2.1 and 2.2, the function h is L h smooth for some L h > 0 which is precised in Appendix A.2.

Another consequence of Assumptions 2.1 and 2.2 is the boundedness of the function v * . Proposition 2.4. Assume that Assumptions 2.1 and 2.2 hold. Then, for R =

L F 0 µ G it holds that for any x ∈ R d , we have ∥v * (x)∥ ≤ R.
We denote Γ the closed ball centered in 0 with radius R and Π Γ the projection onto Γ.

For (z, v, x) ∈ R p × R p × R d , we denote Π(z, v, x) = (z, Π Γ (v), x).

Hypergradient Approximation

The gradient of h given by ( 2) is intractable in practice because it requires the perfect knowledge of z * (x) and v * (x) which are usually costly to compute. As classically done in the stochastic bilevel literature (Ji et al., 2021;[START_REF] Arbel | Amortized Implicit Differentiation for Stochastic Bilevel Optimization[END_REF]Li et al., 2022), z * (x) and v * (x) are replaced by approximate surrogate variables z and v. The variable z is typically the output of one or several steps of an optimization procedure applied to G( . , x). The variable v can be computed by using Neumann approximations or doing some optimization steps on the quadratic v →

1 2 v ⊤ ∇ 2 11 G(z, x)v + ∇ 1 F (z, x) ⊤ v.
We consider the approximate hypergradient given by

D x (z, v, x) = ∇ 2 21 G(z, x)v + ∇ 2 F (z, x) .
The motivation behind this direction is that if we take z = z * (x) and v = v * (x), we recover the true gradient, that is D x (z * (x), v * (x), x) = ∇h(x). Proposition 2.5 from (Dagréou et al., 2022, Lemma 3.4) controls the hypergradient approximation error by the distances between z and z * (x) and between v and v * (x).

Proposition 2.5. Let x ∈ R d . Assume that F is differentiable and L F 1 smooth with bounded gradient, G is twice differentiable with Lipschitz gradient and Hessian and G( . , x) is µ G -strongly convex. Then there exists a constant L x such that

∥D x (z, v, x)-∇h(x)∥ 2 ≤ L 2 x (∥z-z * (x)∥ 2 +∥v-v * (x)∥ 2 ).
Thus, it is natural to make z and v move towards their respective equilibrium values which are given by z * (x) and v * (x). As a consequence, we also introduce the directions D z and D x as follows

D z (z, v, x) = ∇ 1 G(z, x) , D v (z, v, x) = ∇ 2 11 G(z, x)v + ∇ 1 F (z, x) .
The interest of considering the directions D z and D v is expressed in Proposition 2.6.

Proposition 2.6. Assume that G is strongly convex with respect to its first variable. Then for any

x ∈ R d , it holds D z (z * (x), v * (x), x) = 0 and D v (z * (x), v * (x), x) = 0.
The directions D z , D v , and D x can be written as sums over the samples. Hence, following these directions enables to adapt any classical algorithm suited for single-level finite sum minimization to bilevel finite sum minimization. In what follows, for two indices i ∈ [n] and j ∈ [m], we consider the sampled directions D z,i,j , D v,i,j and D x,i,j defined by

D z,i,j (z, v, x) = ∇ 1 G i (z, x) (4) D v,i,j (z, v, x) = ∇ 2 11 G i (z, x)v + ∇ 1 F j (z, x) (5) D x,i,j (z, v, x) = ∇ 2 21 G i (z, x)v + ∇ 2 F j (z, x) . ( 6
)
When i and j are randomly sampled uniformly, these directions are unbiased estimators of the true directions D z , D v , and D x . Yet, as in Nguyen et al. (2017), we use them to recursively build biased estimators of the directions that enable fast convergence.

SRBA: Stochastic Recursive Bilevel Algorithm

In Algorithm 1, we present SRBA, a combination of the idea of recursive gradient coming from (Fang et al.,

Algorithm 1 Stochastic Recursive Bilevel Algorithm Input: initializations z 0 ∈ R p , x 0 ∈ R d , v 0 ∈ R p ,
number of iterations T and q, step sizes ρ and γ. Set ũ0 = (z 0 , v 0 , x 0 ) for t = 0, . . . , T -1 do Reset ∆: ∆ t,0 = (ρD z (ũ t ), ρD v (ũ t ), γD x (ũ t )) Update u: u t,1 = Π(ũ t -∆ t,0 ) , for k = 1, . . . , q -1 do Draw i ∈ {1, . . . , n} and j ∈ {1, . . . , m}

∆ t,k z = ρ(D z,i,j (u t,k ) -D z,i,j (u t,k-1 )) + ∆ t,k-1 z ∆ t,k v = ρ(D v,i,j (u t,k ) -D v,i,j (u t,k-1 )) + ∆ t,k-1 v ∆ t,k x = γ(D x,i,j (u t,k ) -D x,i,j (u t,k-1 )) + ∆ t,k-1
x Update u: u t,k+1 = Π(u t,k -∆ t,k ) end for Set ũt+1 = u t+1,q end for Return (z T , ṽT , xT ) = ũT 2018; Nguyen et al., 2022) and the framework proposed in (Dagréou et al., 2022). It relies on a recursive estimation of each direction D z , D v , D x which is updated following the same strategy as SARAH. Let us denote by ρ the step size of the update for the variables z and v, and γ the step size for the update of the variable x. We use the same step size for z and v because the problems of minimizing the inner function G and solving the linear system (3) have the same conditioning driven by ∇ 2 11 G. For simplicity, we denote the joint variable u = (z, v, x) and the joint directions weighted by the step sizes ∆ = (ρD z , ρD v , γD x ) = (∆ z , ∆ v , ∆ x ).

At iteration t, the estimate direction ∆ is initialized by computing full batch directions:

∆ t,0 = (ρD z (ũ t ), ρD v (ũ t ), γD x (ũ t ))
and a first update is performed by moving from ũt in the direction -∆ t,0 . As done in Hu et al. (2022), we project the variable v onto Γ to leverage the boundedness property of v * . Then, during the kth iteration of an inner loop of size q -1, two indices i ∈ [n] and j ∈ [m] are sampled and the estimate directions are updated according to Equations ( 7) to (9)

∆ t,k z = ρ(D z,i,j (u t,k ) -D z,i,j (u t,k-1 )) + ∆ t,k-1 z (7) ∆ t,k v = ρ(D v,i,j (u t,k ) -D v,i,j (u t,k-1 )) + ∆ t,k-1 v (8) ∆ t,k x = γ(D x,i,j (u t,k ) -D x,i,j (u t,k-1 )) + ∆ t,k-1 x (9)
where the sampled directions D z,i,j , D v,i,j and D x,i,j are defined by Equations ( 4) to (6). Then the joint variable u is updated by

u t,k+1 = Π(u t,k -∆ t,k ) . ( 10
)
Recall that the projection is only performed on the variable v. The other variables z and x remain un-changed after the projection step. At the end of the inner procedure, we set ũt+1 = u t,q .

In Algorithm 1, the variables z, v, and x are updated simultaneously rather than alternatively. From a computational perspective, this enables sharing the common computations between the different oracles and doing the update of each variable in parallel. So there is no sub-procedure to approximate the solution of the inner problem and the solution of the linear system.

Theoretical Analysis of SRBA

In this section we provide the theoretical analysis of Algorithm 1 leading to a final sample complexity in O (n + m)

1 2 ε -1 ∨ (n + m) .
The detailed proofs of the results are deferred to the appendix. In Definition 3.1, we recall what is an ε-stationary point.

Definition 3.1. Let d a positive integer, f : R d → R a differentiable function and ε > 0. We say that a point

x ∈ R d is an ε-stationary point of f if ∥∇f (x)∥ 2 ≤ ϵ. In a stochastic context, we call ε-stationary point a random variable x such that E[∥∇f (x)∥ 2 ] ≤ ε.
In this paper, the theoretical complexity of the algorithms is given in terms of number of calls to oracle, that is to say, the number of times the quantity

[∇F j (z, x), ∇G i (z, x), ∇ 2 11 G i (z, x)v, ∇ 2 21 G i (z, x)v] (11) is queried for i ∈ [n], j ∈ [m], z ∈ R p , v ∈ R p and x ∈ R d .
Note that in practice, although the secondorder derivatives of the inner functions ∇ 2 11 G i (z, x) ∈ R p×p and ∇ 2 21 G i (z, x) ∈ R d×p are involved, they are never computed or stored explicitly. We rather work with Hessian-vector products ∇ 2 11 G i (z, x)v ∈ R p and Jacobian-vector products ∇ 2 21 G i (z, x)v ∈ R d which can be computed efficiently thanks to automatic differentiation with a computational cost similar to the cost of computing the gradients ∇ 1 G i (z, x) and ∇ 2 G i (z, x) Pearlmutter (1994). The cost of one query (11) is therefore of the same order of magnitude as that of computing one stochastic gradient.

Mean Squared Error of the Estimated Directions

A strength of our method is its simple expression of the estimation error of the directions coming from the bias-variance decomposition provided by Nguyen et al. (2017). Let us denote the estimate directions

D t,k z = ∆ t,k z /ρ, D t,k v = ∆ t,k v /ρ and D t,k x = ∆ t,k
x /γ. We also introduce the residuals

S t,k • = k r=1 E[∥D • (u t,r ) -D • (u t,r-1 )∥ 2 ], St,k • = k r=1 E[∥D t,r • -D t,r-1 • ∥ 2 ] .
We provide a link between the mean squared error

E[∥D t,k • -D • (u t,k )∥ 2
] and the residuals. Proposition 3.2 (MSE of the estimate directions). For any t ≥ 0 and k ∈ {1, . . . , q -1}, the estimate D t,k

• of the direction D • (u t,k ) satisfies E[∥D t,k • -D • (u t,k )∥ 2 ] = St,k • -S t,k • .
The above error has two components: the accumulation of the difference between two successive full batch directions and the accumulation of the difference between two successive estimate directions.

Fundamental Lemmas

We establish descent lemmas which are key ingredients to get the final convergence result. Lemma 3.3 characterizes the dynamic of u on the inner problem. To do so, we define the function ϕ z as

ϕ z (z, x) = G(z, x) -G(z * (x), x) .
In the bilevel literature, direct control on the distance to optimum δ

t,k z ≜ E[∥z t,k -z * (x t,k )∥ 2 ] is es- tablished.
Here, the biased nature of the estimate direction D t,k z makes it hard to upper bound appropriately the scalar product ⟨D z (u t,k )-D t,k z , z t,k -z * (x t,k )⟩. Therefore, we rather consider ϕ t,k z . By combining the smoothness property of ϕ z and the bias-variance decomposition provided in Proposition 3.2, we can show some descent property on the sequence ϕ

t,k z defined by ϕ t,k z = E[ϕ z (z t,k , x t,k )]. Before stating Lemma 3.3, let us define G t,k v = 1 ρ v t,k -Π Γ (v t,k -ρD t,k v ) so that v t,k+1 = v t,k -ρG t,k
v . This is the actual update direction of v. If there were no projections, we would have

G t,k v = D t,k v .
Hence, it acts as a surrogate of D t,k v in our analysis. We also define

V t,k z = E[∥D t,k z ∥ 2 ], V t,k v = E[∥G t,k v ∥ 2 ], V t,k x = E[∥D t,k x ∥ 2 ]
the variances and their respective sums over the inner loop

V t,k z = k r=1 E[∥D t,r-1 z ∥ 2 ], V t,k v = k r=1 E[∥G t,r-1 v ∥ 2 ], V t,k x = k r=1 E[∥D t,r-1 x ∥ 2 ] .
Lemma 3.3 (Descent on the inner level). Assume that the step sizes ρ and γ verify γ ≤ C z ρ for some positive constant C z specified in the appendix. Then it holds

ϕ t,k+1 z ≤ 1 - µ G 2 ρ ϕ t,k z - ρ 2 (1 -Λ z ρ) V t,k z (12) + ρ 3 β zz V t,k z + γ 2 ρβ zv V t,k v + γ 2 ρβ zx V t,k x + Λ z 2 γ 2 V t,k x + γ 2 ρ β zx E[∥D x (u t,k )∥ 2 ]
for some positive constants Λ z , β zz , β zx and β zx that are specified in the appendix.

In ( 12) we recover a decrease of ϕ t,k z by a factor (1ρµ G ). But the outer variable's movement and the noise make appear D x (u t,k ) and the variance hindering the convergence of z towards z * (x).

For the variable v, the quantity we consider is

ϕ v (v, x) = Ψ(z * (x), v, x) -Ψ(z * (x), v * (x), x) where Ψ(z, v, x) is defined as Ψ(z, v, x) = 1 2 v ⊤ ∇ 2 11 G(z, x)v + ∇ 1 F (z, x) ⊤ v .
The intuition behind considering this quantity is that solving the linear system (3) is equivalent to minimizing over v the function Ψ(z * (x), v, x).

Lemma 3.4. Assume that the step sizes ρ and γ verify ρ ≤ B v and γ ≤ C v ρ for some positive constants B v and C v specified in the appendix. Then it holds

ϕ t,k+1 v ≤ 1 - ρµ G 16 ϕ t,k v -βvv ρV t v + ρ 3 β vz V t,k z + 2ρ 3 β vv V t,k v + γ 2 ρβ vx V t,k x + ρα vz ϕ t,k z + Λ v 2 γ 2 E[∥D t,k x ∥ 2 ] + γ 2 ρ β vx E[∥D x (u t,k )∥ 2 ]
for some positive constants Λ v , β vz , β vx , βvv and β vx that are specified in the appendix.

Lemma 3.4 is similar to Lemma 3.3 with a term in ϕ t,k z taking into account the error of z * (x)'s approximation. Its proof harnesses the generalization of Polyak-Łojasiewicz inequality for composite functions introduced in Karimi et al. (2016).

The following lemma is a consequence of the smoothness of h.

Let us denote the expected values h t,k = E[h(x t,k )] and expected gradient

g t,k = E[∥∇h(x t,k )∥ 2 ]. Lemma 3.5. There exist constants β hz , β hv , β hx > 0 such that h t,k+1 ≤ h t,k - γ 2 g t,k + γ 2L 2 x µ G (ϕ t,k z + ϕ t,k v ) + γρ 2 + ρα vz ϕ t,k z + γρ 2 β hv V t,k v + γ 3 β hx V t,k x - γ 2 1 -L h γ V t,k x .
This lemma shows that the control of the approximation error ϕ • (Lemma 3.3 and Lemma 3.4) and the sum of variances

V • is crucial to get a decrease of E[h(x t,k )].

Complexity Analysis of SRBA

In Theorem 1, we provide the convergence rate of SRBA towards a stationary point.

Theorem 1. Assume that Assumptions 2.1 and 2.2 hold. Assume that the step sizes verify ρ ≤ ρ and γ ≤ min(γ, ξρ) for some constants ξ, ρ and γ specified in appendix. Then it holds

1 T q T -1 t=0 q-1 k=0 E[∥∇h(x t,k )∥ 2 ] = O 1 qT γ
where O hides regularity constants that are independent from n and m.

The proof combines classical proof techniques from the bilevel literature and elements from SARAH's analysis (Nguyen et al., 2017(Nguyen et al., , 2022)). We introduce the Lyapunov function

L(u t,k ) = h t,k + ψ z ϕ t,k z + ψ v ϕ t,k
v where ψ z and ψ v are non-negative constants chosen so that we have the inequality

L(u t,k+1 ) ≤ L(u t,k ) -γ 4 g t,k
. Summing and telescoping this inequality provides the result.

Note that increasing q allows a faster convergence in terms of iterations but makes each iteration more expensive since the number of oracle calls per iteration is (2n + 3m) + 2 × 5(q -1). Thus, there is a trade-off between the convergence rate and the overall complexity. In Corollary 3.6, we state that the value of q that gives the best sample complexity is O(n + m).

Corollary 3.6. Suppose that Assumptions 2.1 and 2.2 hold.

If we take ρ

= ρ(n + m) -1 2 , γ = min(γ, ξρ)(n + m) -1 2 and q = n + m, then O (n + m) 1 2 ε -1 ∨ (n + m
) calls to oracles are sufficient to find an ε-stationary point with SRBA.

This sample complexity is analogous to the sample complexity of SARAH in the nonconvex finite-sum setting. To the best of our knowledge, such a rate is the best known for bilevel empirical risk minimization problems in terms of dependency on the number of samples n + m and the precision ε. This improves by a factor (n+m) -1 6 the previous result which was achieved by SABA (Dagréou et al., 2022). As a comparison, VRBO (Yang et al., 2021) achieves a sample complexity in Õ(ε -3

2 ). Note that, for large value of n + m we can have actually (n + m)

1 2 ε -1 ≳ ε -2
. This means that, just like single-level SARAH, the complexity of SRBA can be beaten by others when the number of samples is too high with respect to the desired accuracy (actually if n + m = Ω(ε -2 )).

Lower Bound for Bilevel ERM

In this section, we derive a lower bound for bilevel empirical risk minimization problems. This shows that SRBA is a near-optimal algorithm for this class of problems.

Function and Algorithm Classes

We define the function and algorithm classes we consider.

((F j ) 1≤j≤m , (G i ) 1≤i≤n ) defined on R p × R d such that for all j ∈ [m], F j is L F 1 smooth and for all i ∈ [n], G i is twice differentiable and µ G -strongly convex.
Note that we consider a class of nonconvex bilevel problems. This class contains, the functions defining the bilevel formulation of the datacleaning task.

For the algorithmic class, we consider algorithms that use approximate implicit differentiation. Definition 4.2. Given initial points z 0 , v 0 , x 0 , a linear bilevel algorithm A is a measurable mapping such that for any

((F j ) 1≤j≤m , (G i ) 1≤i≤n ) ∈ C L F 1 ,µ G , the output of A((F j ) 1≤j≤m , (G i ) 1≤i≤n ) is a sequence {(z t , v t , x t , i t , j t )} t≥0 of points (z t , v t , x t ) and random variables i t ∈ [n] and j t ∈ [m] such that for all t ≥ 0 z t+1 ∈ z 0 +Span{∇ 1 G i0 (z 0 , x 0 ), . . . , ∇ 1 G it (z t , x t )} v t+1 ∈ v 0 +Span{∇ 2 11 G i0 (z 0 , x 0 )v 0 + ∇ 1 F j0 (z 0 , x 0 ), . . . , ∇ 2 11 G it (z t , x t )v t + ∇ 1 F jt (z t , x t )} x t+1 ∈ x 0 +Span{∇ 2 21 G i0 (z 0 , x 0 )v 0 + ∇ 2 F j0 (z 0 , x 0 ), . . . , ∇ 2 21 G it (z t , x t )v t + ∇ 2 F jt (z t , x t )}.
This algorithm class includes popular stochastic bilevel first-order algorithms, such as AmIGO [START_REF] Arbel | Amortized Implicit Differentiation for Stochastic Bilevel Optimization[END_REF]), FSLA (Li et al., 2022), SOBA, and SABA (Dagréou et al., 2022). Moreover, despite the projection step, SRBA is part of this algorithm class since the projection of a vector onto Γ is actually just a rescaling.

Main Theorem Problem (1) is actually a smooth nonconvex optimization problem. The lower complexity bound for nonconvex finite sum problem has been studied in Fang et al. ( 2018); Zhou and Gu (2019).

In particular, they show that the number of gradient calls needed to get an ε-stationary point for a smooth nonconvex finite sum is at least O(m

1 2 ε -1
), where m is the number of terms in the finite sum.

Intuitively, we expect the lower complexity bound to solve (1) to be larger. Indeed, bilevel problems are harder than single-level problems because a bilevel problem involves the resolution of several subproblems to progress in its resolution. Theorem 2 formalizes this intuition by showing that the classical single-level lower bound is also a lower bound for bilevel problems.

Theorem 2. For any linear bilevel algorithm A, and any

L F , n, ∆, ε, p such that ε ≤ (∆L F m -1 )/10 3 , there exists a dimension d = O(∆ε -1 m 1 2 L F ), an ele- ment ((F j ) 1≤j≤m , (G i ) 1≤i≤n ) ∈ C L F 1 ,µ G such that the value function h defined as in (1) satisfies h(x 0 ) - inf x∈R d h(x) ≤ ∆ and in order to find x ∈ R d such that E[∥∇h(x)∥ 2 ] ≤ ε, A needs at least Ω(m 1 2 ε -1 ) calls to oracles of the form (11).
The proof is an adaptation of the proof of Zhou and Gu (2019, Theorem 4.7). We take as outer function F defined by F (z, x) = m j=1 f (U (j) z) where f is the "worst-case function" used by Carmon et al. (2021), U = [U (j) , . . . , U (m) ] ⊤ is an orthogonal matrix and G(z, x) = 1 2 ∥z -x∥ 2 . We leverage the fact that ∥∇f (y)∥ 2 > K as long as the two last coordinates of y are zero for some known constant K. Then we use the "zero chain property" to bound the number of indices j such the two last components of U (j) x t are zero at a given iteration t, implying ∥∇h(x t )∥ 2 > ϵ when t is smaller than O(m

1 2 ε -1 ).
As a comparison to the existing lower bound for bilevel optimization in Ji and Liang (2023), we consider randomized algorithms and do not assume the value function h to be convex or strongly convex.

Numerical Experiments

Even though our contribution is mostly theoretical, we run several experiments to highlight to compare the proposed algorithm with state-of-the-art stochastic bilevel solvers. We compare our method to AmIGO [START_REF] Arbel | Amortized Implicit Differentiation for Stochastic Bilevel Optimization[END_REF], F 2 SA (Kwon et al., 2023), MRBO (Yang et al., 2021), VRBO (Yang et al., 2021), StocBiO (Ji et al., 2021) and SABA (Dagréou et al., 2022). They are run on a synthetic problem with quadratic functions and on a hyperparameter selection problem for ℓ 2 -regularized logistic regression with the dataset IJCNN11 . A more detailed description of the experiments is available in Appendix C and an additional experiment the datacleaning task is available in Appendix D.

Experiments on quadratics To evaluate the performance of stochastic bilevel optimizers in a controlled setting, we perform a benchmark on quadratic loss functions described in Appendix C. Here F and G are quadratic jointly in (z, x), allowing us to choose freely the conditioning of F , G, and h. We take for the Hessian and cross derivative matrices of each sample, the empirical correlation of random vectors drawn with a prescribed covariance matrix. The generation process is detailed in Appendix C. In Figure 1, we report the norm of the gradient of the value function function with respect to time. Our first observation is that among all the methods, SRBA and SABA converge the fastest. These two solvers share two key ingredients: variance reduction and warm-starting. Variance reduction makes the variance of the gradient estimate go to zero without using decreasing step sizes. The warm-starting strategy in both the approximation of z * (x t ) and the approximation of v * (x t ) enables getting an estimator of ∇h(x t ) which is asymptotically unbiased, without requiring an increasing number of inner iterations or batch-size. Note that solvers using Neumann iterations (VRBO, MRBO, stocBiO) fail to converge because Neumann iterations provide a biased estimate of v * (x). Moreover, AmiIGO and stocBiO evolve slowly after some iterations because they require vanishing step sizes to converge. Finally, SRBA is faster than SABA, which is consistent with the theory.

Hyperparameter selection We also run an experiment on hyperparameter selection problem for ℓ 2regularized logistic regression with the IJCNN1 dataset. SRBA shows good performances in the experiment, both in speed and accuracy. It is competitive with other state-of-the-art methods AmIGO and SABA, while going faster than Amigo and requiring less memory than SABA. VRBO -another extension of SARAH for bilevel problems-is slower in all problems. This is due to the burden of computing the approximate hypergradient at each inner iteration without updating the outer parameter. We can also notice that in the experiment on IJCNN1, the slowest method are method implementing Neumann approximations to approximate v * (x). Note that this last experiment does not include F 2 SA because we find that on this problem, the norm of the iterates of F 2 SA goes towards infinity.

10 -2 10 -1 h(x t ) -h * AmIGO MRBO VRBO SABA StocBiO SRBA
Figure 1: Comparison of the behavior of SRBA with other stochastic bilevel solvers. For each experiment, the solvers are run with 10 different seeds and the median performance over these seeds is reported. The shaded area corresponds to the performances between the 20% and the 80% percentiles. The performances are reported with respect to wall-clock time. Top: Experiments on quadratic functions. We report the gradient norm of the value function. Bottom: Hyperparameter selection with the IJCNN1 dataset.

Conclusion

In this paper, we have introduced SRBA, an algorithm for bilevel empirical risk minimization. We have demonstrated that the sample complexity of SRBA is O((n + m) 1 2 ε -1 ) for any bilevel problem where the inner problem is strongly convex. Then, we have demonstrated that any bilevel empirical risk minimization algorithm has a sample complexity of at least O(m 1 2 ε -1 ) on some problems where the inner problem is strongly convex. This demonstrates that SRBA is optimal, up to constant factors, and that bilevel ERM is as hard as single-level nonconvex ERM.

D v (z * (x), v, x) = ∇ 2 11 G(z * (x), x)v + ∇ 1 F (z * (x), x) = 0 (13)
admits a unique solution given by v * (x).

A.2 Smoothness constant of h

We can find in Ghadimi and Wang (2018, Lemma 2.2) the following value for the smoothness constant of h

L h = L F 1 + 2L F 1 L G 2 + (L F 0 ) 2 L G 2 µ G + L G 11 L G 1 L F 0 + L G 1 L G 2 L F 0 + (L G 1 ) 2 L F 1 µ 2 G + (L G 1 ) 2 L G 2 L F 0 µ 3 G . A.3 Proof of Proposition 3.2 Proof. Let t > 0 and k ∈ [q -1]. For k = 0, we directly have E[∥D t,k • -D • (u t,k )∥ 2 ] = 0.
For k ≥ 1 and r ∈ {1, . . . , k}, the bias/variance decomposition of D t,r

• reads

E t,r [∥D t,r • -D • (u t,r )∥ 2 ] = E t,r [∥D t,r • -D t,r-1 • + D • (u t,r-1 ) -D • (u t,r )∥ 2 ] + ∥D • (u t,r ) + D • (u t,r-1 ) -D t,r-1 • -D • (u t,r )∥ 2 = E t,r [∥D t,r • -D t,r-1 • -(D • (u t,r-1 ) -D • (u t,r ))∥ 2 ] + ∥D t,r-1 • -D • (u t,r-1 )∥ 2
The term

E t,r [∥D t,r • -D t,r-1 • -(D • (u t,r-1 ) -D • (u t,r ))∥ 2 ] is the variance of D t,r • -D t,r-1 •
, and then can written as

E t,r [∥D t,r • -D t,r-1 • -(D • (u t,r-1 ) -D • (u t,r ))∥ 2 ] = E t,r [∥D t,r • -D t,r-1 • ∥ 2 ] -∥D • (u t,r ) -D • (u t,r-1 )∥ 2
Plugging this in the previous inequality and taking the total expectation leads to

E[∥D t,r • -D • (u t,r )∥ 2 ] = E[∥D t,r • -D t,r-1 • ∥ 2 ] -E[∥D • (u t,r ) -D • (u t,r-1 )∥ 2 ] + E[∥D t,r-1 • -D t,r-1 • (u t,r-1 )∥ 2 ]
Summing for r ∈ {1, . . . , k} and telescoping gives the final result (taking into account that D t,0

• = D • (u t,0 )): E[∥D t,k • -D • (u t,k )∥ 2 ] = k r=1 E[∥D t,r • -D t,r-1 • ∥ 2 ] - k r=1 E[∥D • (u t,r ) -D • (u t,r-1 )∥ 2 ] .

A.4 Technical lemmas

Lemma A.1. There exists constant L z * and L v * such that for any x 1 , x 2 ∈ R d , we have

∥z * (x 1 ) -z * (x 2 )∥ ≤ L z * ∥x 1 -x 2 ∥ and ∥v * (x 1 ) -v * (x 2 )∥ ≤ L v * ∥x 1 -x 2 ∥ Proof. The Jacobian of z * reads dz * (x) = [∇ 2 11 G(z * (x), x)] -1 ∇ 2 12 G(z * (x), x). By µ G -strong convexity and L G 1 -smoothness of G, we have ∥dz * (x)∥ ≤ L G 1 µ G which implies that z * is L z * -Lipschtiz with L z * = L G 1 µ G .
For v * we do the computation directly:

∥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 )∥ ≤ ∥[∇ 2 11 G(z * (x 1 ), x 1 )] -1 (∇ 1 F (z * (x 1 ), x 1 ) -∇ 1 F (z * (x 2 ), x 2 ))∥ + ∥([∇ 2 11 G(z * (x 1 ), x 1 ) -[∇ 2 11 G(z * (x 2 ), x 2 )] -1 ] -1 ∇ 1 F (z * (x 2 ), x 2 )∥ ≤ L F 1 µ G + L G 2 L F 0 µ 2 G ∥(z * (x 1 ), x 1 ) -(z * (x 2 ), x 2 )∥ ≤ L F 1 µ G + L G 2 L F 0 µ 2 G (∥z * (x 1 ) -z * (x 2 )∥ + ∥x 1 -x 2 ∥) ≤ 1 + L G 1 µ G L F 1 µ G + L G 2 L F 0 µ 2 G ∥x 1 -x 2 ∥ Then taking L v * = 1 + L G 1 µ G L F 1 µ G + L G 2 L F 0 µ 2 G concludes the proof. Lemma A.2. Let us consider the update directions D t,k z = ∆ t,k z /ρ, D t,k v = ∆ t,k v /ρ and D t,k x = ∆ t,k x /γ where ∆ t,k z , ∆ t,k v and ∆ t,k
x verify Equations (7) to (9). Then it holds

E[∥D t,k z -D z (u t,k )∥ 2 ] ≤ k r=1 L G 1 (ρ 2 E[∥D t,r-1 z ∥ 2 ] + γ 2 E[∥D t,r-1 z ∥ 2 ]) E[∥D t,k v -D v (u t,k )∥ 2 ] ≤ 4ρ 2 (L G 2 R) 2 + (L F 1 ) 2 k r=1 E[∥D t,r-1 z ∥ 2 ] + 4ρ 2 (L G 1 ) 2 k r=1 E[∥G t,r-1 v ∥ 2 ] + 4γ 2 (L G 2 R) 2 + (L F 1 ) 2 k r=1 E[∥D t,r-1 x ∥ 2 ] E[∥D t,k x -D x (u t,k )∥ 2 ] ≤ 4ρ 2 (L G 2 R) 2 + (L F 1 ) 2 k r=1 E[∥D t,r-1 z ∥ 2 ] + 4ρ 2 (L G 1 ) 2 k r=1 E[∥G t,r-1 v ∥ 2 ] + 4γ 2 (L G 2 R) 2 + (L F 1 ) 2 k r=1 E[∥D t,r-1 x ∥ 2 ] .

Proof. Direction D z

We start from Proposition 3.2.

E[∥D t,k z -D z (u t,k )∥ 2 ] = E[∥D t,k z -∇ 1 G(z t,k , x t,k )∥ 2 ] = k r=1 E[∥D t,r z -D t,r-1 z ∥ 2 ] - k r=1 E[∥∇ 1 G(z t,r , x t,r ) -∇ 1 G(z t,r-1 , x t,r-1 )∥ 2 ] ≤ k r=1 E[∥D t,r z -D t,r-1 z ∥ 2 ] ≤ k r=1 L G 1 (ρ 2 E[∥D t,r-1 z ∥ 2 ] + γ 2 E[∥D t,r-1 z ∥ 2 ])
where the last inequality comes from the smoothness of each G i .

Direction D v For D v , the proof is almost the same. Proposition 3.2 gives us

E[∥D t,k v -D v (u t,k )∥ 2 ] ≤ k r=1 E[∥D t,r v -D t,r-1 v ∥ 2 ] .
Then, using the boundedness of v and regularity of each G i and F j , we have

E[∥D t,r v -D t,r-1 v ∥ 2 ] ≤ 2(E[∥∇ 2 11 G i (z t,r , x t,r )v t,r -∇ 2 11 G i (z t,r-1 , x t,r-1 )v t,r-1 ∥ 2 ] + E[∥∇ 2 F j (z t,r , x t,r ) -∇ 2 F j (z t,r-1 , x t,r-1 )∥ 2 ]) ≤ 4(E[∥∇ 2 11 G i (z t,r , x t,r )(v t,r -v t,r-1 )∥ 2 ] + E[∥(∇ 2 11 G i (z t,r , x t,r ) -∇ 2 11 G i (z t,r-1 , x t,r-1 ))v t,r-1 ∥ 2 ] + (L F 1 ) 2 (γ 2 E[∥D t,r-1 z ∥] + ρ 2 E[∥D t,r-1 x ∥ 2 ])) ≤ 4((L G 1 ) 2 ρ 2 E[∥G t,r-1 v ∥ 2 ] + (L G 2 ) 2 R 2 (ρ 2 E[∥D t,r-1 z ∥] + γ 2 E[∥D t,r-1 x ∥ 2 ]) + (L F 1 ) 2 (ρ 2 E[∥D t,r-1 z ∥] + γ 2 E[∥D t,r-1 x ∥ 2 ])) ≤ 4ρ 2 (L G 2 R) 2 + (L F 1 ) 2 E[∥D t,r-1 z ∥ 2 ] + 4ρ 2 (L G 1 ) 2 E[∥G t,r-1 v ∥ 2 ] + 4γ 2 (L G 2 R) 2 + (L F 1 ) 2 E[∥D t,r-1 x ∥ 2 ] .
Direction D x The proof is the same as the proof for D v .

A.5 Proof of Lemma 3.3

Let ϕ z (z, x) = G(z, x) -G(z * (x), x)
the inner suboptimality gap. The proof of Lemma 3.3 is based on the smoothness of ϕ z , which is the object of the following lemma.

Lemma A.3. The function ϕ z has Λ z -Lipschitz continuous gradient on R p × R d , for some constant Λ z .

Proof. For any (z, x) ∈ R p × R d , we have

∇ 1 ϕ z (z, x) = ∇ 1 G(z, x) and ∇ 2 ϕ z (z, x) = ∇ 2 G(z, x) -∇ 2 G(z * (x), x) . Let us consider (z, x) ∈ R p × R d and (z ′ , x ′ ) ∈ R p × R d . Since ∇G is L G 1 -Lipschitz continuous, we have directly ∥∇ 1 ϕ z (z, x) -∇ 1 ϕ z (z ′ , x ′ )∥ ≤ L G 1 ∥(z, x) -(z ′ , x ′ )∥ . Moreover, we have ∥∇ 2 ϕ z (z, x) -∇ 2 ϕ z (z ′ , x ′ )∥ ≤ ∥∇ 2 G(z, x) -∇ 2 G(z ′ , x ′ )∥ + ∥∇ 2 G(z * (x), x) -∇ 2 G(z * (x ′ ), x ′ )∥ ≤ L G 1 ∥(z, x) -(z ′ , x ′ )∥ + L G 1 ∥(z * (x), x) -(z * (x ′ ), x ′ )∥ ≤ L G 1 ∥(z, x) -(z ′ , x ′ )∥ + L G 1 (∥z * (x) -z * (x ′ )∥ + ∥x -x ′ ∥) . From Lemma A.1, z * is L * Lipschitz continuous, so ∥∇ 2 ϕ z (z, x) -∇ 2 ϕ z (z ′ , x ′ )∥ ≤ L G 1 ∥(z, x) -(z ′ , x ′ )∥ + L G 1 (∥z * (x) -z * (x ′ )∥ + ∥x -x ′ ∥) ≤ L G 1 ∥(z, x) -(z ′ , x ′ )∥ + L G 1 (L * + 1)∥x -x ′ ∥ ≤ L G 1 (L z * + 2)∥(z, x) -(z ′ , x ′ )∥ . As a consequence ∥∇ϕ z (z, x) -∇ϕ z (z ′ , x ′ )∥ ≤ ∥∇ 1 ϕ z (z, x) -∇ 1 ϕ z (z ′ , x ′ )∥ + ∥∇ 2 ϕ z (z, x) -∇ 2 ϕ z (z ′ , x ′ )∥ ≤ L G 1 (L z * + 3)∥(z, x) -(z ′ , x ′ )∥ . Hence, ϕ z is Λ z smooth with Λ z = L G 1 (L z * + 3).
We can now turn to the proof of Lemma 3.3.

Proof. The smoothness of ϕ z provides us the following upper bound

ϕ z (z t,k+1 , x t,k+1 ) ≤ ϕ z (z t,k , x t,k ) -ρ⟨D t,k z , ∇ 1 G(z t,k , x t,k )⟩ + Λ z 2 ρ 2 ∥D t,k z ∥ 2 (14) -γ⟨D t,k x , ∇ 2 G(z t,k , x t,k ) -∇ 2 G(z * (x t,k ), x t,k )⟩ + Λ z 2 γ 2 ∥D t,k x ∥ 2 .
Using the equality ⟨a, b⟩ = 1 2 (∥a∥ 2 + ∥b∥ 2 -∥a -b∥ 2 ), we get

-⟨D t,k z , ∇ 1 G(z t,k , x t,k )⟩ + Λ z 2 ρ∥D t,k z ∥ 2 = 1 2 (∥D t,k z -∇ 1 G(z t,k , x t,k )∥ 2 (15) -∥∇ 1 G(z t,k , x t,k )∥ 2 -(1 -Λ z ρ) ∥D t,k z ∥ 2 ) .
Plugging Equation ( 15) into Equation ( 14) and tacking the expectation conditionally to the past iterates yields

E t,k [ϕ t,k+1 z ] ≤ ϕ t,k z + ρ 2 E t,k [∥D t,k z -∇ 1 G(z t,k , x t,k )∥ 2 ] (16) - ρ 2 ∥∇ 1 G(z t,k , x t,k )∥ 2 - ρ 2 (1 -Λ z ρ) E t,k [∥D t,k z ∥ 2 ] -γ⟨E t,k [D t,k x ], ∇ 2 G(z t,k , x t,k ) -∇ 2 G(z * (x t,k ), x t,k )⟩ + Λ z 2 γ 2 E t,k [∥D t,k x ∥ 2 ] .
From Young inequality, we have for any c > 0

⟨E t,k [D t,k x ], ∇ 2 G(z t,k , x t,k ) -∇ 2 G(z * (x t,k ), x t,k )⟩ ≤ 1 2c ∥E t,k [D t,k x ]∥ 2 (17) + c 2 ∥∇ 2 G(z t,k , x t,k ) -∇ 2 G(z * (x t,k ), x t,k )∥ 2
The smoothness of G and strong convexity give us

∥∇ 2 G(z t,k , x t,k ) -∇ 2 G(z * (x t,k ), x t,k )∥ 2 ≤ L G 1 ∥z t,k -z * (x t,k )∥ 2 ≤ 2L G 1 µ G ϕ z (z t,k , x t,k ) (18) Let us denote L ′ = L G 1 µ G .
Plugging Inequalities ( 17) and ( 18) into Equation ( 16) yields

E t,k [ϕ z (z t,k+1 , x t,k+1 )] ≤ (1 + cL ′ γ)ϕ z (z t,k+1 , x t,k+1 ) - ρ 2 E t,k [∥∇ 1 G(z t,k , x t,k )∥ 2 ] (19) + ρ 2 E t,k [∥D t,k z -∇ 1 G(z t,k , x t,k )∥ 2 ] - ρ 2 (1 -Λ z ρ) E t,k [∥D t,k z ∥ 2 ] + γ 2c ∥E t,k [D t,k x ]∥ 2 + Λ z 2 γ 2 E t,k [∥D t,k x ∥ 2 ]
From Lemma A.2, we have

E[∥D t,k z -∇ 1 G(z t,k , x t,k )∥ 2 ] ≤ k r=1 L G 1 (ρ 2 E[∥D t,r-1 z ∥ 2 ] + γ 2 E[∥D t,r-1 z ∥ 2 ]) .
Taking the total expectation and plugging the previous inequality into Equation ( 19) yields

ϕ t,k+1 z ≤ (1 + cL ′ γ)ϕ t,k + L G 1 2 k r=1 (ρ 3 E[∥D t,r-1 z ∥ 2 ] + γ 2 ρE[∥D t,r-1 x ∥ 2 ]) (20) - ρ 2 E[∥∇ 1 G(z t,k , x t,k )∥ 2 ] - ρ 2 (1 -Λ z ρ) E[∥D t,k z ∥ 2 ] + γ 2c E[∥E[D t,k x ]∥ 2 ] + Λ z 2 γ 2 E[∥D t,k x ∥ 2 ]
Since G is µ G -strongly convex with respect to z, Polyak-Łojasiewicz inequality holds:

∥∇ 1 G(z t,k , x t,k )∥ 2 ≥ 2µ G ϕ z (z t,k , x t,k )
As a consequence, Equation (20) becomes

ϕ t,k+1 z ≤ (1 + cL ′ γ -µ G ρ) ϕ t,k + L G 1 2 k r=1 (ρ 3 E[∥D t,r-1 z ∥ 2 ] + γ 2 ρE[∥D t,r-1 x ∥ 2 ]) - ρ 2 (1 -Λ z ρ) E[∥D t,k z ∥ 2 ] + γ 2c E[∥E[D t,k x ]∥ 2 ] + Λ z 2 γ 2 E[∥D t,k x ∥ 2 ] Taking c = µ G ρ 2L ′ γ yields ϕ t,k+1 z ≤ 1 - µ G 2 ρ ϕ t,k + L G 1 2 k r=1 (ρ 3 E[∥D t,r-1 z ∥ 2 ] + γ 2 ρE[∥D t,r-1 x ∥ 2 ]) - ρ 2 (1 -Λ z ρ) E[∥D t,k z ∥ 2 ] + L ′ µ G γ 2 ρ E[∥E[D t,k x ]∥ 2 ] + Λ z 2 γ 2 E[∥D t,k x ∥ 2 ] For the term E[∥E t,k [D t,k z ]∥ 2 ],
we have

E[∥E t,k [D t,k x ]∥ 2 ] = E[∥D x (z t,k , v t,k , x t,k ) -D x (z t,k-1 , v t,k-1 , x t,k-1 ) + D t,k-1 x ∥ 2 ] = E[∥D x (z t,k , v t,k , x t,k ) -D x (z t,k-1 , v t,k-1 , x t,k-1 ) -E[D t,k-1 x ]∥ 2 ] + E[∥D t,k-1 x -E[D t,k-1 x ]∥ 2 ] = E[∥D x (z t,k , v t,k , x t,k )∥ 2 ] (21) + E[∥D t,k-1 x -D x (z t,k-1 , v t,k-1 , x t,k-1 )∥ 2 ] .
Using Lemma A.2, we get

E[∥D t,k-1 x -D x (u t,k-1 )∥ 2 ] ≤ 4ρ 2 (L G 2 R) 2 + (L F 1 ) 2 k-1 r=1 E[∥D t,r-1 z ∥ 2 ] + 4ρ 2 (L G 1 ) 2 k-1 r=1 E[∥G t,r-1 v ∥ 2 ] + 4γ 2 (L G 2 R) 2 + (L F 1 ) 2 k-1 r=1 E[∥D t,r-1 x ∥ 2 ] .
Putting all together yields

ϕ t,k+1 z ≤ 1 - µ G 2 ρ ϕ t,k - ρ 2 (1 -Λ z ρ) E[∥D t,k z ∥ 2 ] + Λ z 2 γ 2 E[∥D t,k x ∥ 2 ] (22) + L ′ µ G γ 2 ρ E[∥D t,k x (u t,k )∥ 2 ] + 4(L G 1 ) 2 L ′ µ G γ 2 ρ k r=1 E[∥G t,r-1 v ∥ 2 ] + ρ ρ 2 L G 1 2 + 4(L G 2 R) 2 L ′ µ G γ 2 + 4(L F 1 ) 2 L ′ µ G γ 2 k r=1 E[∥D t,r-1 z ∥ 2 ] + γ 2 ρ L G 1 2 + 4(L G 2 R) 2 L ′ µ G γ 2 ρ + 4(L F 1 ) 2 L ′ µ G γ 2 ρ k r=1 E[∥D t,r-1 x ∥ 2 ] By assumption, γ ≤ C z ρ, with C z = µ G L G 1 8L ′ ((L G 2 R) 2 +(L F ϕ t,k+1 z ≤ 1 - µ G 2 ρ ϕ t,k z - ρ 2 (1 -Λ z ρ) E[∥D t,k z ∥ 2 ] + Λ z 2 γ 2 E[∥D t,k x ∥ 2 ] + L ′ µ G γ 2 ρ E[∥D t,k x (u t,k )∥ 2 ] + ρ 3 L G 1 k r=1 E[∥D t,r-1 z ∥ 2 ] + 4(L G 1 ) 2 L ′ µ G γ 2 ρ k r=1 E[∥G t,r-1 v ∥ 2 ] + γ 2 ρL G 1 k r=1 E[∥D t,r-1 x ∥ 2 ] ≤ 1 - µ G 2 ρ ϕ t,k - ρ 2 (1 -Λ z ρ) V t,k z + Λ z 2 γ 2 V t,k x + β zx γ 2 ρ E[∥D t,k x (u t,k )∥ 2 ] + ρ 3 β zz V t,k z + γ 2 ρβ zv V t,k v + γ 2 ρβ zx V t,k x with β zz = L G 1 , β zv = 4(L G 1 ) 2 L ′ µ G , β zx = L G 1 and β zx = L ′ µ G .
A.6 Proof of Lemma 3.4

Recall that we denote

Ψ(z, v, x) = 1 2 v ⊤ ∇ 2 11 G(z, x)v + ∇ 1 F (z, x) ⊤ v and ϕ v (v, x) = Ψ(z * (x), v, x) - Ψ(z * (x), v * (x), x).
As for Lemma 3.3, the key property we need is the smoothness of ϕ v . The derivatives of ϕ v involve the third derivative of G. For a tensor T ∈ R p1×p2×p3 and a vector a ∈ R p3 we denote (T |a) the matrix in R p1×p2 defined by:

(T |a) = p3 k=1 T i,j,k a k 1≤i≤p1 1≤j≤p2 . Lemma A.4. The function ϕ v has Λ v -Lipschitz continuous gradient on Γ × R d , for some constant Λ v . Proof. For any (v, x) ∈ Γ × R d , we have ∇ 1 ϕ v (v, x) = D v (z * (x), v, x)
and

∇ 2 ϕ v (v, x) = (dz * (x)) ⊤ 1 2 (∇ 3 111 G(z * (x), x)|v)v - 1 2 (∇ 3 111 G(z * (x), x)|v * (x))v * (x) +∇ 2 11 F (z * (x), x)v -∇ 2 11 F (z * (x), x)v * (x) + 1 2 (∇ 3 211 G(z * (x), x)|v)v - 1 2 (∇ 3 211 G(z * (x), x)|v * (x))v * (x) +∇ 2 21 F (z * (x), x)v -∇ 2 21 F (z * (x), x)v * (x) . Let us consider (v, x) ∈ Γ × R d and (v ′ , x ′ ) ∈ Γ × R d . We have ∥∇ 1 ϕ v (v, x) -∇ 1 ϕ v (v ′ , x ′ )∥ ≤ ∥∇ 2 11 G(z * (x), x)v -∇ 2 11 G(z * (x ′ ), x ′ )v ′ ∥ + ∥∇ 1 F (z * (x), x) -∇ 1 F (z * (x ′ ), x ′ )∥ For the first term, ∥∇ 2 11 G(z * (x), x)v -∇ 2 11 G(z * (x ′ ), x ′ )v ′ ∥ ≤ ∥∇ 2 11 G(z * (x), x)(v -v ′ )∥ + ∥(∇ 2 11 G(z * (x), x) -∇ 2 11 G(z * (x ′ ), x ′ ))v ′ ∥ + ∥∇ 2 11 G(z * (x ′ ), x ′ )(v -v ′ )∥ ≤ 2L G 1 ∥v -v ′ ∥ + L G 2 (L z * + 1)∥v ′ ∥∥x -x ′ ∥ ≤ [2L G 1 + L G 2 (L z * + 1)R]∥(v, x) -(v ′ , x ′ )∥
For the second terms, we use the smoothness of F and the Lipschitz continuity of z * (Lemma A.1):

∥∇ 1 F (z * (x), x) -∇ 1 F (z * (x ′ ), x ′ )∥ ≤ L F 1 ∥(z * (x), x) -(z * (x ′ ), x ′ )∥ ≤ L F 1 (∥z * (x) -z * (x ′ )∥ + ∥x -x ′ ∥) ≤ L F 1 (L z * + 1)∥x -x ′ ∥ ≤ L F 1 (L z * + 1)∥(x, v) -(x ′ , v ′ )∥ . As a consequence ∥∇ 1 ϕ v (v, x) -∇ 1 ϕ v (v ′ , x ′ )∥ ≤ Λ 1 ∥(v, x) -(v ′ , x ′ )∥ (23)
with

Λ 1 = L F 1 (L z * + 1) + 2L G 1 + L G 2 (L z * + 1)R . (24)
To prove the Lipschitz continuity of ∇ 2 ϕ v , we remark that ∇ 3 111 G, ∇ 3 211 G are Lipschitz and bounded by assumption. (v → v) is Lipschitz and bounded on Γ. Also by Lemma A.1, z * and v * are Lipschitz and bounded. Finally, dz * is bounded (Lemma A.1) and Lipschitz according to Chen et al. ( 2021

)[Lemma 9]. As a consequence, ∇ 2 ϕ v is Λ 2 -Lpischitz for some constant Λ 2 > 0. Hence, ∇ϕ v is Λ v -Lipschitz continuous with Λ v = Λ 1 + Λ 2 . Lemma A.5. Let t > 0. For k ∈ [q -1], we have 0 ≤ - 1 ρ (v t,k+1 -v t,k ) + D t,k v , v t,k+1 -v t,k
Proof. The function ι Γ being convex (since Γ is convex), let us consider its sub-differential

∂ι γ (v) = {η ∈ R p , ∀v ′ ∈ R p , ι Γ (v ′ ) ≥ ι Γ (v) + ⟨η, v ′ -v⟩} By definition v t,k+1 = arg min v (ι Γ (v) + 1 2ρ ∥v -(v t,k -ρD t,k v )∥ 2 ) .
Using Fermat's rule, we get

- 1 ρ (v t,k+1 -v t,k ) -D t,k v ∈ ∂ι Γ (v t,k+1 ) .
We can use the definition of the sub-differential with η

= -1 ρ (v t,k+1 -v t,k ) -D t,k v to get ι Γ (v t,k+1 ) =0 ≤ ι Γ (v t,k ) =0 - 1 ρ (v t,k+1 -v t,k ) + D t,k v , v t,k+1 -v t,k .
We can now turn to the proof of Lemma 3.4.

Proof. The smoothness of ϕ v provides us the following upper bound

ϕ v (v t,k+1 , x t,k+1 ) ≤ ϕ v (v t,k , x t,k ) + ⟨Π Γ (v t,k -ρD t,k v ) -v t,k , D v (z * (x t,k ), v t,k , x t,k )⟩ (25) + Λ v 2 ρ 2 ∥Π Γ (v t,k -ρD t,k v ) -v t,k ∥ 2 -γ⟨D t,k x , ∇ 2 ϕ v (v t,k , x t,k )⟩ + Λ v 2 γ 2 ∥D t,k x ∥ 2 . Let us denote ∆ t,k Π = Π Γ (v t,k -ρD t,k v ) -Π Γ (v t,k -ρD v (z * (x t,k ), v t,k , x t,k )). Adding and subtracting ⟨Π Γ (v t,k -ρD v (z * (x t,k ), v t,k , x t,k ) -v t,k , D v (z * (x t,k ), v t,k , x t,k )⟩ + Λ v 2 ∥Π Γ (v t,k -ρD v (z * (x t,k ), v t,k , x t,k ) -v t,k ∥ 2 yields ϕ v (v t,k+1 , x t,k+1 ) ≤ ϕ v (v t,k , x t,k ) + ⟨∆ t,k Π , D v (z * (x t,k ), v t,k , x t,k )⟩ (26) + Λ v 2 ∥Π Γ (v t,k -ρD v (z * (x t ), v t,k , x t,k )) -v t,k ∥ 2 + ⟨Π Γ (v t,k -ρD v (z * (x t,k ), v t,k , x t,k )) -v t,k , D v (z * (x t,k ), v t,k , x t,k )⟩ + Λ v 2 ∥∆ t,k Π ∥ 2 + Λ v ⟨∆ t,k Π , Π Γ (v t,k -ρD v (z * (x t,k ), v t,k , x t,k )) -v t,k ⟩ -γ⟨D t,k x , ∇ 2 ϕ v (v t,k , x t,k )⟩ + Λ v 2 γ 2 ∥D t,k x ∥ 2 . Taking ρ ≤ 1 Γ v gives ϕ v (v t,k+1 , x t,k+1 ) ≤ ϕ v (v t,k , x t,k ) + ⟨∆ t,k Π , D v (z * (x t,k ), v t,k , x t,k )⟩ (27) + 1 2ρ ∥Π Γ (v t,k -ρD v (z * (x t ), v t,k , x t,k )) -v t,k ∥ 2 + ⟨Π Γ (v t,k -ρD v (z * (x t,k ), v t,k , x t,k )) -v t,k , D v (z * (x t,k ), v t,k , x t,k )⟩ + Λ v 2 ∥∆ t,k Π ∥ 2 + Λ v ⟨∆ t,k Π , Π Γ (v t,k -ρD v (z * (x t,k ), v t,k , x t,k )) -v t,k ⟩ -γ⟨D t,k x , ∇ 2 ϕ v (v t,k , x t,k )⟩ + Λ v 2 γ 2 ∥D t,k x ∥ 2 .
Let ι Γ the indicator function of the convex set Γ. Similarly to Karimi et al. (2016, Equation 13) we define for any α > 0 and v ∈ R p

D ιΓ (v, x, α) = -2α min v ′ ∈R p ⟨∇ 1 ϕ v (v, x), v ′ -v⟩ + α 2 ∥v ′ -v∥ 2 + ι Γ (v ′ ) -ι Γ (v) .
Hence, for v ∈ Γ and x ∈ R d , we have

- ρ 2 D ιΓ v, x, 1 ρ = ⟨Π Γ (v -ρD v (z * (x), v, x)) -v, D v (z * (x), v, x)⟩ + 1 2ρ ∥Π Γ (v -ρD v (z * (x), v, x)) -v∥ 2 .
Therefore, Equation ( 27) can be written as

ϕ v (v t,k+1 , x t,k+1 ) ≤ ϕ v (v t,k , x t,k ) - ρ 2 D ιΓ v t,k , x t,k , 1 ρ + ⟨∆ t,k Π , D v (z * (x t,k ), v t,k , x t,k )⟩ + Λ v 2 ∥∆ t,k Π ∥ 2 + Λ v ⟨∆ t,k Π , Π Γ (v t,k -ρD v (z * (x t,k ), v t,k , x t,k )) -v t,k ⟩ -γ⟨D t,k x , ∇ 2 ϕ v (v t,k , x t,k )⟩ + Λ v 2 γ 2 ∥D t,k x ∥ 2 .
By strong convexity of ϕ v with respect top v and smoothness, we have

D ιΓ (v t,k , x t,k , Λ v ) ≥ 2µ G ϕ v (v t,k , x t,k ). According to Karimi et al. (2016, Lemma 1), D ιΓ (v t,k , x t,k , •) is an increasing function. As a consequence, since Λ v ≤ 1 ρ , we have D ιΓ v t,k , x t,k , 1 ρ ≥ 2µ G ϕ v (v t,k , x t,k ). This leads to ϕ v (v t,k+1 , x t,k+1 ) ≤ (1 -ρµ G )ϕ v (v t,k , x t,k ) + ⟨∆ t,k Π , D v (z * (x t,k ), v t,k , x t,k )⟩ (28) + Λ v 2 ∥∆ t,k Π ∥ 2 + Λ v ⟨∆ t,k Π , Π Γ (v t,k -ρD v (z * (x t,k ), v t,k , x t,k )) -v t,k ⟩ -γ⟨D t,k x , ∇ 2 ϕ v (v t,k , x t,k )⟩ + Λ v 2 γ 2 ∥D t,k x ∥ 2 .
The non-expansiveness of Π Γ yields

∥∆ t,k Π ∥ ≤ ρ∥D t,k v -D v (z * (x t,k ), v t,k , x t,k )∥ (29)
and

∥Π Γ (v t,k -ρD v (z * (x t,k ), v t,k , x t,k ) -v t,k ∈Γ ∥ = ∥Π Γ (v t,k -ρD v (z * (x t,k ), v t,k , x t,k )) -Π Γ (v t,k )∥ ≤ ρ∥D v (z * (x t,k ), v t,k , x t,k )∥ . (30)
Moreover, using Equation ( 29) and Young Inequality, we have for any c > 0

⟨∆ t,k Π , D v (z * (x t,k ), v t,k , x t,k )⟩ ≤ c 2 ∥∆ Π ∥ 2 + 1 2c ∥D v (z * (x t,k ), v t,k , x tk, )∥ 2 ≤ cρ 2 2 ∥D t,k v -D v (z * (x t,k ), v t,k , x t,k )∥ 2 + 1 2c ∥D v (z * (x t,k ), v t,k , x tk, ) -D v (z * (x t,k ), v * (x t,k ), x tk, ) =0 ∥ 2 ≤ cρ 2 2 ∥D t,k v -D v (z * (x t,k ), v t,k , x t,k )∥ 2 (31) + L G 1 µ G c ϕ v (v t,k , x t,k )
Plugging Equation (31) into Equation ( 28) with c =

2L G 1 µ 2 G ρ yields ϕ v (v t,k+1 , x t,k+1 ) ≤ 1 - ρµ G 2 ϕ v (v t,k , x t,k ) + L G 1 ρ µ 2 G ∥D t,k v -D v (z * (x t,k ), v t,k , x t,k )∥ 2 (32)
+ Λ v 2 ∥∆ t,k Π ∥ 2 + Λ v ⟨∆ t,k Π , Π Γ (v t,k -ρD v (z * (x t,k ), v t,k , x t,k )) -v t,k ⟩ -γ⟨D t,k x , ∇ 2 ϕ v (v t,k , x t,k )⟩ + Λ v 2 γ 2 ∥D t,k x ∥ 2 .
Using Equation (29), Equation (30) and Young Inequality for d > 0 yields

⟨∆ t,k Π , Π Γ (v t,k -ρD v (z * (x t,k ), v t,k , x t,k )) -v t,k ⟩ ≤ d 2 ∥∆ t,k Π ∥ 2 + 1 2d ∥Π Γ (v t,k -ρD v (z * (x t,k ), v t,k , x t,k )) -v t,k ∥ 2 ≤ dρ 2 2 ∥D t,k v -D v (z * (x t,k ), v t,k , x t,k )∥ 2 (33) + ρ 2 2d ∥D v (z * (x t,k ), v t,k , x t,k )∥ 2 ≤ dρ 2 2 ∥D t,k v -D v (z * (x t,k ), v t,k , x t,k )∥ 2 (34) + L G 1 ρ 2 µ G d ϕ v (v t,k , x t,k ) .
Plugging Equation (34) into Equation ( 32) with d =

4L G 1 Λv)ρ µ 2 G gives ϕ v (v t,k+1 , x t,k+1 ) ≤ 1 - ρµ G 4 ϕ v (v t,k , x t,k ) (35)
+ L G 1 ρ µ 2 G + 2L G 1 Λ 2 v ρ 3 µ 2 G ∥D t,k v -D v (z * (x t,k ), v t,k , x t,k )∥ 2 + Λ v 2 ∥∆ t,k Π ∥ 2 -γ⟨D t,k x , ∇ 2 ϕ v (v t,k , x t,k )⟩ + Λ v 2 γ 2 ∥D t,k x ∥ 2 .
Using once again (29), we get

ϕ v (v t,k+1 , x t,k+1 ) ≤ 1 - ρµ G 4 ϕ v (v t,k , x t,k ) (36) + L G 1 ρ µ 2 G + 2L G 1 Λ 2 v ρ 3 µ 2 G + Λ v ρ 2 2 ∥D t,k v -D v (z * (x t,k ), v t,k , x t,k )∥ 2 -γ⟨D t,k x , ∇ 2 ϕ v (v t,k , x t,k )⟩ + Λ v 2 γ 2 ∥D t,k x ∥ 2 .
By Lemma A.5, we have for any α > 0

0 ≤ -α 1 ρ (v t,k+1 -v t,k ) + D t,k v , v t,k+1 -v t,k .
By adding this to Equation (36), we get

ϕ v (v t,k+1 , x t,k+1 ) ≤ 1 - ρµ G 4 ϕ v (v t,k , x t,k ) (37) - α ρ ∥v t,k+1 -v t,k ∥ 2 -α⟨D t,k v , v t,k+1 -v t,k ⟩ + L G 1 ρ µ 2 G + 2L G 1 Λ 2 v ρ 3 µ 2 G + Λ v ρ 2 2 ∥D t,k v -D v (z * (x t,k ), v t,k , x t,k )∥ 2 -γ⟨D t,k x , ∇ 2 ϕ v (v t,k , x t,k )⟩ + Λ v 2 γ 2 ∥D t,k x ∥ 2 .
We can control -D t,k v , v t,k+1 -v t,k by Cauchy-Schwarz and Young for some c, d, e, f > 0

-D t,k v , v t,k+1 -v t,k = -D v (z * (x t,k ), v t,k , x t,k ), Π Γ (v t,k -ρD v (z * (x t,k ), v t,k , x t,k ) -v t,k -D v (z * (x t,k ), v t,k , x t,k ), ∆ t,k Π -D t,k v -D v (z * (x t,k ), v t,k , x t,k ), Π Γ (v t,k -ρD v (z * (x))) -v t,k -D t,k v -D v (z * (x t,k ), v t,k , x t,k ), ∆ t,k Π ≤ c 2 ∥D v (z * (x t,k ), v t,k , x t,k )∥ 2 + 1 2c ∥Π Γ (v t,k -ρD v (z * (x t,k ), v t,k , x t,k )) -v t,k ∥ 2 + d 2 ∥D v (z * (x t,k ), v t,k , x t,k )∥ 2 + 1 2d ∥∆ t,k Π ∥ 2 + e 2 ∥D t,k v -D v (z * (x t,k ), v t,k , x t,k )∥ 2 + 1 2e ∥Π Γ (v t,k -ρD v (z * (x t,k ), v t,k , x t,k )) -v t,k ∥ 2 + f 2 ∥D t,k v -D v (z * (x t,k ), v t,k , x t,k )∥ 2 + 1 2f ∥∆ t,k Π ∥ 2 ≤ c + d 2 + ρ 2 1 2c + 1 2e ∥D v (z * (x t,k ), v t,k , x t,k )∥ 2 + e + f 2 + ρ 2 1 2d + 1 2f ∥D t,k v -D v (z * (x t,k ), v t,k , x t,k )∥ 2 ≤ c + d 2 + ρ 2 1 2c + 1 2e 2L G 1 µ G ϕ v (v t,k , x t,k ) + e + f 2 + ρ 2 1 2d + 1 2f ∥D t,k v -D v (z * (x t,k ), v t,k , x t,k )∥ 2
Let us take c = d = e = f = ρ. We get

-D t,k v , v t,k+1 -v t,k ≤ 4L G 1 µ G ρϕ v (v t,k , x t,k ) + 2ρ∥D t,k v -D v (z * (x t,k ), v t,k , x t,k )∥ 2 . ( 38
)
Then, by plugging the last Inequality in Equation ( 37) and setting α =

µ 2 G 32L G 1 , we end up with ϕ v (v t,k+1 , x t,k+1 ) ≤ 1 - µ G 8 ρ ϕ v (v t,k , x t,k ) - α ρ ∥v t,k+1 -v t,k ∥ 2 + ρ L G 1 µ 2 G + µ 2 G 16L G 1 + Λ v ρ 2 + 2L G 1 Λ 2 v ρ 2 µ 2 G ∥D t,k v -D v (z * (x t,k ), v t,k , x t,k )∥ 2 -γ⟨D t,k x , ∇ 2 ϕ v (v t,k , x t,k )⟩ + Λ v 2 γ 2 ∥D t,k x ∥ 2 ≤ 1 - µ G 8 ρ ϕ v (v t,k , x t,k ) - µ 2 G 32L G 1 ρ∥G t,k v ∥ 2 + ρ L G 1 µ 2 G + µ 2 G 16L G 1 + Λ v ρ 2 + 2L G 1 Λ 2 v ρ 2 µ 2 G ∥D t,k v -D v (z * (x t,k ), v t,k , x t,k )∥ 2 -γ⟨D t,k x , ∇ 2 ϕ v (v t,k , x t,k )⟩ + Λ v 2 γ 2 ∥D t,k x ∥ 2 . Since ρ ≤ B v ≜ L G 1 µ 2 G + µ 2 G 16L G 1 min 2 Λv , µ G √ 2L G 1 Λv yields ϕ v (v t,k+1 , x t,k+1 ) ≤ 1 - µ G 8 ρ ϕ v (v t,k , x t,k ) - µ 2 G 32L G 1 ρ∥G t,k v ∥ 2 + 3ρ L G 1 µ 2 G + µ 2 G 16L G 1 ∥D t,k v -D v (z * (x t,k ), v t,k , x t,k )∥ 2 -γ⟨D t,k x , ∇ 2 ϕ v (v t,k , x t,k )⟩ + Λ v 2 γ 2 ∥D t,k x ∥ 2 ≤ 1 - µ G 8 ρ ϕ v (v t,k , x t,k ) - µ 2 G 32L G 1 ρ∥G t,k v ∥ 2 (39) + 6ρ L G 1 µ 2 G + µ 2 G 16L G 1 ∥D t,k v -D v (u t,k )∥ 2 + 6ρ L G 1 µ 2 G + µ 2 G 16L G 1 ∥D v (u t,k ) -D v (z * (x t,k ), v t,k , x t,k )∥ 2 -γ⟨D t,k x , ∇ 2 ϕ v (v t,k , x t,k )⟩ + Λ v 2 γ 2 ∥D t,k x ∥ 2 .
Tacking the expectation conditionally to the past iterates yields

E t,k [ϕ v (v t,k+1 , x t,k+1 )] ≤ 1 - µ G 8 ρ ϕ v (v t,k , x t,k ) - µ 2 G 32L G 1 ρE t,k [∥G t,k v ∥ 2 ] (40) + 6ρ L G 1 µ 2 G + µ 2 G 16L G 1 E t,k [∥D t,k v -D v (u t,k )∥ 2 ] + 6ρ L G 1 µ 2 G + µ 2 G 16L G 1 E t,k [∥D v (u t,k ) -D v (z * (x t,k ), v t,k , x t,k )∥ 2 ] -γ⟨E t,k [D t,k x ], ∇ 2 ϕ v (v t,k , x t,k )⟩ + Λ v 2 γ 2 E t,k [∥D t,k x ∥ 2 ] .
From Young inequality, we have for any c > 0

⟨E t,k [D t,k x ], ∇ 2 ϕ v (v t,k , x t,k )⟩ ≤ c -1 ∥E t,k [D t,k x ]∥ 2 + c∥∇ 2 ϕ v (v t,k , x t,k )∥ 2 . ( 41
)
Moreover, using the Lipschitz continuity of z * , of ∇ 2 11 G and ∇ F and the fact that v and v * are bounded, we have

∥∇ 2 ϕ v (v, x)∥ ≤ ∥dz * (x)∥ 1 2 (∇ 3 111 G(z * (x), x)|v)v - 1 2 (∇ 3 111 G(z * (x), x)|v * (x))v * (x) +∥∇ 2 11 F (z * (x), x)v -∇ 2 11 F (z * (x), x)v * (x)∥ + ∥ 1 2 (∇ 3 211 G(z * (x), x)|v)v - 1 2 (∇ 3 211 G(z * (x), x)|v * (x))v * (x)∥ + ∥∇ 2 21 F (z * (x), x)v -∇ 2 21 F (z * (x), x)v * (x)∥ ≤ L * 1 2 (∇ 3 111 G(z * (x), x)|v -v * (x))v - 1 2 (∇ 3 111 G(z * (x), x)|v * (x))(v -v * (x)) +L F 2 ∥v -v * (x)∥ + 1 2 (∇ 3 211 G(z * (x), x)|v -v * (x))v - 1 2 (∇ 3 211 G(z * (x), x)|v * (x))(v -v * (x)) + L F 2 ∥v -v * (x)∥ ≤ L * 1 2 (∇ 3 111 G(z * (x), x)|v -v * (x))v + 1 2 (∇ 3 111 G(z * (x), x)|v * (x))(v -v * (x)) + L F 2 ∥v -v * (x)∥ + 1 2 (∇ 3 211 G(z * (x), x)|v -v * (x))v + 1 2 (∇ 3 211 G(z * (x), x)|v * (x))(v -v * (x)) + L F 2 ∥v -v * (x)∥ ≤ L * L G 2 2 (∥v∥ + ∥v * (x)∥)∥v -v * (x))∥ + L F 2 ∥v -v * (x)∥ + L G 2 2 (∥v∥ + ∥v * (x)∥)∥v -v * (x))∥ + L F 2 ∥v -v * (x)∥ ≤ (L * + 1) L G 2 R + L F 2 ∥v -v * (x)∥ .
On the other hand, we have by strong convexity

∥v -v * (x)∥ 2 ≤ 2 µ G ϕ v (v, x) .
As a consequence, we have

∥∇ 2 ϕ v (v t,k , x t,k )∥ 2 ≤ L ′′ ϕ v (v t,k , x t,k ) (42)
with

L ′′ = 2(L * +1) 2 [L G 2 R+L F 2 ] 2 µ G .
Plugging Inequalities ( 41) and ( 42) into (40) yields

E t,k [ϕ v (v t,k+1 , x t,k+1 )] ≤ 1 - µ G 8 ρ + cL ′′ γ ϕ v (v t,k , x t,k ) - µ 2 G 32L G 1 ρE t,k [∥G t,k v ∥ 2 ] + 6ρ L G 1 µ 2 G + µ 2 G 16L G 1 E t,k [∥D t,k v -D v (u t,k )∥ 2 ] + 6ρ L G 1 µ 2 G + µ 2 G 16L G 1 E t,k [∥D v (u t,k ) -D v (z * (x t,k ), v t,k , x t,k )∥ 2 ] + γ c ∥E t,k [D t,k x ]∥ 2 + Λ v 2 γ 2 E t,k [∥D t,k x ∥ 2 ] .
The Lipschitz continuity of ∇ 2 11 G and ∇ 1 F and the boundedness of v give us

∥D v (u t,k ) -D v (z * (x t,k ), v t,k , x t,k )∥ 2 ≤ ∥∇ 2 11 G(z t,k , x t,k )v t,k -∇ 2 11 G(z * (x t,k ), x t,k )v t,k ∥ +∥∇ 1 F (z t,k , x t,k ) -∇ 1 F (z * (x t,k ), x t,k )∥ 2 ≤ (L G 2 R + L F 1 ) 2 ∥z t,k -z * (x t,k )∥ 2 ≤ 2(L G 2 R + L F 1 ) 2 µ G ϕ z (z t,k , x t,k ) . As a consequence E t,k [ϕ v (v t,k+1 , x t,k+1 )] ≤ 1 - µ G 8 ρ + cL ′′ γ ϕ v (v t,k , x t,k ) - µ 2 G 32L G 1 ρE t,k [∥G t,k v ∥ 2 ] (43) + 6ρ L G 1 µ 2 G + µ 2 G 16L G 1 E t,k [∥D t,k v -D v (u t,k )∥ 2 ] + 6ρ L G 1 µ 2 G + µ 2 G 16L G 1 2(L G 2 R + L F 1 ) 2 µ G ϕ z (z t,k , x t,k ) + γ c ∥E t,k [D t,k x ]∥ 2 + Λ v 2 γ 2 E t,k [∥D t,k x ∥ 2 ] .
From Lemma A.2, we have

E[∥D t,k v -D v (u t,k )∥ 2 ] ≤ 4ρ 2 (L G 2 R) 2 + (L F 1 ) 2 k r=1 E[∥D t,r-1 z ∥ 2 ] + 4ρ 2 (L G 1 ) 2 k r=1 E[∥G t,r-1 v ∥ 2 ] + 4γ 2 (L G 2 R) 2 + (L F 1 ) 2 k r=1 E[∥D t,r-1 x ∥ 2 ]
Taking the total expectation and plugging the previous inequality in Equation ( 43) yields

ϕ t,k+1 v ≤ 1 - µ G 8 ρ + cL ′′ γ ϕ t,k v - µ 2 G 32L G 1 ρE t,k [∥G t,k v ∥ 2 ] + 24ρ 3 (L G 2 R) 2 + (L F 1 ) 2 L G 1 µ 2 G + µ 2 G 16L G 1 k r=1 E[∥D t,r-1 z ∥ 2 ] + 24ρ 3 (L G 1 ) 2 L G 1 µ 2 G + µ 2 G 16L G 1 k r=1 E[∥G t,r-1 v ∥ 2 ] + 24ργ 2 (L G 2 R) 2 + (L F 1 ) 2 L G 1 µ 2 G + µ 2 G 16L G 1 k r=1 E[∥D t,r-1 x ∥ 2 ] + L G 1 µ 2 G + µ 2 G 16L G 1 12(L G 2 R + L F 1 ) 2 µ G ρϕ t,k z + γ c E[∥[E t,k D t,k x ]∥ 2 ] + Λ v 2 γ 2 E[∥D t,k x ∥ 2 ] . Taking c = µ G ρ 16L ′′ γ yields ϕ t,k+1 v ≤ 1 - µ G 16 ρ ϕ t,k v - µ 2 G 32L G 1 ρE t,k [∥G t,k v ∥ 2 ] + 24ρ 3 (L G 2 R) 2 + (L F 1 ) 2 L G 1 µ 2 G + µ 2 G 16L G 1 k r=1 E[∥D t,r-1 z ∥ 2 ] + 24ρ 3 (L G 1 ) 2 L G 1 µ 2 G + µ 2 G 16L G 1 k r=1 E[∥G t,r-1 v ∥ 2 ] + 24ργ 2 (L G 2 R) 2 + (L F 1 ) 2 L G 1 µ 2 G + µ 2 G 16L G 1 k r=1 E[∥D t,r-1 x ∥ 2 ] + L G 1 µ 2 G + µ 2 G 16L G 1 12(L G 2 R + L F 1 ) 2 µ G ρϕ t,k z + 16L ′′ µ G γ 2 ρ E[∥[E t,k D t,k x ]∥ 2 ] + Λ v 2 γ 2 E[∥D t,k x ∥ 2 ] .
Then, we use the identity ⟨a, b⟩ = 1 2 (∥a∥ 2 + ∥b∥ 2 -∥a -b∥ 2 ) to get

h(x t,k+1 ) ≤ h(x t,k ) - γ 2 ∥∇h(x t,k )∥ 2 - γ 2 ∥D t,k x ∥ 2 + γ 2 ∥∇h(x t,k ) -D t,k x ∥ 2 + γ 2 L h 2 ∥D t,k x ∥ 2 ≤ h(x t,k ) - γ 2 ∥∇h(x t,k )∥ 2 - γ 2 ∥D t,k x ∥ 2 + γ∥∇h(x t,k ) -D x (u t,k )∥ 2 + γ∥D x (u t,k ) -D t,k x ∥ 2 + γ 2 L h 2 ∥D t,k x ∥ 2 .
Then taking the expectation gives and using Proposition 2.5 yields

h t,k+1 ≤ h t,k - γ 2 g t,k + γE[∥∇h(x t,k ) -D x (u t,k )∥ 2 ] + γE[∥D x (u t,k ) -D t,k x ∥ 2 ] - γ 2 1 -L h γ E[∥D t,k x ∥ 2 ] ≤ h t,k - γ 2 g t,k + γL 2 x (E[∥z t,k -z * (x t,k )∥ 2 ] + E[∥v t,k -v * (x t,k )∥ 2 ]) + γE[∥D x (u t,k ) -D t,k x ∥ 2 ] - γ 2 1 -L h γ E[∥D t,k x ∥ 2 ] . The µ G -strong convexity of G( . , x) ensures that ∥z -z * (x)∥ 2 ≤ 2 µ G ϕ z (z, x) and ∥v -v * (x)∥ 2 ≤ 2 µ G ϕ v (v, x). As a consequence h t,k+1 ≤ h t,k - γ 2 g t,k + γ 2L 2 x µ G (ϕ t,k z + ϕ t,k v ) + γE[∥D x (z t,k , v t,k , x t,k ) -D t,k x ∥ 2 ] - γ 2 1 -L h γ E[∥D t,k x ∥ 2 ] .
From Lemma A.2, we have

E[∥D t,k x -D x (u t,k )∥ 2 ] ≤ 4ρ 2 (L G 2 R) 2 + (L F 1 ) 2 k r=1 E[∥D t,r-1 z ∥ 2 ] + 4ρ 2 (L G 1 ) 2 k r=1 E[∥G t,r-1 v ∥ 2 ] + 4γ 2 (L G 2 R) 2 + (L F 1 ) 2 k r=1 E[∥D t,r-1 x ∥ 2 ] .
As a consequence

h t,k+1 ≤ h t,k - γ 2 g t,k + γ 2L 2 x µ G (ϕ t,k z + ϕ t,k v ) + 4γρ 2 (L G 2 R) 2 + 2(L F 1 ) 2 k r=1 E[∥D t,r-1 z ∥ 2 ] + 4γρ 2 (L G 1 ) 2 k r=1 E[∥G t,r-1 v ∥ 2 ] + 4γ 3 (L G 2 R) 2 + 2(L F 1 ) 2 k r=1 E[∥D t,r-1 x ∥ 2 ] - γ 2 1 -L h γ E[∥D t,k x ∥ 2 ] ≤ h t,k - γ 2 g t,k + γ 2L 2 x µ G (ϕ t,k z + ϕ t,k v ) + γρ 2 β hz V t,k z + γρ 2 β hv V t,k v + γ 3 β hx V t,k x - γ 2 1 -L h γ E[∥D t,k x ∥ 2 ] with β hz = 4 (L G 2 R) 2 + 2(L F 1 ) 2 , β hv = 4(L G 1 ) 2 and β hx = 4 (L G 2 R) 2 + 2(L F 1 ) 2 .
A.8 Proof of Theorem 1 and Corollary 3.6

The constants involved in Theorem 1 are

ψ z = 1 16β zx , ψ v = min 1 16β vx , α zv µ G 12 ψ z ρ = min ψ z 12q(ψ z β zz + ψ v β zv ) , 1 6Λ z , 1 12qβ vv , sqrt βvv 3Λ v , B v , ξ = min C z , C v , 1, ψ v µ 2 G 16L 2 x , µ G 8β vx , ψ z µ 2 G 24L 2 x , µ G 12β zx , γ = min   1 12q(ψ z β zx + ψ v β vx ) , 1 12qβ hx , 1 6(L h + ψ z Λ z + ψ v Λ v ) , ψ v βvv 6q(β hv + ψ z β vz ) , ψ z 12qβ hz   .
Proof. The proof is a classical Lyapunov analysis. Consider the following Lyapunov function L

t,k = h t,k + ψ z ϕ t,k z + ψ v ϕ t,k
v for some positive constants ψ z and ψ v . We use use Lemmas 3.3 to 3.5 to upper bound L t,k -L t,k+1 . We have

L t,k+1 -L t,k ≤ - γ 2 g t,k + (ψ z β zx + ψ v β vx ) γ 2 ρ E[∥D x (u t,k )∥ 2 ] (44) + 2L 2 x µ G γ -ψ z µ G 2 ρ + ψ v α zv ρ ϕ t,k z + 2L 2 x µ G γ -ψ v µ G 16 ρ ϕ t,k v + ψ z Λ z 2 ρ 2 -ψ z 1 2 ρ V t,k z -ψ v βvv ρV t,k v + L h 2 γ 2 + ψ z Λ z 2 γ 2 + ψ v Λ v 2 γ 2 - γ 2 V t,k x + β hz ργ 2 + ψ z β zz ρ 3 + ψ v β zv ρ 3 V t,k z + β hv ργ 2 + ψ z β vz γ 2 ρ + ψ v β vv ρ 3 V t,k v + β hx γ 3 + ψ z β zx γ 2 ρ + ψ v β vx ρ 3 V t,k x .
We bound E[∥D x (u t,k )∥ 2 ] crudely by using Proposition 2.5

E[∥D x (u t,k )∥ 2 ] ≤ 2E[∥∇h(x t,k )∥ 2 ] + 2E[∥D x (u t,k ) -∇h(x t,k )∥ 2 ] ≤ 2g t,k + 2(E[∥z t,k -z * (x t,k )∥ 2 ] + E[∥v t,k -v * (x t,k )∥ 2 ]) ≤ 2g t,k + 4 µ G (ϕ t,k z + ϕ t,k v ) .
Summing in (44) for k = 0, . . . , q -1 yields

L t,q -L t,0 ≤ - γ 2 -2ψ z β zx γ 2 ρ -2ψ v β vx γ 2 ρ q-1 k=0 g t,k (45)
+ 2L 2 x µ G γ -ψ z µ G 2 ρ + ψ v α zv ρ + ψ z β zx γ 2 ρ q-1 k=0 ϕ t,k z + 2L 2 x µ G γ -ψ v µ G 16 ρ + ψ v β vx γ 2 ρ q-1 k=0 ϕ t,k v -ψ z βvv ρ q-1 k=0 V t,k z + ψ v Λ v 2 ρ 2 -ψ v 1 2 ρ q-1 k=0 V t,k v + L h 2 γ 2 + ψ z Λ z 2 γ 2 + ψ v Λ v 2 γ 2 - γ 2 q-1 k=0 V t,k x + β hz ργ 2 + ψ z β zz ρ 3 + ψ v β zv ρ 3 q-1 k=0 V t,k z + β hv ργ 2 + ψ z β vz γ 2 ρ + ψ v β vv ρ 3 q-1 k=0 V t,k v + β hx γ 3 + ψ z β zx γ 2 ρ + ψ v β vx ρ 3 q-1 k=0 V t,k
x .

Since we have

q-1 k=0 V t,k • = q-1 k=0 k r=1 E[∥D t,r-1 • ∥ 2 ] = q-1 r=1 q-1 k=r E[∥D t,r-1 • ∥ 2 ] = q-1 r=1 (q -r)E[∥D t,r-1 • ∥ 2 ] ≤ q q-1 k=1 E[∥D t,k-1 • ∥ 2 ] we get L t,q -L t,0 ≤ - γ 2 -2ψ z β zx γ 2 ρ -2ψ v β vx γ 2 ρ q-1 k=0 g t,k (46)
+ 2L 2 x µ G γ -ψ z µ G 2 ρ + ψ v α zv ρ + ψ z β zx γ 2 ρ q-1 k=0 ϕ t,k z + 2L 2 x µ G γ -ψ v µ G 2 ρ + ψ v β vx γ 2 ρ q-1 k=0 ϕ t,k v + ψ z Λ z 2 ρ 2 -ψ z 1 2 ρ + q β hz ργ 2 + ψ z β zz ρ 3 + ψ v β zv ρ 3 q-1 k=0 V t,k z + ψ v Λ v 2 ρ 2 -ψ v βvv ρ + q β hv ργ 2 + ψ z β vz γ 2 ρ + ψ v β vv ρ 3 q-1 k=0 V t,k v + L h 2 γ 2 + ψ z Λ z 2 γ 2 + ψ v Λ v 2 γ 2 - γ 2 + q β hx γ 3 + ψ z β zx γ 2 ρ + ψ v β vx ργ 2 q-1 k=0 V t,k x . Since ρ ≤ ρ ≤ min ψz 12q(ψzβzz+ψvβzv) , 1 6Λz and γ ≤ γ ≤ ψz 12qβ hz , we have ψ z Λ z 2 ρ 2 -ψ z 1 2 ρ + q β hz ργ 2 + ψ z β zz ρ 3 + ψ v β zv ρ 3 < 0 . (47)
Moreover, the conditions ρ ≤ ρ ≤ min βvv 6qβvv , βvv 3Λv and γ ≤ γ ≤ ψv βvv 6q(β hv +ψzβvz) , ensure that

ψ v Λ v 2 ρ 2 -ψ v βvv ρ + q β hv ργ 2 + ψ z β vz γ 2 ρ + ψ v β vv ρ 3 < 0 . (48)
The conditions ρ ≤ ρ ≤ 1 12q(ψzβzx+ψvβvx) and γ ≤ γ ≤ min 1 12q(ψzβzx+ψvβvx) ,

1 12qβ hx , 1 6(L h +ψzΛz+ψvΛv) yield L h 2 γ 2 + ψ z Λ z 2 γ 2 + ψ v Λ v 2 γ 2 - γ 2 + q β hx γ 3 + ψ z β zx γ 2 ρ + ψ v β vx ργ 2 < 0 . ( 49
)
The condition γ ≤ ξρ ≤ min

ψvµ 2 G 16L 2 x , µ G 8β vx ρ ensures 2L 2 x µ G γ -ψ v µ G 2 ρ + ψ v β vx γ 2 ρ ≤ 0 (50)
By definition, we have ψ v ≤ αzvµ G 12 ψ z and by assumptions γ ≤ ξρ ≤ min

ψzµ 2 G 24L 2 x , µ G 12β zx ρ. As a consequence 2L 2 x µ G γ -ψ z µ G 2 ρ + ψ v α zv ρ + ψ z β zx γ 2 ρ < 0 . (51)
Plugging Inequalities (47) to (51) into Equation ( 46) gives L t,q -L t,0 ≤ -γ 2 -2ψ z β zx γ 2 ρ -2ψ v β vx γ 2 ρ q-1 k=0 g t,k .

Since ψ z = ρ 16β zx and ψ v ≤ ρ 16β vx and γ 2 ρ ≤ ξ ≤ 1, we get L t,q -L t,0 L t+1,0 -L t,0 ≤ -γ 4

q-1 k=0 g t,k .

Summing, telescoping and dividing by T q gives 1 T q

T -1 t=0 q-1 k=0 g t,k ≤ 4 T qγ h 0,0 -h * + ψ z ϕ 0,0 + ψ v ϕ 0,0 ] Γ 0

.

From Theorem 1 we deduce Corollary 3.6.

Proof. Let us take ρ = ρ(n + m) -1 2 , γ = min(ξρ, γ) and q = n + m. Then Theorem 1 holds:

1 T q T -1 t=0 q-1 k=0 g t,k ≤ 4 T qγ Γ 0 .
with Γ 0 = O(1). To get an ε-stationary solution, we set T ≥ 4 qγ Γ 0 ε -1 ∨ 1 = O 1 qγε ∨ 1 . One iteration has Θ(q) = Θ(n + m) oracle complexity. As a consequence, the sample complexity to get an ε-stationary point is

O (n + m) 1 2 ϵ -1 ∨ (n + m) .

B Lower bound for bilevel problems (proof of Theorem 2)

The proof of Theorem 2 is an adaptation of the proof of (Zhou and Gu, 2019, Theorem 4.7) from single-level to bilevel problems. We build the outer function from the worst-case instance of (Zhou and Gu, 2019, Theorem 4.7) and we add a bilevel component by using as inner function the function G defined by G(z, x) = µ G 2 ∥z -x∥ 2 . We start by introducing the different tools used in this proof.

B.1 Preliminary results

In what follows, we provide the building blocks of our worst-case instance. The proof uses the following quadratic function presented by (Nesterov, 2018). Proof. Let x ∈ R q such that U x = U (≤k) x. For 0 ≥ k ≤ d, we denote The essential properties of f nc come from (Carmon et al., 2021, Lemmas 2, 3, 4). The first part provides the regularity properties of f nc . The second part bounds the distance between f nc ( • ; α, d) and the optimal value of the function. The third part will be key to the overall proof. In words, it states that as long x ∈ R d+1 has its two last components equal to zero, the norm of the gradient of f nc ( • ; α, d) is higher than a constant controlled by α.

R k,d = {v ∈ R d , v k+1 = • • • = v d = 0} . Let us write Q(x; ξ, d) = 1 2 x ⊤ Ax + b ⊤ x + c with A =          1 + ξ -1 0 • • • 0 -1 2
As a consequence, if α is properly chosen, as soon as x d = x d+1 = 0, we are ensured that ∥∇f nc (x; α, d)∥ ≥ ϵ.

Proposition B.6. For α ∈ [0 , 1], it holds 1. -αc ⪯ ∇ 2 f nc ⪯ 4 + αc.

B.2 Main proof

Now we are ready to prove Theorem 2.

Proof. We consider U ∈ O((T + 1)m, (T + 1)m) and we denote 1) . . .

U =    U (
U (m)   
with U (j) = (u (j) 1 , . . . , u

T +1 ) ⊤ ∈ O(T + 1, (T + 1)m). For j ∈ [m], we choose F j : R (T +1)m+(T +1)m → R defined by F j (z, x) = f nc (U (j) z; α, T ) and we set F = 1 m m j=1 F j . We also define for i ∈ [n] G i (z, x) = 1 2 ∥z -x∥ 2 , G = 1 n n i=1 G i , z * (x) = arg min z G(z, x) and h(x) = F (z * (x), x) = f nc (U (j) x; α, T ). By Proposition B.6, F j is 4 + αc m smooth, and G i is 1-smooth and 1-strongly convex.

We have h(0) -inf

x h(x) ≤ √ α + 10αT .

  Proceedings of the 27 th International Conference on Artificial Intelligence and Statistics (AISTATS) 2024, Valencia, Spain. PMLR: Volume 238. Copyright 2024 by the author(s).

  one variable minimizes another function. This can be formalized as followsmin x∈R d h(x) = F (z * (x), x) , subject to z * (x) ∈ arg min z∈R p G(z, x) .

  Definition 4.1. Let n, m two positive integers, L F 1 and µ G two positive constants. The class of the smooth empirical risk minimization problems denoted by C L F 1 ,µ G is the set of pairs of real-valued function families

Definition B. 1 .

 1 Let d ∈ N >0 , ξ ∈ [0 , +∞) and ζ ≤ 1. We define Q(.; ξ, d) : R d → R by Q(x; ξ, d) -x k ) 2 . Proposition B.2 proposition comes directly from (Zhou and Gu, 2019, Proposition 3.5). The first part of the proposition gives us the regularity of Q. In the second part shows that a function defined as Q(U × • ; ξ, d) + q p=1 g(⟨u p , • ⟩) verifies the so-called "zero-chain property" Carmon et al. (2020): if U x ∈ Span(u 1 , . . . , u k ), then we gain a non zero coordinate by calling the gradient ∇ [Q (U × • ; ξ, d) + q p=1 g(⟨u p , • ⟩)](x). In other words, that makes us progress in the problem resolution. Proposition B.2. For d ∈ N >0 , ξ ∈ [0 , +∞) and ζ ≤ 1. The following holds:1. Q( • ; ξ, d) is convex and 4-smooth. 2. Let q ∈ N >0 , U = [u 1 , . . . , u d ] ⊤ ∈ R d×q such that U U ⊤ = I and for k ≤ d, U (≤k) = [u 1 , . . . , u k , 0, . . . , 0] ⊤ ∈ R d×q .Let g : R → R differentiable such that g ′ (0) = 0. Then for any x ∈ R q such that U x = U(≤k) x, then∇ Q(U × • ; ξ, d) + d p=1g(⟨u p , • ⟩) (x) ∈ Span(u 1 , . . . , u k , u k+1 ) .

  0, . . . , 0) ⊤ and c = ξ 2 (1, 0, . . . , 0) ⊤ . On the one hand it is known from(Nesterov, 2018, Lemma 2.5.1) that if v ∈ R k,d , ∇Q(v; ξ, d) ∈ R k+1,d As a consequence, ∇Q(U x; ξ, d) = ∇Q(U (≤k) x ∈R k,d ; ξ, d) ∈ R k+1,d and ∇[Q(U × • ; ξ, d)](x) = U ⊤ ∇Q(U x; ξ, d) ∈ Span(u 1 , . . . , u k+1 ) . ⟨u p , x⟩)u p = k p=1 g ′ (⟨u p , x⟩)u p ∈ Span(u 1 , . . . , u k+1 ) . Thus ∇ Q(U × • ; ξ, d) + d p=1 g(⟨u p , • ⟩) (x) ∈ Span(u 1 , . . . , u k , u k+1 ) .However, the function Q is convex. That is why we also use the function Γ introduced in Carmon et al. (2021).As explained inCarmon et al. (2021), this function is essential to lower bound the gradient of our worst wase instance.Definition B.3. Let d ∈ N >0 . We define Γ( • ; d) : R d+1 → R by Γ(x; d) t -1) 1 + t 2 dt .An important property of Γ shown in(Carmon et al., 2021, Lemma 2) is the smoothness of the function Γ. Proposition B.4. There exists a constant c > 0 such that Γ( • ; d) is c-smooth. Now we introduce the function f nc we use to build our worst-case instance. This function comes from (Zhou and Gu, 2019, Definition 3.5). It is the sum of the quadratic function defined by B.1 and the nonconvex component given by Definition B.3. Definition B.5. For α > 0 and d ∈ N >0 , f nc ( • ; α, d) : R d+1 → R is defined a f nc (x; α, d) = Q(x; √ α, d + 1) + αΓ(x) .

2

  . f nc (0; α, d) -inf x f nc (x; α, d) ≤ x ∈ R d+1 such that x d = x d+1 = 0, ∥∇f nc (x; α, d)∥ ≥ α b) = {U ∈ R a×b , U U ⊤ = I a } .The following Lemma adapted from Zhou and Gu (2019) is fundamental for our lower bound proof.Lemma B.7. Let d, m ∈ N >0 and U ∈ O((d + 1)m, (d + 1)m). We denote U = U (i) ∈ O(d + 1, (d + 1)m). Let {h j } j∈[m] with h j (x) = f nc (U (j) x; α, d) and h = 1 m m j=1 h j . Let x ∈ R (d+1)m and y (j) = U (j) x ∈ R d+1 . Let I = {j ∈ [m], y (i) d = y (i) d+1 = 0}. Then it holds ∥∇h(x)∥ 2 ≥ α

  Cong Fang, Chris Junchi Li, Zhouchen Lin, and Tong Zhang. SPIDER: Near-Optimal Non-Convex Optimization via Stochastic Path Integrated Differential Estimator. In Advances in Neural Information Processing Systems (NeurIPS), 2018. Prashant Khanduri, Siliang Zeng, Mingyi Hong, Hoi-To Wai, Zhaoran Wang, and Zhuoran Yang. A Near-Optimal Algorithm for Stochastic Bilevel Optimization via Double-Momentum. In Advances in Neural Information Processing Systems (NeurIPS), 2021.

	Mao Ye, Bo Liu, Stephen Wright, Peter Stone, and	In International Conference on Machine Learning (b) All the training details (e.g., data splits, hy-
	Qiang Liu. BOME! Bilevel Optimization Made Easy:	(ICML), 2017. perparameters, how they were chosen). Yes
	Jeongyeol Kwon, Dohyun Kwon, Stephen Wright, and Robert Nowak. A Fully First-Order Method for A Simple First-Order Approach. In Advances in Neural Information Processing Systems (NeurIPS), 2022. Miao Zhang, Steven W. Su, Shirui Pan, Xiaojun Chang, Ehsan M Abbasnejad, and Reza Haffari. iDARTS:	Luca Franceschi, Michele Donini, Paolo Frasconi, and Massimiliano Pontil. Forward and Reverse Gradient-[Appendix C and Appendix D] Lam M. Nguyen, Marten van Dijk, Dzung T. Phan, (c) A clear definition of the specific measure or Phuong Ha Nguyen, Tsui-Wei Weng, and Jayant R. statistics and error bars (e.g., with respect to Kalagnanam. Finite-sum smooth optimization with the random seed after running experiments SARAH. Computational Optimization and Applica-tions, 82(3):561-593, 2022. ISSN 0926-6003, 1573-multiple times). Yes [Figure 1]
	Stochastic Bilevel Optimization. In International Differentiable Architecture Search with Stochastic	Based Hyperparameter Optimization. In Interna-2894. doi: 10.1007/s10589-022-00375-x. (d) A description of the computing infrastructure
	Conference on Machine Leaning (ICML), 2023. Junyi Li, Bin Gu, and Heng Huang. A Fully Single Loop Implicit Gradients. In International Conference on Machine Learning (ICML), 2021. Algorithm for Bilevel Optimization without Hessian Dongruo Zhou and Quanquan Gu. Lower Bounds Inverse. In Proceedings of the Thirty-sixth AAAI for Smooth Nonconvex Finite-Sum Optimization. Conference on Artificial Intelligence, AAAI'22, 2022. Yonggang Li, Guosheng Hu, Yongtao Wang, Timothy In International Conference on Machine Learning (ICML), 2019. Hospedales, Neil M. Robertson, and Yongxin Yang. DADA: Differentiable Automatic Data Augmenta-tion. arXiv preprint arXiv:2003.03780, 2020.	tional Conference on Machine Learning (ICML), 2017. used. (e.g., type of GPUs, internal cluster, or Barak A. Pearlmutter. Fast Exact Multiplication by the cloud provider). Yes [Appendix C] Hessian. Neural Computation, 6(1):147-160, 1994. Luca Franceschi, Paolo Frasconi, Saverio Salzo, Ric-ISSN 0899-7667, 1530-888X. doi: 10.1162/neco.1994. 4. If you are using existing assets (e.g., code, data, cardo Grazzi, and Massimilano Pontil. Bilevel Pro-6.1.147. models) or curating/releasing new assets, check if gramming for Hyperparameter Optimization and Meta-Learning. In International Conference on Ma-chine Learning (ICML), 2018. Saeed Ghadimi and Guanghui Lan. Stochastic first-and zeroth-order methods for nonconvex stochastic programming. SIAM Journal on Optimization, 23 you include: Fabian Pedregosa. Hyperparameter optimization with approximate gradient. In International Conference (a) Citations of the creator If your work uses on Machine Learning (ICML), 2016. existing assets. Yes Aravind Rajeswaran, Chelsea Finn, Sham Kakade, and (b) The license information of the assets, if appli-Sergey Levine. Meta-Learning with Implicit Gradi-cable. Not applicable
	Hanxiao Liu, Karen Simonyan, and Yiming Yang. DARTS: Differentiable Architecture Search. In In-ternational Conference on Learning Representations (ICLR), 2019.	(4):2341-2368, 2013. Saeed Ghadimi and Mengdi Wang. Approximation Methods for Bilevel Programming. arXiv preprint arXiv:1802.02246, 2018. ents. In Advances in Neural Information Processing (c) New assets either in the supplemental material Systems (NeurIPS), 2019. or as a URL, if applicable. Yes (d) Information about consent from data providers/curators. Not applicable
	Yair Carmon, John C. Duchi, Oliver Hinder, and Aaron Sidford. Lower bounds for finding stationary points II: First-order methods. Mathematical Programming, 185(1-2):315-355, 2021. Tianyi Chen, Yuejiao Sun, and Wotao Yin. Closing the Gap: Tighter Analysis of Alternating Stochastic Gra-dient Methods for Bilevel Problems. In Advances in Neural Information Processing Systems (NeurIPS), 2021. Ashok Cutkosky and Francesco Orabona. Momentum-based variance reduction in non-convex SGD. In Advances in Neural Information Processing Systems (NeurIPS), 2019. Mathieu Dagréou, Pierre Ablin, Samuel Vaiter, and Thomas Moreau. A framework for bilevel optimiza-tion that enables stochastic and global variance re-duction algorithms. In Advances in Neural Informa-Jonathan Lorraine, Paul Vicol, and David Duvenaud. Optimizing Millions of Hyperparameters by Implicit Differentiation. In International Conference on Ar-tificial Intelligence and Statistics (AISTAT), pages Checklist Arkadiȋ Semenovich Nemirovsky and David Berkovich Yudin. Problem complexity and method efficiency in optimization. Wiley-Interscience series in discrete ISBN 978-0-471-10345-5. [Section 2.1] mathematics. Wiley, Chichester ; New York, 1983. (c) Clear explanations of any assumptions. Yes [Appendix A and Appendix B] (b) Complete proofs of all theoretical results. Yes Advances in Neural Information Processing Systems all theoretical results. Yes [Section 2.1] (NeurIPS), 2022. (a) Statements of the full set of assumptions of cient and collaborative optimization benchmarks. In 2. For any theoretical claim, check if you include: and Samuel Vaiter. Benchopt: Reproducible, effi-Rakotomamonjy, Zaccharie Ramzi, Joseph Salmon, external libraries. Yes [Appendix C] Malézieux, Badr Moufad, Binh T. Nguyen, Alain specification of all dependencies, including 1540-1552, 2020. 1. For all models and algorithms presented, check if Dougal Maclaurin, David Duvenaud, and Ryan P. Adams. Gradient-based Hyperparameter Optimiza-tion through Reversible Learning. In International Conference on Machine Learning (ICML), 2015. Thomas Moreau, Mathurin Massias, Alexandre Gram-fort, Pierre Ablin, Pierre-Antoine Bannier Benjamin Charlier, Mathieu Dagréou, Tom Dupré la Tour, Ghislain Durif, Cassio F. Dantas, Quentin Klopfen-stein, Johan Larsson, En Lai, Tanguy Lefort, Benoit you include: (a) A clear description of the mathematical set-ting, assumptions, algorithm, and/or model. Yes [Section 2] (b) An analysis of the properties and complexity (time, space, sample size) of any algorithm. Yes [Section 3] (c) (Optional) Anonymized source code, with	(e) Discussion of sensible content if applicable, Riccardo Grazzi, Luca Franceschi, Massimiliano Pontil, e.g., personally identifiable information or of-and Saverio Salzo. On the iteration complexity of hypergradient computation. In Hal Daumé III and fensive content. Not Applicable Aarti Singh, editors, International Conference on Machine Learning (ICML), 2020. Mingyi Hong, Hoi-To Wai, Zhaoran Wang, and Zhuo-ran Yang. A Two-Timescale Stochastic Algorithm Framework for Bilevel Optimization: Complexity Analysis and Application to Actor-Critic. SIAM Journal on Optimization, 33(1):147-180, 2023. Quanqi Hu, Yongjian Zhong, and Tianbao Yang. Multi-block Min-max Bilevel Optimization with Applica-tions in Multi-task Deep AUC Maximization. In Processing Systems (NeurIPS), 2022. erated Algorithms for Bilevel Optimization. Journal of Machine Learning Research, 24(22):1-56, 2023. Kaiyi Ji, Junjie Yang, and Yingbin Liang. Bilevel Optimization: Convergence Analysis and Enhanced Design. In International Conference on Machine Bounds for Optimizing Composite Objectives. In Blake E Woodworth and Nati Srebro. Tight Complexity tems (NeurIPS), 2022. In Advances in Neural Information Processing Sys-vergence Theory for Hessian-Free Bilevel Algorithms. Kaiyi Ji and Yingbin Liang. Lower Bounds and Accel-Daouda Sow, Kaiyi Ji, and Yingbin Liang. On the Con-tations (ICLR), 2022. man Distance. In Advances in Neural Information In International Conference on Learning Represen-Huang. Enhanced Bilevel Optimization via Breg-Differentiable Data Augmentation for EEG Signals. Feihu Huang, Junyi Li, Shangqian Gao, and Heng Alexandre Gramfort. CADDA: Class-wise Automatic (NeurIPS), 2022. Cédric Rommel, Thomas Moreau, Joseph Paillard, and Advances in Neural Information Processing Systems 5. If you used crowdsourcing or conducted research Sashank J. Reddi, Suvrit Sra, Barnabas Poczos, and with human subjects, check if you include: Alex Smola. Fast Incremental Method for Nonconvex (a) The full text of instructions given to partici-Optimization. In 2016 IEEE 55th Conference on De-pants and screenshots. Not Applicable cision and Control (CDC), IEEE, pages 1971-1977, (b) Descriptions of potential participant risks, 2016. with links to Institutional Review Board (IRB) Herbert Robbins and Sutton Monro. A stochastic approvals if applicable. Not Applicable approximation method. The Annals of Mathematical (c) The estimated hourly wage paid to partici-Statistics, 22(3):400-407, 1951. doi: 10.1214/aoms/ pants and the total amount spent on partici-1177729586. pant compensation. Not Applicable
	tion Processing Systems (NeurIPS), 2022.	Learning (ICML), 2021. Advances in Neural Information Systems Processing
		(NeurIPS), 2016.
		Junjie Yang, Kaiyi Ji, and Yingbin Liang. Provably
		Faster Algorithms for Bilevel Optimization. In Ad-
		vances in Neural Information Processing Systems
	Yes [Appendix C]	(NeurIPS), 2021.

Aaron Defazio, Francis Bach, and Simon Lacoste-Julien. SAGA: A Fast Incremental Gradient Method With Support for Non-Strongly Convex Composite Objectives. In Advances in Neural Information Processing Systems (NeurIPS), 2014. Hamed Karimi, Julie Nutini, and Mark Schmidt. Linear Convergence of Gradient and Proximal-Gradient Methods Under the Polyak-\Lojasiewicz Condition. In European Conference on Machine Learning (ECML), 2016. Yurii Nesterov. Lectures on Convex Optimization. Springer Berlin Heidelberg, New York, NY, 2018. ISBN 978-3-319-91577-7. Lam M. Nguyen, Jie Liu, Katya Scheinberg, and Martin Takáč. SARAH: A Novel Method for Machine Learning Problems Using Stochastic Recursive Gradient. Zaccharie Ramzi, Florian Mannel, Shaojie Bai, Jean-Luc Starck, Philippe Ciuciu, and Thomas Moreau. SHINE: SHaring the INverse Estimate from the forward pass for bi-level optimization and implicit models. In International Conference on Learning Representations (ICLR), 2022. 3. For all figures and tables that present empirical results, check if you include:

(a) The code, data, and instructions needed to reproduce the main experimental results (either in the supplemental material or as a URL).

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

)

) therefore

The code of the benchmark is available at https://github.com/benchopt/benchmark_bilevel and the results are displayed in https://benchopt.github.io/results/benchmark_bilevel.html.

Acknowledgements

SV acknowledges the support of the ANR GraVa ANR-18-CE40-0005. This work is supported by a public grant overseen by the French National Research Agency (ANR) through the program UDOPIA, project funded by the ANR-20-THIA-0013-01 and DATAIA convergence institute (ANR-17-CONV-0003).

Appendix A contains the necessary lemmas and proofs of Section 3. Appendix B contains the proof of the lower bound for stochastic bilevel optimization. Appendix C details the setting of the numerical experiments. Finally, Appendix D contains two more experiments on hyperparameter selection and datacleaning tasks.

A Convergence analysis of SRBA

A.1 Proof of Proposition 2.6

Proof. Let x ∈ R d . Since G( . , x) is differentiable and z * (x) minimizes G( . , x), the first order optimality condition ensures ∇ 1 G(z * (x), x) = 0 = D z (z * (x), v * (x), x). Since G is strongly convex with respect to z, the Hessian ∇ 2 11 G(z * (x), x) is invertible. As a consequence, the equation in v

Combining Equation ( 21) and Lemma A.2 yields

We get finally

Proof. The smoothness of h (Proposition 2.3) gives us

We finally consider

As a consequence, we have z * (x) = arg min G = z * (x) and h(x) = F (z * (x), x). We also consider a fixed indices sequence (i t , j t ). We set

We can check that each F j is L F 1 -smooth, and each

Let us assume without loss of generality that the algorithm at initialization we have z 0 = v 0 = x 0 = 0 and consider (z t , v t , x t ) the output of an algorithm with the known sequence (i t , j t ).

Given our inner function and the fact that ∇ 2 F (z, x) = 0 for any (z, x) ∈ R (m+1)d+(m+1)d , we have

Since v 0 = 0, we have by Equation ( 53) v 1 ∈ Span(∇ 1 F j0 (z 0 , x 0 )) and by induction

Therefore, using Equation ( 54), we have

Since z 0 = 0, by Equation ( 52), z 1 ∈ Span(x 0 ) and by induction

As a consequence,

Let us denote y (j,t) = U (j) x t . Since x 0 = 0, y (j0,0) = 0 and by the second part of Proposition B.2, x 1 ∈ Span(u

1 ). Now we assume that for all s ≤ t we have

s

) .

There exist scalars α 1 , . . . , α r , β 1,1 , β 2,1 , β 2,2 , . . . , β t,1 , . . . , β t,t such that

Let X r = r s=1 β r,s x s . For r ∈ {1, . . . , t}, we have by induction hypothesis

By orthogonality, we have

T +1 ) .

As a consequence U (jr) X r = (⟨u

1 , X r ⟩, . . . , ⟨u (jr) r , X r ⟩, 0, . . . , 0) . We can use Proposition B.2 to say

r+1 , . . . , u

r+1 ) .

And we get finally x t+1 ∈ Span(u

By induction, for any t, we have and so x t ⊥ Span((u

T +1 ) j∈[m]\{j0,...,jt} , (u

T +1 ) j∈{j0,...,jt} ) .

As a consequence, for t ≤ m 2 T , let I = {j, y

and by Lemma B.7, we have ∥∇h(x t )∥ ≥ ϵ .

If we define T ((x t ) t , h) = inf{t ∈ N, ∥∇h(x t )∥ 2 ≤ ϵ}, we just showed that for the fixed sequence (i t , j t ), we have

The right-hand side being independent from the sequence (i t , j t ), for t ≤ m 2 T , we have

where the expectation is taken over the random choice of i 0 , . . . , i t-1 , j 0 , . . . , j t-1 .

C Details on the experiments

We performed the experiments with the Python package Benchopt (Moreau et al., 2022) 2 . For each experiment, we use minibatches instead of single samples to estimate oracles because it is more efficient in practice. We use a batch size of 64 for the stochastic inner and outer oracles. All the experiments were performed on processors AMD EPYC 7742 (4 CPUs/experiment).

C.1 Benchmark on quadratics

For this benchmark, we consider

The functions F j and G i are defined as To get B Gi , we generate x i ∼ N (0, I p ) and set B Gi = (V Σx i )(U x i ) ⊤ . We proceed similarly for B Fj . In our experiment, we take n = 32768 and m = 1024. To select the parameters of the solvers, we perform a grid search. More precisely, for each solver, we take the inner step size in the form of αt -a where a is the theoretical decrease rate of each solver and α is chosen in {0.01, 0.1}. The outer step size is taken as α r t -b where b is the theoretical decrease rate and r is chosen in {0.1, 1, 10, 100}. For the two-loops algorithms (i.e. StocBiO, VRBO, AmIGO), the number of inner steps is set to 10 after a manual search. In the methods implementing Neumann approximations (MRBO, VRBO, StocBiO), the number of terms in the Neumann series is also set to 10 and the scaling parameter η is set to 0.5. To get the fastest convergence, we keep for each solver the set of parameters that give the best decrease of h on the 100 first epochs. The period of full batch computation of VRBO and SRBA q is parametrized as q = a n+m b where b = 64 is the batch size and a is chosen in {2 -6 , 2 -3 , 2 -1 , 2 3 , 2 6 }. For F 2 SA, we take λ 0 = 1 and δ t = αt -1 7 with α chosen in {0.01, 0.1, 1}.

C.2 Hyperparameter selection with IJCNN1

We solve a regularization selection problem for an ℓ 2 -regularized logistic regression problem. Here, we assume that we have a regularization parameter per feature.

where φ is the logistic loss defined by φ(u) = log(1 + e -u ). The inner and outer step sizes are set to 0.05.

To make our comparison, we select the parameters of each solver with an extensive grid search. More precisely, for each solver, we take the inner step size in the form of αt -a where a is the theoretical decrease rate of each solver and α is chosen in {2 -5 , 2 -4 , 2 -3 , 2 -2 }. The outer step size is taken as α r t -b where b is the theoretical decrease rate and r is chosen in {10 -2 , 10 -1.5 , 10 -1 , 10 -0.5 , 10 0 }. For the two-loops algorithms (i.e. StocBiO, VRBO, AmIGO), the number of inner steps is set to 10 after a manual search. In the methods implementing Neumann approximations (MRBO, VRBO, StocBiO), the number of terms in the Neumann series is also set to 10 and the scaling parameter η is set to 0.5. To get the fastest convergence, we keep for each solver the set of parameters that give the best decrease of h on the 100 first epochs. The period of full batch computation of VRBO and SRBA q is parametrized as q = a n+m b where b = 64 is the batch size and a is chosen in {2 -6 , 2 -3 , 2 -1 , 2 3 , 2 6 , 2 9 }. For F 2 SA, we take λ 0 = 1 and δ t = αt -1 7 with α chosen in {0.01, 0.1, 1}.

D Additional experiment: Datacleaning task

We run an additional experiment. For each experiment, the parameters of the solvers are chosen by an extensive grid search. Then we select the curve that gives the best validation accuracy for each solver and finally plot the corresponding test error on where θ ∈ R C×p , λ ∈ R ntrain , ℓ is the cross entropy loss and σ is the sigmoid function defined by σ(λ) = 1 1+e -λ ∈ (0 , 1].

We run this experiment on the MNIST dataset. We used 20, 000 training samples, 5, 000 validation samples, and 10, 000 test samples. The parameter C r is set to 0.2 after a manual search to get the best performance. For the tuning of the step sizes of each method, we set (ρ t , γ t ) = (αt -a , βt -b ) where (a, b) are the rate provided by the analysis of each method, α is chosen among 4 values between 10 -3 and 10 0 spaced on a logarithmic scale. The scaling parameter β is set to β r where r is chosen among 6 values between 10 -5 and 10 0 spaced on a logarithmic scale. The other parameters are chosen in the same way as the IJCNN1 experiments (see Appendix C.2).

We plot the test error on the Figure D.1 (right). On the one hand, SRBA reaches the best final value. On the other hand, in terms of speed, it is the second fastest after SABA. The other methods are slower and reach a worse final accuracy.