AMORTIZED IMPLICIT DIFFERENTIATION FOR STOCHASTIC BILEVEL OPTIMIZATION
Year: 2021

Abstract:
We study a class of algorithms for solving bilevel optimization problems in both stochastic and deterministic settings when the inner-level objective is strongly convex. Specifically, we consider algorithms based on inexact implicit differentiation and we exploit a warm-start strategy to amortize the estimation of the exact gradient. We then introduce a unified theoretical framework inspired by the study of singularly perturbed systems (Habets, 1974) to analyze such amortized algorithms. By using this framework, our analysis shows these algorithms to match the computational complexity of oracle methods that have access to an unbiased estimate of the gradient, thus outperforming many existing results for bilevel optimization. We illustrate these findings on synthetic experiments and demonstrate the efficiency of these algorithms on hyper-parameter optimization experiments involving several thousands of variables.

heading: 1 INTRODUCTION
text: Bilevel optimization refers to a class of algorithms for solving problems with a hierarchical structure involving two levels: an inner and an outer level. The inner-level problem seeks a solution y?(x) minimizing a cost g(x, y) over a set Y given a fixed outer variable x in a set X . The outer-level problem minimizes an objective of the form L(x)=f(x, y?(x)) over X for some upper-level cost f . When the solution y?(x) is unique, the bilevel optimization problem takes the following form:
min x∈X L(x) := f(x, y?(x)), such that y?(x) = arg min y∈Y g(x, y). (1)
First introduced in the field of economic game theory by Stackelberg (1934), this problem has recently received increasing attention in the machine learning community (Domke, 2012; Gould et al., 2016; Liao et al., 2018; Blondel et al., 2021; Liu et al., 2021; Shaban et al., 2019). Indeed, many machine learning applications can be reduced to (1) including hyper-parameter optimization (Feurer and Hutter, 2019), meta-learning (Bertinetto et al., 2018), reinforcement learning (Hong et al., 2020b; Liu et al., 2021) or dictionary learning (Mairal et al., 2011; Lecouat et al., 2020a;b).
The hierarchical nature of (1) introduces additional challenges compared to standard optimization problems, such as finding a suitable trade-off between the computational budget for approximating the inner and outer level problems (Ghadimi and Wang, 2018; Dempe and Zemkoho, 2020). These considerations are exacerbated in machine learning applications, where the costs f and g often come as an average of functions over a large or infinite number of data points (Franceschi et al., 2018). All these challenges highlight the need for methods that are able to control the computational costs inherent to (1) while dealing with the large-scale setting encountered in machine learning.
Gradient-based bilevel optimization methods appear to be viable approaches for solving (1) in largescale settings (Lorraine et al., 2020). They can be divided into two categories: Iterative differentiation (ITD) and Approximate implicit differentiation (AID). ITD approaches approximate the map y?(x) by a differentiable optimization algorithm A(x) viewed as a function of x. The resulting surrogate loss L̃(x) = f(x,A(x)) is optimized instead of L(x) using reverse-mode automatic differentiation (see Baydin et al., 2018). AID approaches (Pedregosa, 2016) rely on an expression of the gradient ∇L resulting from the implicit function theorem (Lang, 2012, Theorem 5.9). Unlike ITD, AID avoids differentiating the algorithm approximating y?(x) and, instead, approximately
solves a linear system using only Hessian and Jacobian-vector products to estimate the gradient∇L (Rajeswaran et al., 2019). These methods can also rely on stochastic approximation to increase scalability (Franceschi et al., 2018; Grazzi et al., 2020; 2021).
In the context of machine-learning, Ghadimi and Wang (2018) provided one of the first comprehensive studies of the computational complexity for a class of bilevel algorithms based on AID approaches. Subsequently, Hong et al. (2020b); Ji et al. (2021); Ji and Liang (2021); Yang et al. (2021) proposed different algorithms for solving (1) and obtained improved overall complexity by achieving a better trade-off between the cost of the inner and outer level problems. Still, the question of whether these complexities can be improved by better exploiting the structure of (1) through heuristics such as warm-start remains open (Grazzi et al., 2020). Moreover, these studies proposed separate analysis of their algorithms depending on the convexity of the loss L and whether a stochastic or deterministic setting is considered. This points out to a lack of unified and systematic theoretical framework for analyzing bilevel problems, which is what the present work addresses.
We consider the Amortized Implicit Gradient Optimization (AmIGO) algorithm, a bilevel optimization algorithm based on Approximate Implicit Differentiation (AID) approaches that exploits a warm-start strategy when estimating the gradient of L. We then propose a unified theoretical framework for analyzing the convergence of AmIGO when the inner-level problem is strongly convex in both stochastic and deterministic settings. The proposed framework is inspired from the early work of Habets (1974) on singularly perturbed systems and analyzes the effect of warm start by viewing the iterates of AmIGO algorithm as a dynamical system. The evolution of such system is described by a total energy function which allows to recover the convergence rates of unbiased oracle methods which have access to an unbiased estimate of ∇L (c.f. Table 1). To the best of our knowledge, this is the first time a bilevel optimization algorithm based on a warm-start strategy provably recovers the rates of unbiased oracle methods across a wide range of settings including the stochastic ones.

heading: 2 RELATED WORK
text: Singularly perturbed systems (SPS) are continuous-time deterministic dynamical systems of coupled variables (x(t), y(t)) with two time-scales where y(t) evolves much faster than x(t). As such, they exhibit a hierarchical structure similar to (1). The early work of Habets (1974); Saberi and Khalil (1984) provided convergence rates for SPS towards equilibria by studying the evolution of a single scalar energy function summarizing these systems. The present work takes inspiration from these works to analyze the convergence of AmIGO which involves three time-scales.
Two time-scale Stochastic Approximation (TTSA) can be viewed as a discrete-time stochastic version of SPS. (Kaledin et al., 2020) showed that TTSA achieves a finite-time complexity of O
( −1 )
for linear systems while Doan (2020) obtained a complexity of O ( −3/2 ) for general non-linear systems by extending the analysis for SPS. Hong et al. (2020b) further adapted the non-linear TTSA for solving (1). In the present work, we obtain faster rates by taking into account the dynamics of a third variable zk appearing in AmIGO, thus resulting in a three time-scale dynamics.
Warm-start in bilevel optimization. Ji et al. (2021); Ji and Liang (2021) used a warm-start for the inner-level algorithm to obtain an improved computational complexity over algorithms without warm-start. In the deterministic non-convex setting, Ji et al. (2021) used a warm-start strategy when solving the linear system appearing in AID approaches to obtain improved convergence rates. However, it remained open whether using a warm-start when solving both inner-level problem and linear system arising in AID approaches can yield faster algorithms in the more challenging stochastic setting (Grazzi et al., 2020). In the present work, we provide a positive answer to this question.

heading: 3 AMORTIZED IMPLICIT GRADIENT OPTIMIZATION

heading: 3.1 GENERAL SETTING AND MAIN ASSUMPTIONS
text: Notations. In all what follows, X and Y are Euclidean spaces. For a differentiable function h(x, y) : X ×Y → R, we denote by∇h its gradient w.r.t. (x, y), by ∂xh and ∂yh its partial derivatives w.r.t. x and y and by ∂xyh and ∂yyh the partial derivatives of ∂yh w.r.t x and y, respectively.
Geometries Setting Algorithms Complexity BA (Ghadimi and Wang, 2018) O(κ2L ∨ κ2g log 2 −1)
AccBio (Ji and Liang, 2021) O(κ1/2L κ 1/2 g log 2 −1)(D) AmIGO (Corollary 1) O(κLκg log −1)
BSA (Ghadimi and Wang, 2018) O(κ4L −2)
TTSA (Hong et al., 2020b) O(κ0.5L (κ 8.5 g + κ 3 L) −3/2 log −1)
(SC)
(S) AmIGO (Corollary 2) O(κ2Lκ 3 g −1 log −1)
BA (Ghadimi and Wang, 2018) O(κ5g −5/4)
AID-BiO (Ji et al., 2021) O(κ4g −1)(D)
AmIGO (Corollary 3) O(κ4g −1)
BSA (Ghadimi and Wang, 2018) O(κ9g −3 + κ6g −2)
TTSA (Hong et al., 2020b) O(κ16g −5/2 log −1)
stocBiO (Ji et al., 2021) O(κ9g −2+κ6g −2 log −1)
MRBO/VRBO? (Yang et al., 2021) O(poly(κg) −3/2 log −1)
(NC)
(S)
AmIGO (Corollary 4) O(κ9g −2)
Table 1: Cost of finding an -accurate solution as measured by E[L(xk)−L?]∧2−1µE
[
‖xk−x?‖2 ]
when L is µ-strongly-convex (SC) and 1k
∑k
i=1 E
[
‖∇L(xi)‖2
]
when L is non-convex (NC). The settings (D) and (S) stand for the deterministic and stochastic settings. The cost corresponds to the total number of gradients, Jacobian and Hessian-vector products used by the algorithm. κL and κg are the conditioning numbers of L and g whenever applicable. The dependence on κL and κg for TTSA and AccBio are derived in Proposition 11 of Appendix A.4. The rate of MRBO/VRBO is obtained under the additional mean-squared smoothness assumption (Arjevani et al., 2019).
To ensure that (1) is well-defined, we consider the setting where the inner-level problem is strongly convex so that the solution y?(x) is unique as stated by the following assumption: Assumption 1. For any x ∈ X , the function y 7→ g(x, y) is Lg-smooth and µg-strongly convex.
Assumption 1 holds in the context of hyper-parameter selection when the inner-level is a kernel regression problem (Franceschi et al., 2018), or when the variable y represents the last linear layer of a neural network as in many meta-learning tasks (Ji et al., 2021). Under Assumption 1 and additional smoothness assumptions on f and g, the next proposition shows that L is differentiable: Proposition 1. Let g be a twice differentiable function satisfying Assumption 1. Assume that f is differentiable and consider the quadratic problem:
min z∈Rdy
Q(x, y, z) := 1
2 z>(∂yyg(x, y))z + z >∂yf(x, y). (2)
Then, (2) admits a unique minimizer z?(x, y) for any (x, y) in X × Y . Moreover, y?(x) is unique and well-defined for any x in X and L is differentiable with gradient given by:
∇L(x) = ∂xf(x, y?(x)) + ∂xyg(x, y?(x))z?(x, y?(x)). (3)
Proposition 1 follows by application of the implicit function theorem (Lang, 2012, Theorem 5.9) and provides an expression for ∇L solely in terms of partial derivatives of f and g evaluated at (x, y?(x)). Following Ghadimi and Wang (2018), we further make two smoothness assumptions on f and g: Assumption 2. There exist positive constants Lf and B such that for all x, x′ ∈ X and y, y′ ∈ Y:
‖∇f(x, y)−∇f(x′, y′)‖ ≤ Lf‖(x, y)− (x′, y′)‖, ‖∂yf(x, y)‖ ≤ B.
Assumption 3. There exit positive constants L′g , Mg such that for any x, x′ ∈ X and y, y′ ∈ Y:
max {‖∂xyg(x, y)− ∂xyg(x′, y′)‖, ‖∂yyg(x, y)− ∂yyg(x′, y′)‖} ≤Mg‖(x, y)− (x′, y′)‖ ‖∂yg(x, y)− ∂yg(x′, y)‖ ≤ L′g‖x− x′‖.
Assumptions 1 to 3 allow a control of the variations of y? and z? and ensure L is L-smooth for some positive constant L as shown in Proposition 6 of Appendix B.2. As an L-smooth function, L is necessarily weakly convex (Davis et al., 2018), meaning thatL satisfies the inequalityL(x)−L(y) ≤ ∇L(x)>(x− y)− µ2 ‖x− y‖
2 for some fixed µ ∈ R with |µ| ≤ L. In particular, L is convex when µ ≥ 0, strongly convex when µ > 0 and generally non-convex when µ < 0. We thus consider two cases for L, the strongly convex case (µ > 0) and the non-convex case (µ < 0). When L is convex, we denote by L? its minimum value achieved at a point x? and define κL=L/µ when µ > 0. Stochastic/deterministic settings. We consider the general setting where f(x, y) and g(x, y) are expressed as an expectation of stochastic functions f̂(x, y, ξ) and ĝ(x, y, ξ) over a noise variable ξ. We recover the deterministic setting as a particular case when the variable ξ has zero variance, thus allowing us to treat both stochastic (S) and deterministic (D) settings in a unified framework. As often in machine-learning, we assume we can always draw a new batch D of i.i.d. samples of the noise variable ξ with size |D| ≥ 1 and use it to compute stochastic approximations of f and g defined by abuse of notation as f̂(x, y,D) := 1|D| ∑ ξ∈D f̂(x, y, ξ) and ĝ(x, y,D) := 1 |D| ∑ ξ∈D ĝ(x, y, ξ). We make the following noise assumptions which are implied by those in Ghadimi and Wang (2018):
Assumption 4. For any batchD,∇f̂(x, y,D) and ∂y ĝ(x, y,D) are unbiased estimator of∇f(x, y) and ∂yg(x, y) with a uniformly bounded variance, i.e. for all x, y ∈ X × Y:
E [∥∥∥∇f̂(x, y,D)−∇f(x, y)∥∥∥2] ≤ σ̃2f |D|−1, E[‖∂y ĝ(x, y,D)− ∂yg(x, y)‖2] ≤ σ̃2g |D|−1.
Assumption 5. For any batch D, the matrices F1(x, y,D) := ∂xy ĝ(x, y,D) − ∂xyg(x, y) and F2(x, y,D) := ∂yy ĝ(x, y,D)− ∂yyg(x, y) have zero mean and satisfy for all x, y ∈ X × Y:∥∥E[F1(x, y,D)>F1(x, y,D)]∥∥op ≤ σ̃2gxy |D|−1, ∥∥E[F2(x, y,D)>F2(x, y,D)]∥∥op ≤ σ̃2gyy |D|−1. For conciseness, we will use the notations σ2f :=σ̃ 2 f |D| −1, σ2g :=σ̃ 2 g |D| −1, σ2gxy :=σ̃ 2 gxy |D|
−1 and σ2gyy :=σ̃ 2 gyy |D| −1, without explicit reference to the batch D. Next, we describe the algorithm.

heading: 3.2 ALGORITHMS
text: Amortized Implicit Gradient Optimization (AmIGO) is an iterative algorithm for solving (1). It constructs iterates xk, yk and zk such that xk approaches a stationary point of L while yk and zk track the quantities y?(xk) and z?(xk, yk). AmIGO computes the iterate xk+1 using an update equation xk+1 = xk − γkψ̂k for some given step-size γk and a stochastic estimate ψ̂k of ∇L(xk) based on (3) and defined according to (4) below for some new batches of samples Df and Dgxy .
ψ̂k := ∂xf̂(xk, yk,Df ) + ∂x,y ĝ(xk, yk,Dgxy )>zk. (4)
Algorithm 1 AmIGO 1: Inputs: x0, y−1, z−1. 2: Parameters: γk, K. 3: for k ∈ {0, ...,K} do 4: yk ← Ak(xk, yk−1) 5: Sample batches Df , Dg . 6: (uk, vk)← ∇f̂(xk, yk,Df ). 7: zk ← Bk(xk, yk, vk, zk−1) 8: wk ← ∂xy ĝ(xk, yk,Dgxy )zk 9: ψ̂k−1 ← uk + wk 10: xk ← xk−1 − γkψ̂k−1 11: end for 12: Return xK .
AmIGO computes ψ̂k in 4 steps given iterates xk, yk−1 and zk−1. A first step computes an approximation yk to y?(xk) using a stochastic algorithm Ak initialized at yk−1. A second step computes unbiased estimates uk = ∂xf̂(xk, yk,Df ) and vk = ∂y f̂(xk, yk,Df ) of the partial derivatives of f w.r.t. x and y. A third step computes an approximation zk to z?(xk, yk) using a second stochastic algorithm Bk for solving (2) initialized at zk−1. To increase efficiency, algorithm Bk uses the pre-computed vector vk for approximating the partial derivative ∂yf in (2). Finally, the stochastic estimate ψ̂k is computed using (4) by summing the pre-computed vector uk with the jacobian-vector product wk = ∂xy ĝ ( xk, yk,Dgxy ) zk. AmIGO is summarized in Algorithm 1.
AlgorithmsAk and Bk. While various choices forAk and Bk are possible, such as adaptive algorithms (Kingma and Ba, 2015), or accelerated stochastic algorithms (Ghadimi and Lan, 2012), we
focus on simple stochastic gradient descent algorithms with a pre-defined number of iterations T and N . These algorithms compute intermediate iterates yt and zn optimizing the functions y 7→ g(xk, y) and z 7→ Q(xk, yk, z) starting from some initial values y0 and z0 and returning the last iterates yT and zN as described in Algorithms 2 and 3. Algorithm Ak updates the current iterate yt−1 using a stochastic gradient ∂y ĝ(xk, yt−1,Dg) for some new batch of samples Dg and a fixed step-size αk. Algorithm Bk updates the current iterate zt−1 using a stochastic estimate of ∂zQ(xk, yk, zt−1) with step-size βk. The stochastic gradient is computed by evaluating the Hessian-vector product ∂yy ĝ(xk, yk,Dgyy )zt−1 for some new batch of samples Dgyy and summing it with a vector vk approximating ∂yf(xk, yk) provided as input to algorithm Bk.
Warm-start for y0 and z0. Following the intuition that y?(xk) remains close to y?(xk−1) when xk ' xk−1, and assuming that yk−1 is an accurate approximation to y?(xk−1), it is natural to initialize Ak with the iterate yk−1. The same intuition applies when initializing Bk with zk−1. Next, we introduce a framework for analyzing the effect of warm-start on the convergence speed of AmIGO.
Algorithm 2 Ak(x, y0) 1: Parameters: αk, T 2: for t ∈ {1, ..., T} do 3: Sample batch Dgt,k. 4: yt ← yt−1 − αk∂y ĝ ( x, yt−1,Dgt,k ) .
5: end for 6: Return yT .
Algorithm 3 Bk(x, y, v, z0) 1: Parameters: βk, N . 2: for n ∈ {1, ..., N} do 3: Sample batch Dgyyn,k . 4: zn ← zn−1−βk ( ∂yy ĝ ( x, y,Dgyyn,k ) zn−1 + v ) .
5: end for 6: Return zN .

heading: 4 ANALYSIS OF AMORTIZED IMPLICIT GRADIENT OPTIMIZATION

heading: 4.1 GENERAL APPROACH AND MAIN RESULT
text: The proposed approach consists in three main steps: (1) Analysis of the outer-level problem , (2) Analysis of the inner-level problem and (3) Analysis of the joint dynamics of both levels.
Outer-level problem. We consider a quantity Exk describing the evolution of xk defined as follows:
Exk :=  δk 2γk E [ ‖xk − x?‖2 ] + (1− u)E[L(xk)− L?], µ ≥ 0
δk 2γkL2
E [ ‖∇L(xk)‖2 ] , µ < 0.
(5)
where u ∈ {0, 1} is set to 1 in the stochastic setting and to 0 in the deterministic one and δk is a positive sequence that determines the convergence rate of the outer-level problem and is defined by:
δk := ηkγk, ηk+1 := ηk(1 + γk+1(ηk − µ))1µ≥0 + L1µ<0.
with η0 such that γ−10 ≥η0≥µ if µ≥0 and η0=L if µ<0 and where we choose the step-size γk to be a non-increasing sequence with γ0 ≤ 1L . With this choice for δk and by setting u = 1 in (5), E x k recovers the quantity considered in the stochastic estimate sequences framework of Kulunchakov and Mairal (2020) to analyze the convergence of stochastic optimization algorithms when L is convex. When L is non-convex, Exk recovers a standard measure of stationarity (Davis and Drusvyatskiy, 2018). In Section 4.3, we control Exk using bias and variance error E ψ k−1 and V ψ k−1 of ψ̂k given by (6) below where Ek denotes expectation conditioned on (xk, yk, zk−1).
Eψk := E [∥∥∥Ek[ψ̂k]−∇L(xk)∥∥∥2], V ψk := E[∥∥∥ψ̂k − Ek[ψ̂k]∥∥∥2]. (6)
Inner-level problems. We consider the mean-squared errors Eyk and E z k between initializations (y0=yk−1 and z0=zk−1) and stationary values (y?(xk) and z?(xk, yk)) of algorithms Ak and Bk:
Eyk := E [∥∥y0k − y?(xk)∥∥2], Ezk := E[∥∥z0k − z?(xk, yk)∥∥2].
In Section 4.3, we show that the warm-start strategy allows to control Eyk and E z k in terms of previous iterates Eyk−1 and E z k−1 as well as the bias and variance errors in (6). We further prove that such bias and variance errors are, in turn, controlled by Eyk and E z k .
Joint dynamics. Following Habets (1974), we consider an aggregate error Etotk defined as a linear combination of Exk , E y k and E z k with carefully selected coefficients ak and bk:
Etotk = E x k + akE y k + bkE z k . (7)
As suchEtotk represents the dynamics of the whole system. The following theorem provides an error bound for Etotk in both convex and non-convex settings for a suitable choice of the coefficients ak and bk provided that T and N are large enough:
Theorem 1. Choose a batch-size ∣∣Dgyy ∣∣≥1∨ σ̃2gyyµgLg and the step-sizes αk=L−1g , βk=(2Lg)−1, γk=L −1. Set the coefficients ak and bk to be ak:=δ0(1−αkµg)1/2 and bk:=δ0 ( 1− 12βkµg
)1/2 and set the number of iterations T and N of Algorithms 2 and 3 to be of order T=O(κg) and N=O(κg) up to a logarithmic dependence on κg . Let x̂k=u(1−δk)x̂k−1+(1−u(1−δk))xk, with x̂0=x0. Then, under Assumptions 1 to 5, Etotk satisfies: E[L(x̂k)− L?] + Etotk ≤ ( 1− (2κL)−1 )k[ Etot0 + E[L(x0)− L?] ] + 2W2 L , µ ≥ 0 1
k k∑ t=1 Etott ≤ 2 k (E[L(x0)− L?] + Ey0 + Ez0 ) + 2W2 L , µ < 0,
whereW2, defined in (20) of Appendix A.2, is the effective variance of the problem withW2=0 in the deterministic setting and, in the stochastic setting,W2>0 is of the following order:
W2 =O ( δ−10 κ 5 g|Dg| −1 σ̃2g + κ 3 g ∣∣Dgyy ∣∣−1σ̃2gyy + κ2g∣∣Dgxy ∣∣−1σ̃2gxy + κ2g|Df |−1σ̃2f), We describe the strategy of the proof in Section 4.3 and provide a proof outline in Appendix A.1 with exact expressions for all variables including the expressions of T , N andW2. The full proof is provided in Appendix A.2. The choice of ak and bk ensures that E y k and E z k contribute less to Etotk as the algorithms Ak and Bk become more accurate. The effective varianceW2 accounts for interactions between both levels in the presence of noise and becomes proportional to the outer-level variance σ2f when the inner-level problem is solved exactly. In the deterministic setting, all variances σ̃2f , σ̃ 2 g , σ̃ 2 gxy and σ̃ 2 gyy vanish so thatW
2=0. Hence, we characterize such setting byW2=0 and the stochastic one byW2>0. Next, we apply Theorem 1 to obtain the complexity of AmIGO.

heading: 4.2 COMPLEXITY ANALYSIS
text: We define the complexity C( ) of a bilevel algorithm to be the total number of queries to the gradients of f and g, Jacobian/hessian-vector products needed by the algorithm to achieve an error according to some pre-defined criterion. Let the number of iterations k, T andN and sizes of the batches |Dg|, |Df |,
∣∣Dgxy ∣∣ and ∣∣Dgyy ∣∣, be such that AmIGO achieves a precision . Then C( ) is given by: C( ) = k ( T |Dg|+N
∣∣Dgyy ∣∣+ ∣∣Dgxy ∣∣+ |Df |), (8) We provide the complexity of AmIGO in the 4 settings of Table 1 in the form of Corrolaries 1 to 4 .
Corollary 1 (Case µ>0 andW2=0). Use batches of size 1. AchievingL(xk)−L?+µ2 ‖xk − x ?‖2≤ requires C( )=O ( κLκg log ( Etot0 )) .
Corollary 1 outperforms the complexities in Table 1 in terms of the dependence on . It is possible to improve the dependence on κg to κ 1/2 g using acceleration in Ak and Bk as discussed in Appendix A.5.1, or using generic acceleration methods such as Catalyst (Lin et al., 2018).
Corollary 2 (Case µW2>0). Choose |Dg|=Θ ( −1κLκ 2 gσ̃ 2 g ) , ∣∣Dgxy ∣∣=Θ( −1σ̃2gxy), |Df |=Θ( σ̃2f )
and ∣∣Dgyy ∣∣=Θ(σ̃2gyy( 1 ∨ κg)). Achieving E[L(x̂k)−L?]+µ2E[‖xk − x?‖2]≤ requires: C( ) = O ( κL ( κLκ 3 gσ̃ 2 g + κg(1 ∨ κg)σ̃2gyy + σ̃ 2 gxy + σ̃ 2 f )1 log ( Etot0 + E[L(x0)− L?] )) .
Corollary 2 improves over the results in Table 1 in the stochastic strongly-convex setting and recovers the dependence on of stochastic gradient descent for smooth and strongly convex functions up to a logarithmic factor.
Corollary 3 (Case µ<0 andW2=0). Choose batches of size 1. Achieving 1k ∑k i=1 ‖∇L(xi)‖
2 ≤ requires C( ) = O ( κ4g ((L(x0)− L ?) + Ey0 + E z 0 ) ) .
Corollary 3 recovers the complexity of AID-BiO (Ji et al., 2021) in the deterministic non-convex setting. This is expected since AID-BiO also exploits warm-start for both Ak and Bk.
Corollary 4 (Case µ<0 andW>0). Choose ∣∣Dgxy ∣∣=Θ(κ2g σ̃2gxy), ∣∣Dgyy ∣∣=Θ(κ3g σ̃2gyy(1 ∨ µ2g)),
|Df |=Θ ( κ2gσ̃ 2 f ) and |Dg|=Θ ( κ5gσ̃ 2 g ) . Achieving an error 1k ∑k i=1 E [ ‖∇L(xi)‖2 ] ≤ requires:
C( ) = O ( κ5g 2 ( κ4gσ̃ 2 g + κ 2 g ( 1 ∨ µ2g ) σ̃2gyy + σ̃ 2 gxy + σ̃ 2 f ) (E[L(xk)− L?] + Ey0 + Ez0 ) )
Corollary 4 recovers the optimal dependence on of O( 1 2 ) achieved by stochastic gradient descent in the smooth non-convex case (Arjevani et al., 2019, Theorem 1). It also improves over the results in (Ji et al., 2021) which involve an additional logarithmic factor log( −1) as N is required to be

heading: O(κg log(
text: −1)). In our case,N remains constant since Bk benefits from warm-start. The faster rates of MRBO/VRBO? (Yang et al., 2021) are obtained under the additional mean-squared smoothness assumption (Arjevani et al., 2019), which we do not investigate in the present work. Such assumption allows to achieve the improved complexity of O( −3/2 log( −1)). However, these algorithms still require N=O(log( −1)), indicating that the use of warm-start in Bk could further reduce the complexity to O( −3/2) which would be an interesting direction for future work.

heading: 4.3 OUTLINE OF THE PROOF
text: The proof of Theorem 1 proceeds by deriving a recursion for both outer-level error Exk and innerlevel errors Eyk and E z k and then combining those to obtain an error bound on the total error E tot k .
Outer-level recursion. To allow a unified analysis of the behavior of Exk in both convex and nonconvex settings, we define Fk as follows:
Fk :=uδkE[L(xk)− L?]1µ>0 + ( E[L(xk)− L(xk−1)] + Exk−1 − Exk ) 1µ<0.
The following proposition, with a proof in Appendix C.1, provides a recursive inequality on Exk involving the errors in (6) due to the inexact gradient ψ̂k: Proposition 2. Let ρk be a non-increasing sequence with 0<ρk<2. Assumptions 1 to 3 ensure that:
Fk + E x k ≤
( 1− ( 1− 2−1ρk ) δk ) Exk−1 + γkskV ψ k−1 + γk ( sk + ρ −1 k ) Eψk−1, (10)
with sk defined as sk := 12δk + ( u 2 δk + (1− u) ) 1µ>0.
In the ideal case where yk = y?(xk) and zk = z?(xk, yk), the bias E ψ k vanishes and (10) simplifies to (Kulunchakov and Mairal, 2019, Proposition 1) which recovers the convergence rates for stochastic gradient methods in the convex case. However, yk and zk are generally inexact solutions and introduce a positive bias Eψk . Therefore, controlling the inner-level iterates is required to control the bias Eψk which, in turn, impacts the convergence of the outer-level as we discuss next.
Controlling the inner-level iterates yk and zk. Proposition 3 below controls the expected mean squared errors between iterates yk and zk and their limiting values y?(xk) and z?(xk, yk): Proposition 3. Let the step-sizes αk and βk be such that αk≤L−1g and βk≤ 12Lg∧ µg µ2g+σ 2 gyy . Let
Λk:=(1−αkµg)T and Πk:= ( 1−βkµg2 )N . Under Assumptions 1, 4 and 5, it holds that:
E [ ‖yk − y?(xk)‖2 ] ≤ ΛkEyk +R y k, E [ ‖zk − z?(xk, yk)‖2 ] ≤ ΠkEzk +Rzk, (11)
where Ryk=O ( κgσ 2 g ) and Rzk=O ( κ3gσ 2 gyy+κ 2 gσ 2 f ) are defined in (14) of Appendix A.2.
While Proposition 3 is specific to the choice of the algorithmsAk and Bk in Algorithms 2 and 3, our analysis directly extends to other algorithms satisfying inequalities similar to (11) such as accelerated or variance reduced algorithms discussed in Appendices A.5.1 and A.5.2. Proposition 4 below controls the bias and variance terms V ψk and E ψ k in terms of the warm-start error E y k and E z k . Proposition 4. Under Assumptions 1 to 5, the following inequalities hold: Eψk ≤ 2L 2 ψ(ΛkE y k + ΠkE z k +R y k), V ψ k ≤ w 2 x + σ 2 xΠkE z k ,
where w2x=O ( κ2g ( σ2f + σ 2 gxy ) + κ3gσ 2 gyy ) , σ2x=O ( σ2gxy + κ 2 gσ 2 gyy ) and Lψ=O(κ2g) are positive
constants defined in (13) and (16) of Appendix A.2 with Lψ controlling the variations of Ek[ψ̂k]. Proposition 4 highlights the dependence of Eψk and V ψ k on the inner-level errors. It suggests analyzing the evolution of Eyk and E z k to quantify how large the bias and variances can get: Proposition 5. Let ζk > 0, a 2×2 matrix Pk, two vectors Uk and Vk in R2 all independent of xk, yk and zk be as defined in Proposition 8 of Appendix A.2. Under Assumptions 1 to 5, it holds that:(
Eyk Ezk
) ≤ Pk ( Λk−1E y k−1 +R y k−1
Πk−1E z k−1 +R z k−1
) + 2γk ( Eψk−1 + V ψ k−1 + ζkE x k−1 ) Uk + Vk. (12)
Proposition 5 describes the evolution of the inner-level errors as the number of iterations k increases. The matrix Pk and vectors Uk and Vk arise from discretization errors and depend on the step-sizes and constants of the problem. The second term of (12) represents interactions with the outer-level throughExk−1, V ψ k−1 andE ψ k−1. Propositions 2, 4 and 5 describe the joint dynamics of (E x k , E y k , E z k) from which the evolution of Etotk can be deduced as shown in Appendices A.1 and A.2.

heading: 5 EXPERIMENTS
text: We run three sets of experiments described in Sections 5.1 to 5.3. In all cases, we consider AmIGO with either gradient descent (AmIGO-GD) or conjugate gradient (AmIGO-CG) for algorithm Bk. We then compare these algorithms with AID methods without warm-start for Bk which we refer to as (AID-GD) and (AID-CG) and with(AID-CG-WS) which uses warm-start for Bk but not for Ak. We also consider other variants using either a fixed-point algorithm (AID-FP) (Grazzi et al., 2020) or Neumann series expansion (AID-N) (Lorraine et al., 2020) for Bk. Finally, we consider two algorithms based on iterative differentiation which we refer to as (ITD) (Grazzi et al., 2020) and (Reverse) (Franceschi et al., 2017). For all methods except (AID-CG-WS), we use warm-start in algorithm Ak, however only AmIGO, AmIGO-CG and AID-CG-WS exploits warm-start in Bk the other AID based methods initializing Bk with z0=0. In Sections 5.2 and 5.3, we also compare with BSA algorithm (Ghadimi and Wang, 2018), TTSA algorithm (Hong et al., 2020a) and stocBiO (Ji et al., 2021).

heading: 5.1 SYNTHETIC PROBLEM
text: To study the behavior of AmIGO in a controlled setting, we consider a synthetic problem where both inner and outer level losses are quadratic functions with thousands of variables as described in details in Appendix F.1. Figure 1(a) shows the complexity C( ) needed to reach 10−6 relative error amongst the best choice for T and M over a grid as the conditioning number κg increases. AmIGO-CG achieves the lowest time and is followed by AID-CG thus showing a favorable effect of warm-start for Bk. The same conclusion holds for AmIGO-GD compared to AID-GD. Note that AID-CG is still faster than AmIGO-CG for larger values of κg highlighting the advantage of using algorithms Bk with O( √ κg) complexity such as (CG) instead non-accelerated ones with O(κ−1g ) such as (GD). Figure 1(b) shows the relative error after 10s and maintains the same conclusions. For moderate values of κg , only AmIGO and AID-CG reach an error of 10−20 as shown in Figure 1(c). We refer to Figures 2 and 3 of Appendix F for additional results on the effect of the choice of T and M showing that AmIGO consistently performs well for a wide range of values of T and M .

heading: 5.2 HYPER-PARAMETER OPTIMIZATION
text: We consider a classification task on the 20Newsgroup dataset using a logistic loss and a linear model. Each dimension of the linear model is regularized using a different hyper-parameter. The
collection of those hyper-parameters form a vector x of dimension d=101631 optimized using an unregularized regression loss over the validation set while the model is learned using the training set. We consider two evaluations settings: A default setting based on Grazzi et al. (2020); Ji et al. (2021) and a grid-seach search setting near the default values of βk, T and N as detailed in Appendix F.2. We also vary the batch-size from 103 ∗ {0.1, 1, 2, 4} and report the best performing choices for each method. Figure 1(d,e,f) show AmIGO-CG to be the fastest, achieving the lowest error and highest validation and test accuracies. The test accuracy of AmIGO-CG decreases after exceeding 80% indicating a potential overfitting as also observed in Franceschi et al. (2018). Similarly, AmIGO-GD outperformed all other methods that uses an algorithm Bk with O(κg) complexity. Moreover, all remaining methods achieved comparable performance matching those reported in Ji et al. (2021), thus indicating that the warm-start in Bk and acceleration in Bk were the determining factors for obtaining an improved performance. Additionally, Figure 4 of Appendix F report similar results for each choice of the batch-size indicating robustness to the choice of the batch-size.

heading: 5.3 DATASET DISTILLATION
text: Dataset distillation (Wang et al., 2018) consists in learning a synthetic dataset so that a model trained on this dataset achieves a small error on the training set. Figure 5 of Appendix F.3 shows the training loss (outer loss), the training and test accuracies of a model trained on MNIST by dataset distillation. Similarly to Figure 1, AmIGO-CG achieves the best performance followed by AID-CG. AmIGO obtains the best performance by far among methods without acceleration for Bk while all the remaining ones fail to improve. This finding is indicative of an ill-conditioned inner-level problem as confirmed when computing the conditioning number of the hessian ∂yyg(x, y) which we found to be of order 7×104. Indeed, when compared to the synthetic example for κg=104 as shown in Figure 2, we also observe that only AmIGO-CG, AmIGO and AID-CG could successfully optimize the loss. Hence, these results confirm the importance of warm-start for an improved performance.

heading: 6 CONCLUSION
text: We studied AmIGO, an algorithm for bilevel optimization based on amortized implicit differentiation and introduced a unified framework for analyzing its convergence. Our analysis showed that AmIGO achieves the same complexity as unbiased oracle methods, thus achieving improved rates compared to methods without warm-start in various settings. We then illustrated empirically such improved convergence in both synthetic and a hyper-optimization experiments. A future research direction consists in extending the proposed framework to non-smooth objectives and analyzing acceleration in both inner and outer level problems as well as variance reduction techniques.

heading: A CONVERGENCE OF AMIGO ALGORITHM
text: In this section, we provide a proof of Theorem 1 as well as its corollaries Corrolaries 1 to 4. In Appendix A.1, we provide an outline of the proof Theorem 1 that states the main intermediary results needed for the proof and provide explicit expressions for the quantities needed throughout the rest of the paper. Appendices A.2 and A.3 provide the proofs of Theorem 1 and Corrolaries 1 to 4. The proofs of the intermediary results are deferred to Appendices B and C.

heading: A.1 PROOF OUTLINE OF THEOREM 1
text: The proof of Theorem 1 proceeds in 8 steps as discussed bellow.
Step 1: Smoothness properties. This step consists in characterizing the smoothness of ∇L, y?, z? as well as the conditional expectation E[ψ̂k|xk, yk, zk] knowing xk, yk and zk. For this purpose, we consider the function Ψ : X × Y × Y → X defined as follows:
Ψ(x, y, z) := ∂xf(x, y) + ∂xyg(x, y)z.
Hence, by definition of ψ̂k, it is easy to see that E[ψ̂k|xk, yk, zk] = Ψ(xk, yk, zk). The following proposition controls the smoothness of ∇L, y?, z? and Ψ and is adapted from (Ghadimi and Wang, 2018, Lemma 2.2). We provide a proof in Appendix B.2 for completeness.
Proposition 6. Under Assumptions 1 to 3, L, Ψ, y? and z? satisfy:
‖y?(x)− y?(x′)‖ ≤ Ly‖x− x′‖, ‖z?(x, y)− z?(x′, y′)‖ ≤ Lz(‖x− x′‖+ ‖y − y′‖) ‖∇L(x)−∇L(x′)‖ ≤ L‖x− x′‖, ‖Ψ(x, y, z)−∇L(x)‖ ≤ Lψ[‖y − y?(x)‖+ ‖z − z?(x, y)‖]
‖z?(x, y?(x))‖ ≤ µ−1g B,
where Ly , Lz , Lψ and L are given by:
Ly := µ −1 g L ′ g, Lz := µ −2 g MgB + µ −1 g Lf , (13) Lψ := max (( Lf +Mgµ −1 g B + L ′ gLz ) , L′g ) L := ( Lf + µ −2 g L ′ gMgB + µ −1 g ( L′gLf +MgB ))( 1 + µ−1g L ′ g
) The expressions of Ly , Lz , Lψ and L suggests the following dependence on the conditioning κg of the inner-level problem which will be useful for the complexity analysis: Ly = O(κg), Lψ = O(κ2g), Lz = O(κ 2 g) and L = O(κ 3 g).
Step 2: Convergence of the inner-level iterates. In this step, we control the mean squared errors E [ ‖yk − y?(xk)‖2 ] and E [ ‖zk − z?(xk, yk)‖2 ] as stated in Proposition 3. In fact we prove a
slightly stronger version stated below:
Proposition 7. Let αk and βk be two positive sequences with αk≤L−1g and βk≤ 12Lg min ( 1, 2Lg
µg(1+µ −2 g σ2gyy )
) and define Λk:=(1−αkµg)T and Πk:= ( 1−βkµg2 )N . De-
note by z̄k the conditional expectation of zk knowing xk, yk and z0k. Let R y k and R z k be defined as:
Ryk := 2αkµ −1 g σ 2 g , R z k := βkB 2µ−3g σ 2 gyy + 3µ −2 g σ 2 f . (14)
Then, under Assumptions 1, 4 and 5 the iterates zk and z̄k satisfy: E [ ‖yk − y?(xk)‖2 ] ≤ ΛkEyk +R y k
E [ ‖zk − z?(xk, yk)‖2 ] ≤ ΠkE [∥∥z0k − z?(xk, yk)∥∥2]+Rzk, E [ ‖z̄k − z?(xk, yk)‖2 ] ≤ ΠkE
[∥∥z0k − z?(xk, yk)∥∥2] E [ ‖zk − z̄k‖2 ] ≤ 4σ2gµ−2g ΠkE [∥∥z0k − z?(xk, yk)∥∥2]+ 2Rzk. (15a)
It is easy to see from the above expressions that Ryk = O ( κgσ 2 g ) while Rzk = O ( κ3gσ 2 gyy + κ 2 gσ 2 f ) as stated in Proposition 3. Controlling E[‖yk − y?(xk)‖]2 follows by standard results on SGD (Kulunchakov and Mairal, 2020, Corollary 31) since the iterates of Algorithm 2 uses i.i.d. samples. The error terms E [ ‖zk − z?(xk, yk)‖2 ] is more delicate since Algorithm 3 uses the same sample
∂y f̂(xk, yk) for updating the iterates, therefore introducing additional correlations between these iterates. We defer the proof of Proposition 7 to Appendix B.3 which relies on a general result for stochastic linear systems with correlated noise provided in Appendix E.
Step 3: Controlling the bias and variance errors V ψk and E ψ k is achieved by Proposition 4. The bias Eψk is controlled simply by using the smoothness of the potential Ψ near the point (xk, y ?(xk), z ?(xk, y ?(xk)) as shown in Proposition 6. The variance term V ψ k is more delicate to control due to the multiplicative noise resulting from the Jacobian-vector product between ∂xy ĝ(xk, yk,Dgxy )zk. We defer the proof to Appendix B.4 and provide below explicit expressions for the constants σ2x and w 2 x:
w2x := ( 1 + 2L′gµ −1 g + 6 ( σ2gxy + (L ′ g) 2 ) µ−2g ) σ2f (16a)
+ 2 ( σ2gxy + (L ′ g) 2 ) B2L−1g µ −3 g σ 2 gyy + 2B 2µ−2g σ 2 gxy .
σ2x :=2σ 2 gxy + 2(L ′ g) 2µ−2g σ 2 gyy .
Note that w2x=O ( κ2g ( σ2f + σ 2 gxy ) +κ3gσ 2 gyy ) and σ2x=O ( σ2gxy+κ 2 gσ 2 gyy ) , as stated in Proposition 4.
Step 4: Outer-level error bound. This step consists in obtaining the inequality in Proposition 2 which extends the result of (Kulunchakov and Mairal, 2020, Proposition 1) to biased gradients and to the non-convex case. We defer the proof of such result to Appendix C.1.
Step 5: Inner-level error bound. This step consists in proving Proposition 5. For clarity, we provide a second statement with explicit expressions for the quantities of interest:
Proposition 8. Let rk and θk be two positive non-increasing sequences no greater than 1. For any 0 ≤ v ≤ 1, denote by φk and R̃k the following non-negative scalars:
φk := (1− v)L2gT 2α2k + 2v, R̃ y k := 2(1− v)Lgµ −1 g T 2σ2gα 2 k + vR y k ζk := 2L ( min ( (1− u)−1, Lη−1k−1 ) 1µ≥0 + 1µ<0 ) Finally, consider the following matrices and vectors:
Pk :=
( 1 + rk, 0
16L2z φk θk , 1 + θk
) , Uk := ( 2L2y γk rk
4L2z γk θk ( 1 + 4L2y φk rk )) , Vk := R̃yk ( 08L2zθ−1k ) .
(17)
Then, under Assumptions 1 to 5, the following holds:( Eyk Ezk ) ≤ Pk ( Λk−1E y k−1 +R y k−1 Πk−1E z k−1 +R z k−1 ) + 2γk ( Eψk−1 + V ψ k−1 + ζkE x k−1 ) Uk + Vk.
We defer the proof of the above result to Appendix C.2.
Step 6: General error bound. By combining the inequalities in Propositions 2 and 8 resulting from the analysis of both outer and inner levels, we obtain a general error bound on Etotk in Proposition 9 with a proof in Appendix C.4.
Proposition 9. Choose the step-sizes αk and βk such that they are non-increasing in k and choose rk and θk such that δkr−1k and δkθ −1 k are non-increasing sequences. Choose the coefficients ak and bk defining Etotk in (7) to be of the form ak = δkr −1 k Λ s k and bk = δkθ −1 k Π s k for some 0 < s < 1. and fix a non-increasing sequence 0 < ρk < 1. Then, under Assumptions 1 to 5 Etotk satisfies:
Fk + E tot k ≤ ‖Ak‖∞E tot k−1 + V tot k .
where Vk and Ak are given by:
Ak :=  1− (1− 12ρk)δk + 2ζkλku I k Λ1−sk ( 1 + rk ( 1 + 16L2zΠ s kθ −2 k φk + 2L 2 ψη −1 k ( 2uIk + sk + ρ −1 k
))) Π1−sk ( 1 + θk ( 1 + η−1k ( 2L2ψ + σ 2 x )( 2uIk + sk + ρ −1 k
)))  (18)
V totk :=δk ( Λsk(1 + r −1 k ) + 16L 2 zφkθ −2 k Π s k + 2L 2 ψη −1 k ( 2uIk + sk + ρ −1 k )) Ryk−1
+ δk ( (1 + θ−1k )Π s kR z k−1 + 8L 2 zθ −2 k Π s kR̃ y k ) + γk ( sk + u I k ) w2x
where we introduced uIk = akU (1) k + bkU (2) k for conciseness with U (1) k and U (2) k being the components of the vector Uk defined in Proposition 8.
Proposition 9 holds without conditions on the error made by Algorithms 2 and 3. The general form of ak and bk allows to account for potentially decreasing step-sizes γk, αk and βk. However, in the present work, we will restrict to the constant step-size for ease of presentation as we discuss next.
Step 7: Controlling the precision of the inner-level algorithms. In this step, we provide conditions on T andN in Proposition 10 bellow so that ‖Ak‖∞ ≤ 1−(1−ρk)δk in the constant step-size case. These conditions are expressed in terms of the following constants:
C1 :=1 + 2 log ( 6 + 24L2ψη −1 0 ) (19a)
C2 :=2 log ( 1 + 4L2yL −2 max (η0, 8ζ0) )
(19b)
C3 := max
( 0,−2 log ( 5L2ψη −1 0 ) ,−2 log ( L
4L2y
)) (19c)
C ′1 :=1 + 2 log ( 4 + 12η−10 (2L 2 ψ + σ 2 x) )
(19d) C ′2 :=2 log ( 1 + 2L2zL −2 y ( 1 + 16L2y )) , (19e)
C ′3 := max
( 0,−2 log ( L2ψ
4L2zη0
) ,−2 log ( γ ( σ2gxy + (L ′ g) 2 )))
(19f)
Proposition 10. Let Assumptions 1 to 5 hold. Choose ρk = 12 , the step-sizes to be constant:
αk=α≤ 1Lg , βk=β≤ 1 2Lg and γk=γ≤ 1L and choose the batch-size ∣∣Dgyy ∣∣ = Θ( σ̃2gyyµgLg). Moreover, set s = 12 , rk = θk = 1. Finally, choose T and N as follows:
T = bα−1µ−1g max (C1, C2, C3)c+ 1, N = b2β−1µ−1g (max (C1, C2, C3) + max (C ′1, C ′2, C ′3))c+ 1
with C1, C2, C3, C ′1, C ′ 2 and C ′ 3 defined in (19a) to (19f). Then, ‖Ak‖∞ ≤ 1 − (1 − ρk)δk and V totk ≤ γδ0W2, withW2 given by:
W2 := ( δ−10
1− u 2 1µ>0 + 3
) w2x +
60L2ψ δ0µgLg σ2g (20)
We provide a proof of Proposition 10 in Appendix D. It is easy to see from Proposition 10 that that T=O(κg) and N=O(κg) when α= 1Lg and β= 1 2Lg
, where the big-O notation hides a logarithmic dependence in κg coming from the constants {Ci, C ′i|i ∈ {1, 2, 3}}. Step 8: Proving the main inequalities. The final step combines Propositions 9 and 10 to get the desired inequality. We provide a full proof in Appendix A.2 assuming Propositions 9 and 10 hold.

heading: A.2 PROOF OF THEOREM 1
text: In order to prove Theorem 1 in the convex case, we need the following averaging strategy lemma, a generalization of (Kulunchakov and Mairal, 2020, Lemma 30):
Lemma 1. Let L be a convex function onX . Let xk be a (potentially stochastic) sequence of iterates in X . Let (Ek)k≥0 , (Vk)k≥0 and (δk)k≥0 be non-negative sequences such that δk ∈ (0, 1). Fix some non-negative number u ∈ [0, 1] and define the following averaged iterates x̂k recursively by x̂k = u(1−δk)x̂k−1+(1− (1− δk)u)xk and starting from any initial point x̂0. Assume the iterates (xk)k≥1 satisfy the following relation for all k ≥ 1:
(1− u(1− δk))E[L(xk)− L?] + Ek ≤ (1− δk)(Ek−1 + (1− u)E[L(xk−1)− L?]) + Vk. (21)
Let Γk := ∏k t=1(1− δk). Then the averaged iterates (x̂k)k≥1 satisfy the following:
E[L(x̂k)− L?] + Ek ≤ Γk ( E0 + E [ L(x0)− L? + uk(L(x̂0)− L?) ]) + Γk ∑ 1≤t≤k Γ−1t Vt.
Proof. For simplicity, we write Fk = E[L(xk)− L?] and F̂k = E[L(x̂k)− L?]. We first multiply (21) by Γ−1k and sum the resulting inequalities for all 1 ≤ k ≤ K to get:
K∑ k=1 Γ−1k (1− u(1− δk))Fk + Γ −1 k Tk ≤ K∑ k=1 Γ−1k (1− δk)(Tk−1 + (1− u)Fk−1) + K∑ k=1 Γ−1k Vk.
Grouping the terms in Fk together and recalling that Γ−1k (1− δk) = Γ −1 k−1 yields:
K∑ k=1 Γ−1k Fk − Γ −1 k−1Fk−1 + u K∑ k=1 Γ−1k−1(Fk−1 − Fk) ≤ K∑ k=1 ( Γ−1k−1Tk−1 − Γ −1 k Tk ) + K∑ k=1 Γ−1k Vk.
Simplifying the telescoping sums and multiplying by Γk, we get:
FK + uΓK K∑ k=1 Γ−1k−1(Fk−1 − Fk) ≤ −TK + ΓK
( F0 + T0 +
K∑ k=1 Γ−1k Vk
) . (22)
Consider now the quantity F̂k. Recalling that L is convex and by definition of the iterates x̂k we apply Jensen’s inequality to write:
F̂k ≤ u(1− δk)F̂k−1 + (1− u(1− δk))Fk.
By iteratively applying the above inequality, we get that:
F̂K ≤uKΓK F̂0 + ΓK K∑ k=1 uK−kΓ−1k (1− u(1− δk))Fk
=uKΓK F̂0 + ΓK K∑ k=1 uK−kΓ−1k Fk − K∑ k=1 uK−k+1Γk−1Fk
=uKΓK F̂0 + FK + ΓK K−1∑ k=1 uK−kΓ−1k (Fk − Fk+1)
We can therefore apply (22) to the above inequality to get the desired result.
We now proceed to prove Theorem 1.
Proof of Theorem 1. By application of Proposition 9 and using the choice of T and N given by Proposition 10, the following inequality holds:
Fk + E tot k ≤ (1− (1− ρk)δk)Etotk−1 + γδ0W2. (23)
withW defined in Proposition 10. We then distinguish two cases depending on the sign of µ: Case µ ≥ 0. Recall that Fk and Exk are given by:
Fk = uδkE[L(xk)− L?], Exk = δk
2γk ‖xk − x?‖2 + (1− u)E[L(xk)− L?]
Since µ > 0, L is a convex function and we can apply Lemma 1 with Vk=γδ0W2 and Ek=
δk 2γk ‖xk − x?‖2+akEyk+bkEzk . The result follows by noting that Γk ∑k t=1 Γ −1 t ≤ δ−10 .
Case µ < 0. In this case, we recall that Fk and Exk are given by: Fk = E[L(xk)− L(xk−1)] + Exk−1 − Exk , Exk = L−1E [ ‖∇L(xk)‖2 ] .
We then sum (23) for all iterations 0 ≤ t ≤ k which, by telescoping, yields: E[L(xk)− L(x0)] + Ex0 − Exk + Etotk − Etot0 + ∑ 1≤t≤k (1− ρt)δtEtott ≤ kγδ0W2.
Using that E[L(xk)− L?] +Etotk −Exk is non-negative since Etotk −Exk = akE y k + bkE z k , we get:∑
1≤t≤k
(1− ρt)δtExt ≤ (E[L(x0)− L?] + a0E y 0 + b0E z 0 ) + kγδ0W2.
Finally, since ρt = 12 , δt = Lγ, the result follows after dividing both sides by 1 2kLγ.

heading: A.3 PROOF OF CORROLARIES 1 TO 4
text: Proof of Corollary 1. Choosing u = 0 implies that Exk = µ 2 ‖xk − x ?‖2 +L(xk)−L? ≤ Etotk . We can then apply Theorem 1 for µ > 0 which yields the following:
µ 2 ‖xk − x?‖+ L(xk)− L? ≤
( 1− (2κL)−1 )k Etot0 + 2W2
L . (24)
In the deterministic setting, it holds all variances vanish : σ2f=σ 2 g=σ 2 gyy=σ 2 gxy=0. Hence,W 2 = 0 by definition of W2. Therefore, to achieve an error L(xk)−L?≤ for some >0, (24) suggests choosing k=O ( κL log ( Etot0 )) . Additionally, T=Θ(κg) and N=Θ(κg) as required by Theorem 1
and since σ2gyy=0, it holds that N=O(κg). Using batches of size 1, yields the desired complexity.
Proof of Corollary 2. Here we choose u = 1 and apply Theorem 1 for µ > 0 which yields: E[L(x̂k)− L?] + Etotk ≤ ( 1− (2κL)−1 )k( Etot0 + E[L(x0)− L?] ) + 2L−1W2
Hence, to achieve an error E[L(x̂k)− L?] ≤ , we need k = O ( κL log ( Etot0 +E[L(x0)−L ?] )) to
guarantee that the first term in the l.h.s. of the above inequality is O( ). Moreover, we recall that L−1 = O(κ−3g ) from Proposition 6 and that Theorem 1 ensure the varianceW satisfies:
W2 =O ( κ5gσ 2 g + κ 3 gσ 2 gyy + κ 2 gσ 2 gxy + κ 2 gσ 2 f ) .
Hence, ensuring the variance term 2L−1W2 is of order is achieved by choosing the size of the batches as follows:
|Df | = Θ
( σ̃2f ) , |Dg| = Θ ( κLκ 2 gσ̃ 2 g ) , ∣∣Dgxy ∣∣ = Θ ( σ̃2gxy ) ,
∣∣Dgyy ∣∣ = Θ(σ̃2gyy(1 ∨ κg ))
Recall that T=Θ(κg) andN=Θ(κg) as required by, Theorem 1, thus yielding the desired result.
Proof of Corollary 3. In the non-convex deterministic case, recall thatExk = 1 L‖∇L(xk)‖ 2 ≤ Etotk . We thus apply Theorem 1 for µ < 0, multiply by L to get:
1
k k∑ t=1 ‖∇L(xt)‖2 ≤ 2L k (L(x0)− L? + (Ey0 + Ez0 )) + 2W2.
The setting being deterministic, it holds that W2=0. Moreover, recall that L = O(κ3g) from Proposition 6. Hence, to achieve an error of order min1≤t≤k ‖∇L(xt)‖2 ≤ , it suffice to choose k=O ( κ3g (L(x0)− L ? + (Ey0 + E z 0 )) ) . Thus using batches of size 1 and T and N of order κg .
Proof of Corollary 4. In the non-convex stochastic case, Exk = 1 LE [ ‖∇L(xk)‖2 ] ≤ Etotk . We thus apply Theorem 1 for µ < 0, multiply by L to get:
1
k k∑ t=1 E [ ‖∇L(xt)‖2 ] ≤ 2L k (E[L(x0)− L?] + (Ey0 + Ez0 )) + 2W2.
to achieve an error of order , we need to ensure each term in the l.h.s. of the above inequality is of order . For the first term, similarly to the deterministic setting Corollary 3, we simply need k = O ( κ3g (L(x0)− L ? + (Ey0 + E z 0 )) )
. For the second term, we need to haveW2=O( ), which is achieved using the following choice for the sizes of the batches:
|Df | = O
( κ2gσ̃ 2 f ) , |Dg| = O ( κ5gσ̃ 2 g ) , ∣∣Dgxy ∣∣ = O ( κ2gσ̃ 2 gxy ) ,
∣∣Dgyy ∣∣ = O ( κ3gσ̃ 2 gyy ( 1 ∨ µ2g )) .
Finally, as required by Theorem 1, we set T = Θ(κg) and N = Θ(κg) thus yielding the desired complexity.

heading: A.4 COMPARAISONS WITH OTHER METHODS
text: In this subsection we derive and discuss the complexities of methods presented in Table 1.

heading: A.4.1 COMPARAISON WITH TTSA (HONG ET AL., 2020B)
text: Proposition 11. strongly-convex case µ > 0. The complexity of the TTSA algorithm in Hong et al. (2020b) to achieve an error µ2E [ ‖xk − x?‖2 ] ≤ is given by:
C( ) := O (( κ 17 2 g κ 1/2 L + κ 7/2 L ) 1
3 2
log 1 ) non-convex case µ < 0. The complexity of the TTSA algorithm in Hong et al. (2020b) to achieve an error 1k ∑ 1≤i≤k E [ ‖∇L(xi)‖2 ] ≤ is given by:
C( ) := k(1 +N) = O (( κ11g + κ 16 g ) 1
5 2
log
( 1 )) .
Proof. strongly-convex case µ > 0 Using the choice of step-sizes in Hong et al. (2020b), the following bound holds:
E [ ‖xk − x?‖2 ] . k∏ i=0 ( 1− 8 3(k + kα) )( ∆0x + L 2 ψµ −2 g ∆ y 0 ) + L2ψ
µgµ 4 3
( µ−1g +
µg µ2 L2y
)( 1
k + kα
) 2 3
where ∆0x = E [ ‖x0 − x?‖2 ] , E [ ‖y1 − y?(x0)‖2 ] and kα given by:
kα = max
( 35 ( L3g µ3g (1 + σ2g) 3 2 ) , (512) 3 2L2ψL 2 y µ2 ) .
By a simple calculation, it is easy to see that ∏k i=0 ( 1− 83(k+kα) ) ≤ (kα−1) 2 (k−1+kα)2 . Moreover, using that Lψ = O(µ−2g ), Ly = O(µ −1 g ), we get that
µ 2 E [ ‖xk − x?‖2 ] . (kα − 1)2 (k + kα − 1)2 ( µ∆0x + µµ −6 g ∆ y 0 ) + µ−6g µ − 13 ( 1 + µ−2 )( 1 k + kα ) 2 3
Using that L=O(µ−3g ), we get µ −6 g µ
− 13 ( 1 + µ−2 ) =O ( µ−5g κ 1 3 L + µgκ 7 3 L ) . Hence, to reach an error
, we need to control both terms in the above inequality. This suggests the following condition on k to control the second term which dominates the error:
k ≥ ( µ 3 2 g κ 7/2 L + κ 1/2 L µ − 152 g ) 1
3 2
.
Moreover, the result in (Hong et al., 2020b, Theorem 1) requires N = Θ ( κg log 1 ) , where N is the number of terms in the Neumann series used to approximate the hessian inverse (∂yyg(x, y))−1 in the expression of the gradient∇L. Hence, the total complexity is given by the following expression:
C( ) := k(1 +N) = O (( κ 17 2 g κ 1/2 L + κ 7/2 L ) 1
3 2
log 1 ) Smooth Non-convex case µ < 0. Following Hong et al. (2020b), consider the proximal map of L for a fixed ρ > 0:
x̂(z) := arg min x∈X
{ L(x) + ρ
2 ‖x− z‖2 } and define the quantity ∆̃kx := E [ ‖x̂(xk)− xk‖2 ] , where xk are the iterates produced by the TTSA
algorithm. Let K be a random variable uniformly distributed on {0, ...,K − 1} and independent from the remaining r.v. used in the TTSA algorithm. The result in (Hong et al., 2020b, Theorem 2) provide the following error bound on ∆̃kx
1
k ∑ 1≤i≤k ∆̃ix . ( L2ψ ( ∆0 + σ2g µ2g ) + µg ) k− 2 5 L2 .
where ρ is set to 2L and ∆0≤max ( E[L(x0)−L?],E [ ‖y1−y?(x0)‖2 ]) . Now, recall that by defi-
nition of the proximal map, we have the following identity:
1
k ∑ 1≤i≤k E [ ‖∇L(xi)‖2 ] ≤ 2(ρ2 + L2) 1 k ∑ 1≤i≤k ∆̃ix . L 2 1 k ∑ 1≤i≤k ∆̃ix.
Hence, we obtain the following error bound:
1
k ∑ 1≤i≤k E [ ‖∇L(xi)‖2 ] ≤ ( κ4g ( ∆0 + κ2gσ 2 g ) + µg ) k− 2 5
where we used that Lψ = O(κ2g). Therefore, to reach an error of order , TTSA requires:
k ≥ 1
5 2
( κ10g ∆ 5 2 0 + κ 15 g σ 5 g ) Moreover, controlling the bias in the estimation of the gradient requires N = O(κg log 1 ) terms in the Neumann series approximating the hessian. Hence, the total complexity of the algorithm is:
C( ) := k(1 +N) = O (( κ11g + κ 16 g ) 1
5 2
log
( 1 )) .

heading: A.4.2 COMPARAISON WITH ACCBIO (JI AND LIANG, 2021)
text: Complexity of AccBio. The bilevel algorithm AccBio introduced in Ji and Liang (2021) uses acceleration for both the inner and outer loops. This allows to obtain the following conditions on k, T and N to achieve an accuracy:
k = O ( κ 1 2 L log 1 ) , T = O ( κ 1 2 g ) , N = O ( κ 1 2 g log 1 ) .
Note that, since AccBio do not use a warm-start strategy when solving the linear system, N is required to grow as log 1 in order to achieve an accuracy. This contributes an additional logarithmic factor to the total complexity so that C( ) = O(κ 1 2
Lκ 1 2 g
( log 1 )2 ). This is by contrast with AmIGO
which exploits warm start when solving the linear system and thus only needs a constant number of iterations N = O(κg) although the dependence on κg is worse compared to AccBio. However, it is possible to improve such dependence by using acceleration in the inner-level algorithms Ak and Bk as we discuss in Appendix A.5.1.
Complexity of AccBio as a function of µ and µg . The authors choose to report the complexity as a function of µ and µg instead of the conditioning numbers κL and κg . To achieve this, they observe that, under the additional assumption that the hessian ∂yyg(x, y) is constant w.r.t. y, the Lipschitz constant L has an improved dependence on µg: L = O(µ−2g ) instead of L = O(µ −3 g ) in the general case where ∂yyg(x, y) is only Lipschitz in y. This allows them to express κL = Lµ = O(µ −2 g µ −1) and to report the following complexities in terms of µ and µg:
C( ) = O ( µ− 1 2µ − 32 g ( log 1 )2) .
Note that, in the general case where L = O(µ−3g ), the complexity as a function of µ and µg becomes O ( µ−
1 2µ−2g
( log 1 )2) , while still maintaining the same expression in terms of κL and κg . Hence,
using the expression in terms of conditioning allows a more general expression for the complexity that is less sensitive to the specific assumptions on g and is therefore more suitable for comparaison with other results in the literature.
A.5 CHOICE OF THE INNER-LEVEL ALGORITHMS Ak AND Bk .
The choice of Ak and Bk has an impact on the total complexity of the algorithm. We discuss two choices for Ak and Bk which improve the total complexity of AmIGO: Accelerated algorithms (in Appendix A.5.1) and variance reduced algorithms (in Appendix A.5.2).

heading: A.5.1 ACCELERATION OF THE INNER-LEVEL FOR AMIGO
text: AmIGO could benefit from acceleration in the inner-loop by using standard acceleration schemes Nesterov (2003) for Ak and Bk. As a consequence, and using analysis of accelerated algorithms (Nesterov, 2003) in the deterministic setting, the error of the inner-level iterates would satisfy:
E [ ‖yk − y?(xk)‖2 ] ≤ Λ̃kEyk , E [ ‖yk − y?(xk)‖2 ] ≤ Π̃kEzk
where Λ̃k and Π̃k are accelerated rates of the form Λ̃k = O((1− √ κg) T ) and Π̃k = O((1− √ κg)
N ). The rest of the proofs are similar provided that Λk and Πk are replaced by their accelerated rates Λ̃k and Π̃k. This implies that T and N need to be only of order T = O( √ κg) and N = O( √ κg) so that the final complexity becomes:
C( ) := O ( κLκ 1/2 g log 1 ) .
Note that using conjugate gradient for Bk also enjoys an accelerated convergence rate Shewchuk et al. (1994). This is confirmed in our experiments of Figure 1 where AmIGO-CG enjoys the fastest convergence.
In order to further improve the dependence on κL to κ 1/2 L , one would need to use an accelerated scheme when updating the iterates xk. The analysis of such scheme along with warm-start would be an interesting direction for future work.
A.5.2 VARIANCE REDUCED ALGORITHMS FOR Ak AND Bk When the inner-level cost function g is a finite average of functions g(x, y) = 1n ∑
1≤i≤n gi(x, y) empirical average, it is possible to use variance reduced algorithms such as SAG (Schmidt et al.,
2017). If every function gi is Lg-smooth, then by (Schmidt et al., 2017, Proposition 1), the inner level error becomes:
E [ ‖yk − y?(xk)‖2 ] . Λ̃k(3E y k + 9
4 L−2g σ 2 g),
with Λ̃k = ( 1− κg8n )T . This has the advantage that the error due to the variance decays exponentially with the number of iterations T . As a consequence, the dependence of the effective varianceW2 on the conditioning numbers κL and κg can be improved to:
W2 =O ( |Dg|−1σ̃2g + κ3g ∣∣Dgyy ∣∣−1σ̃2gyy + κ2g∣∣Dgxy ∣∣−1σ̃2gxy + κ2g|Df |−1σ̃2f). This can be achieved by taking T = O(nκg) up to a logarithmic dependence on the condition numbers. As a consequence, the complexity in the strongly convex stochastic setting becomes:
C( ) = O ( κL ( nσ̃2g + κg(1 ∨ κg)σ̃2gyy + σ̃ 2 gxy + σ̃ 2 f )1 log ( Etot0 + E[L(x0)− L?] )) .
In the non-convex setting, the complexity becomes:
C( ) = O ( κ4g 2 ( nσ̃2g + κ 3 g ( 1 ∨ µ2g ) σ̃2gyy + κgσ̃ 2 gxy + κgσ̃ 2 f ) (E[L(xk)− L?] + Ey0 + Ez0 ) ) .
The downside of this approach is the dependence on the number n of functions gi in the total complexity.

heading: B PRELIMINARY RESULTS

heading: B.1 EXPRESSION OF THE GRADIENT
text: We provide a proof of Proposition 1 which shows that L is differentiable and provides an expression of its gradient.
Proof. Assumption 1 ensures that y 7→ g(x, y) admits a unique minimizer y?(x) defined as the unique solution to the implicit equation ∂yg(x, y?(x)) = 0. Moreover, since g is twice continuously differentiable and strongly convex, it follows that ∂yyg(x, y?(x)) is invertible for any x ∈ X . Therefore the implicit function theorem (Lang, 2012, Theorem 5.9), ensures that x 7→ y?(x) is continuously differentiable with Jacobian given by ∇y?(x) = −∂xyg(x, y?(x))∂yyg(x, y?(x))−1. Hence, by composition of differentiable functions, L is also differentiable with gradient given by:
∇L(x) = ∂xf(x, y?(x))− ∂xyg(x, y?(x))∂yyg(x, y?(x))−1∂yf(x, y?(x)).
We can thus define z?(x, y) = −∂yyg(x, y)−1∂yf(x, y) to get the desired expression for ∇L(x) and note that z? is the solution to (2).
B.2 SMOOTHNESS PROPERTIES OF L, Y? , z? AND Ψ
Proof of Proposition 6. Lipschitz continuity of x 7→ y?(x). By Assumptions 1 and 3, the implicit function theorem (Lang, 2012, Theorem 5.9) ensures y?(x) is differentiable with Jacobian given by:
∇y?(x) = −∂xyg(x, y?(x))(∂yyg(x, y?(x)))−1.
Moreover, by Assumption 3, we know that ∂yg(x, y) is L′g-Lipchitz in x for any y ∈ Y , hence, ‖∂xyg(x, y?(x))‖op is upper-bounded by L ′ g . Moreover, by Assumption 1, g is µg-strongly convex
in y uniformly on X . Therefore, it holds that ∥∥∂yyg(x, y?(x))−1∥∥op ≤ µ−1g . This allows to deduce that ‖∇y?(x)‖op ≤ µ−1g L′g , and by application of the fundamental theorem of calculus that:
‖y?(x)− y?(x′)‖ ≤ µ−1g L′g‖x− x′‖.
This shows that y? is Ly-Lipschitz continuous with Ly := µ−1g L ′ g .
Lipchitz continuity of x 7→ z?(x, y). Let (x, y) and (x′, y′) be two points in X × Y . Recalling the definition of z?(x, y) in Proposition 1, it is easy to see that z?(x, y) admits the following expression:
z?(x, y) = −∂yyg(x, y)−1∂yf(x, y). (25) Recalling the expression of z?(x, y), the following holds:
z?(x, y)− z?(x′, y′) =∂yyg(x′, y′)−1∂yf(x′, y′)− ∂yyg(x, y)−1∂yf(x, y) = ( ∂yyg(x ′, y′)−1 − ∂yyg(x, y)−1 ) ∂yf(x ′, y′)
+ ∂yyg(x, y) −1(∂yf(x ′, y′)− ∂yf(x, y)) =∂yyg(x
′, y′)−1(∂yyg(x, y)− ∂yyg(x′, y′))∂yyg(x, y)−1∂yf(x′, y′) + ∂yyg(x, y) −1(∂yf(x ′, y′)− ∂yf(x, y))
Hence, by taking the norm of the above quantity a triangular inequality followed by operator inequalities, it follows that:
‖z?(x, y)− z?(x′, y′)‖ ≤ ∥∥H−12 ∥∥op‖H1 −H2‖op∥∥H−11 ∥∥op‖∂yf(x′, y′)‖ + ∥∥H−11 ∥∥‖∂yf(x′, y′)− ∂yf(x, y)‖.
where we introduced H1 = ∂yyg(x, y) and H2 = ∂yyg(x′, y′) for conciseness. Using Assumption 1, we can upper-bound ∥∥H−11 ∥∥op and ∥∥H−12 ∥∥op by µ−1g . By Assumption 3, we know that ‖H1 −H2‖op ≤ Mg‖(x, y)− (x′, y′)‖. Finally by Assumption 2, we also have that ‖∂yf(x′, y′)− ∂yf(x, y)‖ ≤ Lf‖(x, y)− (x′, y′)‖ and that ‖∂yf(x′, y′)‖ ≤ B ensuring that:
‖z?(x, y)− z?(x′, y′)‖ ≤ ( µ−2g MgB + µ −1 g Lf ) ‖(x, y)− (x′, y′)‖.
Hence, we conclude that z? is Lz-Lipchitz continuous with Lz defined as in (13).
boundedness of z?(x, y) Recalling the expression of z? in (25), it is easy to see that ‖z?(x, y)‖ is upper-bounded by µ−1g B since ∂yyg(x, y) is µg-strongly convex in y by Assumption 1 and ∂yf(x, y) is bounded by B by Assumption 2.

heading: Regularity of Ψ.
text: Ψ(x, y, z)−Ψ(x′, y′, z′) =∂xf(x, y)− ∂xf(x′, y′) + ∂xyg(x, y)z − ∂xyg(x′, y′)z′
=∂xf(x, y)− ∂xf(x′, y′) + ∂xyg(x, y)(z − z′) + (∂xyg(x, y)− ∂xyg(x′, y′))z′.
By taking the norm of the above expression and applying a triangular inequality followed by operator inequalities, it follows that:
‖Ψ(x, y, z)−Ψ(x′, y′, z′)‖ ≤‖∂xf(x, y)− ∂xf(x′, y′)‖+ ‖∂xyg(x, y)‖op‖z − z ′‖ (26)
+ ‖∂xyg(x, y)− ∂xyg(x′, y′)‖op‖z ′‖
≤Lf (‖x− x′‖+ ‖y − y′‖) + L′g‖z − z′‖ +Mg‖z′‖(‖x− x′‖+ ‖y − y′‖).
To get the first term of the last inequality above, we used that ∂xf is Lf -Lipschitz by Assumption 2. To get the second term, we used that ∂xyg(x, y) is bounded since ∂yg(x, y) is L′g-Lipschitz by Assumption 3. Finally, for the last term, we used that ∂xyg(x, y) is Mg-Lipschitz by Assumption 3.
By choosing x′ = x, y′ = y?(x) and z′ = z?(x, y?(x)), it is easy to see from Proposition 1 that Ψ(x, y?(x), z?(x, y?(x))) = ∇L(x). Hence, applying the above inequality yields: ‖Ψ(x, y, z)−∇L(x)‖ ≤(Lf +Mg‖z?(x, y?(x))‖)‖y − y?(x)‖+ L′g‖z − z?(x, y?(x))‖
≤(Lf +Mg‖z?(x, y?(x))‖)‖y − y?(x)‖+ L′g‖z − z?(x, y)‖ + L′g‖z?(x, y)− z?(x, y?(x))‖
As shown earlier, ‖z?(x, y?(x))‖ is upper-bounded by µ−1g B, while ‖z?(x, y)− z?(x, y?(x))‖ is bounded by Lz‖y − y?(x)‖. This allows to conclude that ‖Ψ(x, y, z)−∇L(x)‖ ≤ Lψ with Lψ defined in (13).
Lipschitz continuity of x 7→ ∇L(x). We apply (26) with (y, z) = (y?(x), z?(x, y?(x))) and (y′, z′) = (y?(x′), z?(x, y?(x′))) which yields:
‖∇L(x)−∇L(x′)‖ ≤(Lf +Mg‖z?(x′, y?(x′))‖)(‖x− x′‖+ ‖y?(x)− y?(x′)‖) + L′g‖z?(x, y?(x))− z?(x′, y?(x′))‖ ≤ ( Lf +Mgµ −1 g B + L ′ gLz ) (1 + Ly)‖x− x′‖,
where we used that ‖z?(x′, y?(x′))‖ is upper-bounded by µ−1g B, z? is Lz-Lipschitz and y? is LyLipschitz. Hence,∇L is L-Lipschitz continuous, with L as given by (13).
B.3 CONVERGENCE OF THE ITERATES OF ALGORITHMS Ak AND Bk
Proof. Controlling the iterates yt of Ak. Consider a new batch Dg of samples ξ. We have by definition of the update equation of yt that:∥∥yt − y?(xk)∥∥2 =∥∥yt−1 − y?(xk)∥∥2 + α2k∥∥∂y ĝ(xk, yt−1,Dg)∥∥2
− 2αk∂y ĝ ( xk, y t−1,Dg )>( yt−1 − y?(xk) )
Taking the expectation conditionally on xk and yt−1, we get:
E [∥∥yt − y?(xk)∥∥2∣∣∣xk, yt−1] =(1− αkµg)∥∥yt−1 − y?(xk)∥∥2
+ α2kE [∥∥∂y ĝ(xk, yt−1,Dg)− ∂yg(xk, yt−1)∥∥2∣∣∣xk, yt−1]
− 2αk∂yg ( xk, y t−1)>(yt−1 − y?(xk)− αk 2 ∂yg(xk, y t−1) )
≤(1− αkµg) ∥∥yt−1 − y?(xk)∥∥2 + 2α2kσ2g
The first line uses that ∂y ĝ ( xk, y t−1,Dg )
is an unbiased estimator of ∂yg(xk, yt−1). For the second line, we use Assumption 4 which allows to upper-bound the variance of ∂y ĝ by σ2g . Moreover, since g is convex and Lg-smooth and since αk ≤ L−1g , it follows that the last term in the above inequality is non-positive and can thus be upper-bounded by 0. By unrolling the resulting inequality recursively for 1 < t ≤ k, we obtain the desired result. Controlling the iterates zn of Bk. The poof follows by direct application of Proposition 15 with β = βk and the following choices for An, A, b̂, b:
An =∂yy ĝ(xk, yk,Dgyy ), b̂ =∂y f̂(xk, yk,Df ) A =∂yyg(xk, yk) b =∂yf(xk, yk).
This directly yields the following inequalities: E [ ‖zk − z?(xk, yk)‖2 ] ≤Π̃kE [∥∥z0k − z?(xk, yk)∥∥2]+ R̃zk, E [ ‖z̄k − z?(xk, yk)‖2 ] ≤Π̃kE
[∥∥z0k − z?(xk, yk)∥∥2]. where Π̃k and R̃zk are given by:
Π̃k := (1− βkµg)N , R̃zk := β2k ( σ2gE [ ‖z?(xk, yk)‖2 ] + 3 ( N ∧ 1
βkµg
) σ2f )( N ∧ 1
βkµg
) .
First we have that Π̃k ≤ Πk. Moreover, Proposition 6, we have that ‖z?(xk, yk)‖ ≤ Bµ−1g hence, R̃zk ≤ Rzk thus yielding the desired inequalities. Finally (15a) also follows similarly using (45) from Proposition 15.
B.4 CONTROLLING THE BIAS AND VARIANCE Eψk AND V ψ k
Proof of Proposition 4 . Recall that the expressions of Eψk and V ψ k in (6) involves the conditional expectation Ek [ ψ̂k ] knowing xk, yk and zk−1. This can be also expressed using Ψ as follows:
Ek [ ψ̂k ] =E [ E [ ψ̂k|xk, yk, zk ] |xk, yk, zk−1 ] =E[Ψ(xk, yk, zk)|xk, yk, zk−1] =E[Ψ(xk, yk,E[zk||xk, yk, zk−1])]
where we used the tower property for conditional expectations in the first line, then the fact that the expectation of ψk conditionally on xk, yk and zk is simply Ψ(xk, yk, zk). Finally, for the last line, we use the independence of the noise and the linearity of Ψ w.r.t. the last variable. In all what follows, we write z̄k = E[zk|xk, yk, zk−1] which is the same object as defined in Proposition 7. We then treat Eψk and V ψ k separately.
Bounding Eψk . Using Propositions 6 and 7 we directly get the desired inequality:
Eψk ≤2L 2 ψ ( E [ ‖yk − y?(xk)‖2 ] + E [ ‖z̄k − z?(xk, yk)‖2 ]) ≤2L2ψ(ΛkE y k + ΠkE z k +R y k)
Bound on V ψk . We decompose V ψ k into a sum of three terms Wk, W ′ k and W ′′ k given by:
Wk := E [∥∥∥∂xf̂(xk, yk,Df )− ∂xf(xk, yk)∥∥∥2],
W ′k := E [∥∥∥∂xy ĝ(xk, yk, ξ̃N+1,k)zk − ∂xyg(xk, yk)z̄k∥∥∥2]
W ′′k := E [( ∂xf̂(xk, yk,Df )− ∂xf(xk, yk) )> ∂xyg(xk, yk)(zk − z̄k) ] .
where we used that ξ̃N+1,k is independent from zk and Df to get the last term. Hence, using Assumption 4 to bound the first term of the above relation, we get V ψk ≤ σ2f +W ′k + 2W ′′k . Thus, it remains to control each of W ′k and W ′′ k .
Bound on W ′′k . Using that Df is independent from ξ̃n,k, we can apply Proposition 14 to write:
W ′′k = βkE [ ∂x(f̂ − f)(xk, yk,Df )>∂xyg(xk, yk) ( N∑ t=1 (I − βkA)N−t ) ∂y(f̂ − f)(xk, yk,Df ) ]
where we used the simplifying notion (f̂ − f)(xk, yk,Df ) = f̂(xk, yk,Df ) − f(xk, yk). Using Assumption 3 to bound ‖∂xyg(xk, yk)‖op by L
′ g , Assumption 1 to upper-bound∥∥∥(∑Nt=1 (I − βkA)N−t)∥∥∥ op by (∑N t=1 (1− βkµg) N−t ) we get
|W ′′k | ≤βkL′g N−1∑ t=0 (1− βkµg)tE [∣∣∣∂x(f̂ − f)(xk, yk,Df )>∂y(f̂ − f)(xk, yk,Df )∣∣∣]
≤L′gµ−1g E [∣∣∣∂x(f̂ − f)(xk, yk,Df )>∂y(f̂ − f)(xk, yk,Df )∣∣∣]
≤L′gµ−1g E [∥∥∥∂x(f̂ − f)(xk, yk,Df )2∥∥∥] 12E[∥∥∥∂y(f̂ − f)(xk, yk,Df )∥∥∥2] 12 ≤ L′gµ−1g σ2f
where we used that ∑N−1 t=0 (1− βkµg)≤ 1 βkµg
for the second line, Cauchy-Schwarz inequality to get the third line and Assumption 4 to get the last line.
Bound on W ′k Using that ξ̃N+1,k is independent from zk, we write:
W ′k =E [∥∥∥∂xy(ĝ − g)(xk, yk, ξ̃N+1,k)zk∥∥∥2]+ E[‖∂xyg(xk, yk)(zk − z̄k)‖2]
(i) ≤σ2gxyE [ ‖zk‖2 ] + (L′g) 2E [ ‖zk − z̄k‖2 ] (ii)
≤ 2σ2gxy ( E [ ‖zk − z?(xk, yk)‖2 ] + E [ ‖z?(xk, yk)‖2 ]) + (L′g) 2E [ ‖zk − z̄k‖2 ] (iii)
≤ 2σ2gxyE [ ‖zk − z?(xk, yk)‖2 ] + (L′g) 2E [ ‖zk − z̄k‖2 ] + 2σ2gxyB 2µ−2g
(iv) ≤ 2σ2gxy ( ΠkE [∥∥z0k − z?(xk, yk)∥∥2 +Rzk])
+ (L′g) 2 ( 4σ2gyyµ −2 g ΠkE [∥∥z0k − z?(xk, yk)∥∥2]+ 2Rzk)+ 2σ2gxyB2µ−2g ≤2 ( σ2gxy + 2(L ′ g) 2µ−2g σ 2 gyy ) ΠkE
[∥∥z0k − z?(xk, yk)∥∥2] + 2 ( σ2gxy + (L ′ g) 2 ) Rzk + 2σ 2 gxyB 2µ−2g
(i) follows from Assumptions 3 and 5, (ii) uses that ‖zk‖2 ≤ 2 ( ‖zk − z‖2 + ‖z‖2 ) , (iii) uses that
‖z?(xk, yk)‖ ≤ Bµ−1g by Proposition 6. Finally (iv) follows by application of Proposition 7. We further have by definition of Rzk that:
Rzk ≤ B2L−1g µ−3g σ2gyy + 3µ −2 g σ 2 f (29)
Combining the inequalities on W ′k, W ′′ k and (29), we get that V ψ k ≤w2x+σ2xΠkEzk , with w2x and σ2x given by (16).

heading: C GENERAL ANALYSIS OF AMIGO

heading: C.1 ANALYSIS OF THE OUTER-LOOP
text: Proof of Proposition 2. We treat both cases µ ≥ 0 and µ < 0 separately. For simplicity we denote by Ek the conditional expectation knowing the iterates xk, yk and zk−1 and write ψk = Ek [ ψ̂k ] .
Case µ ≥ 0. Recall that Exk is given by:
Exk = ηk 2 E [ ‖xk − x?‖2 ] + (1− u)E[L(xk)− L?].
For simplicity define k = uδk + (1 − u), ek = (1 − u)(L(xk)− L?) + ηk2 ‖xk − x ?‖2 and e′k = uδk(L(xk)− L?). It is then easy to see that E[ek] is equal to the l.h.s of (10), i.e. E[ek] = Exk . We
will start by bounding the difference between two successive iterates of ek:
e′k + ek − ek−1 ≤uδk(L(xk)− L?) + (1− u)(L(xk)− L(xk−1))
+ ηk 2 ‖xk − x?‖2 − ηk−1 2 ‖xk−1 − x?‖2
(i) ≤uδk(L(xk)− L?) + (1− u)(L(xk)− L(xk−1))
− δkηk−1 2 ‖xk−1 − x?‖2 + δk (µ 2 ‖xk−1 − x?‖2 −∇L(xk−1)>(xk−1 − x?) ) + ηk 2 ( γ2k
∥∥∥ψ̂k−1∥∥∥2 − 2γk(ψ̂k−1 −∇L(xk−1))>(xk−1 − x?)) (ii)
≤ uδk(L(xk)− L?) + (1− u)(L(xk)− L(xk−1))
− δkηk−1 2 ‖xk−1 − x?‖2 − δk(L(xk−1)− L?)
+ ηk 2
( γ2k ∥∥∥ψ̂k−1∥∥∥2 − 2γk(ψ̂k−1 −∇L(xk−1))>(xk−1 − x?)) ≤(uδk + (1− u))(L(xk)− L(xk−1))
− δkηk−1 2 ‖xk−1 − x?‖2 − δk(1− u)(L(xk−1)− L?)
+ ηk 2
( γ2k ∥∥∥ψ̂k−1∥∥∥2 − 2γk(ψ̂k−1 −∇L(xk−1))>(xk−1 − x?)) ≤ k(L(xk)− L(xk−1))− δkek−1
+ ηk 2
( γ2k ∥∥∥ψ̂k−1∥∥∥2 − 2γk(ψ̂k−1 −∇L(xk−1))>(xk−1 − x?)) (iii)
≤ − δkek−1 + k∇L(xk−1)>(xk − xk−1) + kL
2 ‖xk − xk−1‖2
+ ηk 2
( γ2k ∥∥∥ψ̂k−1∥∥∥2 − 2γk(ψ̂k−1 −∇L(xk−1))>(xk−1 − x?)) =− δkek−1 − γk k∇L(xk−1)>ψ̂k−1 + kγ 2 kL
2 ∥∥∥ψ̂k−1∥∥∥2 + ηk 2 ( γ2k
∥∥∥ψ̂k−1∥∥∥2 − 2γk(ψ̂k−1 −∇L(xk−1))>(xk−1 − x?)). (i) follows from the update expression xk = xk−1 − γkψ̂k−1, (ii) follows from the convexity of L and (iii) follows by L-smoothness of L. Taking the expectation conditionally on the randomness at iteration k − 1 and using that Ek−1 [ ψ̂k−1 ] = ψk−1, we therefore get
Ek−1[e′k + ek − ek−1] ≤− δkek−1 − γk k∇L(xk−1)>ψk−1 + γk 2
(δk + k)Ek−1 [∥∥∥ψ̂k−1∥∥∥2]
− δk(ψk−1 −∇L(xk−1))>(xk−1 − x?) =− δkek−1 + γksk ( Ek−1 [∥∥∥ψ̂k−1 − ψk−1∥∥∥]2 + ‖ψk−1 −∇L(xk−1)‖2) − δk(ψk−1 −∇L(xk−1))>(xk−1 − γk∇L(xk−1)− x?)
− γk 2 ( k − δk)‖∇L(xk−1)‖2
(i) ≤ − δkek−1 + γksk ( Ek−1 [∥∥∥ψ̂k−1 − ψk−1∥∥∥]2 + ‖ψk−1 −∇L(xk−1)‖2) − δk(ψk−1 −∇L(xk−1))>(xk−1 − γk∇L(xk−1)− x?)
where (i) follows from δk ≤ k since by construction δk ≤ 1. Taking the expectation w.r.t. all the randomness and applying Cauchy-Schwarz inequality to the last term yields the following inequality:
uδk(L(xk)− L?) + Exk ≤(1− δk)Exk−1 + γksk ( V ψk−1 + E ψ k−1 ) (30)
+ δk ( Eψk−1 ) 1 2E [ ‖xk−1 − γk∇L(xk−1)− x?‖2 ] 1 2 .
SinceL is convex, we have the inequality: ‖xk−1 − γk∇L(xk−1)− x?‖2 ≤ ‖xk−1 − x?‖2. Hence, we can deduce that:
δk‖xk−1 − γk∇L(xk−1)− x?‖2 ≤ δk‖xk−1 − x?‖2 ≤ 2γkηkη−1k−1E x k−1 ≤ 2γkExk−1,
where we used that ηk is non-increasing by construction. Combining the above inequality with (30) yields:
Fk + E x k ≤(1− δk)Exk−1 + γksk ( V ψk−1 + E ψ k−1 ) + √ 2γ 1 2 k δ 1 2 k ( Eψk−1 ) 1 2 ( Exk−1 ) 1 2 .
Case µ < 0. Recall that for µ < 0, we set Exk = 1 LE [ ‖∇L(xk)‖2 ] . Using that L is L-smooth, we have that:
L(xk)− L(xk−1) ≤∇L(xk−1)>(xk − xk−1) + L
2 ‖xk − xk−1‖2
≤− γk∇L(xk−1)>ψ̂k−1 + Lγ2k
2 ∥∥∥ψ̂k−1∥∥∥2 ≤− γk‖∇L(xk−1)‖2 − γk∇L(xk−1)> ( ψ̂k−1 −∇L(xk−1)
) + Lγ2k
2 (∥∥∥ψ̂k−1 − ψk−1∥∥∥2 + 2(ψ̂k−1 − ψk−1)>ψk−1 + ‖ψk−1‖2). Taking the expectation w.r.t. all randomness in the algorithm in the above inequality, we get:
E[L(xk)− L(xk−1)] ≤− γkE [ ‖∇L(xk−1)‖2 ] − γkE [ ∇L(xk−1)>(ψk−1 −∇L(xk−1)) ] + Lγ2k
2
( E [(∥∥∥ψ̂k−1 − ψk−1∥∥∥2)]+ E[‖ψk−1‖2])
=− γk(1− Lγk 2 )E [ ‖∇L(xk−1)‖2 ] + Lγ2k 2 ( V ψk−1 + E ψ k−1 ) − γk(1− Lγk)E [ ∇L(xk−1)>(ψk−1 −∇L(xk−1))
] ≤− γk 2 E [ ‖∇L(xk−1)‖2 ] + Lγ2k 2 ( V ψk−1 + E ψ k−1
) + γkE [ ‖∇L(xk−1)‖2 ] 1 2 ( Eψk−1 ) 1 2 .
=− δkExk−1 + δkγk
2
( V ψk−1 + E ψ k−1 ) + √ 2δ 1 2 k γ 1 2 k ( Exk−1 ) 1 2 ( Eψk−1 ) 1 2 .
where we used that 1− Lγk2 ≥ 1 2 and 0 ≤ 1−Lγk ≤ 1 to get the last inequality. Using the definition of Fk yields an inequality of the form:
Fk + E x k ≤(1− δk)Exk−1 + γksk ( V ψk−1 + E ψ k−1 ) + √ 2γ 1 2 k δ 1 2 k ( Eψk−1 ) 1 2 ( Exk−1 ) 1 2 .
Hence, in both cases µ ≥ 0 and µ < 0 we get an inequality of the of the same form, but with different expressions for Fk and sk. We get the desired result using Young’s inequality, to upper-bound the last term in the r.h.s. of the above inequality. More precisely, we use that for any 0 < ρk < 1:
√ 2γ 1 2
k δ 1 2 k ( Eψk−1 ) 1 2 ( Exk−1 ) 1 2 ≤ 1
2 ρkδkE
x k−1 + ρ −1 k γkE ψ k−1.
C.2 INNER-LEVEL ERROR BOUND
In this section we prove Proposition 5 which controls the evolutions of the warm-start errors Eyk and Ezk . As a first step, in Proposition 12, we provide a result controlling the mean squared error between two successive iterates xk−1, xk and yk−1, yk which will be used in the proof of Proposition 5.
Proposition 12 (Control of the increments of xk and yk). Consider ζk, φk and R̃yk as defined in Proposition 8 for some fixed 0 ≤ v ≤ 1. Then, the following holds:
γ2kE [∥∥∥ψ̂k−1∥∥∥2] = E[‖xk − xk−1‖2] ≤γ2k(V ψk−1 + 2Eψk−1 + 2ζkExk−1)
E [ ‖yk − yk−1‖2 ] ≤2φkEyk + 2R̃ y k
Proof. Proof of Proposition 12 We prove each inequality separately.
Increments of xk. By the update equation, we have that xk = xk−1− γkψ̂k−1, hence we only need to control E [∥∥∥ψ̂k−1∥∥∥2]. We have the following:
E [∥∥∥ψ̂k−1∥∥∥2] ≤E[∥∥∥ψ̂k−1 − ψk−1∥∥∥2]+ 2E[‖ψk−1 −∇L(xk−1)‖2]+ 2E[‖∇L(xk−1)‖2]
=V ψk−1 + 2E ψ k−1 + 2E [ ‖∇L(xk−1)‖2 ] .
In the case (µ < 0), we have Exk−1 = 1 2LE [ ‖∇L(xk−1)‖2 ] , hence by setting ζk := 2L, we get the desired inequality. In the convex case (µ ≥ 0), since L is L-smooth, we have that:
‖∇L(xk−1)‖2 ≤ 2L(L(xk−1)− L?) ≤ 2L(1− u)−1Exk−1,
provided that u < 1. We also have that (L(xk−1)− L?) ≤ L2 ‖xk−1 − x ?‖2 ≤ Lη−1k−1Exk−1 which yields ‖∇L(xk−1)‖2 ≤ 2L2η−1k−1Exk−1. Hence, we can set ζk = 2Lmin ( (1− u)−1, Lη−1k−1 ) .
Increments of yk. Denoting by Dtg a batch of samples at time iteration t of algorithm Ak and using the update equation of yt we get the following inequality by application of the triangular inequality: E [ ‖yk − yk−1‖2 ] 1 2 ≤αk
T−1∑ t=0 E [∥∥∂y ĝ(xk, yt,Dtg)∥∥2] 12 ≤ αk T−1∑ t=0 ( σ2g + L 2 gE [∥∥yt − y?(xk)∥∥2]) 12
≤αk T−1∑ t=0 ( σ2g + L 2 gRk + L 2 gΛt,kE y k ) 1 2 ≤ αkT ( σ2g ( 1 + 2L2gαkµ −1 g ) + L2gE y k ) 1 2
where we applied Proposition 7 for every 0<t≤T−1 to get the second line with Λt,k:=(1− αkµg)t. This directly implies the following bound:
E [ ‖ỹk − yk−1‖2 ] ≤ 2α2kT 2 ( σ2g ( 1 + 2L2gµ −1 g αk ) + L2gE y k ) . (31)
On the other hand, using a triangular inequality and applying Proposition 7, we also have that: E [ ‖yk − yk−1‖2 ] ≤ 2E [ ‖yk − y?(xk)‖2 ] + 2E [ ‖yk−1 − y?(xk)‖2 ] ≤ (4Eyk + 2R y k). (32)
The result follows by combining (31) and (32) using coefficients 1−v and v.

heading: C.3 PROOF OF PROPOSITION 5
text: Proof of Proposition 5 . We will control each of Eyk and E z k separately.
Upper-bound on Eyk . Let rk be a non-increasing sequences between 0 and 1. The following holds:
Eyk =E [ ‖yk−1 − y?(xk−1)‖2 ] + E [ ‖y?(xk − y?(xk−1)‖2 ] + 2E [ (yk−1 − y?(xk−1))>(y?(xk − y?(xk−1))
] (i)
≤(1 + rk)E [ ‖yk−1 − y?(xk−1)‖2 ] + ( 1 + r−1k ) E [ ‖y?(xk − y?(xk−1)‖2 ] (ii)
≤ (1 + rk) ( Λk−1E y k−1 +R y k−1 ) + 2r−1k E [ ‖y?(xk)− y?(xk−1)‖2 ] (iii)
≤ (1 + rk) ( Λk−1E y k−1 +R y k−1 ) + 2L2yr −1 k E [ ‖xk − xk−1‖2 ] (iv)
≤ (1 + rk) ( Λk−1E y k−1 +R y k−1 ) + 2L2yr −1 k γ 2 kE [∥∥∥ψ̂k−1∥∥∥2] (33a)
(i) follows by Young’s inequality, (ii) uses Proposition 7 to bound the first term and that (1+r−1k ) ≤ 2r−1k for the second term, (iii) uses that y
? is Ly-Lipschitz by Proposition 6 and (iv) uses the update equation xk = xk−1 − γkψ̂k−1. Upper-bound on Ezk . Similarly, for a non-increasing sequence 0 < θk ≤ 1, we have that:
Ezk =E [ ‖zk−1 − z?(xk−1, yk−1)‖2 ] + E [∥∥∥z?(xk, yk)− z?(xk−1, yk−1)2∥∥∥] + 2E [ (zk−1 − z?(xk−1, yk−1))>(z?(xk, yk)− z?(xk−1, yk−1))
] (i)
≤(1 + θk)E [ ‖zk−1 − z?(xk−1, yk−1)‖2 ] + (1 + θ−1k )E [ ‖z?(xk, yk)− z?(xk−1, yk−1)‖2 ] (ii)
≤ (1 + θk) ( Πk−1E z k−1 +R z k−1 ) + 2θ−1k E [ ‖z?(xk, yk)− z?(xk−1, yk−1)‖2 ] (iii)
≤ (1 + θk) ( Πk−1E z k−1 +R z k−1 ) + 4L2zθ −1 k ( E [ ‖xk − xk−1‖2 + ‖yk − yk−1‖2 ]) (iv)
≤ (1 + θk) ( Πk−1E z k−1 +R z k−1 ) + 4L2zθ −1 k ( γ2kE [∥∥∥ψ̂k−1∥∥∥2]+ 2φkEyk + 2R̃yk) (34a) (i) follows by Young’s inequality, (ii) uses Proposition 7 to bound the first term and that (1+θ−1k ) ≤ 2θ−1k for the second term, (iii) uses that z ?(x, y) is Lz-Lipschitz in x and y by Proposition 6 and,
finally, (iv) uses the update equation xk = xk−1 − γkψ̂k−1 for the term E [ ‖xk − xk−1‖2 ] and
Proposition 12 to control the increments E [ ‖yk − yk−1‖2 ] .
In order to express the upper-bound on Ezk in terms of E y k−1 instead of E y k , we substitute E y k in (34a) by its upper-bound in (33a) and use that (1 + rk) ≤ 2 to write:
Ezk ≤(1 + θk) ( Πk−1E z k−1 +R z k−1 ) + 4L2zθ −1 k γ 2 k ( 1 + 4L2yφkr −1 k ) E [∥∥∥ψ̂k−1∥∥∥2]
+ 8L2zθ −1 k ( 2φk ( Λk−1E y k−1 +R y k−1 ) + R̃yk ) (35)
We can then express (33a) and (35) jointly in matrix form as follows:( Eyk Ezk ) ≤ Pk ( Λk−1E y k−1 +R y k−1 Πk−1E z k−1 +R z k−1 ) + γkE [∥∥∥ψ̂k−1∥∥∥2]Uk + Vk where the Pk is a 2× 2 matrix and Uk and Vk are 2-dimensional vectors given by (17). The desired
result follows directly by substituting E [∥∥∥ψ̂k−1∥∥∥2] by its upper-bound from Proposition 12 in the
above inequality.

heading: C.4 GENERAL ERROR BOUND
text: Proof of Proposition 9. First note that, by assumption, we have that δkr−1k ≤ δk−1r −1 k−1 and δkθ −1 k ≤ δk−1θ −1 k−1. Moreover, since αk and βk are non-increasing, we also have that Λk−1 ≤ Λk and Πk−1 ≤ Πk. This implies the following inequalities which will be used in the rest of the proof:
a−1k−1Λk−1 ≤ a −1 k Λk, b −1 k−1Πk−1 ≤ b −1 k Πk. (36)
Now, let ρk be a non-increasing sequence with 0 < ρk < 1. By Proposition 2, it follows that Exk satisfies the inequality:
Fk + E x k ≤
( 1− ( 1− 1
2 ρk
) δk ) Exk−1 + γkskV ψ k−1 + γk ( sk + ρ −1 k ) Eψk−1 (37)
On the other hand, by Proposition 5 we know that Eyk and E z k satisfy:(
Eyk Ezk
) ≤ Pk ( Λk−1E y k−1 +R y k−1
Πk−1E z k−1 +R z k−1
) + γk ( V ψk−1 + 2E ψ k−1 + 2ζkE x k−1 ) Uk + Vk (38)
where the Pk, Uk and Vk are defined in (17). For conciseness, define Sk and EIk to be:
Sk := ( ak 0 0 bk ) , EIk = Sk ( Eyk Ezk ) By (36), we directly have that:
SkPkS −1 k−1
( Λk−1 0
0 Πk−1
) ≤ SkPkS−1k ( Λk 0 0 Πk ) := P̃k, (39)
where the inequality in (39) holds component-wise. Therefore, multiplying (38) by Sk and using (39) yields:
EIk ≤P̃kE I k−1 + Sk ( Pk ( Ryk−1 Rzk−1 ) + γk ( V ψk−1 + 2E ψ k−1 + 2ζkE x k−1 ) Uk + Vk ) (40)
Furthermore, by Proposition 4 we can bound Eψk−1 and V ψ k−1 as follows:
Eψk−1 ≤2L 2 ψ ( Λk(ak) −1 ak−1E y k−1 + Πk(bk) −1 bk−1E z k−1 ) + 2L2ψR y k−1
V ψk−1 ≤w 2 x + σ 2 xΠk(bk) −1 bk−1E z k−1,
where we used (36) a second time to replace Λk−1(ak−1) −1 and Πk−1(bk−1) −1 by Λk(ak) −1 and Πk(bk) −1. By summing both inequalities (37) and (40) and substituting all terms Eψk−1 and V ψ k−1 by their upper-bounds we obtain an inequality of the form:
Fk + E tot k ≤AxkExk−1 +A y kak−1E y k−1 +A z kbk−1E z k−1 + V tot k
where Axk , A y k , A z k are the components of the vector Ak defined in (18) and V tot k is the variance term also defined in (18). The desired inequality follows by upper-bounding Axk , A y k , A z k by their maximum value ‖Ak‖∞ .

heading: D CONTROLLING THE PRECISION OF THE INNER-LEVEL ALGORITHMS.
text: In this section, we prove Proposition 10. To achieve this, we first provide general conditions on Λk and Πk for controlling the rate ‖Ak‖∞ and which hold regardless of the choice of step-sizes. This is achieved in Proposition 13 of Appendix D.1. Then we prove Proposition 10 in Appendix D.2.
D.1 CONTROLLING Πk AND Λk .
We introduce the following quantities:
D (1) k :=
1
1− s
log  1− (1− ρk)δk
1 + 2rk [ 1 + 2L2ψγkδ −1 k [ sk + 1 + [2ρk]
−1 ]]  (41a)
D (2) k :=
1 s log ( (16L2y) −1ρkζ −1 k γ −2 k r 2 k ) (41b)
D (3) k :=−
1 s log ( 4L2yδkγkr −2 k ) (41c)
D (4) k :=
1
1− s log  1− (1− ρk)δk 1 + θk [ 1 + 2γkδ −1 k ( 2L2ψ + σ 2 x )[ sk + 1 + (2ρk) −1 ]]  (41d)
D (5) k :=−
1 s log ( 16L2zθ −2 k φk ) (41e)
D (6) k :=−
1 s log ( 2L2zL −2 y θ −2 k r 2 k ( 1 + 4L2yr −1 k φk )) (41f)
Proposition 13. Let ρk be a non-increasing sequence of positive numbers smaller than 1. Consider Λk and Πk so that:
log Λk ≤min ( D (1) k , D (2) k , D (3) k ) , (42a)
log Πk ≤min ( D (1) k , D (2) k , log Λk +D (3) k ) . (42b)
Then, the following inequalities holds:
‖Ak‖∞ ≤ (1− (1− ρk)δk), V tot k ≤ Ṽ totk .
where Ak and V totk are defined in (18) of Proposition 9 and Ṽ tot k is defined as: Ṽ totk :=δk ( 3r−1k + 2L 2 ψγkδ −1 k ( 2 + (sk + ρ −1 k ) )) Ryk−1
+ δkΠ s k ( 2θ−1k R z k−1 + 8L 2 zθ −2 k R̃ y k ) + γk ( sk + 4L 2 yΛ s kδkγkr −2 k ) w2x
Proof. We first prove that uIk≤u + k ≤1 and p + k ≤1 under (42a) and (42b) with u + k , p + k given by:
u+k := 4L 2 yδkγkr −2 k Λ s k, p + k = 16L 2 zΠ s kθ −2 k φk.
A direct calculation shows u+k ≤1 whenever (42a) holds. Moreover, recall that uIk=akU (1) k +bkU (2) k with U (1)k and U (2) k being the components of the vector Uk defined in (17). Thus by direct substitution, we get the following expression for uIk:
uIk = 2L 2 yδkγkr −2 k Λ s k ( 1 + 2L2zL −2 y θ −2 k r 2 k ( 1 + 4L2yr −1 k φk )Πsk Λsk ) .
Therefore, (42b) suffices to ensure that uIk ≤ u + k . Finally, (42b) implies directly that p + k ≤ 1.
We will control each component Axk , A y k and A z k of the vector Ak separately.
Controlling Axk . Recalling the expression of Axk , the first component of Ak in (18), it holds that: Axk = 1− (
1− 1 2 ρk
) δk + 2ζkλku I k (i) ≤1− (
1− 1 2 ρk
) δk + 2ζkγku + k
(ii) ≤ 1− (
1− 1 2 ρk
) δk + 1
2 ρkδk = 1− (1− ρk)δk.
(i) holds since uIk≤u + k while (ii) follows from (42a) which ensures that 2ζkγku + k ≤ 1 2ρkδk.
Controlling Ayk. Recall the expression of second component of A y k in (18), we have:
Ayk =Λ 1−s k ( 1 + rk ( 1 + 16L2zΠ s kθ −2 k φk + 2L 2 ψγkδ −1 k ( 2uIk + (sk + ρ −1 k ) )))
(i) ≤ (
1 + 2rk ( 1 + 2L2ψγkδ −1 k (sk + 1 + (2ρk) −1 ) )) Λ1−sk (ii) ≤ 1− (1− ρk)δk.
(i) holds since p+k := 16L 2 zΠ s kθ −2 k φk ≤ 1 and uIk ≤ 1 while (ii) is a consequence of (42a).
Controlling Azk. Similarly, recalling the expression of the third component of Ak we get that: Azk =Π 1−s k ( 1 + θk ( 1 + γkδ −1 k ( 2L2ψ + σ 2 x )( 2uIk + (sk + ρ −1 k ) )))
(i) ≤Π1−sk ( 1 + θk ( 1 + 2γkδ −1 k ( 2L2ψ + σ 2 x )( 1 + (sk + (2ρk) −1) ))) (ii) ≤ 1− (1− ρk)δk,
where (i) uses that uIk ≤ 1 and (ii) follows from (42b). Controlling V totk . Recalling the expression of V totk from (18), we have that:
V totk :=δk ( Λsk(1 + r −1 k ) + 16L 2 zφkθ −2 k Π s k + 2L 2 ψγkδ −1 k ( 2uIk + (sk + ρ −1 k ) )) Ryk−1
+ δk ( (1 + θ−1k )Π s kR z k−1 + 8L 2 zθ −2 k Π s kR̃ y k ) + γk ( sk + u I k ) w2x
≤δk ( 3r−1k + 2L 2 ψγkδ −1 k ( 2 + (sk + ρ −1 k ) )) Ryk−1
+ δk ( 2θ−1k Π s kR z k−1 + 8L 2 zθ −2 k Π s kR̃ y k ) + γk ( sk + u + k ) w2x = Ṽ tot k .
where we use 16L2zφkθ −2 k Π s k ≤ 1 and uIk ≤ 1 for the first line and uIk ≤ u + k for the last line.

heading: D.2 CONTROLLING THE NUMBER OF INNER-LEVEL ITERATIONS
text: We provide now a proof of Proposition 10 which is a consequence of Proposition 13.
Proof of Proposition 10. We first provide conditions on the number of iterations T and N of algorithms Ak and Bk to control the rate ‖Ak‖∞ and then provide an upper-bound on V totk . Conditions on T and N . We consider the setting with constant step-size γk=γ, αk=α and βk=β and choose rk=θk = 1 and δk=δ0 for some 0<δ0<1. We also take v=1 so that φk=2 and R̃ y k=R y k. By direct substitution of the parameters rk, θk, φk, γk, δk, ρk and ζk, in the expressions of D (1) k , D (2) k , D (3) k , D (4) k , D (5) k and D (6) k defined in (41a) to (41f), we verify that:
−C1 ≤ D(1)k , −C ′ 1 ≤ D (4) k , −C2 ≤ min
( D
(2) k , D (3) k ) , −C ′2 ≤ min ( D (5) k , D (6) k ) .
Hence, we can ensure the conditions of Proposition 13 hold by choosing T and N so that: log Λk ≤ −max (1, C1, C2), log Πk ≤ − log Λk −max (1, C ′1, C ′2).
This is achieved by for the following choice: T = bα−1µ−1g max (C1, C2, C3)c+ 1, (43) N = b2β−1µ−1g (max (C1, C2, C3) + max (C ′1, C ′2, C ′3))c+ 1 Hence, for such choice, we are guaranteed by Proposition 13 that ‖Ak‖∞ ≤ 1− (1− ρk)δk. Bound on the variance V totk . By choosing T and N as in (43), we know that Λk and Πk satisfy (42), so that the variance term V totk is upper-bounded by Ṽ tot k . Moreover, by direct substitution of the sequences appearing in the expression of Ṽ totk by their values, we get:
V totk ≤ Ṽ totk =δk ( Λsk(1 + r −1 k ) + 16L 2 zφkθ −2 k Π s k + 2L 2 ψη −1 k ( 2uIk + (sk + ρ −1 k ) )) Ryk−1
+ δk ( (1 + θ−1k )Π s kR z k−1 + 8L 2 zθ −2 k Π s kR̃ y k ) + γk ( sk + u I k ) w2x
≤δ0 ( η−10
1− u 2 +
( 1 + u
2 γ + 4L2yΛ s kγ 2
)) w2x
+ δ0 ( 2Λsk + 32L 2 zΠ s k + 10L 2 ψη −1 0 ) Ryk−1 + δ0 ( 2ΠskR z k−1 + 8L 2 zΠ s kR̃ y k )
Furthermore, by definition of Ryk−1 and R z k−1, we have that:
Ryk−1 ≤ 2µ −1 g L −1 g σ 2 g , R z k−1 ≤ µ−3g ( B2L−1g σ 2 gyy + 3µgσ 2 f ) .
Moreover, recall that R̃yk=R y k since we chose v=1. Thus R̃ y k ≤ 2µ−1g ασ2g . This implies that:
V totk δ −1 0 γ −1 ≤ ( δ−10
1− u 2
+ ( 1 + 4L2yΛ s kγ )) w2x + 2Π s kγ −1µ−3g ( B2βσ2gyy + 3µgσ 2 f ) (44)
+ ( 4Λsk + 80L 2 zΠ s k + 20L 2 ψη −1 0 ) µ−1g αγ −1σ2g
By choosing T and N as in (43), the following conditions hold:
Λsk ≤ L 4L2y , Λsk ≤ 5L2ψη−10 , Πsk ≤ γ
( σ2gxy + (L ′ g) 2 ) , Πsk ≤ 1
4L2z L2ψη −1 0 .
By applying these inequalities in (44), we get:
V totk δ −1 0 γ −1 ≤ ( δ−10
1− u 2 + 2
) w2x + 2µ −3 g ( σ2gxy + (L ′ g) 2 )( B2βσ2gyy + 3µgσ 2 f ) + 60L2ψη −1 0 µ −1 g L −1 g γ −1σ2g
≤ ( δ−10
1− u 2 + 3
) w2x +
60L2ψ δ0µgLg σ2g =W2
where we used that 2µ−3g ( σ2gxy + (L ′ g) 2 )( B2βσ2gyy + 3µgσ 2 f ) ≤ w2x by definition of w2x in (16a). Therefore, we have shown that V totk ≤ γδ0W2, withW2 given by (20).

heading: E STOCHASTIC LINEAR DYNAMICAL SYSTEM WITH CORRELATED NOISE
text: Let A be a positive definite matrix in Rd × Rd satisfying 0 < µg ≥ |σi(A)| ≤ Lg and b a vector in Rd. We denote by z? = −A−1b. Consider Am be a sequence of i.i.d. positive symmetric matrices in Rd × Rd such that E[Am] = A, and b̂ a random vector in Rd such that E [ b̂ ] = b, with Am and
b̂ being mutually independent. Define ΣA = E [ (An −A)>(An −A) ] and denote by σA and LA
the largest singular values of ΣA and A−1ΣAA−1. Let β be such that β ≤ 1Lg . Finally let σ 2 c be an upper-bound on E [∥∥∥b̂− b∥∥∥2]. Let z and z′ be two vectors in Rd and define the iterates zn and z̄n
such that z0 = z and z̄0 = z′ and using the recursion:
zn = (I − βAn)zn−1 − βb̂, z̄n = (I − βA)z̄n−1 − βb.
Hence, from the definition of zn and z̄n we directly have that:
zn = n∏ t=1 (I − βAn)z − n∑ t=1 n∏ j=t+1 β(I − βAj)b̂, z̄n = n∏ t=1 (I − βA)z′ − n∑ t=1 n∏ j=t+1 β(I − βA)b.
The next proposition computes the bias E [ zn − z̄n ∣∣∣b̂]. Proposition 14. The following identities hold:
E [ zn − z̄n ∣∣∣b̂] =(I − βA)n(z − z′)− β( n∑ t=1 (I − βA)n−t )( b̂− b )
E[zn − z̄n] =(I − βA)n(z − z′)
Proof. The proof is a consequence of An and b̂ being i.i.d. and unbiased estimates of A and b.
The next proposition controls the mean squared errors E [ ‖zn − z?‖2 ] and E [ ‖zn − z̄n‖2 ] .
Proposition 15. Define n?(n, β) = min(n, 1βµg ). Let β such that:
β ≤ 1 2Lg min
( 1,
2Lg µg ( 1 + µ−2g σ2A )) Then, the following inequalities holds:
E [ ‖zn − z?‖2 ] ≤(1− βµg)n‖z − z?‖2 + β2 ( σ2A‖z?‖ 2 + 3 ( n ∧ 1
βµg
) σ2c )( n ∧ 1
βµg
) ,
‖z̄n − z?‖2 ≤(1− βµg)n‖z′ − z?‖ 2 .
Moreover, if z′ = z, then we have:
E [ ‖zn − z̄n‖2 ] ≤4µ−2g σ2A(1−
βµg 2 )n‖z − z?‖2 (45)
+ 2β2 ( σ2A‖z?‖ 2 + 3 ( n ∧ 1
βµg
) σ2c )( n ∧ 1
βµg
) .
Proof. It is straightforward to see that:
‖z̄n − z?‖2 ≤ (1− βµg)‖z′ − z?‖ 2 . Now, let’s control E [ ‖zn − z̄n‖2 ] . The following identity holds by definition of zn and z̄n:
E [ ‖zn − z̄n‖2 ] =E [( zn−1 − z̄n−1 )>( (I − βA)2 + β2ΣA )( zn−1 − z̄n−1 )] − 2βE [ (zn−1 − z̄n−1)>(I − βA)(b̂− b)
] + β2 ( E [∥∥∥b̂− b∥∥∥2]+ E[(z̄n−1)>ΣAz̄n−1])
=E [( zn−1 − z̄n−1 )>( (I − βA)2 + β2ΣA )( zn−1 − z̄n−1 )] + β2 (( z̄n−1 )> ΣAz̄ n−1 + E [( b̂− b )> (I + 2(I − βA)Dn) ( b̂− b
)]) where ΣA = E [ (An −A)>(An −A) ] and Dn = ∑n−1 t=0 (I − βA)t. By simple calculation we can upper-bound the last term by:
β2E [( b̂− b )> (I + 2(I − βA)Dn) ( b̂− b )] ≤ 3 ( n ∧ 1
βµg
) β2E [∥∥∥b̂− b∥∥∥2]. Moreover, provided that β ≤ 1Lg(1+LA) , where LA is the highest eigenvalue of A −1ΣAA −1, then we have the following:
E [ ‖zn − z̄n‖2 ] ≤(1− βµg)E [∥∥zn−1 − z̄n−1∥∥2] + β2 (( z̄n−1 )> ΣAz̄ n−1 + 3 ( n ∧ 1
βµg
) E [∥∥∥b̂− b∥∥∥2])
≤(1− βµg)E [∥∥zn−1 − z̄n−1∥∥2]+ β2(σ2A∥∥z̄n−1∥∥2 + 3(n ∧ 1βµg ) σ2c ) ,
Unrolling the recursion, it follows that: E [ ‖zn − z̄n‖2 ] ≤(1− βµg)n‖z − z′‖ 2 + β2 n∑ t=1 (1− βµg)n−t ( σ2A ∥∥z̄n−1∥∥2 + 3(n ∧ 1 βµg ) σ2c ) (46)
In particular, if z ′ = z?, then z̄n = z? and we get:
E [ ‖zn − z?‖2 ] ≤(1− βµg)n‖z − z?‖2 + β2 ( σ2A‖z?‖ 2 + 3 ( n ∧ 1
βµg
) σ2c )( n ∧ 1
βµg
) .
To get the last inequality, we simply choose z′ = z and recall that:∥∥z̄t−1∥∥ ≤ 2((1− βµg)t−1‖z′ − z?‖+ ‖z?‖). Using the above in in (46) yields:
E [ ‖zn − z̄n‖2 ] ≤2nβ2σ2A(1− βµg)n−1‖z − z?‖ 2 + 2β2 ( σ2A‖z?‖ 2 + 3 ( n ∧ 1
βµg
) σ2c )( n ∧ 1
βµg
) .
Moreover, by Lemma 3 we know that nβ2µ2g(1−βµg)n−1 ≤ (1− βµg 2 ) n−1 and since βµg ≤ 1, we have that (1− βµg2 ) −1 ≤ 2 so that (1− βµg2 ) n−1 ≤ 2(1− βµg2 ) n. Hence, we can write:
E [ ‖zn − z̄n‖2 ] ≤4µ−2g σ2A(1−
βµg 2
)n‖z − z?‖2 + 2β2 ( σ2A‖z?‖ 2 + 3 ( n ∧ 1
βµg
) σ2c )( n ∧ 1
βµg
) .
Lemma 2. Let A and ΣA be symmetric positive matrix in Rd × Rd with σ2A its largest singular value of ΣA and 0 < µg ≤ σi(A) ≤ Lg . Let β be a positive number such that:
β ≤ 1 2Lg min
( 1,
2Lg µg ( 1 + µ−2g σ2A )) Then the following holds: ∥∥(I − βA)2 + β2ΣA∥∥op ≤ 1− βµg. Proof. First note that β ≤ 1Lg , so that I − βA is positive. Now, we observe that∥∥(I − βA)2 + β2ΣA∥∥op ≤ (1 − βµg)2 + β2σ2A which holds since I − βA is positive. And since β ≤ µg
µ2g+σ 2 A
, we further have (1− βµg)2 + β2σ2A ≤ 1− βµg , which yields the desired result.
Lemma 3. Let 0 ≤ b < 1 and n ≥ 1, then the following inequality holds:
nb2(1− b)n−1 ≤ (1− b 2 )n−1.
Proof. We consider the function h(n, b) defined by:
h(n, b) := (n− 1) log ( 1− b2 1− b ) − log(nb2).
We need to show that h(n, b) is non-negative for any n ≥ 1 and 0 ≤ b < 1. For this purpose, we fix b and consider the variations of h(n, b) in n:
∂nh(n, b) = log ( 1− b2 1− b ) − 1 n .
∂nh(n, b) is non-negative for n ≥ n? = log ( 1− b2 1−b )−1 and non-positive for all n ≤ n?. Hence, h(n, b) achieves its minimum value in n? over the (0,+∞). We distinguish two case depending on whether n? is greater of smaller than 1.
Case n? ≤ 1. In this case n 7→ h(n, b) is increasing on the interval [1,+∞) since ∂nh(n, b) ≥ 0 for n ≥ n?. Hence, h(n, b) ≥ h(1, b) for all n ≥ 1. Moreover, since h(1, b) = − log(b2) ≥ 0 the result follows directly.
Case n? > 1. In this case we still have h(n, b) ≥ h(n?, b) for all n ≥ 1, since n? achieves the minimum value of h. Thus we only need to show that h(n?, b) ≥ 0. Using the expression of n?, we have:
h(n?, b) =1− log ( 1− b2 1− b ) − log ( n?b2
) =1− 1
n? − log
( n?b2 ) Since n? > 1, the first term 1− 1n? is non-negative, thus we only need to show that n
?b2 ≤ 1 so that the last term is also non-negative. It is easy to see that n?b2 ≤ 1 is equivalent to having h̃(b) ≥ 0, where we define the function h̃(b) as:
h̃(b) = log 1− b2 1− b − b2.
We can analyze the variations of b̃ be computing its derivative which is given by:
∂bh̃(b) = 1
(1− b)(2− b) − 2b.
Hence, we have the following equivalence:
∂bh̃(b) ≥ 0 ⇐⇒ 2b(1− b)(2− b) ≤ 1
This is always true for 0 ≤ b < 1 since b(1 − b) ≤ 14 so that 2b(1 − b)(2 − b) ≤ 2−b 2 ≤ 1. Thus we have shown that h̃ is increasing over [0, 1) so that h̃(b) ≥ h̃(0) = 0. As discussed above, this is equivalent to having n?b2 ≤ 1, so that h(n?, b) ≥ 0 which concludes the proof.

heading: F EXPERIMENTS

heading: F.1 DETAILS OF THE SYNTHETIC EXAMPLE
text: We choose the functions f and g to be of the form: f(x, y) := 12x >Afx + y >Cf and g(x, y) := 1 2y >Agy + y
>Bgx where Af and Ag are symmetric definite positive matrices of size dx × dx and dy × dy , Bg is a dy × dx matrix and Cf is a dy vector with dx = 2000 and dy = 1000. We generate the parameters of the problem so that the smoothness constants L and Lg are fixed to 1, κL=10 and κg taking values in {10i, i ∈ {0, .., 7}}. We then solve each problem using different methods and perform a grid-search on the number of iterations T and M of algorithms Ak and Bk . We fix the step-sizes to γk=1/L and αk=βk=1/Lg and perform a grid-search on the number of iterations T and M of algorithms Ak and Bk from {10i, i ∈ 0, 1, 2, 3}. For AID methods without warm-start in Bk, we consider an additional setting where M increases logarithmically with k, as suggested in Ji et al. (2021), with M=b103 log(k)c. Similarly, for (ITD) and (Reverse), we additionally use an increasing T of the same form.

heading: F.2 EXPERIMENTAL DETAILS FOR LOGISTIC REGRESSION
text: The inner-level and outer-level cost functions for such task take the following form:
f(x, y) = 1 |Dval| ∑
ξ∈Dval
L(y, ξ), g(x, y) = 1 |Dtr| ∑ ξ∈Dtr L(y, ξ) + 1 pd ∑ i=1 exp(xi)‖y.,i‖2
For the default setting, we use the well-chosen parameters reported in Grazzi et al. (2020); Ji et al. (2021) where αk=γk=100, βk=0.5, and T=N=10. For the grid-search setting, we select the best performing parameters T , M and βk from a grid {10, 20} × {5, 10} × {0.5, 10}, while the batchsize (chosen to be the same for all steps of the algorithms) varies from 10 ∗ {0.1, 1, 2, 4}. We also compared with VRBO (Yang et al., 2021) using the implementation available online and noticed instabilities for large values of T and N , as reported by the authors, but also a drop in performance compared to stocBiO for smaller T and N due to inexact estimates of the gradient.

heading: F.3 DATASET DISTILLATION
text: Dataset distillation (Wang et al., 2018; Lorraine et al., 2020) consists in learning a small synthetic dataset such that a model trained on this dataset achieves a small error on the training set. Specifically, we consider a classification problem of C classes using a linear model and a training dataset Dtr where each training point ξ ∈ Dtr is a d-dimensional vector with a class cξ ∈ {1, ..., C}. The linear model is represented by a matrix y ∈ Rc×d multiplying a data point yξ and providing the logits of each class. The dataset distillation can be cas as a bilevel problem of the form:
min x∈Rc×d,λ∈Rd
1 |Dtr| ∑ ξ∈Dtr CE(y?(x, λ)ξ, cξ),
y?(x, λ) ∈ arg min y∈Rc×d
1
C C∑ c=1 CE(yxc, c) + 1 Cd d∑ i=1 exp(λi)‖y.,i‖2.
where λ ∈ Rd is a vector of hyper-parameter for regularizing the inner-level problem which we found beneficial to add.
Experimental setup. We perform the distillation task on MNIST dataset. We set the stepsizes αk=βk=0.1 and T=N=10. We perform a grid-search on the outer-level step-size γk∈{0.01, 0.001, 0.0001} and run the algorithms for k = 40000 iterations.

Review summaries:
The paper studies the problem of bi-level optimization. In particular, it considers algorithms based on inexact implicit differentiation, where the inner problem and the implicit differentiation are not solved exactly. Warm-up is used when solving the inner problem and implicit differentiation. The convergence of the proposed method is analyzed by viewing the iterates of the proposed method as a dynamical system using the idea of singularly perturbed systems. 

This work focuses on bilevel optimization. The key innovation here is the warm start, which enables improved complexity bounds (under some settings). The analysis nicely builds on singularly perturbed systems (SPS). Numerical tests are also provided on both synthetic and real problems, where the merit of warm start is demonstrated.