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

Mathieu Dagréou111Corresponding author<EMAIL_ADDRESS>222Université
Paris-Saclay, Inria, CEA, Palaiseau, 91120, France , Thomas Moreau†, Samuel
Vaiter333CNRS & Université Côte d’Azur, LJAD, Nice, France , Pierre
Ablin444Apple

###### Abstract

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

##  1 Introduction

In the last few years, bilevel optimization has become an essential tool for
the machine learning community thanks to its numerous applications. Among
them, we can cite hyperparameter selection (Bengio, 2000; Pedregosa, 2016;
Franceschi et al., 2017; Lorraine et al., 2020), implicit deep learning (Bai
et al., 2019), neural architecture search (Liu et al., 2019; Zhang et al.,
2021), data augmentation (Li et al., 2020; Rommel et al., 2022) and meta-
learning (Franceschi et al., 2018; Rajeswaran et al., 2019) to name a few. In
bilevel optimization, we are interested in minimizing a function under the
constraint that one variable minimizes another function. This can be
formalized as follows

$\min_{x\in\mathbb{R}^{d}}h(x)=F(z^{*}(x),x),\quad\text{subject to
}z^{*}(x)\in\operatorname*{arg\,min}_{z\in\mathbb{R}^{p}}G(z,x)\enspace.$ (1)

The function $F$ is called the outer function and the function $G$ is the
inner function. Likewise, we refer to $z$ as the inner variable and $x$ as the
outer variable. A strategy to solve bilevel problems consists in using
implicit differentiation that provides the following expression for the
gradient of $h$

$\nabla h(x)=\nabla_{2}F(z^{*}(x),x)+\nabla_{21}^{2}G(z^{*}(x),x)v^{*}(x)$ (2)

where $v^{*}(x)$ is the solution of a linear system

$v^{*}(x)=-\left[\nabla^{2}_{11}G(z^{*}(x),x)\right]^{-1}\nabla_{1}F(z^{*}(x),x)\enspace.$
(3)

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

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

$F(z,x)=\frac{1}{m}\sum_{j=1}^{m}F_{j}(z,x),\quad
G(z,x)=\frac{1}{n}\sum_{i=1}^{n}G_{i}(z,x)\enspace.$

For instance, in hyperparameter selection, $F$ is the validation loss which is
an average on the validation set and $G$ is the training loss which is an
average on the training set. In single-level optimization, the finite sum
structure has been widely leveraged to produce fast first-order algorithms
that provably converge faster than gradient descent. These methods are the
cornerstone of many successful machine learning applications. Among these
algorithms, we can cite stochastic methods such as stochastic gradient descent
(Robbins and Monro, 1951; Bottou, 2010) and its variance-reduced variants such
as SAGA (Defazio et al., 2014), STORM (Cutkosky and Orabona, 2019) or
SPIDER/SARAH (Fang et al., 2018; Nguyen et al., 2017) that use only a handful
of samples at a time to make progress. In order to get faster methods than
full-batch approaches, it is natural to extend these methods to the bi-level
setting. The main obstacle comes from the difficulty of obtaining stochastic
approximations of $\nabla h(x)$ because of its structure (2). In the
literature, several strategies have been proposed to overcome this obstacle,
and some works demonstrate that stochastic implicit differentiation based
algorithms for solving (1) have the same complexity as single-level analogous
algorithms. For instance, ALSET from (Chen et al., 2021) and SOBA from Dagréou
et al. (2022) have the same convergence rate as stochastic gradient descent
for smooth nonconvex single-level problems Ghadimi and Lan (2013); Bottou et
al. (2018). Furthermore, Dagréou et al. (2022) show that SABA, an adaptation
of SAGA algorithm (Defazio et al., 2014), has a sample complexity in
$\mathcal{O}((n+m)^{\frac{2}{3}}\varepsilon^{-1})$ which is analogous to the
sample complexity of SAGA for nonconvex single-level problems (Reddi et al.,
2016).

However, in classical single-level optimization, it is known that neither of
these algorithms is optimal: the SARAH algorithm (Nguyen et al., 2017)
achieves a better sample complexity of
$\mathcal{O}(m^{\frac{1}{2}}\varepsilon^{-1})$ with $m$ the number of samples.
Furthermore, this algorithm is _near-optimal_ (i.e. optimal up to constant
factors), because the lower bound for single-level non-convex optimization is
also $\mathcal{O}(m^{\frac{1}{2}}\varepsilon^{-1})$ as proved by Zhou and Gu
(2019). It is natural to ask if we can extend these results to bilevel
optimization:

_Are the optimal complexity bounds for solving bilevel optimization the same
as in single-level optimization?_

#### Contributions

In Section 2, we introduce SRBA, an adaptation of the SARAH algorithm to the
bilevel setting. We then demonstrate in Section 3 that, similarly to the
single-level setting, it requires
$\mathcal{O}\left((n+m)^{\frac{1}{2}}\varepsilon^{-1}\vee(n+m)\right)$ calls
to oracles to reach an $\varepsilon$-stationary point. This is therefore an
upper bound on the complexity of solving bilevel empirical risk minimization
(ERM) problems. As shown in Table 1, it achieves the best-known complexity in
the regime $n+m\lesssim\mathcal{O}(\varepsilon^{-2})$. In Section 4, we
analyze the lower bounds for such problems. We demonstrate that the number of
iterations required to reach an $\varepsilon$-stationary point (see Definition
3.1) is at least $\Omega(m^{\frac{1}{2}}\varepsilon^{-1})$, hereby matching
the previous upper-bound in the case where $n\asymp m$ and $\varepsilon\leq
m^{-\frac{1}{2}}$. SRBA is therefore near-optimal in that regime. Even though
our main contribution is theoretical, we illustrate the numerical performances
of the algorithm in Section 5.

| Sample complexity | Stochastic setting | $F$ | $G$
---|---|---|---|---
StocBiO Ji et al. (2021) | $\tilde{\mathcal{O}}(\varepsilon^{-2})$ | General expectation | $\mathcal{C}^{1,1}_{L}$ | SC and $\mathcal{C}^{2,2}_{L}$
AmIGO Arbel and Mairal (2022) | $\mathcal{O}(\varepsilon^{-2})$ | General expectation | $\mathcal{C}^{1,1}_{L}$ | SC and $\mathcal{C}^{2,2}_{L}$
MRBO Yang et al. (2021) | $\tilde{\mathcal{O}}(\varepsilon^{-\frac{3}{2}})$ | General expectation | $\mathcal{C}^{1,1}_{L}$ | SC and $\mathcal{C}^{2,2}_{L}$
VRBO Yang et al. (2021) | $\tilde{\mathcal{O}}(\varepsilon^{-\frac{3}{2}})$ | General expectation | $\mathcal{C}^{1,1}_{L}$ | SC and $\mathcal{C}^{2,2}_{L}$
SABA (Dagréou et al., 2022) | $\mathcal{O}((n+m)^{\frac{2}{3}}\varepsilon^{-1})$ | Finite sum | $\mathcal{C}^{2,2}_{L}$ | SC and $\mathcal{C}^{3,3}_{L}$
SRBA | $\mathcal{O}((n+m)^{\frac{1}{2}}\varepsilon^{-1})$ | Finite sum | $\mathcal{C}^{2,2}_{L}$ | SC and $\mathcal{C}^{3,3}_{L}$
Table 1: Comparison between the sample complexities and the Assumptions of
some stochastic algorithms for bilevel optimization. It corresponds to the
number of calls to gradient, Hessian-vector products and Jacobian-vector
product sufficient to get an $\varepsilon$-stationary point. The tilde on the
$\tilde{\mathcal{O}}$ hide a factor $\log(\varepsilon^{-1})$. ”SC” means
”strongly-convex”. $\mathcal{C}^{p,p}_{L}$ means $p$-times differentiable with
Lipschitz $k$th order derivatives for $k\leq p$.

#### Related work

There are several strategies to solve (1) in a stochastic fashion. They can be
separated into two groups: iterative differentiation algorithms (ITD) and
approximate implicit differentiation algorithms (AID). On the one hand, in ITD
algorithms, the Jacobian of $z^{*}$ is estimated by differentiating the
different steps used to compute an approximation of $z^{*}$. On the other
hand, AID algorithms leverage the implicit gradient given by (2) replacing
$z^{*}$ and $v^{*}$ by some approximations $z$ and $v$. In the class of ITD
algorithms, Maclaurin et al. (2015) propose to approximate the Jacobian of the
solution of the inner problem by differentiating through the iterations of SGD
with momentum. The complexity of the hypergradient computation in ITD solvers
is studied in Franceschi et al. (2017); Grazzi et al. (2020); Ablin et al.
(2020). For AID algorithms, Ghadimi and Wang (2018); Chen et al. (2021); Ji et
al. (2021) propose to perform several SGD steps in the inner problem and then
use Neumann approximations to approximate $v^{*}(x)$ defined in (3). A method
consisting of alternating steps in the inner and outer variables was proposed
in Hong et al. (2021). These methods can be improved by using a warm start
strategy for the inner problem Ji et al. (2021); Chen et al. (2021) and for
the linear system Arbel and Mairal (2022). Some works elaborate on these ideas
by adapting variance reduction methods like STORM Khanduri et al. (2021); Yang
et al. (2021) or SAGA Dagréou et al. (2022). We take a similar approach and
extend the SARAH variance reduction method to the bilevel setting. Finally,
recent works propose to approximate the Jacobian of $z^{*}$ by stochastic
finite difference Sow et al. (2022) or Bregman divergence-based methods Huang
et al. (2022).

In single-level optimization, the problem of finding complexity lower bound
for optimization problems has been widely studied since the seminal work of
Nemirovsky and Yudin (1983). On the one hand, Agarwal and Bottou (2015)
provided a lower bound to minimize strongly convex and smooth finite sum with
deterministic algorithms that have access to individual gradients. These
results were extended to randomized algorithms for (strongly) convex and
eventually nonsmooth finite sum objective by Woodworth and Srebro (2016). On
the other hand, Carmon et al. (2017b) provided a lower bound for minimizing
nonconvex functions with deterministic and randomized algorithms. The
nonconvex finite sum case is treated in Fang et al. (2018); Zhou and Gu
(2019). In the bilevel case, Ji and Liang (2023) showed a lower bound for
deterministic, full-batch algorithms. However, this result is restricted to
the case where the value function $h$ is convex or strongly convex, which is
not the case with most ML-related bilevel problems. Our results are instead in
a non-convex setting.

#### Notation

The quantity $A_{\bullet}$ refers to $A_{z}$, $A_{v}$, or $A_{x}$, depending
on the context. If $f:\mathbb{R}^{p}\times\mathbb{R}^{d}\to\mathbb{R}$ is a
twice differentiable function, we denote $\nabla_{i}f(z,x)$ its gradient
w.r.t. its $i^{\text{th}}$ variable. Its Hessian with respect to $z$ is
denoted $\nabla_{11}^{2}f(z,x)~{}\in~{}\mathbb{R}^{p\times p}$ and its cross
derivative matrix $\left(\frac{\partial^{2}f}{\partial z_{i}\partial
x_{j}}\right)_{\begin{subarray}{c}i\in[p]\\\ j\in[d]\end{subarray}}$ is
denoted $\nabla^{2}_{12}f(z,x)~{}\in~{}\mathbb{R}^{p\times d}$. We denote
$\Pi_{\mathcal{C}}$ the projection on a closed convex set $\mathcal{C}$.

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

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

###  2.1 Assumptions

Before presenting our algorithm, we formulate several Assumptions on the
functions $F$ and $G$. As for SARAH, the regularity assumptions are made on
the individual functions $(G_{i})_{1\leq i\leq n}$ and $(F_{j})_{1\leq j\leq
m}$ rather than on the empirical means $G$ and $F$. In 2.1 and 2.2, we state
the regularity needed on the outer function $F$ and inner function $G$
respectively.

######  Assumption 2.1 (Regularity of $F$).

For all $j\in[m]$, the function $F_{j}$ is twice differentiable. The function
$F_{j}$ is $L^{F}_{0}$-Lipschitz continuous, its gradient $\nabla F_{j}$ is
$L^{F}_{1}$-Lipschitz continuous and the Hessian $\nabla^{2}F_{j}$ is
$L^{F}_{2}$-Lipschitz continuous.

######  Assumption 2.2 (Regularity of $G$).

For all $i\in[n]$, The function $G_{i}$ is three times differentiable. Its
first, second, and third order derivatives are respectively
$L^{G}_{1}$-Lipschitz continuous, $L^{G}_{2}$-Lipschitz continuous, and
$L^{G}_{3}$-Lipschitz continuous. For any $x\in\mathbb{R}^{d}$, the function
$G_{i}(\,.\,,x)$ is $\mu_{G}$-strongly convex.

The strong convexity and the smoothness with respect to $z$ hold for instance
when we consider an $\ell^{2}$-regularized logistic regression problem with
non-separable data. These regularity assumptions up to first-order for $F$ and
second-order for $G$ are standard in the stochastic bilevel literature (Arbel
and Mairal, 2022; Ji et al., 2021; Yang et al., 2021). The second-order
regularity for $F$ and third-order regularity for $G$ are necessary for the
analysis of the dynamics of $v$, as it is the case in Dagréou et al. (2022).

As shown in Ghadimi and Wang (2018, Lemma 2.2), these assumptions are
sufficient to get the smoothness of $h$, which is a fundamental property to
get a descent on $h$.

######  Proposition 2.3 (Smoothness of the value function).

Under Assumptions 2.1 and 2.2, the function $h$ is $L^{h}$ smooth for some
$L^{h}>0$ which is precised in Section A.2.

Another consequence of Assumptions 2.1 and 2.2 is the boundedness of the
function $v^{*}$.

######  Proposition 2.4 (Boundedness of $v^{*}$).

Assume that Assumptions 2.1 and 2.2 hold. Then, for
$R=\frac{L^{F}_{0}}{\mu_{G}}$ it holds that for any $x\in\mathbb{R}^{d}$, we
have $\|v^{*}(x)\|\leq R$.

In what follows we denote $\Gamma$ the closed ball centered in $0$ with radius
$R$ and $\Pi_{\Gamma}$ the projection onto $\Gamma$. Moreover, for
$(z,v,x)\in\mathbb{R}^{p}\times\mathbb{R}^{p}\times\mathbb{R}^{d}$, we denote
$\Pi(z,v,x)=(z,\Pi_{\Gamma}(v),x)$.

###  2.2 Hypergradient Approximation

The gradient of $h$ given by (2) is intractable in practice because it
requires the perfect knowledge of $z^{*}(x)$ and $v^{*}(x)$ which are usually
costly to compute, for instance when the inner problem is ill-conditioned. As
classically done in the stochastic bilevel literature (Ji et al., 2021; Arbel
and Mairal, 2022; Li et al., 2022), $z^{*}(x)$ and $v^{*}(x)$ are replaced by
approximate surrogate variables $z$ and $v$. The variable $z$ is typically the
output of one or several steps of an optimization procedure applied to
$G(\,.\,,x)$. The variable $v$ can be computed by using Neumann approximations
or doing some optimization steps on the quadratic
$v\mapsto\frac{1}{2}v^{\top}\nabla^{2}_{11}G(z,x)v+\nabla_{1}F(z,x)^{\top}v$.
We consider the approximate hypergradient given by

$D_{x}(z,v,x)=\nabla^{2}_{21}G(z,x)v+\nabla_{2}F(z,x)\enspace.$

The motivation behind this direction is that if we take $z=z^{*}(x)$ and
$v=v^{*}(x)$, we recover the true gradient, that is
$D_{x}(z^{*}(x),v^{*}(x),x)=~{}\nabla h(x)$. Proposition 2.5 from (Dagréou et
al., 2022, Lemma 3.4) controls the hypergradient approximation error with the
distances between $z$ and $z^{*}(x)$ and between $v$ and $v^{*}(x)$.

######  Proposition 2.5 (Hypergradient approximation error).

Let $x~{}\in~{}\mathbb{R}^{d}$. Assume that $F$ is differentiable and
$L^{F}_{1}$ smooth with bounded gradient, $G$ is twice differentiable with
Lipschitz gradient and Hessian and $G(\,.\,,x)$ is $\mu_{G}$-strongly convex.
Then there exists a constant $L_{x}$ such that

$\|D_{x}(z,v,x)-\nabla
h(x)\|^{2}\\!\leq\\!L_{x}^{2}(\|z-z^{*}(x)\|^{2}+\|v-v^{*}(x)\|^{2}).$

Thus, it is natural to make $z$ and $v$ move towards their respective
equilibrium values which are given by $z^{*}(x)$ and $v^{*}(x)$. As a
consequence, we also introduce the directions $D_{z}$ and $D_{x}$ as follows

$\displaystyle D_{z}(z,v,x)$ $\displaystyle=\nabla_{1}G(z,x)\enspace,$
$\displaystyle D_{v}(z,v,x)$
$\displaystyle=\nabla_{11}^{2}G(z,x)v+\nabla_{1}F(z,x)\enspace.$

The interest of considering the directions $D_{z}$ and $D_{v}$ is expressed in
Proposition 2.6.

######  Proposition 2.6 (First-order conditions).

Assume that $G$ is strongly convex with respect to its first variable. Then
for any $x\in\mathbb{R}^{d}$, it holds $D_{z}(z^{*}(x),v^{*}(x),x)=0$ and
$D_{v}(z^{*}(x),v^{*}(x),x)=0$.

The directions $D_{z}$, $D_{v}$, and $D_{x}$ can be written as sums over the
samples. Hence, as mentioned by Dagréou et al. (2022), following these
directions enables us to adapt any classical algorithm suited for single-level
finite sum minimization to bilevel finite sum minimization. In what follows,
for two indices $i\in[n]$ and $j\in[m]$, we consider the sampled directions
$D_{z,i,j}$, $D_{v,i,j}$ and $D_{x,i,j}$ defined by

$\displaystyle D_{z,i,j}(z,v,x)$ $\displaystyle=\nabla_{1}G_{i}(z,x)$ (4)
$\displaystyle D_{v,i,j}(z,v,x)$
$\displaystyle=\nabla_{11}^{2}G_{i}(z,x)v+\nabla_{1}F_{j}(z,x)$ (5)
$\displaystyle D_{x,i,j}(z,v,x)$
$\displaystyle=\nabla_{21}^{2}G_{i}(z,x)v+\nabla_{2}F_{j}(z,x)\enspace.$ (6)

When $i$ and $j$ are randomly sampled uniformly, these directions are unbiased
estimators of the true directions $D_{z}$, $D_{v}$, and $D_{x}$. Nevertheless,
as in Nguyen et al. (2017), we use them to recursively build biased estimators
of the directions that enable fast convergence.

###  2.3 SRBA: Stochastic Recursive Bilevel Algorithm

We propose SRBA which is a combination of the idea of recursive gradient
coming from (Fang et al., 2018; Nguyen et al., 2022) and the framework
proposed in (Dagréou et al., 2022). The SRBA algorithm relies on a recursive
estimation of each direction $D_{z}$, $D_{v}$, $D_{x}$ which is updated
following the same strategy as SARAH. Let us denote by $\rho$ the step size of
the update for the variables $z$ and $v$ and $\gamma$ the step size for the
update of the variable $x$. We use the same step size for $z$ and $v$ because
the problems of minimizing the inner function $G$ and solving the linear
system (3) have the same conditioning driven by $\nabla^{2}_{11}G$. For
simplicity, we denote the joint variable $\mathbf{u}=(z,v,x)$ and the joint
directions weighted by the step sizes $\mathbf{\Delta}~{}=~{}(\rho D_{z},\rho
D_{v},\gamma
D_{x})=(\mathbf{\Delta}_{z},\mathbf{\Delta}_{v},\mathbf{\Delta}_{x})$.

At iteration $t$, the estimate direction $\mathbf{\Delta}$ is initialized by
computing full batch directions:

$\displaystyle\mathbf{\Delta}^{t,0}=(\rho D_{z}(\mathbf{\tilde{u}}^{t}),\rho
D_{v}(\mathbf{\tilde{u}}^{t}),\gamma D_{x}(\mathbf{\tilde{u}}^{t}))$

and a first update is performed by moving from $\mathbf{\tilde{u}}^{t}$ in the
direction $-\mathbf{\Delta}^{t,0}$. As done in Hu et al. (2022), we project
the variable $v$ onto $\Gamma$ to leverage the boundedness property of
$v^{*}$. Then, during the $k$th iteration of an inner loop of size $q-1$, two
indices $i\in[n]$ and $j\in[m]$ are sampled and the estimate directions are
updated according to Equations (7) to (9)

$\displaystyle\mathbf{\Delta}_{z}^{t,k}$
$\displaystyle=\rho(D_{z,i,j}(\mathbf{u}^{t,k})-D_{z,i,j}(\mathbf{u}^{t,k-1}))+\mathbf{\Delta}_{z}^{t,k-1}$
(7) $\displaystyle\mathbf{\Delta}_{v}^{t,k}$
$\displaystyle=\rho(D_{v,i,j}(\mathbf{u}^{t,k})-D_{v,i,j}(\mathbf{u}^{t,k-1}))+\mathbf{\Delta}_{v}^{t,k-1}$
(8) $\displaystyle\mathbf{\Delta}_{x}^{t,k}$
$\displaystyle=\gamma(D_{x,i,j}(\mathbf{u}^{t,k})-D_{x,i,j}(\mathbf{u}^{t,k-1}))+\mathbf{\Delta}_{x}^{t,k-1\\!}$
(9)

where the sampled directions $D_{z,i,j}$, $D_{v,i,j}$ and $D_{x,i,j}$ are
defined by the Equations (4) to (6).

Input: initializations $z_{0}\in\mathbb{R}^{p}$, $x_{0}\in\mathbb{R}^{d}$,
$v_{0}\in\mathbb{R}^{p}$, number of iterations $T$ and $q$, step sizes $\rho$
and $\gamma$. Set $\mathbf{\tilde{u}^{0}}=(z_{0},v_{0},x_{0})$

for $t=0,\dots,T-1$ do

Reset $\mathbf{\Delta}$: $\mathbf{\Delta}^{t,0}=(\rho
D_{z}(\mathbf{\tilde{u}}^{t}),\rho D_{v}(\mathbf{\tilde{u}}^{t}),\gamma
D_{x}(\mathbf{\tilde{u}}^{t}))$ Update $\mathbf{u}$:
$\mathbf{u}^{t,1}=\Pi(\mathbf{\tilde{u}}^{t}-\mathbf{\Delta^{t,0}})\enspace,$

for $k=1,\dots,q-1$ do

Draw $i\in\\{1,\dots,n\\}$ and $j\in\\{1,\dots,m\\}$
$\mathbf{\Delta}_{z}^{t,k}=\rho(D_{z,i,j}(\mathbf{u}^{t,k})-D_{z,i,j}(\mathbf{u}^{t,k-1}))+\mathbf{\Delta}_{z}^{t,k-1}$
$\mathbf{\Delta}_{v}^{t,k}=\rho(D_{v,i,j}(\mathbf{u}^{t,k})-D_{v,i,j}(\mathbf{u}^{t,k-1}))+\mathbf{\Delta}_{v}^{t,k-1}$
$\mathbf{\Delta}_{x}^{t,k}=\gamma(D_{x,i,j}(\mathbf{u}^{t,k})-D_{x,i,j}(\mathbf{u}^{t,k-1}))+\mathbf{\Delta}_{x}^{t,k-1}$
Update $\mathbf{u}$:
$\mathbf{u}^{t,k+1}=\Pi(\mathbf{u}^{t,k}-\mathbf{\Delta}^{t,k})$

end forSet $\mathbf{\tilde{u}}^{t+1}=\mathbf{u}^{t+1,q}$

end forReturn
$(\tilde{z}^{T},\tilde{v}^{T},\tilde{x}^{T})=\mathbf{\tilde{u}}^{T}$

Algorithm 1 Stochastic Recursive Bilevel Algorithm

Then the joint variable $\mathbf{u}$ is updated by

$\mathbf{u}^{t,k+1}=\Pi(\mathbf{u}^{t,k}-\mathbf{\Delta}^{t,k})\enspace.$ (10)

Recall that the projection is only performed on the variable $v$. The other
variables $z$ and $x$ keep unchanged after the projection step. At the end of
the inner procedure, we set $\mathbf{\tilde{u}}^{t+1}=\mathbf{u}^{t,q}$. The
method is summarized in Algorithm 1. Note that this two loops structure with
periodic full batch computations is similar to the structure of SVRG. Unlike
SVRG, there is no reference point and the directions are updated recursively.

In Algorithm 1, the three variables $z$, $v$, and $x$ are updated
simultaneously rather than alternatively. From a computational perspective,
this allows us to share the common computations between the different oracles
and to do the update of each variable in parallel. As a consequence, there is
no sub-procedure to approximate the solution of the inner problem and the
solution of the linear system.

Note that in Yang et al. (2021), the authors propose VRBO, another adaptation
of SPIDER/SARAH for bilevel problems. VRBO has a double loop structure where
the inner variable is updated by several steps in an inner loop. In this inner
loop, the estimate of the gradient of $G$ and the gradient of $h$ are also
updated using SARAH’s update rules. SRBA has a different structure. First, in
SRBA, the inner variable $z$ is updated only once between two updates of the
outer variable instead of several times. Second, the solution of the linear
system evolves following optimization steps whereas in VRBO a Neumann
approximation is used. Finally, in Yang et al. (2021), the algorithm VRBO is
analyzed in the case where the functions $F$ and $G$ are general expectations
but not in the specific case of empirical risk minimization, as we do in
Section 3, and achieves a worse sample complexity (see Table 1).

##  3 Theoretical Analysis of SRBA

In this section we provide the theoretical analysis of Algorithm 1 leading to
a final sample complexity in
$\mathcal{O}\left((n+m)^{\frac{1}{2}}\varepsilon^{-1}\vee(n+m)\right)$. The
detailed proofs of the results are deferred to the appendix.

Before diving into the details, let us define a few concepts. In Definition
3.1, we recall the definition of $\varepsilon$-stationary point.

######  Definition 3.1 ($\varepsilon$-stationary point).

Let $d$ a positive integer, $f:\mathbb{R}^{d}\to\mathbb{R}$ a differentiable
function and $\varepsilon>0$. We say that a point $x\in\mathbb{R}^{d}$ is an
$\varepsilon$-stationary point of $f$ if $\|\nabla f(x)\|^{2}\leq\epsilon$.

With a slight abuse of language, in a stochastic context, we also call
$\varepsilon$-stationary point a random variable $x$ such that
$\mathbb{E}[\|\nabla f(x)\|^{2}]\leq\varepsilon$. This notion of
$\varepsilon$-stationary point is necessary since we are dealing with
nonconvex objectives.

In this paper, the theoretical complexity of the algorithms is given in terms
of number of calls to oracle, that is to say the number of times the quantity

$\hskip-6.45831pt[\nabla F_{j}(z,x),\nabla
G_{i}(z,x),\nabla^{2}_{11}G_{i}(z,x)v,\nabla^{2}_{21}G_{i}(z,x)v]$ (11)

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

###  3.1 Mean Squared Error of the Estimated Directions

One strength of our method is its simple expression of the estimation error of
the directions which comes from the bias-variance decomposition of the mean
squared error provided by Nguyen et al. (2017). Let us denote the estimate
directions $D^{t,k}_{z}~{}=~{}\mathbf{\Delta}_{z}^{t,k}/\rho$,
$D^{t,k}_{v}=\mathbf{\Delta}_{v}^{t,k}/\rho$ and
$D^{t,k}_{x}=\mathbf{\Delta}_{x}^{t,k}/\gamma$. We also introduce the
residuals

$\displaystyle S_{\bullet}^{t,k}$
$\displaystyle=\sum_{r=1}^{k}\mathbb{E}[\|D_{\bullet}(\mathbf{u}^{t,r})-D_{\bullet}(\mathbf{u}^{t,r-1})\|^{2}]$
$\displaystyle\tilde{S}_{\bullet}^{t,k}$
$\displaystyle=\sum_{r=1}^{k}\mathbb{E}[\|D^{t,r}_{\bullet}-D^{t,r-1}_{\bullet}\|^{2}]\enspace.$

Proposition 3.2 provides a simple link between the mean squared error
$\mathbb{E}[\|D_{\bullet}^{t,k}-D_{\bullet}(\mathbf{u}^{t,k})\|^{2}]$ and the
residuals $S_{\bullet}^{t,k}$ and $\tilde{S}_{\bullet}^{t,k}$.

######  Proposition 3.2 (MSE of the estimate directions).

For any $t\geq 0$ and $k\in\\{1,\dots,q-1\\}$, the estimate
$D_{\bullet}^{t,k}$ of the direction $D_{\bullet}(\mathbf{u}^{t,k})$ satisfies

$\displaystyle\mathbb{E}[\|D^{t,k}_{\bullet}-D_{\bullet}(\mathbf{u}^{t,k})\|^{2}]$
$\displaystyle=\tilde{S}_{\bullet}^{t,k}-S_{\bullet}^{t,k}\enspace.$

We observe that the above error has two components: the accumulation of the
difference between two successive full batch directions and the accumulation
of the difference between two successive estimate directions. Proposition 3.2
will play a critical role in the analysis of SRBA.

###  3.2 Fundamental Lemmas

As usually done in optimization, we start by establishing descent lemmas which
are key ingredients to get the final convergence result. Lemma 3.3 aims at
characterizing the joint dynamic of $\mathbf{u}$ on the inner problem. To do
so, we introduce the function $\phi_{z}$ defined as

$\phi_{z}(z,x)=G(z,x)-G(z^{*}(x),x)\enspace.$

In the bilevel literature, direct control on the distance to optimum
$\delta_{z}^{t,k}\triangleq\mathbb{E}[\|z^{t,k}-z^{*}(x^{t,k})\|^{2}]$ is
established. Here, the biased nature of the estimate direction $D^{t,k}_{z}$
makes it hard to upper bound appropriately the scalar product $\langle
D_{z}(\mathbf{u}^{t,k})-D_{z}^{t,k},z^{t,k}-z^{*}(x^{t,k})\rangle$. This
explains the choice of considering $\phi_{z}^{t,k}$ instead of
$\delta^{t,k}_{z}$. By combining the smoothness property of $\phi_{z}$ and the
bias-variance decomposition provided in Proposition 3.2, we can show some
descent property on the sequence $\phi_{z}^{t,k}$ defined by
$\phi_{z}^{t,k}=\mathbb{E}[\phi_{z}(z^{t,k},x^{t,k})]$. Before stating Lemma
3.3, let us define
$\mathcal{G}_{v}^{t,k}=\frac{1}{\rho}\left(v^{t,k}-\Pi_{\Gamma}(v^{t,k}-\rho
D^{t,k}_{v})\right)$ so that $v^{t,k+1}=v^{t,k}-\rho\mathcal{G}^{t,k}_{v}$.
This is the actual update direction of $v$. Note that if there were no
projections, we would have $\mathcal{G}^{t,k}_{v}=D^{t,k}_{v}$. As a
consequence, it acts as a surrogate of $D^{t,k}_{v}$ in our analysis. We also
define

$\displaystyle V^{t,k}_{z}=$
$\displaystyle\mathbb{E}[\|D^{t,k}_{z}\|^{2}],\quad
V^{t,k}_{v}=\mathbb{E}[\|\mathcal{G}^{t,k}_{v}\|^{2}],\quad
V^{t,k}_{x}=\mathbb{E}[\|D^{t,k}_{x}\|^{2}]$

the variances and their respective sums over the inner loop

$\displaystyle\mathcal{V}^{t,k}_{z}=\sum_{r=1}^{k}$
$\displaystyle\mathbb{E}[\|D^{t,r-1}_{z}\|^{2}],\quad\mathcal{V}^{t,k}_{v}=\sum_{r=1}^{k}\mathbb{E}[\|\mathcal{G}^{t,r-1}_{v}\|^{2}],\quad\mathcal{V}^{t,k}_{x}=\sum_{r=1}^{k}\mathbb{E}[\|D^{t,r-1}_{x}\|^{2}]\enspace.$

######  Lemma 3.3 (Descent on the inner level).

Assume that the step sizes $\rho$ and $\gamma$ verify
$\gamma~{}\leq~{}C_{z}\rho$ for some positive constant $C_{z}$ specified in
the appendix. Then it holds

$\displaystyle\phi_{z}^{t,k+1}$
$\displaystyle\leq\left(1-\frac{\mu_{G}}{2}\rho\right)\phi_{z}^{t,k}-\frac{\rho}{2}\left(1-\Lambda_{z}\rho\right)V_{z}^{t,k}+\rho^{3}\beta_{zz}\mathcal{V}_{z}^{t,k}+\gamma^{2}\rho\beta_{zv}\mathcal{V}_{v}^{t,k}$
(12)
$\displaystyle\qquad+\gamma^{2}\rho\beta_{zx}\mathcal{V}_{x}^{t,k}+\frac{\Lambda_{z}}{2}\gamma^{2}V_{x}^{t,k}+\frac{\gamma^{2}}{\rho}\overline{\beta}_{zx}\mathbb{E}[\|D^{t,k}_{x}(\mathbf{u}^{t,k})\|^{2}]$

for some positive constants $\Lambda_{z},\beta_{zz}$, $\beta_{zx}$ and
$\overline{\beta}_{zx}$ that are specified in the appendix.

In (12) we recover a linear decrease of $\phi_{z}^{t,k}$ by a factor
$(1-\rho\mu_{G})$ but the outer variable’s movement and the stochasticity make
appear the direction $D_{x}(\mathbf{u}^{t,k})$ and the noise that slow down
the convergence of $z$ towards $z^{*}(x)$.

For the variable $v$, the quantity we consider is

$\phi_{v}(v,x)~{}=~{}\Psi(z^{*}(x),v,x)-\Psi(z^{*}(x),v^{*}(x),x)$

where $\Psi(z,v,x)$ is defined as

$\Psi(z,v,x)=\frac{1}{2}v^{\top}\nabla^{2}_{11}G(z,x)v+\nabla_{1}F(z,x)^{\top}v\enspace.$

The intuition behind considering this quantity is that solving the linear
system (3) is equivalent to minimizing over $v$ the function
$\Psi(z^{*}(x),v,x)$.

######  Lemma 3.4 (Descent on the linear system).

Assume that the step sizes $\rho$ and $\gamma$ verify $\rho\leq B_{v}$ and
$\gamma~{}\leq~{}C_{v}\rho$ for some positive constants $B_{v}$ and $C_{v}$
specified in the appendix. Then it holds

$\displaystyle\phi_{v}^{t,k+1}$
$\displaystyle\leq\left(1-\frac{\rho\mu_{G}}{16}\right)\phi_{v}^{t,k}-\tilde{\beta}_{vv}\rho
V_{v}^{t}+\rho^{3}\beta_{vz}\mathcal{V}_{z}^{t,k}+2\rho^{3}\beta_{vv}\mathcal{V}_{v}^{t,k}+\gamma^{2}\rho\beta_{vx}\mathcal{V}_{x}^{t,k}$
$\displaystyle\qquad+\rho\alpha_{vz}\phi_{z}^{t,k}+\frac{\Lambda_{v}}{2}\gamma^{2}\mathbb{E}[\|D_{x}^{t,k}\|^{2}]+\frac{\gamma^{2}}{\rho}\overline{\beta}_{vx}\mathbb{E}[\|D_{x}(\mathbf{u}^{t,k})\|^{2}]$

for some positive constants $\Lambda_{v},\beta_{vz}$, $\beta_{vx}$,
$\tilde{\beta}_{vv}$ and $\overline{\beta}_{vx}$ that are specified in the
appendix.

Lemma 3.4 is similar to Lemma 3.3. The appearance of $\phi_{z}^{t,k}$ is a
consequence of the fact $D_{v}(z,v,x)$ is a step towards
$-[\nabla^{2}_{11}G(z,x)]^{-1}\nabla_{1}F(z,x)$ instead of
$-[\nabla^{2}_{11}G(z^{*}(x),x)]^{-1}\nabla_{1}F(z^{*}(x),x)$. The proof of
this lemma harnesses the generalization of Polyak-Łojasiewicz inequality for
composite functions introduced in Karimi et al. (2016).

The following lemma is a consequence of the smoothness of $h$. Let us denote
the expected values $h^{t,k}=\mathbb{E}[h(x^{t,k})]$ and expected gradient
$g^{t,k}~{}=~{}\mathbb{E}[\|\nabla h(x^{t,k})\|^{2}]$.

######  Lemma 3.5 (Descent on the value function $h$).

There exist constants $\beta_{hz}$, $\beta_{hv}$, $\beta_{hx}>0$ such that

$\displaystyle h^{t,k+1}$ $\displaystyle\leq
h^{t,k}-\frac{\gamma}{2}g^{t,k}+\gamma\frac{2L_{x}^{2}}{\mu_{G}}(\phi_{z}^{t,k}+\phi_{v}^{t,k})+\gamma\rho^{2}\beta_{hz}\mathcal{V}^{t,k}_{z}$
$\displaystyle\qquad+\gamma\rho^{2}\beta_{hv}\mathcal{V}^{t,k}_{v}+\gamma^{3}\beta_{hx}\mathcal{V}^{t,k}_{x}-\frac{\gamma}{2}\left(1-L^{h}\gamma\right)V^{t,k}_{x}\enspace.$

This lemma shows that the control of the approximation error $\phi_{\bullet}$
(Lemma 3.3 and Lemma 3.4) and the sum of variances $\mathcal{V}_{\bullet}$ is
crucial to get a decrease of $\mathbb{E}[h(x^{t,k})]$.

###  3.3 Complexity Analysis of SRBA

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

######  Theorem 1 (Convergence rate of SRBA).

Assume that Assumptions 2.1 and 2.2 hold. Assume that the step sizes verify
$\rho\leq\overline{\rho}$ and $\gamma\leq\min(\overline{\gamma},\xi\rho)$ for
some constants $\xi$, $\overline{\rho}$ and $\overline{\gamma}$ specified in
appendix. Then it holds

$\frac{1}{Tq}\sum_{t=0}^{T-1}\sum_{k=0}^{q-1}\mathbb{E}[\|\nabla
h(x^{t,k})\|^{2}]=\mathcal{O}\left(\frac{1}{qT\gamma}\right)$

where $\mathcal{O}$ hides regularity constants that are independent from $n$
and $m$.

The proof combines classical proof techniques from the bilevel literature and
elements from SARAH’s analysis (Nguyen et al., 2017, 2022). We introduce the
Lyapunov function
$\mathcal{L}(\mathbf{u}^{t,k})=h^{t,k}+\psi_{z}\phi_{z}^{t,k}+\psi_{v}\phi_{v}^{t,k}$
where $\psi_{z}$ and $\psi_{v}$ are non-negative constants chosen so that we
have the inequality
$\mathcal{L}(\mathbf{u}^{t,k+1})\leq\mathcal{L}(\mathbf{u}^{t,k})-\frac{\gamma}{4}g^{t,k}$.
Summing and telescoping this inequality provides the convergence rate.

Note that if we set $q=1$, we are actually in a nonstochastic regime and we
can observe that we recover the convergence rate of Gradient Descent for
nonconvex single-level problems (Nesterov, 2018) since the step size $\gamma$
depends neither on the current iteration $t$ nor the horizon $T$. Increasing
$q$ allows a faster convergence in terms of iterations but makes each
iteration more expensive since the number of oracle calls per iteration is
$(2n+3m)+2\times 5(q-1)$. Thus, there is a trade-off between the convergence
rate and the overall complexity. In Corollary 3.6, we state that the value of
$q$ that gives the best sample complexity is $\mathcal{O}(n+m)$.

######  Corollary 3.6 (Sample complexity of SRBA).

Suppose that Assumptions 2.1 and 2.2 hold. If we take
$\rho~{}=~{}\overline{\rho}(n+m)^{-\frac{1}{2}}$,
$\gamma~{}=~{}\min(\overline{\gamma},\xi\rho)(n+m)^{-\frac{1}{2}}$ and
$q~{}=~{}n+m$, then
$\mathcal{O}\left((n+m)^{\frac{1}{2}}\varepsilon^{-1}\vee(n+m)\right)$ calls
to oracles are sufficient to find an $\varepsilon$-stationary point with SRBA.

This sample complexity is analogous to the sample complexity of SARAH in the
nonconvex finite-sum setting. To the best of our knowledge, such a rate is the
best known for bilevel empirical risk minimization problems in terms of
dependency on the number of samples $n+m$ and the precision $\varepsilon$.
This improve by a factor $(n+m)^{-\frac{1}{6}}$ the previous result which was
achieved by SABA (Dagréou et al., 2022). As a comparison, VRBO (Yang et al.,
2021) achieves a sample complexity in
$\tilde{\mathcal{O}}(\varepsilon^{-\frac{3}{2}})$. Note that, for large value
of $n+m$ we can have actually
$(n+m)^{\frac{1}{2}}\varepsilon^{-1}\gtrsim\varepsilon^{-2}$. This means that,
just like single-level SARAH, the complexity of SRBA can be beaten by others
when the number of samples is too high with respect to the desired accuracy
(actually, if $n+m=\Omega(\varepsilon^{-2})$).

##  4 Lower Bound for Bilevel ERM

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

###  4.1 Functions and Algorithms Classes

We start by defining the function class and the algorithm class we consider.

######  Definition 4.1 (Function class).

Let $n,m$ two positive integers, $L^{F}_{1}$ and $\mu_{G}$ two positive
constants. The class of the smooth empirical risk minimization problems
denoted by $\mathcal{C}^{L^{F}_{1},\mu_{G}}$ is the set of pairs of real-
valued function families $((F_{j})_{1\leq j\leq m},(G_{i})_{1\leq i\leq n})$
defined on $\mathbb{R}^{p}\times\mathbb{R}^{d}$ such that for all $j\in[m]$,
$F_{j}$ is $L^{F}_{1}$ smooth and for all $i\in[n]$, $G_{i}$ is twice
differentiable and $\mu_{G}$-strongly convex.

Note that we do not require $F$ to be convex in the class
$\mathcal{C}^{L^{F}_{1},\mu_{G}}$. In particular, the class of bilevel
problems that we consider is nonconvex. This class contains, for instance, the
functions defining the bilevel formulation of the datacleaning task (see
Section 5).

For the algorithmic class, we consider algorithms that implement approximate
implicit differentiation, using oracles of the form (11).

######  Definition 4.2 (Algorithmic class).

Given initial points $z^{0},v^{0},x^{0}$, a _linear bilevel algorithm
$\mathcal{A}$_ is a measurable mapping such that for any $((F_{j})_{1\leq
j\leq m},(G_{i})_{1\leq i\leq n})~{}\in~{}\mathcal{C}^{L^{F}_{1},\mu_{G}}$,
the output of $\mathcal{A}((F_{j})_{1\leq j\leq m},(G_{i})_{1\leq i\leq n})$
is a sequence $\\{(z^{t},v^{t},x^{t},i_{t},j_{t})\\}_{t\geq 0}$ of points
$(z^{t},v^{t},x^{t})$ and random variables $i_{t}\in[n]$ and $j_{t}\in[m]$
such that for all $t\geq 0$

$\displaystyle z^{t+1}\in z^{0}+$
$\displaystyle\mathrm{Span}\\{\nabla_{1}G_{i_{0}}(z^{0},x^{0}),\dots,\nabla_{1}G_{i_{t}}(z^{t},x^{t})\\}$
$\displaystyle v^{t+1}\in v^{0}+$
$\displaystyle\mathrm{Span}\\{\nabla^{2}_{11}G_{i_{0}}(z^{0},x^{0})v^{0}+\nabla_{1}F_{j_{0}}(z^{0},x^{0}),$
$\displaystyle\dots,\nabla^{2}_{11}G_{i_{t}}(z^{t},x^{t})v^{t}+\nabla_{1}F_{j_{t}}(z^{t},x^{t})\\}$
$\displaystyle x^{t+1}\in x^{0}+$
$\displaystyle\mathrm{Span}\\{\nabla^{2}_{21}G_{i_{0}}(z^{0},x^{0})v^{0}+\nabla_{2}F_{j_{0}}(z^{0},x^{0}),$
$\displaystyle\dots,\nabla^{2}_{21}G_{i_{t}}(z^{t},x^{t})v^{t}+\nabla_{2}F_{j_{t}}(z^{t},x^{t})\\}.$

This algorithmic class includes popular stochastic bilevel first-order
algorithms, such as AmIGO (Arbel and Mairal, 2022), FSLA (Li et al., 2022),
SOBA, and SABA (Dagréou et al., 2022). Note that despite the projection step,
SRBA is part of this algorithmic class since the projection of a vector onto
$\Gamma$ is actually just a rescaling.

###  4.2 Main Theorem

Problem (1) is actually a smooth nonconvex optimization problem. The lower
complexity bound for nonconvex finite sum problem has been studied in Fang et
al. (2018); Zhou and Gu (2019). In particular, they show that the number of
gradient calls needed to get an $\varepsilon$-stationary point for a smooth
nonconvex finite sum is at least
$\mathcal{O}(m^{\frac{1}{2}}\varepsilon^{-1})$, where $m$ is the number of
terms in the finite sum.

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

######  Theorem 2 (Lower bound for bilevel ERM).

For any linear bilevel algorithm $\mathcal{A}$, and any $L^{F}$, $n$,
$\Delta$, $\varepsilon$, $p$ such that $\varepsilon^{2}\leq(\Delta
L^{F}m^{-1})/10^{3}$, there exists a dimension
$d=\mathcal{O}(\Delta\varepsilon^{-1}m^{\frac{1}{2}}L^{F})$, an element
$((F_{j})_{1\leq j\leq m},(G_{i})_{1\leq i\leq
n})\in\mathcal{C}^{L^{F}_{1},\mu_{G}}$ such that the value function $h$
defined as in (1) satisfies
$h(x^{0})~{}-~{}\inf_{x\in\mathbb{R}^{d}}h(x)\leq\Delta$ and in order to find
$\hat{x}\in\mathbb{R}^{d}$ such that $\mathbb{E}[\|\nabla
h(\hat{x})\|^{2}]\leq\varepsilon$, $\mathcal{A}$ needs at least
$\Omega(m^{\frac{1}{2}}\varepsilon^{-1})$ calls to oracles of the form (11).

The proof is an adaptation of the proof of Zhou and Gu (2019, Theorem 4.7). It
consists in taking as outer function $F$ defined by
$F(z,x)=\sum_{j=1}^{m}f(U^{(j)}z)$ where $f$ is the “worst-case function” used
by Carmon et al. (2017a), $U=[U^{(j)},\dots,U^{(m)}]^{\top}$ is an orthogonal
matrix and $G(z,x)=\frac{1}{2}\|z-x\|^{2}$. We harness the fact that $\|\nabla
f(y)\|^{2}>K$ as long as the two last coordinates of $y$ are zero for some
known constant $K$. Then we use the “zero chain property” to upper bound the
number of indices $j$ such the two last components of $U^{(j)}x^{t}$ are zero
at a given iteration $t$, implying $\|\nabla h(x^{t})\|^{2}>\epsilon$ when $t$
is smaller than $\mathcal{O}(m^{\frac{1}{2}}\varepsilon^{-1})$.

Note that here, the function class considered is less restrictive than the
function class that verifies the upper complexity bound achieved by SRBA in
Corollary 3.6. Considering a class as restrictive as the class needed for the
analysis of SRBA could lead to a smaller lower bound.

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

Figure 1: Left: Comparison of the performances of SRBA for different values of
$q$ on a hyperparameter selection task for $\ell^{2}$-regularized logistic
regression with IJCNN1 dataset with respect to time and iterations. Right:
Comparison of SRBA with other stochastic bilevel methods on the datacleaning
task with the MNIST dataset. The solvers are run with 10 different seeds and
the median performance over these seeds is reported. The shaded area
corresponds to the performances between the 20% and the 80% percentiles. We
report the test error with respect to wall clock time. We notice that SRBA
achieves the best final accuracy even though it is slower than the others at
the beginning.

##  5 Numerical Experiments

Even though our contribution is mostly theoretical, we run several experiments
to highlight the influence of the inner loop size on the performances of SRBA
and to compare the proposed algorithm with state-of-the-art stochastic bilevel
solvers. A more detailed description of the experiments is available in
Appendix C.

###  5.1 Influence of the Period $q$

We are interested in the impact of the period $q$ on the algorithm’s
performance. We consider the hyperparameter selection problem for
$\ell^{2}$-regularized logistic regression with the dataset
IJCNN1555https://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/binary.html.
In this case, $F$ is the validation loss and $G$ the training loss. We run
SRBA for several values of $q$. In Figure 1, we display the suboptimality
$h(x^{t})-h^{*}$ where $h^{*}$ is the minimum value reached among all the
runs. The performances are reported both relatively to wall clock time and
iterations. An iteration corresponds to one update of the variables $z$, $v$,
and $x$ with the full batch or stochastic directions.

The first observation is that the parameter $q$ impacts dramatically the
convergence speed of SRBA in practice. While all the runs converge, the
variance and the speed of the suboptimality $h(x^{t})-h^{*}$ differ. The
figure shows that reducing the period $q$ gives performances with less
variance, due to improved gradient estimators. However, we notice a difference
between the performances with respect to time and with respect to iteration.
For instance, the curve corresponding to $q=n+m$ is among the best curves when
looking in terms of iterations while it becomes the second slowest with
respect to time. This suggests that there is a trade-off between computing too
many times the full batch quantities versus improving the gradient estimates.
In the presented experiment, $q=4(n+m)$ gives the best performances.

###  5.2 Comparison of SRBA with Competitors

We compare the performances of SRBA with stochastic bilevel solvers on the
datacleaning problem (Franceschi et al., 2017) with MNIST
dataset666http://yann.lecun.com/exdb/mnist/. In this task, the training set is
composed of $n_{\mathrm{train}}$ labeled samples
$(d_{i}^{\mathrm{train}},y_{i}^{\mathrm{train}})_{i\in[n_{\mathrm{train}}]}$
that possibly have corrupted labels with a probability $p_{c}$. The validation
set has $n_{{\mathrm{val}}}$ clean labels. The datacleaning task consists in
learning simultaneously a classifier and a weighting on the training samples.
This problem can be cast as a bilevel optimization problem where the inner
loss is the training loss where the training samples are weighted and the
outer loss is the validation loss. The inner variable is the parameter of the
classifier and belongs to $\mathbb{R}^{p}$ and the outer variable is the
weighting of the training set and belongs to
$\mathbb{R}^{n_{\mathrm{train}}}$. A more formal formulation of the problem is
provided in the appendix.

In Figure 1 (right), we compare the final test error of SRBA with AmIGO (Arbel
and Mairal, 2022), MRBO (Yang et al., 2021), StocBiO (Ji et al., 2021) and
SABA (Dagréou et al., 2022). In this experiment, we have $p_{c}=0.5$. The
parameters of each algorithm have been selected by a grid search. We observe
that SRBA reaches the lowest plateau. Nevertheless, it is the slowest at the
beginning due to the full batch computations.

##  6 Conclusion

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

## Acknowledgements

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

## References

  * Ablin et al. (2020) P. Ablin, G. Peyré, and T. Moreau.  Super-efficiency of automatic differentiation for functions defined as a minimum.  In _International Conference on Machine Learning (ICML)_ , 2020.
  * Agarwal and Bottou (2015) A. Agarwal and L. Bottou.  A Lower Bound for the Optimization of Finite Sums.  In _International Conference on Machine Learning (ICML)_ , 2015.
  * Arbel and Mairal (2022) M. Arbel and J. Mairal.  Amortized Implicit Differentiation for Stochastic Bilevel Optimization.  In _International Conference on Learning Representations (ICLR)_ , 2022.
  * Bai et al. (2019) S. Bai, J. Z. Kolter, and V. Koltun.  Deep Equilibrium Models.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , 2019.
  * Bengio (2000) Y. Bengio.  Gradient-Based Optimization of Hyperparameters.  _Neural Computation_ , 12(8):1889–1900, 2000\.  ISSN 0899-7667, 1530-888X.  doi: 10.1162/089976600300015187.
  * Bottou (2010) L. Bottou.  Large-Scale Machine Learning with Stochastic Gradient Descent.  In Y. Lechevallier and G. Saporta, editors, _Proceedings of COMPSTAT’2010_ , pages 177–186. Physica-Verlag HD, Heidelberg, 2010.  ISBN 978-3-7908-2603-6 978-3-7908-2604-3.  doi: 10.1007/978-3-7908-2604-3˙16.
  * Bottou et al. (2018) L. Bottou, F. E. Curtis, and J. Nocedal.  Optimization methods for large-scale machine learning.  _Siam Reviews_ , 60(2):223–311, 2018.
  * Carmon et al. (2017a) Y. Carmon, J. C. Duchi, O. Hinder, and A. Sidford.  Lower Bounds for Finding Stationary Points II: First-Order Methods, 2017a.
  * Carmon et al. (2017b) Y. Carmon, J. C. Duchi, O. Hinder, and A. Sidford.  Lower Bounds for Finding Stationary Points I.  _arXiv preprint arXiv:1710.11606_ , 2017b.
  * Chen et al. (2021) T. Chen, Y. Sun, and W. Yin.  Closing the Gap: Tighter Analysis of Alternating Stochastic Gradient Methods for Bilevel Problems.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , 2021.
  * Cutkosky and Orabona (2019) A. Cutkosky and F. Orabona.  Momentum-based variance reduction in non-convex SGD.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , 2019.
  * Dagréou et al. (2022) M. Dagréou, P. Ablin, S. Vaiter, and T. Moreau.  A framework for bilevel optimization that enables stochastic and global variance reduction algorithms.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , 2022.
  * Defazio et al. (2014) A. Defazio, F. Bach, and S. Lacoste-Julien.  SAGA: A Fast Incremental Gradient Method With Support for Non-Strongly Convex Composite Objectives.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , 2014.
  * Fang et al. (2018) C. Fang, C. J. Li, Z. Lin, and T. Zhang.  SPIDER: Near-Optimal Non-Convex Optimization via Stochastic Path Integrated Differential Estimator.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , 2018.
  * Franceschi et al. (2017) L. Franceschi, M. Donini, P. Frasconi, and M. Pontil.  Forward and Reverse Gradient-Based Hyperparameter Optimization.  In _International Conference on Machine Learning (ICML)_ , 2017.
  * Franceschi et al. (2018) L. Franceschi, P. Frasconi, S. Salzo, R. Grazzi, and M. Pontil.  Bilevel Programming for Hyperparameter Optimization and Meta-Learning.  In _International Conference on Machine Learning (ICML)_ , 2018.
  * Ghadimi and Lan (2013) S. Ghadimi and G. Lan.  Stochastic first- and zeroth-order methods for nonconvex stochastic programming.  _SIAM Journal on Optimization_ , 23(4):2341–2368, 2013.  doi: 10.1137/120880811.
  * Ghadimi and Wang (2018) S. Ghadimi and M. Wang.  Approximation Methods for Bilevel Programming.  _arXiv preprint arXiv:1802.02246_ , 2018.
  * Grazzi et al. (2020) R. Grazzi, L. Franceschi, M. Pontil, and S. Salzo.  On the iteration complexity of hypergradient computation.  In _Internation Conference on Machine Learning (ICML)_ , 2020.
  * Hong et al. (2021) M. Hong, H.-T. Wai, Z. Wang, and Z. Yang.  A Two-Timescale Framework for Bilevel Optimization: Complexity Analysis and Application to Actor-Critic.  _arXiv preprint arXiv:2007.05170_ , 2021.
  * Hu et al. (2022) Q. Hu, Y. Zhong, and T. Yang.  Multi-block Min-max Bilevel Optimization with Applications in Multi-task Deep AUC Maximization.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , 2022.
  * Huang et al. (2022) F. Huang, J. Li, and H. Huang.  Enhanced Bilevel Optimization via Bregman Distance.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , 2022.
  * Ji and Liang (2023) K. Ji and Y. Liang.  Lower Bounds and Accelerated Algorithms for Bilevel Optimization.  _Journal of Machine Learning Research_ , 24(22):1–56, 2023.
  * Ji et al. (2021) K. Ji, J. Yang, and Y. Liang.  Bilevel Optimization: Convergence Analysis and Enhanced Design.  In _International Conference on Machine Learning (ICML)_ , 2021.
  * Karimi et al. (2016) H. Karimi, J. Nutini, and M. Schmidt.  Linear Convergence of Gradient and Proximal-Gradient Methods Under the Polyak-\Lojasiewicz Condition.  In _European Conference on Machine Learning (ECML)_ , 2016\.
  * Khanduri et al. (2021) P. Khanduri, S. Zeng, M. Hong, H.-T. Wai, Z. Wang, and Z. Yang.  A Near-Optimal Algorithm for Stochastic Bilevel Optimization via Double-Momentum.  _arXiv preprint arXiv:2102.07367_ , 2021.
  * Li et al. (2022) J. Li, B. Gu, and H. Huang.  A Fully Single Loop Algorithm for Bilevel Optimization without Hessian Inverse.  In _Proceedings of the Thirty-sixth AAAI Conference on Artificial Intelligence, AAAI’22_ , 2022.
  * Li et al. (2020) Y. Li, G. Hu, Y. Wang, T. Hospedales, N. M. Robertson, and Y. Yang.  DADA: Differentiable Automatic Data Augmentation.  _arXiv preprint arXiv:2003.03780_ , 2020.
  * Liu et al. (2019) H. Liu, K. Simonyan, and Y. Yang.  DARTS: Differentiable Architecture Search.  _arXiv preprint arXiv:1806.09055_ , 2019.
  * Lorraine et al. (2020) J. Lorraine, P. Vicol, and D. Duvenaud.  Optimizing Millions of Hyperparameters by Implicit Differentiation.  In _International Conference on Artificial Intelligence and Statistics (AISTATS)_ , 2020.
  * Maclaurin et al. (2015) D. Maclaurin, D. Duvenaud, and R. P. Adams.  Gradient-based Hyperparameter Optimization through Reversible Learning.  In _International Conference on Mahcine Learning (ICML)_ , 2015.
  * Moreau et al. (2022) T. Moreau, M. Massias, A. Gramfort, P. Ablin, P.-A. B. B. Charlier, M. Dagréou, T. D. la Tour, G. Durif, C. F. Dantas, Q. Klopfenstein, J. Larsson, E. Lai, T. Lefort, B. Malézieux, B. Moufad, B. T. Nguyen, A. Rakotomamonjy, Z. Ramzi, J. Salmon, and S. Vaiter.  Benchopt: Reproducible, efficient and collaborative optimization benchmarks.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , 2022.
  * Nemirovsky and Yudin (1983) A. S. Nemirovsky and D. B. Yudin.  _Problem complexity and method efficiency in optimization_.  Wiley-Interscience series in discrete mathematics. Wiley, 1983.  ISBN 978-0-471-10345-5.
  * Nesterov (2018) Y. Nesterov.  _Lectures on Convex Optimization_.  Springer Berlin Heidelberg, New York, NY, 2018.  ISBN 978-3-319-91577-7.
  * Nguyen et al. (2017) L. M. Nguyen, J. Liu, K. Scheinberg, and M. Takáč.  SARAH: A Novel Method for Machine Learning Problems Using Stochastic Recursive Gradient.  In _International Conference on Machine Learning (ICML)_ , 2017.
  * Nguyen et al. (2022) L. M. Nguyen, M. van Dijk, D. T. Phan, P. H. Nguyen, T.-W. Weng, and J. R. Kalagnanam.  Finite-sum smooth optimization with SARAH.  _Computational Optimization and Applications_ , 82(3):561–593, 2022.  ISSN 0926-6003, 1573-2894.  doi: 10.1007/s10589-022-00375-x.
  * Pearlmutter (1994) B. A. Pearlmutter.  Fast exact multiplication by the hessian.  _Neural computation_ , 6(1):147–160, 1994.
  * Pedregosa (2016) F. Pedregosa.  Hyperparameter optimization with approximate gradient.  In _International Conference on Machine Learning (ICML)_ , 2016.
  * Rajeswaran et al. (2019) A. Rajeswaran, C. Finn, S. Kakade, and S. Levine.  Meta-Learning with Implicit Gradients.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , 2019.
  * Ramzi et al. (2022) Z. Ramzi, F. Mannel, S. Bai, J.-L. Starck, P. Ciuciu, and T. Moreau.  SHINE: SHaring the INverse Estimate from the forward pass for bi-level optimization and implicit models.  In _International Conference on Learning Representations (ICLR)_ , 2022.
  * Reddi et al. (2016) S. J. Reddi, S. Sra, B. Poczos, and A. Smola.  Fast Incremental Method for Nonconvex Optimization.  In _2016 IEEE 55th Conference on Decision and Control (CDC)_ , IEEE, pages 1971–1977, 2016.
  * Robbins and Monro (1951) H. Robbins and S. Monro.  A stochastic approximation method.  _The Annals of Mathematical Statistics_ , 22(3):400–407, 1951.  doi: 10.1214/aoms/1177729586.
  * Rommel et al. (2022) C. Rommel, T. Moreau, J. Paillard, and A. Gramfort.  CADDA: Class-wise Automatic Differentiable Data Augmentation for EEG Signals.  In _International Conference on Learning Representations (ICLR)_ , 2022.
  * Sow et al. (2022) D. Sow, K. Ji, and Y. Liang.  On the Convergence Theory for Hessian-Free Bilevel Algorithms.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , 2022.
  * Woodworth and Srebro (2016) B. E. Woodworth and N. Srebro.  Tight Complexity Bounds for Optimizing Composite Objectives.  In _Advances in Neural Information Systems Processing (NeurIPS)_ , 2016.
  * Yang et al. (2021) J. Yang, K. Ji, and Y. Liang.  Provably Faster Algorithms for Bilevel Optimization.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , 2021.
  * Zhang et al. (2021) M. Zhang, S. W. Su, S. Pan, X. Chang, E. M. Abbasnejad, and R. Haffari.  iDARTS: Differentiable Architecture Search with Stochastic Implicit Gradients.  In _International Conference on Machine Learning (ICML)_ , 2021.
  * Zhou and Gu (2019) D. Zhou and Q. Gu.  Lower Bounds for Smooth Nonconvex Finite-Sum Optimization.  In _International Conference on Machine Learning (ICML)_ , 2019.

Appendix A contains the necessary lemmas and proofs of Section 3. Appendix B
contains the proof of the lower bound for stochastic bilevel optimization.
Finally, Appendix C details the setting of the numerical experiments.

##  Appendix A Convergence analysis of SRBA

###  A.1 Proof of Proposition 2.6

###### Proof.

Let $x\in\mathbb{R}^{d}$. Since $G(\,.\,,x)$ is differentiable and $z^{*}(x)$
minimizes $G(\,.\,,x)$, the first order optimality condition ensures
$\nabla_{1}G(z^{*}(x),x)=0=D_{z}(z^{*}(x),v^{*}(x),x)$. Since $G$ is strongly
convex with respect to $z$, the Hessian $\nabla^{2}_{11}G(z^{*}(x),x)$ is
invertible. As a consequence, the equation in $v$

$\displaystyle
D_{v}(z^{*}(x),v,x)=\nabla^{2}_{11}G(z^{*}(x),x)v+\nabla_{1}F(z^{*}(x),x)=0$
(13)

admits a unique solution given by $v^{*}(x)$. ∎

###  A.2 Smoothness constant of $h$

We can find in Ghadimi and Wang [2018, Lemma 2.2] the following value for the
smoothness constant of $h$

$L^{h}=L^{F}_{1}+\frac{2L^{F}_{1}L^{G}_{2}+(L^{F}_{0})^{2}L^{G}_{2}}{\mu_{G}}+\frac{L_{11}^{G}L^{G}_{1}L_{0}^{F}+L^{G}_{1}L^{G}_{2}L^{F}_{0}+(L^{G}_{1})^{2}L^{F}_{1}}{\mu_{G}^{2}}+\frac{(L^{G}_{1})^{2}L^{G}_{2}L^{F}_{0}}{\mu_{G}^{3}}\enspace.$

###  A.3 Proof of Proposition 3.2

###### Proof.

Let $t>0$ and $k\in[q-1]$. For $k=0$, we directly have
$\mathbb{E}[\|D_{\bullet}^{t,k}-D_{\bullet}(\mathbf{u}^{t,k})\|^{2}]=0$. For
$k\geq 1$ and $r\in\\{1,\dots,k\\}$, the bias/variance decomposition of
$D_{\bullet}^{t,r}$ reads

$\displaystyle\mathbb{E}_{t,r}[\|D_{\bullet}^{t,r}-D_{\bullet}(\mathbf{u}^{t,r})\|^{2}]$
$\displaystyle=\mathbb{E}_{t,r}[\|D_{\bullet}^{t,r}-D_{\bullet}^{t,r-1}+D_{\bullet}(\mathbf{u}^{t,r-1})-D_{\bullet}(\mathbf{u}^{t,r})\|^{2}]$
$\displaystyle\qquad+\|D_{\bullet}(\mathbf{u}^{t,r})+D_{\bullet}(\mathbf{u}^{t,r-1})-D_{\bullet}^{t,r-1}-D_{\bullet}(\mathbf{u}^{t,r})\|^{2}$
$\displaystyle=\mathbb{E}_{t,r}[\|D_{\bullet}^{t,r}-D_{\bullet}^{t,r-1}-(D_{\bullet}(\mathbf{u}^{t,r-1})-D_{\bullet}(\mathbf{u}^{t,r}))\|^{2}]$
$\displaystyle\qquad+\|D_{\bullet}^{t,r-1}-D_{\bullet}(\mathbf{u}^{t,r-1})\|^{2}$

The term
$\mathbb{E}_{t,r}[\|D_{\bullet}^{t,r}-D_{\bullet}^{t,r-1}-(D_{\bullet}(\mathbf{u}^{t,r-1})-D_{\bullet}(\mathbf{u}^{t,r}))\|^{2}]$
is the variance of $D_{\bullet}^{t,r}-D_{\bullet}^{t,r-1}$, and then can
written as

$\displaystyle\mathbb{E}_{t,r}[\|D_{\bullet}^{t,r}-D_{\bullet}^{t,r-1}-(D_{\bullet}(\mathbf{u}^{t,r-1})-D_{\bullet}(\mathbf{u}^{t,r}))\|^{2}]$
$\displaystyle=\mathbb{E}_{t,r}[\|D_{\bullet}^{t,r}-D_{\bullet}^{t,r-1}\|^{2}]$
$\displaystyle\qquad-\|D_{\bullet}(\mathbf{u}^{t,r})-D_{\bullet}(\mathbf{u}^{t,r-1})\|^{2}$

Plugging this in the previous inequality and taking the total expectation
leads to

$\displaystyle\mathbb{E}[\|D_{\bullet}^{t,r}-D_{\bullet}(\mathbf{u}^{t,r})\|^{2}]$
$\displaystyle=\mathbb{E}[\|D_{\bullet}^{t,r}-D_{\bullet}^{t,r-1}\|^{2}]-\mathbb{E}[\|D_{\bullet}(\mathbf{u}^{t,r})-D_{\bullet}(\mathbf{u}^{t,r-1})\|^{2}]$
$\displaystyle\qquad+\mathbb{E}[\|D^{t,r-1}_{\bullet}-D^{t,r-1}_{\bullet}(\mathbf{u}^{t,r-1})\|^{2}]$

Summing for $r\in\\{1,\dots,k\\}$ and telescoping gives the final result
(taking into account that
$D^{t,0}_{\bullet}~{}=~{}D_{\bullet}(\mathbf{u}^{t,0})$):

$\mathbb{E}[\|D_{\bullet}^{t,k}-D_{\bullet}(\mathbf{u}^{t,k})\|^{2}]=\sum_{r=1}^{k}\mathbb{E}[\|D_{\bullet}^{t,r}-D_{\bullet}^{t,r-1}\|^{2}]-\sum_{r=1}^{k}\mathbb{E}[\|D_{\bullet}(\mathbf{u}^{t,r})-D_{\bullet}(\mathbf{u}^{t,r-1})\|^{2}]\enspace.$

∎

###  A.4 Technical lemmas

######  Lemma A.1.

There exists a constant $L_{*}$ such that for any $x_{1}$,
$x_{2}\in\mathbb{R}^{d}$, we have

$\displaystyle\|z^{*}(x_{1})-z^{*}(x_{2})\|\leq
L_{*}\|x_{1}-x_{2}\|\quad\text{and}\quad\|v^{*}(x_{1})-v^{*}(x_{2})\|\leq
L_{*}\|x_{1}-x_{2}\|$

###### Proof.

The Jacobian of $z^{*}$ reads
$\mathrm{d}z^{*}(x)=[\nabla^{2}_{11}G(z^{*}(x),x)]^{-1}\nabla^{2}_{12}G(z^{*}(x),x)$.
By $\mu_{G}$-strong convexity and $L^{G}_{1}$-smoothness of $G$, we have
$\|\mathrm{d}z^{*}(x)\|\leq\frac{L^{G}_{1}}{\mu_{G}}$ which implies that
$z^{*}$ is $\frac{L^{G}_{1}}{\mu_{G}}$ Lipschtiz.

For $v^{*}$ we do the computation directly:

$\displaystyle\|v^{*}(x_{1})-v^{*}(x_{2})\|$
$\displaystyle=\|[\nabla^{2}_{11}G(z^{*}(x_{1}),x_{1})]^{-1}\nabla_{1}F(z^{*}(x_{1}),x_{1})$
$\displaystyle\qquad-[\nabla^{2}_{11}G(z^{*}(x_{2}),x_{2})]^{-1}\nabla_{1}F(z^{*}(x_{2}),x_{2})\|$
$\displaystyle\leq\|[\nabla^{2}_{11}G(z^{*}(x_{1}),x_{1})]^{-1}(\nabla_{1}F(z^{*}(x_{1}),x_{1})-\nabla_{1}F(z^{*}(x_{2}),x_{2}))\|$
$\displaystyle\qquad+\|([\nabla^{2}_{11}G(z^{*}(x_{1}),x_{1})-[\nabla^{2}_{11}G(z^{*}(x_{2}),x_{2})]^{-1}]^{-1}\nabla_{1}F(z^{*}(x_{2}),x_{2})\|$
$\displaystyle\leq\left(\frac{1}{\mu_{G}}L^{F}_{1}+L^{G}_{2}L^{F}_{0}\right)\|(z^{*}(x_{1}),x_{1})-(z^{*}(x_{2}),x_{2})\|$
$\displaystyle\leq\left(\frac{1}{\mu_{G}}L^{F}_{1}+L^{G}_{2}L^{F}_{0}\right)(\|z^{*}(x_{1})-z^{*}(x_{2})\|+\|x_{1}-x_{2}\|)$
$\displaystyle\leq\frac{L^{G}_{1}}{\mu_{G}}\left(\frac{1}{\mu_{G}}L^{F}_{1}+L^{G}_{2}L^{F}_{0}\right)\|x_{1}-x_{2}\|$

Then taking
$L_{*}=\max\left[\frac{L^{G}_{1}}{\mu_{G}}\left(\frac{1}{\mu_{G}}L^{F}_{1}+L^{G}_{2}L^{F}_{0}\right),\frac{L^{G}_{1}}{\mu_{G}}\right]$
concludes the proof. ∎

######  Lemma A.2.

Let us consider the update directions
$D_{z}^{t,k}=\mathbf{\Delta}_{z}^{t,k}/\rho$,
$D_{v}^{t,k}=\mathbf{\Delta}_{v}^{t,k}/\rho$ and
$D_{x}^{t,k}=\mathbf{\Delta}_{x}^{t,k}/\gamma$ where
$\mathbf{\Delta}_{z}^{t,k}$, $\mathbf{\Delta}_{v}^{t,k}$ and
$\mathbf{\Delta}_{x}^{t,k}$ verify Equations (7) to (9). Then it holds

$\displaystyle\mathbb{E}[\|D_{z}^{t,k}-D_{z}(\mathbf{u}^{t,k})\|^{2}]$
$\displaystyle\leq\sum_{r=1}^{k}L^{G}_{1}(\rho^{2}\mathbb{E}[\|D^{t,r-1}_{z}\|^{2}]+\gamma^{2}\mathbb{E}[\|D^{t,r-1}_{z}\|^{2}])$
$\displaystyle\mathbb{E}[\|D_{v}^{t,k}-D_{v}(\mathbf{u}^{t,k})\|^{2}]$
$\displaystyle\leq
4\rho^{2}\left((L^{G}_{2}R)^{2}+(L^{F}_{1})^{2}\right)\sum_{r=1}^{k}\mathbb{E}[\|D^{t,r-1}_{z}\|^{2}]+4\rho^{2}(L^{G}_{1})^{2}\sum_{r=1}^{k}\mathbb{E}[\|\mathcal{G}^{t,r-1}_{v}\|^{2}]$
$\displaystyle\qquad+4\gamma^{2}\left((L^{G}_{2}R)^{2}+(L^{F}_{1})^{2}\right)\sum_{r=1}^{k}\mathbb{E}[\|D_{x}^{t,r-1}\|^{2}]$
$\displaystyle\mathbb{E}[\|D_{x}^{t,k}-D_{x}(\mathbf{u}^{t,k})\|^{2}]$
$\displaystyle\leq
4\rho^{2}\left((L^{G}_{2}R)^{2}+(L^{F}_{1})^{2}\right)\sum_{r=1}^{k}\mathbb{E}[\|D^{t,r-1}_{z}\|^{2}]+4\rho^{2}(L^{G}_{1})^{2}\sum_{r=1}^{k}\mathbb{E}[\|\mathcal{G}^{t,r-1}_{v}\|^{2}]$
$\displaystyle\qquad+4\gamma^{2}\left((L^{G}_{2}R)^{2}+(L^{F}_{1})^{2}\right)\sum_{r=1}^{k}\mathbb{E}[\|D_{x}^{t,r-1}\|^{2}]\enspace.$

###### Proof.

Direction $D_{z}$

We start from Proposition 3.2.

$\displaystyle\mathbb{E}[\|D^{t,k}_{z}-D_{z}(\mathbf{u}^{t,k})\|^{2}]$
$\displaystyle=\mathbb{E}[\|D^{t,k}_{z}-\nabla_{1}G(z^{t,k},x^{t,k})\|^{2}]$
$\displaystyle=\sum_{r=1}^{k}\mathbb{E}[\|D^{t,r}_{z}-D^{t,r-1}_{z}\|^{2}]-\sum_{r=1}^{k}\mathbb{E}[\|\nabla_{1}G(z^{t,r},x^{t,r})-\nabla_{1}G(z^{t,r-1},x^{t,r-1})\|^{2}]$
$\displaystyle\leq\sum_{r=1}^{k}\mathbb{E}[\|D^{t,r}_{z}-D^{t,r-1}_{z}\|^{2}]$
$\displaystyle\leq\sum_{r=1}^{k}L^{G}_{1}(\rho^{2}\mathbb{E}[\|D^{t,r-1}_{z}\|^{2}]+\gamma^{2}\mathbb{E}[\|D^{t,r-1}_{z}\|^{2}])$

where the last inequality comes from the smoothness of each $G_{i}$.

Direction $D_{v}$ For $D_{v}$, the proof is almost the same. Proposition 3.2
gives us

$\displaystyle\mathbb{E}[\|D^{t,k}_{v}-D_{v}(\mathbf{u}^{t,k})\|^{2}]\leq\sum_{r=1}^{k}\mathbb{E}[\|D^{t,r}_{v}-D^{t,r-1}_{v}\|^{2}]\enspace.$

Then, using the boundedness of $v$ and regularity of each $G_{i}$ and $F_{j}$,
we have

$\displaystyle\mathbb{E}[\|D^{t,r}_{v}-D^{t,r-1}_{v}\|^{2}]$
$\displaystyle\leq
2(\mathbb{E}[\|\nabla^{2}_{11}G_{i}(z^{t,r},x^{t,r})v^{t,r}-\nabla^{2}_{11}G_{i}(z^{t,r-1},x^{t,r-1})v^{t,r-1}\|^{2}]$
$\displaystyle\quad+\mathbb{E}[\|\nabla_{2}F_{j}(z^{t,r},x^{t,r})-\nabla_{2}F_{j}(z^{t,r-1},x^{t,r-1})\|^{2}])$
$\displaystyle\leq
4(\mathbb{E}[\|\nabla^{2}_{11}G_{i}(z^{t,r},x^{t,r})(v^{t,r}-v^{t,r-1})\|^{2}]$
$\displaystyle\quad+\mathbb{E}[\|(\nabla^{2}_{11}G_{i}(z^{t,r},x^{t,r})-\nabla^{2}_{11}G_{i}(z^{t,r-1},x^{t,r-1}))v^{t,r-1}\|^{2}]$
$\displaystyle\quad+(L^{F}_{1})^{2}(\gamma^{2}\mathbb{E}[\|D^{t,r-1}_{z}\|]+\rho^{2}\mathbb{E}[\|D^{t,r-1}_{x}\|^{2}]))$
$\displaystyle\leq
4((L^{G}_{1})^{2}\rho^{2}\mathbb{E}[\|\mathcal{G}^{t,r-1}_{v}\|^{2}]$
$\displaystyle\quad+(L^{G}_{2})^{2}R^{2}(\rho^{2}\mathbb{E}[\|D^{t,r-1}_{z}\|]+\gamma^{2}\mathbb{E}[\|D^{t,r-1}_{x}\|^{2}])$
$\displaystyle\quad+(L^{F}_{1})^{2}(\rho^{2}\mathbb{E}[\|D^{t,r-1}_{z}\|]+\gamma^{2}\mathbb{E}[\|D^{t,r-1}_{x}\|^{2}]))$
$\displaystyle\leq
4\rho^{2}\left((L^{G}_{2}R)^{2}+(L^{F}_{1})^{2}\right)\mathbb{E}[\|D^{t,r-1}_{z}\|^{2}]+4\rho^{2}(L^{G}_{1})^{2}\mathbb{E}[\|\mathcal{G}^{t,r-1}_{v}\|^{2}]$
$\displaystyle\qquad+4\gamma^{2}\left((L^{G}_{2}R)^{2}+(L^{F}_{1})^{2}\right)\mathbb{E}[\|D_{x}^{t,r-1}\|^{2}]\enspace.$

Direction $D_{x}$ The proof is the same as the proof for $D_{v}$. ∎

###  A.5 Proof of Lemma 3.3

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

######  Lemma A.3.

The function $\phi_{z}$ has $\Lambda_{z}$-Lipschitz continuous gradient on
$\mathbb{R}^{p}\times\mathbb{R}^{d}$, for some constant $\Lambda_{z}$.

###### Proof.

For any $(z,x)\in\mathbb{R}^{p}\times\mathbb{R}^{d}$, we have

$\nabla_{1}\phi_{z}(z,x)=\nabla_{1}G(z,x)\text{ and
}\nabla_{2}\phi_{z}(z,x)=\nabla_{2}G(z,x)-\nabla_{2}G(z^{*}(x),x)\enspace.$

Let us consider $(z,x)\in\mathbb{R}^{p}\times\mathbb{R}^{d}$ and
$(z^{\prime},x^{\prime})\in\mathbb{R}^{p}\times\mathbb{R}^{d}$. Since $\nabla
G$ is $L^{G}_{1}$-Lipschitz continuous, we have directly

$\|\nabla_{1}\phi_{z}(z,x)-\nabla_{1}\phi_{z}(z^{\prime},x^{\prime})\|\leq
L^{G}_{1}\|(z,x)-(z^{\prime},x^{\prime})\|\enspace.$

Moreover, we have

$\displaystyle\|\nabla_{2}\phi_{z}(z,x)-\nabla_{2}\phi_{z}(z^{\prime},x^{\prime})\|$
$\displaystyle\leq\|\nabla_{2}G(z,x)-\nabla_{2}G(z^{\prime},x^{\prime})\|$
$\displaystyle\qquad+\|\nabla_{2}G(z^{*}(x),x)-\nabla_{2}G(z^{*}(x^{\prime}),x^{\prime})\|$
$\displaystyle\leq
L^{G}_{1}\|(z,x)-(z^{\prime},x^{\prime})\|+L^{G}_{1}\|(z^{*}(x),x)-(z^{*}(x^{\prime}),x^{\prime})\|$
$\displaystyle\leq
L^{G}_{1}\|(z,x)-(z^{\prime},x^{\prime})\|+L^{G}_{1}\|(z^{*}(x),x)-(z^{*}(x^{\prime}),x^{\prime})\|$
$\displaystyle\leq
L^{G}_{1}\|(z,x)-(z^{\prime},x^{\prime})\|+L^{G}_{1}(\|z^{*}(x)-z^{*}(x^{\prime})\|+\|x-x^{\prime}\|)\enspace.$

From Lemma A.1, $z^{*}$ is $L_{*}$ Lipschitz continuous, so

$\displaystyle\|\nabla_{2}\phi_{z}(z,x)-\nabla_{2}\phi_{z}(z^{\prime},x^{\prime})\|$
$\displaystyle\leq
L^{G}_{1}\|(z,x)-(z^{\prime},x^{\prime})\|+L^{G}_{1}(\|z^{*}(x)-z^{*}(x^{\prime})\|+\|x-x^{\prime}\|)$
$\displaystyle\leq
L^{G}_{1}\|(z,x)-(z^{\prime},x^{\prime})\|+L^{G}_{1}(L_{*}+1)\|x-x^{\prime}\|$
$\displaystyle\leq
L^{G}_{1}(L_{*}+2)\|(z,x)-(z^{\prime},x^{\prime})\|\enspace.$

As a consequence

$\displaystyle\|\nabla\phi_{z}(z,x)-\nabla\phi_{z}(z^{\prime},x^{\prime})\|$
$\displaystyle\leq\|\nabla_{1}\phi_{z}(z,x)-\nabla_{1}\phi_{z}(z^{\prime},x^{\prime})\|+\|\nabla_{2}\phi_{z}(z,x)-\nabla_{2}\phi_{z}(z^{\prime},x^{\prime})\|$
$\displaystyle\leq
L^{G}_{1}(L_{*}+3)\|(z,x)-(z^{\prime},x^{\prime})\|\enspace.$

Hence, $\phi_{z}$ is $\Lambda_{z}$ smooth with
$\Lambda_{z}=L^{G}_{1}(L_{*}+3)$. ∎

We can now turn to the proof of Lemma 3.3.

###### Proof.

The smoothness of $\phi_{z}$ provides us the following upper bound

$\displaystyle\phi_{z}(z^{t,k+1},x^{t,k+1})$
$\displaystyle\leq\phi_{z}(z^{t,k},x^{t,k})-\rho\langle
D_{z}^{t,k},\nabla_{1}G(z^{t,k},x^{t,k})\rangle+\frac{\Lambda_{z}}{2}\rho^{2}\|D_{z}^{t,k}\|^{2}$
(14) $\displaystyle\qquad-\gamma\langle
D_{x}^{t,k},\nabla_{2}G(z^{t,k},x^{t,k})-\nabla_{2}G(z^{*}(x^{t,k}),x^{t,k})\rangle+\frac{\Lambda_{z}}{2}\gamma^{2}\|D_{x}^{t,k}\|^{2}\enspace.$

Using the equality $\langle
a,b\rangle=\frac{1}{2}(\|a\|^{2}+\|b\|^{2}-\|a-b\|^{2})$, we get

$\displaystyle-\langle
D_{z}^{t,k},\nabla_{1}G(z^{t,k},x^{t,k})\rangle+\frac{\Lambda_{z}}{2}\rho\|D_{z}^{t,k}\|^{2}$
$\displaystyle=\frac{1}{2}(\|D^{t,k}_{z}-\nabla_{1}G(z^{t,k},x^{t,k})\|^{2}$
(15)
$\displaystyle\qquad-\|\nabla_{1}G(z^{t,k},x^{t,k})\|^{2}-\left(1-\Lambda_{z}\rho\right)\|D^{t,k}_{z}\|^{2})\enspace.$

Plugging Equation 15 into Equation 14 and tacking the expectation
conditionally to the past iterates yields

$\displaystyle\mathbb{E}_{t,k}[\phi_{z}^{t,k+1}]$
$\displaystyle\leq\phi_{z}^{t,k}+\frac{\rho}{2}\mathbb{E}_{t,k}[\|D^{t,k}_{z}-\nabla_{1}G(z^{t,k},x^{t,k})\|^{2}]$
(16)
$\displaystyle\qquad-\frac{\rho}{2}\|\nabla_{1}G(z^{t,k},x^{t,k})\|^{2}-\frac{\rho}{2}\left(1-\Lambda_{z}\rho\right)\mathbb{E}_{t,k}[\|D^{t,k}_{z}\|^{2}]$
$\displaystyle\qquad-\gamma\langle\mathbb{E}_{t,k}[D_{x}^{t,k}],\nabla_{2}G(z^{t,k},x^{t,k})-\nabla_{2}G(z^{*}(x^{t,k}),x^{t,k})\rangle+\frac{\Lambda_{z}}{2}\gamma^{2}\mathbb{E}_{t,k}[\|D_{x}^{t,k}\|^{2}]\enspace.$

From Young inequality, we have for any $c>0$

$\displaystyle\langle\mathbb{E}_{t,k}[D_{x}^{t,k}],\nabla_{2}G(z^{t,k},x^{t,k})-\nabla_{2}G(z^{*}(x^{t,k}),x^{t,k})\rangle$
$\displaystyle\leq\frac{1}{2c}\|\mathbb{E}_{t,k}[D_{x}^{t,k}]\|^{2}$ (17)
$\displaystyle\qquad+\frac{c}{2}\|\nabla_{2}G(z^{t,k},x^{t,k})-\nabla_{2}G(z^{*}(x^{t,k}),x^{t,k})\|^{2}$

The smoothness of $G$ and strong convexity give us

$\displaystyle\|\nabla_{2}G(z^{t,k},x^{t,k})-\nabla_{2}G(z^{*}(x^{t,k}),x^{t,k})\|^{2}\leq
L^{G}_{1}\|z^{t,k}-z^{*}(x^{t,k})\|^{2}\leq\frac{2L^{G}_{1}}{\mu_{G}}\phi_{z}(z^{t,k},x^{t,k})$
(18)

Let us denote $L^{\prime}=\frac{L^{G}_{1}}{\mu_{G}}$. Plugging Inequalities
(17) and (18) into Equation 16 yields

$\displaystyle\mathbb{E}_{t,k}[\phi_{z}(z^{t,k+1},x^{t,k+1})]$
$\displaystyle\leq(1+cL^{\prime}\gamma)\phi_{z}(z^{t,k+1},x^{t,k+1})-\frac{\rho}{2}\mathbb{E}_{t,k}[\|\nabla_{1}G(z^{t,k},x^{t,k})\|^{2}]$
(19)
$\displaystyle\qquad+\frac{\rho}{2}\mathbb{E}_{t,k}[\|D^{t,k}_{z}-\nabla_{1}G(z^{t,k},x^{t,k})\|^{2}]-\frac{\rho}{2}\left(1-\Lambda_{z}\rho\right)\mathbb{E}_{t,k}[\|D^{t,k}_{z}\|^{2}]$
$\displaystyle\qquad+\frac{\gamma}{2c}\|\mathbb{E}_{t,k}[D^{t,k}_{x}]\|^{2}+\frac{\Lambda_{z}}{2}\gamma^{2}\mathbb{E}_{t,k}[\|D^{t,k}_{x}\|^{2}]$

From Lemma A.2, we have

$\displaystyle\mathbb{E}[\|D^{t,k}_{z}-\nabla_{1}G(z^{t,k},x^{t,k})\|^{2}]$
$\displaystyle\leq\sum_{r=1}^{k}L^{G}_{1}(\rho^{2}\mathbb{E}[\|D^{t,r-1}_{z}\|^{2}]+\gamma^{2}\mathbb{E}[\|D^{t,r-1}_{z}\|^{2}])\enspace.$

Taking the total expectation and plugging the previous inequality into
Equation 19 yields

$\displaystyle\phi_{z}^{t,k+1}$
$\displaystyle\leq(1+cL^{\prime}\gamma)\phi^{t,k}+\frac{L^{G}_{1}}{2}\sum_{r=1}^{k}(\rho^{3}\mathbb{E}[\|D^{t,r-1}_{z}\|^{2}]+\gamma^{2}\rho\mathbb{E}[\|D^{t,r-1}_{x}\|^{2}])$
(20)
$\displaystyle\qquad-\frac{\rho}{2}\mathbb{E}[\|\nabla_{1}G(z^{t,k},x^{t,k})\|^{2}]-\frac{\rho}{2}\left(1-\Lambda_{z}\rho\right)\mathbb{E}[\|D^{t,k}_{z}\|^{2}]$
$\displaystyle\qquad+\frac{\gamma}{2c}\mathbb{E}[\|\mathbb{E}[D^{t,k}_{x}]\|^{2}]+\frac{\Lambda_{z}}{2}\gamma^{2}\mathbb{E}[\|D^{t,k}_{x}\|^{2}]$

Since $G$ is $\mu_{G}$-strongly convex with respect to $z$, Polyak-Łojasiewicz
inequality holds:

$\|\nabla_{1}G(z^{t,k},x^{t,k})\|^{2}\geq 2\mu_{G}\phi_{z}(z^{t,k},x^{t,k})$

As a consequence, Equation 20 becomes

$\displaystyle\phi_{z}^{t,k+1}$
$\displaystyle\leq\left(1+cL^{\prime}\gamma-\mu_{G}\rho\right)\phi^{t,k}+\frac{L^{G}_{1}}{2}\sum_{r=1}^{k}(\rho^{3}\mathbb{E}[\|D^{t,r-1}_{z}\|^{2}]+\gamma^{2}\rho\mathbb{E}[\|D^{t,r-1}_{x}\|^{2}])$
$\displaystyle\qquad-\frac{\rho}{2}\left(1-\Lambda_{z}\rho\right)\mathbb{E}[\|D^{t,k}_{z}\|^{2}]+\frac{\gamma}{2c}\mathbb{E}[\|\mathbb{E}[D^{t,k}_{x}]\|^{2}]+\frac{\Lambda_{z}}{2}\gamma^{2}\mathbb{E}[\|D^{t,k}_{x}\|^{2}]$

Taking $c=\frac{\mu_{G}\rho}{2L^{\prime}\gamma}$ yields

$\displaystyle\phi_{z}^{t,k+1}$
$\displaystyle\leq\left(1-\frac{\mu_{G}}{2}\rho\right)\phi^{t,k}+\frac{L^{G}_{1}}{2}\sum_{r=1}^{k}(\rho^{3}\mathbb{E}[\|D^{t,r-1}_{z}\|^{2}]+\gamma^{2}\rho\mathbb{E}[\|D^{t,r-1}_{x}\|^{2}])$
$\displaystyle\qquad-\frac{\rho}{2}\left(1-\Lambda_{z}\rho\right)\mathbb{E}[\|D^{t,k}_{z}\|^{2}]+\frac{L^{\prime}}{\mu_{G}}\frac{\gamma^{2}}{\rho}\mathbb{E}[\|\mathbb{E}[D^{t,k}_{x}]\|^{2}]+\frac{\Lambda_{z}}{2}\gamma^{2}\mathbb{E}[\|D^{t,k}_{x}\|^{2}]$

For the term $\mathbb{E}[\|\mathbb{E}_{t,k}[D^{t,k}_{z}]\|^{2}]$, we have

$\displaystyle\mathbb{E}[\|\mathbb{E}_{t,k}[D^{t,k}_{x}]\|^{2}]$
$\displaystyle=\mathbb{E}[\|D_{x}(z^{t,k},v^{t,k},x^{t,k})-D_{x}(z^{t,k-1},v^{t,k-1},x^{t,k-1})+D^{t,k-1}_{x}\|^{2}]$
$\displaystyle=\mathbb{E}[\|D_{x}(z^{t,k},v^{t,k},x^{t,k})-D_{x}(z^{t,k-1},v^{t,k-1},x^{t,k-1})-\mathbb{E}[D^{t,k-1}_{x}]\|^{2}]$
$\displaystyle\quad+\mathbb{E}[\|D^{t,k-1}_{x}-\mathbb{E}[D^{t,k-1}_{x}]\|^{2}]$
$\displaystyle=\mathbb{E}[\|D_{x}(z^{t,k},v^{t,k},x^{t,k})\|^{2}]$ (21)
$\displaystyle\qquad+\mathbb{E}[\|D^{t,k-1}_{x}-D_{x}(z^{t,k-1},v^{t,k-1},x^{t,k-1})\|^{2}]\enspace.$

Using Lemma A.2, we get

$\displaystyle\mathbb{E}[\|D^{t,k-1}_{x}-D_{x}(\mathbf{u}^{t,k-1})\|^{2}]$
$\displaystyle\leq
4\rho^{2}\left((L^{G}_{2}R)^{2}+(L^{F}_{1})^{2}\right)\sum_{r=1}^{k-1}\mathbb{E}[\|D^{t,r-1}_{z}\|^{2}]$
$\displaystyle\qquad+4\rho^{2}(L^{G}_{1})^{2}\sum_{r=1}^{k-1}\mathbb{E}[\|\mathcal{G}^{t,r-1}_{v}\|^{2}]$
$\displaystyle\qquad+4\gamma^{2}\left((L^{G}_{2}R)^{2}+(L^{F}_{1})^{2}\right)\sum_{r=1}^{k-1}\mathbb{E}[\|D_{x}^{t,r-1}\|^{2}]\enspace.$

Putting all together yields

$\displaystyle\phi_{z}^{t,k+1}$
$\displaystyle\leq\left(1-\frac{\mu_{G}}{2}\rho\right)\phi^{t,k}-\frac{\rho}{2}\left(1-\Lambda_{z}\rho\right)\mathbb{E}[\|D^{t,k}_{z}\|^{2}]+\frac{\Lambda_{z}}{2}\gamma^{2}\mathbb{E}[\|D^{t,k}_{x}\|^{2}]$
(22)
$\displaystyle\qquad+\frac{L^{\prime}}{\mu_{G}}\frac{\gamma^{2}}{\rho}\mathbb{E}[\|D^{t,k}_{x}(\mathbf{u}^{t,k})\|^{2}]+4(L^{G}_{1})^{2}\frac{L^{\prime}}{\mu_{G}}\gamma^{2}\rho\sum_{r=1}^{k}\mathbb{E}[\|\mathcal{G}^{t,r-1}_{v}\|^{2}]$
$\displaystyle\qquad+\rho\left[\rho^{2}\frac{L^{G}_{1}}{2}+\frac{4(L^{G}_{2}R)^{2}L^{\prime}}{\mu_{G}}\gamma^{2}+\frac{4(L^{F}_{1})^{2}L^{\prime}}{\mu_{G}}\gamma^{2}\right]\sum_{r=1}^{k}\mathbb{E}[\|D^{t,r-1}_{z}\|^{2}]$
$\displaystyle\qquad+\gamma^{2}\left[\rho\frac{L^{G}_{1}}{2}+4(L^{G}_{2}R)^{2}\frac{L^{\prime}}{\mu_{G}}\frac{\gamma^{2}}{\rho}+4(L^{F}_{1})^{2}\frac{L^{\prime}}{\mu_{G}}\frac{\gamma^{2}}{\rho}\right]\sum_{r=1}^{k}\mathbb{E}[\|D^{t,r-1}_{x}\|^{2}]$

By assumption, $\gamma\leq C_{z}\rho$, with
$C_{z}=\sqrt{\frac{\mu_{G}L^{G}_{1}}{8L^{\prime}((L^{G}_{2}R)^{2}+(L^{F}_{1})^{2})}}$
therefore

$\displaystyle\phi_{z}^{t,k+1}$
$\displaystyle\leq\left(1-\frac{\mu_{G}}{2}\rho\right)\phi^{t,k}_{z}-\frac{\rho}{2}\left(1-\Lambda_{z}\rho\right)\mathbb{E}[\|D^{t,k}_{z}\|^{2}]+\frac{\Lambda_{z}}{2}\gamma^{2}\mathbb{E}[\|D^{t,k}_{x}\|^{2}]$
$\displaystyle\qquad+\frac{L^{\prime}}{\mu_{G}}\frac{\gamma^{2}}{\rho}\mathbb{E}[\|D^{t,k}_{x}(\mathbf{u}^{t,k})\|^{2}]+\rho^{3}L^{G}_{1}\sum_{r=1}^{k}\mathbb{E}[\|D^{t,r-1}_{z}\|^{2}]$
$\displaystyle\qquad+4(L^{G}_{1})^{2}\frac{L^{\prime}}{\mu_{G}}\gamma^{2}\rho\sum_{r=1}^{k}\mathbb{E}[\|\mathcal{G}^{t,r-1}_{v}\|^{2}]+\gamma^{2}\rho
L^{G}_{1}\sum_{r=1}^{k}\mathbb{E}[\|D^{t,r-1}_{x}\|^{2}]$
$\displaystyle\leq\left(1-\frac{\mu_{G}}{2}\rho\right)\phi^{t,k}-\frac{\rho}{2}\left(1-\Lambda_{z}\rho\right)V_{z}^{t,k}+\frac{\Lambda_{z}}{2}\gamma^{2}V_{x}^{t,k}+\overline{\beta}_{zx}\frac{\gamma^{2}}{\rho}\mathbb{E}[\|D^{t,k}_{x}(\mathbf{u}^{t,k})\|^{2}]$
$\displaystyle\qquad+\rho^{3}\beta_{zz}\mathcal{V}_{z}^{t,k}+\gamma^{2}\rho\beta_{zv}\mathcal{V}_{v}^{t,k}+\gamma^{2}\rho\beta_{zx}\mathcal{V}_{x}^{t,k}$

with $\beta_{zz}=L^{G}_{1}$,
$\beta_{zv}=\frac{4(L^{G}_{1})^{2}L^{\prime}}{\mu_{G}}$,
$\beta_{zx}=L^{G}_{1}$ and $\overline{\beta}_{zx}=\frac{L^{\prime}}{\mu_{G}}$.
∎

###  A.6 Proof of Lemma 3.4

Recall that we denote
$\Psi(z,v,x)=\frac{1}{2}v^{\top}\nabla^{2}_{11}G(z,x)v+\nabla_{1}F(z,x)^{\top}v$
and $\phi_{v}(v,x)=\Psi(z^{*}(x),v,x)-\Psi(z^{*}(x),v^{*}(x),x)$. As for Lemma
3.3, the key property we need is the smoothness of $\phi_{v}$. The derivatives
of $\phi_{v}$ involve the third derivative of $G$. For a tensor
$T\in\mathbb{R}^{p_{1}\times p_{2}\times p_{3}}$ and a vector
$a\in\mathbb{R}^{p_{3}}$ we denote $(T|a)$ the matrix in
$\mathbb{R}^{p_{1}\times p_{2}}$ defined by:

$\displaystyle(T|a)=\left[\sum_{k=1}^{p_{3}}T_{i,j,k}a_{k}\right]_{\begin{subarray}{c}1\leq
i\leq p_{1}\\\ 1\leq j\leq p_{2}\end{subarray}}\enspace.$

######  Lemma A.4.

The function $\phi_{v}$ has $\Lambda_{v}$-Lipschitz continuous gradient on
$\Gamma\times\mathbb{R}^{d}$, for some constant $\Lambda_{v}$.

###### Proof.

For any $(v,x)\in\Gamma\times\mathbb{R}^{d}$, we have

$\nabla_{1}\phi_{v}(v,x)=D_{v}(z^{*}(x),v,x)$

and

$\displaystyle\nabla_{2}\phi_{v}(v,x)$
$\displaystyle=(\mathrm{d}z^{*}(x))^{\top}\left[\frac{1}{2}(\nabla^{3}_{111}G(z^{*}(x),x)|v)v-\frac{1}{2}(\nabla^{3}_{111}G(z^{*}(x),x)|v^{*}(x))v^{*}(x)\right.$
$\displaystyle\qquad\left.+\nabla_{11}^{2}F(z^{*}(x),x)v-\nabla_{11}^{2}F(z^{*}(x),x)v^{*}(x)\right]$
$\displaystyle\qquad+\left[\frac{1}{2}(\nabla^{3}_{211}G(z^{*}(x),x)|v)v-\frac{1}{2}(\nabla^{3}_{211}G(z^{*}(x),x)|v^{*}(x))v^{*}(x)\right.$
$\displaystyle\qquad\left.+\nabla_{21}^{2}F(z^{*}(x),x)v-\nabla_{21}^{2}F(z^{*}(x),x)v^{*}(x)\right]\enspace.$

Let us consider $(v,x)\in\Gamma\times\mathbb{R}^{d}$ and
$(v^{\prime},x^{\prime})\in\Gamma\times\mathbb{R}^{d}$. We have

$\displaystyle\|\nabla_{1}\phi_{v}(v,x)-\nabla_{1}\phi_{v}(v^{\prime},x^{\prime})\|$
$\displaystyle\leq\|\nabla_{11}^{2}G(z^{*}(x),x)v-\nabla_{11}^{2}G(z^{*}(x^{\prime}),x^{\prime})v^{\prime}\|$
$\displaystyle\qquad+\|\nabla_{1}F(z^{*}(x),x)-\nabla_{1}F(z^{*}(x^{\prime}),x^{\prime})\|$

For the first term,

$\displaystyle\|\nabla_{11}^{2}G(z^{*}(x),x)v-\nabla_{11}^{2}G(z^{*}(x^{\prime}),x^{\prime})v^{\prime}\|$
$\displaystyle\leq\|\nabla_{11}^{2}G(z^{*}(x),x)(v-v^{\prime})\|$
$\displaystyle\qquad+\|(\nabla_{11}^{2}G(z^{*}(x),x)-\nabla_{11}^{2}G(z^{*}(x^{\prime}),x^{\prime}))v^{\prime}\|$
$\displaystyle\qquad+\|\nabla_{11}^{2}G(z^{*}(x^{\prime}),x^{\prime})(v-v^{\prime})\|$
$\displaystyle\leq
2L^{G}_{1}\|v-v^{\prime}\|+L^{G}_{2}(L_{*}+1)\|v^{\prime}\|\|x-x^{\prime}\|$
$\displaystyle\leq[2L^{G}_{1}+L^{G}_{2}(L_{*}+1)R]\|(v,x)-(v^{\prime},x^{\prime})\|$

For the second terms, we use the smoothness of $F$ and the Lipschitz
continuity of $z^{*}$ (Lemma A.1):

$\displaystyle\|\nabla_{1}F(z^{*}(x),x)-\nabla_{1}F(z^{*}(x^{\prime}),x^{\prime})\|$
$\displaystyle\leq L^{F}_{1}\|(z^{*}(x),x)-(z^{*}(x^{\prime}),x^{\prime})\|$
$\displaystyle\leq L^{F}_{1}(\|z^{*}(x)-z^{*}(x^{\prime})\|+\|x-x^{\prime}\|)$
$\displaystyle\leq L^{F}_{1}(L_{*}+1)\|x-x^{\prime}\|$ $\displaystyle\leq
L^{F}_{1}(L_{*}+1)\|(x,v)-(x^{\prime},v^{\prime})\|\enspace.$

As a consequence

$\|\nabla_{1}\phi_{v}(v,x)-\nabla_{1}\phi_{v}(v^{\prime},x^{\prime})\|\leq\Lambda_{1}\|(v,x)-(v^{\prime},x^{\prime})\|$
(23)

with

$\Lambda_{1}=L^{F}_{1}(L_{*}+1)+2L^{G}_{1}+L^{G}_{2}(L_{*}+1)R\enspace.$ (24)

To prove the Lipschitz continuity of $\nabla_{2}\phi_{v}$, we remark that
$\nabla^{3}_{111}G$, $\nabla^{3}_{211}G$ are Lipschitz and bounded by
assumption. $(v\mapsto v)$ is Lipschitz and bounded on $\Gamma$. Also by Lemma
A.1, $z^{*}$ and $v^{*}$ are Lipschitz and bounded. Finally, $\mathrm{d}z^{*}$
is bounded (Lemma A.1) and Lipschitz according to Chen et al. [2021][Lemma 9].
As a consequence, $\nabla_{2}\phi_{v}$ is $\Lambda_{2}$-Lpischitz for some
constant $\Lambda_{2}>0$. Hence, $\nabla\phi_{v}$ is $\Lambda_{v}$-Lipschitz
continuous with $\Lambda_{v}=\Lambda_{1}+\Lambda_{2}$.

∎

######  Lemma A.5.

Let $t>0$. For $k\in[q-1]$, we have

$0\leq-\left\langle\frac{1}{\rho}(v^{t,k+1}-v^{t,k})+D^{t,k}_{v},v^{t,k+1}-v^{t,k}\right\rangle$

###### Proof.

The function $\iota_{\Gamma}$ beging convex (since $\Gamma$ is convex), let us
consider its sub-differential

$\partial\iota_{\gamma}(v)=\\{\eta\in\mathbb{R}^{p},\forall
v^{\prime}\in\mathbb{R}^{p},\iota_{\Gamma}(v^{\prime})\geq\iota_{\Gamma}(v)+\langle\eta,v^{\prime}-v\rangle\\}$

By definition

$\displaystyle
v^{t,k+1}=\operatorname*{arg\,min}_{v}(\iota_{\Gamma}(v)+\frac{1}{2\rho}\|v-(v^{t,k}-\rho
D^{t,k}_{v})\|^{2})\enspace.$

Using Fermat’s rule, we get

$\displaystyle-\frac{1}{\rho}(v^{t,k+1}-v^{t,k})-D^{t,k}_{v}\in\partial\iota_{\Gamma}(v^{t,k+1})\enspace.$

We can use the definition of the sub-differential with
$\eta=-\frac{1}{\rho}(v^{t,k+1}-v^{t,k})-D^{t,k}_{v}$ to get

$\displaystyle\underbrace{\iota_{\Gamma}(v^{t,k+1})}_{=0}\leq\underbrace{\iota_{\Gamma}(v^{t,k})}_{=0}-\left\langle\frac{1}{\rho}(v^{t,k+1}-v^{t,k})+D^{t,k}_{v},v^{t,k+1}-v^{t,k}\right\rangle\enspace.$

∎

We can now turn to the proof of Lemma 3.4.

###### Proof.

The smoothness of $\phi_{v}$ provides us the following upper bound

$\displaystyle\phi_{v}(v^{t,k+1},x^{t,k+1})$
$\displaystyle\leq\phi_{v}(v^{t,k},x^{t,k})+\langle\Pi_{\Gamma}(v^{t,k}-\rho
D_{v}^{t,k})-v^{t,k},D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k})\rangle$ (25)
$\displaystyle\qquad+\frac{\Lambda_{v}}{2}\rho^{2}\|\Pi_{\Gamma}(v^{t,k}-\rho
D_{v}^{t,k})-v^{t,k}\|^{2}$ $\displaystyle\qquad-\gamma\langle
D_{x}^{t,k},\nabla_{2}\phi_{v}(v^{t,k},x^{t,k})\rangle+\frac{\Lambda_{v}}{2}\gamma^{2}\|D_{x}^{t,k}\|^{2}\enspace.$

Let us denote $\Delta_{\Pi}^{t,k}=\Pi_{\Gamma}(v^{t,k}-\rho
D_{v}^{t,k})-\Pi_{\Gamma}(v^{t,k}-\rho
D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k}))$. Adding and subtracting

$\displaystyle\langle\Pi_{\Gamma}(v^{t,k}-\rho
D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k})-v^{t,k},D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k})\rangle$
$\displaystyle\qquad+\frac{\Lambda_{v}}{2}\|\Pi_{\Gamma}(v^{t,k}-\rho
D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k})-v^{t,k}\|^{2}$

yields

$\displaystyle\phi_{v}(v^{t,k+1},x^{t,k+1})$
$\displaystyle\leq\phi_{v}(v^{t,k},x^{t,k})+\langle\Delta_{\Pi}^{t,k},D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k})\rangle$
(26) $\displaystyle\qquad+\frac{\Lambda_{v}}{2}\|\Pi_{\Gamma}(v^{t,k}-\rho
D_{v}(z^{*}(x^{t}),v^{t,k},x^{t,k}))-v^{t,k}\|^{2}$
$\displaystyle\qquad+\langle\Pi_{\Gamma}(v^{t,k}-\rho
D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k}))-v^{t,k},D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k})\rangle$
$\displaystyle\qquad+\frac{\Lambda_{v}}{2}\|\Delta_{\Pi}^{t,k}\|^{2}+\Lambda_{v}\langle\Delta_{\Pi}^{t,k},\Pi_{\Gamma}(v^{t,k}-\rho
D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k}))-v^{t,k}\rangle$
$\displaystyle\qquad-\gamma\langle
D_{x}^{t,k},\nabla_{2}\phi_{v}(v^{t,k},x^{t,k})\rangle+\frac{\Lambda_{v}}{2}\gamma^{2}\|D_{x}^{t,k}\|^{2}\enspace.$

Taking $\rho\leq\frac{1}{\Gamma}_{v}$ gives

$\displaystyle\phi_{v}(v^{t,k+1},x^{t,k+1})$
$\displaystyle\leq\phi_{v}(v^{t,k},x^{t,k})+\langle\Delta_{\Pi}^{t,k},D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k})\rangle$
(27) $\displaystyle\qquad+\frac{1}{2\rho}\|\Pi_{\Gamma}(v^{t,k}-\rho
D_{v}(z^{*}(x^{t}),v^{t,k},x^{t,k}))-v^{t,k}\|^{2}$
$\displaystyle\qquad+\langle\Pi_{\Gamma}(v^{t,k}-\rho
D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k}))-v^{t,k},D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k})\rangle$
$\displaystyle\qquad+\frac{\Lambda_{v}}{2}\|\Delta_{\Pi}^{t,k}\|^{2}+\Lambda_{v}\langle\Delta_{\Pi}^{t,k},\Pi_{\Gamma}(v^{t,k}-\rho
D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k}))-v^{t,k}\rangle$
$\displaystyle\qquad-\gamma\langle
D_{x}^{t,k},\nabla_{2}\phi_{v}(v^{t,k},x^{t,k})\rangle+\frac{\Lambda_{v}}{2}\gamma^{2}\|D_{x}^{t,k}\|^{2}\enspace.$

Let $\iota_{\Gamma}$ the indicator function of the convex set $\Gamma$.
Similarly to Karimi et al. [2016, Equation 13] we define for any $\alpha>0$
and $v\in\mathbb{R}^{p}$

$\mathcal{D}_{\iota_{\Gamma}}(v,x,\alpha)=-2\alpha\min_{v^{\prime}\in\mathbb{R}^{p}}\left[\langle\nabla_{1}\phi_{v}(v,x),v^{\prime}-v\rangle+\frac{\alpha}{2}\|v^{\prime}-v\|^{2}+\iota_{\Gamma}(v^{\prime})-\iota_{\Gamma}(v)\right]\enspace.$

Hence, for $v\in\Gamma$ and $x\in\mathbb{R}^{d}$, we have

$\displaystyle-\frac{\rho}{2}\mathcal{D}_{\iota_{\Gamma}}\left(v,x,\frac{1}{\rho}\right)$
$\displaystyle=\langle\Pi_{\Gamma}(v-\rho
D_{v}(z^{*}(x),v,x))-v,D_{v}(z^{*}(x),v,x)\rangle$
$\displaystyle\qquad+\frac{1}{2\rho}\|\Pi_{\Gamma}(v-\rho
D_{v}(z^{*}(x),v,x))-v\|^{2}\enspace.$

Therefore, Equation 27 can be written as

$\displaystyle\phi_{v}(v^{t,k+1},x^{t,k+1})$
$\displaystyle\leq\phi_{v}(v^{t,k},x^{t,k})-\frac{\rho}{2}\mathcal{D}_{\iota_{\Gamma}}\left(v^{t,k},x^{t,k},\frac{1}{\rho}\right)$
$\displaystyle\qquad+\langle\Delta_{\Pi}^{t,k},D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k})\rangle$
$\displaystyle\qquad+\frac{\Lambda_{v}}{2}\|\Delta_{\Pi}^{t,k}\|^{2}+\Lambda_{v}\langle\Delta_{\Pi}^{t,k},\Pi_{\Gamma}(v^{t,k}-\rho
D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k}))-v^{t,k}\rangle$
$\displaystyle\qquad-\gamma\langle
D_{x}^{t,k},\nabla_{2}\phi_{v}(v^{t,k},x^{t,k})\rangle+\frac{\Lambda_{v}}{2}\gamma^{2}\|D_{x}^{t,k}\|^{2}\enspace.$

By strong convexity of $\phi_{v}$ with respect top $v$ and smoothness, we have
$\mathcal{D}_{\iota_{\Gamma}}(v^{t,k},x^{t,k},\Lambda_{v})~{}\geq~{}2\mu_{G}\phi_{v}(v^{t,k},x^{t,k})$.
According to Karimi et al. [2016, Lemma 1],
$\mathcal{D}_{\iota_{\Gamma}}(v^{t,k},x^{t,k},\bullet)$ is an increasing
function. As a consequence, since $\Lambda_{v}\leq\frac{1}{\rho}$, we have
$\mathcal{D}_{\iota_{\Gamma}}\left(v^{t,k},x^{t,k},\frac{1}{\rho}\right)\geq
2\mu_{G}\phi_{v}(v^{t,k},x^{t,k})$. This leads to

$\displaystyle\phi_{v}(v^{t,k+1},x^{t,k+1})$
$\displaystyle\leq(1-\rho\mu_{G})\phi_{v}(v^{t,k},x^{t,k})+\langle\Delta_{\Pi}^{t,k},D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k})\rangle$
(28)
$\displaystyle\qquad+\frac{\Lambda_{v}}{2}\|\Delta_{\Pi}^{t,k}\|^{2}+\Lambda_{v}\langle\Delta_{\Pi}^{t,k},\Pi_{\Gamma}(v^{t,k}-\rho
D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k}))-v^{t,k}\rangle$
$\displaystyle\qquad-\gamma\langle
D_{x}^{t,k},\nabla_{2}\phi_{v}(v^{t,k},x^{t,k})\rangle+\frac{\Lambda_{v}}{2}\gamma^{2}\|D_{x}^{t,k}\|^{2}\enspace.$

The non-expansiveness of $\Pi_{\Gamma}$ yields

$\displaystyle\|\Delta_{\Pi}^{t,k}\|\leq\rho\|D^{t,k}_{v}-D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k})\|$
(29)

and

$\displaystyle\|\Pi_{\Gamma}(v^{t,k}-\rho
D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k})-\underbrace{v^{t,k}}_{\in\Gamma}\|$
$\displaystyle=\|\Pi_{\Gamma}(v^{t,k}-\rho
D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k}))-\Pi_{\Gamma}(v^{t,k})\|$
$\displaystyle\leq\rho\|D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k})\|\enspace.$ (30)

Moreover, using Equation 29 and Young Inequality, we have for any $c>0$

$\displaystyle\langle\Delta_{\Pi}^{t,k},D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k})\rangle$
$\displaystyle\leq\frac{c}{2}\|\Delta_{\Pi}\|^{2}+\frac{1}{2c}\|D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{tk,})\|^{2}$
$\displaystyle\leq\frac{c\rho^{2}}{2}\|D^{t,k}_{v}-D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k})\|^{2}$
$\displaystyle\qquad+\frac{1}{2c}\|D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{tk,})-\underbrace{D_{v}(z^{*}(x^{t,k}),v^{*}(x^{t,k}),x^{tk,})}_{=0}\|^{2}$
$\displaystyle\leq\frac{c\rho^{2}}{2}\|D^{t,k}_{v}-D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k})\|^{2}$
(31) $\displaystyle\qquad+\frac{L^{G}_{1}}{\mu_{G}c}\phi_{v}(v^{t,k},x^{t,k})$

Plugging Equation 31 into Equation 28 with
$c=\frac{2L^{G}_{1}}{\mu_{G}^{2}\rho}$ yields

$\displaystyle\phi_{v}(v^{t,k+1},x^{t,k+1})$
$\displaystyle\leq\left(1-\frac{\rho\mu_{G}}{2}\right)\phi_{v}(v^{t,k},x^{t,k})+\frac{L^{G}_{1}\rho}{\mu_{G}^{2}}\|D^{t,k}_{v}-D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k})\|^{2}$
(32)
$\displaystyle\qquad+\frac{\Lambda_{v}}{2}\|\Delta_{\Pi}^{t,k}\|^{2}+\Lambda_{v}\langle\Delta_{\Pi}^{t,k},\Pi_{\Gamma}(v^{t,k}-\rho
D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k}))-v^{t,k}\rangle$
$\displaystyle\qquad-\gamma\langle
D_{x}^{t,k},\nabla_{2}\phi_{v}(v^{t,k},x^{t,k})\rangle+\frac{\Lambda_{v}}{2}\gamma^{2}\|D_{x}^{t,k}\|^{2}\enspace.$

Using Equation 29, Equation 30 and Young Inequality for $d>0$ yields

$\displaystyle\langle\Delta_{\Pi}^{t,k},\Pi_{\Gamma}(v^{t,k}-\rho
D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k}))-v^{t,k}\rangle$
$\displaystyle\leq\frac{d}{2}\|\Delta_{\Pi}^{t,k}\|^{2}$
$\displaystyle\qquad+\frac{1}{2d}\|\Pi_{\Gamma}(v^{t,k}-\rho
D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k}))-v^{t,k}\|^{2}$
$\displaystyle\leq\frac{d\rho^{2}}{2}\|D^{t,k}_{v}-D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k})\|^{2}$
(33)
$\displaystyle\qquad+\frac{\rho^{2}}{2d}\|D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k})\|^{2}$
$\displaystyle\leq\frac{d\rho^{2}}{2}\|D^{t,k}_{v}-D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k})\|^{2}$
(34)
$\displaystyle\qquad+\frac{L_{1}^{G}\rho^{2}}{\mu_{G}d}\phi_{v}(v^{t,k},x^{t,k})\enspace.$

Plugging Equation 34 into Equation 32 with
$d=\frac{4L^{G}_{1}\Lambda_{v})\rho}{\mu_{G}^{2}}$ gives

$\displaystyle\phi_{v}(v^{t,k+1},x^{t,k+1})$
$\displaystyle\leq\left(1-\frac{\rho\mu_{G}}{4}\right)\phi_{v}(v^{t,k},x^{t,k})$
(35)
$\displaystyle\qquad+\left[\frac{L^{G}_{1}\rho}{\mu_{G}^{2}}+\frac{2L^{G}_{1}\Lambda_{v}^{2}\rho^{3}}{\mu_{G}^{2}}\right]\|D^{t,k}_{v}-D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k})\|^{2}$
$\displaystyle\qquad+\frac{\Lambda_{v}}{2}\|\Delta_{\Pi}^{t,k}\|^{2}$
$\displaystyle\qquad-\gamma\langle
D_{x}^{t,k},\nabla_{2}\phi_{v}(v^{t,k},x^{t,k})\rangle+\frac{\Lambda_{v}}{2}\gamma^{2}\|D_{x}^{t,k}\|^{2}\enspace.$

Using once again (29), we get

$\displaystyle\phi_{v}(v^{t,k+1},x^{t,k+1})$
$\displaystyle\leq\left(1-\frac{\rho\mu_{G}}{4}\right)\phi_{v}(v^{t,k},x^{t,k})$
(36)
$\displaystyle\qquad+\left[\frac{L^{G}_{1}\rho}{\mu_{G}^{2}}+\frac{2L^{G}_{1}\Lambda_{v}^{2}\rho^{3}}{\mu_{G}^{2}}+\frac{\Lambda_{v}\rho^{2}}{2}\right]\|D^{t,k}_{v}-D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k})\|^{2}$
$\displaystyle\qquad-\gamma\langle
D_{x}^{t,k},\nabla_{2}\phi_{v}(v^{t,k},x^{t,k})\rangle+\frac{\Lambda_{v}}{2}\gamma^{2}\|D_{x}^{t,k}\|^{2}\enspace.$

By Lemma A.5, we have for any $\alpha>0$

$\displaystyle
0\leq-\alpha\left\langle\frac{1}{\rho}(v^{t,k+1}-v^{t,k})+D^{t,k}_{v},v^{t,k+1}-v^{t,k}\right\rangle\enspace.$

By adding this to Equation 36, we get

$\displaystyle\phi_{v}(v^{t,k+1},x^{t,k+1})$
$\displaystyle\leq\left(1-\frac{\rho\mu_{G}}{4}\right)\phi_{v}(v^{t,k},x^{t,k})$
(37)
$\displaystyle\qquad-\frac{\alpha}{\rho}\|v^{t,k+1}-v^{t,k}\|^{2}-\alpha\langle
D^{t,k}_{v},v^{t,k+1}-v^{t,k}\rangle$
$\displaystyle\qquad+\left[\frac{L^{G}_{1}\rho}{\mu_{G}^{2}}+\frac{2L^{G}_{1}\Lambda_{v}^{2}\rho^{3}}{\mu_{G}^{2}}+\frac{\Lambda_{v}\rho^{2}}{2}\right]\|D^{t,k}_{v}-D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k})\|^{2}$
$\displaystyle\qquad-\gamma\langle
D_{x}^{t,k},\nabla_{2}\phi_{v}(v^{t,k},x^{t,k})\rangle+\frac{\Lambda_{v}}{2}\gamma^{2}\|D_{x}^{t,k}\|^{2}\enspace.$

We can control $-\left\langle D^{t,k}_{v},v^{t,k+1}-v^{t,k}\right\rangle$ by
Cauchy-Schwarz and Young for some $c,d,e,f>0$

$\displaystyle-\left\langle D^{t,k}_{v},v^{t,k+1}-v^{t,k}\right\rangle$
$\displaystyle=-\left\langle
D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k}),\Pi_{\Gamma}(v^{t,k}-\rho
D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k})-v^{t,k}\right\rangle$
$\displaystyle\qquad-\left\langle
D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k}),\Delta_{\Pi}^{t,k}\right\rangle$
$\displaystyle\qquad-\left\langle
D^{t,k}_{v}-D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k}),\Pi_{\Gamma}(v^{t,k}-\rho
D_{v}(z^{*}(x)))-v^{t,k}\right\rangle$ $\displaystyle\qquad-\left\langle
D^{t,k}_{v}-D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k}),\Delta_{\Pi}^{t,k}\right\rangle$
$\displaystyle\leq\frac{c}{2}\|D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k})\|^{2}$
$\displaystyle\qquad+\frac{1}{2c}\|\Pi_{\Gamma}(v^{t,k}-\rho
D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k}))-v^{t,k}\|^{2}$
$\displaystyle\qquad+\frac{d}{2}\|D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k})\|^{2}+\frac{1}{2d}\|\Delta_{\Pi}^{t,k}\|^{2}$
$\displaystyle\qquad+\frac{e}{2}\|D^{t,k}_{v}-D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k})\|^{2}$
$\displaystyle\qquad+\frac{1}{2e}\|\Pi_{\Gamma}(v^{t,k}-\rho
D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k}))-v^{t,k}\|^{2}$
$\displaystyle\qquad+\frac{f}{2}\|D^{t,k}_{v}-D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k})\|^{2}+\frac{1}{2f}\|\Delta_{\Pi}^{t,k}\|^{2}$
$\displaystyle\leq\left(\frac{c+d}{2}+\rho^{2}\left(\frac{1}{2c}+\frac{1}{2e}\right)\right)\|D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k})\|^{2}$
$\displaystyle\qquad+\left(\frac{e+f}{2}+\rho^{2}\left(\frac{1}{2d}+\frac{1}{2f}\right)\right)\|D^{t,k}_{v}-D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k})\|^{2}$
$\displaystyle\leq\left(\frac{c+d}{2}+\rho^{2}\left(\frac{1}{2c}+\frac{1}{2e}\right)\right)\frac{2L^{G}_{1}}{\mu_{G}}\phi_{v}(v^{t,k},x^{t,k})$
$\displaystyle\qquad+\left(\frac{e+f}{2}+\rho^{2}\left(\frac{1}{2d}+\frac{1}{2f}\right)\right)\|D^{t,k}_{v}-D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k})\|^{2}$

Let us take $c=d=e=f=\rho$. We get

$\displaystyle-\left\langle D^{t,k}_{v},v^{t,k+1}-v^{t,k}\right\rangle$
$\displaystyle\leq\frac{4L^{G}_{1}}{\mu_{G}}\rho\phi_{v}(v^{t,k},x^{t,k})+2\rho\|D^{t,k}_{v}-D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k})\|^{2}\enspace.$
(38)

Then, by plugging the last Inequality in Equation 37 and setting
$\alpha=\frac{\mu_{G}^{2}}{32L^{G}_{1}}$, we end up with

$\displaystyle\phi_{v}(v^{t,k+1},x^{t,k+1})$
$\displaystyle\leq\left(1-\frac{\mu_{G}}{8}\rho\right)\phi_{v}(v^{t,k},x^{t,k})-\frac{\alpha}{\rho}\|v^{t,k+1}-v^{t,k}\|^{2}$
$\displaystyle\qquad+\rho\left[\frac{L^{G}_{1}}{\mu_{G}^{2}}+\frac{\mu_{G}^{2}}{16L^{G}_{1}}+\frac{\Lambda_{v}\rho}{2}+\frac{2L^{G}_{1}\Lambda_{v}^{2}\rho^{2}}{\mu_{G}^{2}}\right]\|D^{t,k}_{v}-D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k})\|^{2}$
$\displaystyle\qquad-\gamma\langle
D_{x}^{t,k},\nabla_{2}\phi_{v}(v^{t,k},x^{t,k})\rangle+\frac{\Lambda_{v}}{2}\gamma^{2}\|D_{x}^{t,k}\|^{2}$
$\displaystyle\leq\left(1-\frac{\mu_{G}}{8}\rho\right)\phi_{v}(v^{t,k},x^{t,k})-\frac{\mu_{G}^{2}}{32L^{G}_{1}}\rho\|\mathcal{G}^{t,k}_{v}\|^{2}$
$\displaystyle\qquad+\rho\left[\frac{L^{G}_{1}}{\mu_{G}^{2}}+\frac{\mu_{G}^{2}}{16L^{G}_{1}}+\frac{\Lambda_{v}\rho}{2}+\frac{2L^{G}_{1}\Lambda_{v}^{2}\rho^{2}}{\mu_{G}^{2}}\right]\|D^{t,k}_{v}-D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k})\|^{2}$
$\displaystyle\qquad-\gamma\langle
D_{x}^{t,k},\nabla_{2}\phi_{v}(v^{t,k},x^{t,k})\rangle+\frac{\Lambda_{v}}{2}\gamma^{2}\|D_{x}^{t,k}\|^{2}\enspace.$

Since $\rho\leq
B_{v}\triangleq\left[\frac{L^{G}_{1}}{\mu_{G}^{2}}+\frac{\mu_{G}^{2}}{16L^{G}_{1}}\right]\min\left(\frac{2}{\Lambda_{v}},\frac{\mu_{G}}{\sqrt{2L^{G}_{1}}\Lambda_{v}}\right)$
yields

$\displaystyle\phi_{v}(v^{t,k+1},x^{t,k+1})$
$\displaystyle\leq\left(1-\frac{\mu_{G}}{8}\rho\right)\phi_{v}(v^{t,k},x^{t,k})-\frac{\mu_{G}^{2}}{32L^{G}_{1}}\rho\|\mathcal{G}^{t,k}_{v}\|^{2}$
$\displaystyle\qquad+3\rho\left[\frac{L^{G}_{1}}{\mu_{G}^{2}}+\frac{\mu_{G}^{2}}{16L^{G}_{1}}\right]\|D^{t,k}_{v}-D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k})\|^{2}$
$\displaystyle\qquad-\gamma\langle
D_{x}^{t,k},\nabla_{2}\phi_{v}(v^{t,k},x^{t,k})\rangle+\frac{\Lambda_{v}}{2}\gamma^{2}\|D_{x}^{t,k}\|^{2}$
$\displaystyle\leq\left(1-\frac{\mu_{G}}{8}\rho\right)\phi_{v}(v^{t,k},x^{t,k})-\frac{\mu_{G}^{2}}{32L^{G}_{1}}\rho\|\mathcal{G}^{t,k}_{v}\|^{2}$
(39)
$\displaystyle\qquad+6\rho\left[\frac{L^{G}_{1}}{\mu_{G}^{2}}+\frac{\mu_{G}^{2}}{16L^{G}_{1}}\right]\|D^{t,k}_{v}-D_{v}(\mathbf{u}^{t,k})\|^{2}$
$\displaystyle\qquad+6\rho\left[\frac{L^{G}_{1}}{\mu_{G}^{2}}+\frac{\mu_{G}^{2}}{16L^{G}_{1}}\right]\|D_{v}(\mathbf{u}^{t,k})-D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k})\|^{2}$
$\displaystyle\qquad-\gamma\langle
D_{x}^{t,k},\nabla_{2}\phi_{v}(v^{t,k},x^{t,k})\rangle+\frac{\Lambda_{v}}{2}\gamma^{2}\|D_{x}^{t,k}\|^{2}\enspace.$

Tacking the expectation conditionally to the past iterates yields

$\displaystyle\mathbb{E}_{t,k}[\phi_{v}(v^{t,k+1},x^{t,k+1})]$
$\displaystyle\leq\left(1-\frac{\mu_{G}}{8}\rho\right)\phi_{v}(v^{t,k},x^{t,k})-\frac{\mu_{G}^{2}}{32L^{G}_{1}}\rho\mathbb{E}_{t,k}[\|\mathcal{G}^{t,k}_{v}\|^{2}]$
(40)
$\displaystyle\qquad+6\rho\left[\frac{L^{G}_{1}}{\mu_{G}^{2}}+\frac{\mu_{G}^{2}}{16L^{G}_{1}}\right]\mathbb{E}_{t,k}[\|D^{t,k}_{v}-D_{v}(\mathbf{u}^{t,k})\|^{2}]$
$\displaystyle\qquad+6\rho\left[\frac{L^{G}_{1}}{\mu_{G}^{2}}+\frac{\mu_{G}^{2}}{16L^{G}_{1}}\right]\mathbb{E}_{t,k}[\|D_{v}(\mathbf{u}^{t,k})-D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k})\|^{2}]$
$\displaystyle\qquad-\gamma\langle\mathbb{E}_{t,k}[D_{x}^{t,k}],\nabla_{2}\phi_{v}(v^{t,k},x^{t,k})\rangle+\frac{\Lambda_{v}}{2}\gamma^{2}\mathbb{E}_{t,k}[\|D_{x}^{t,k}\|^{2}]\enspace.$

From Young inequality, we have for any $c>0$

$\displaystyle\langle\mathbb{E}_{t,k}[D_{x}^{t,k}],\nabla_{2}\phi_{v}(v^{t,k},x^{t,k})\rangle\leq
c^{-1}\|\mathbb{E}_{t,k}[D_{x}^{t,k}]\|^{2}+c\|\nabla_{2}\phi_{v}(v^{t,k},x^{t,k})\|^{2}\enspace.$
(41)

Moreover, using the Lipschitz continuity of $z^{*}$, of $\nabla^{2}_{11}G$ and
$\nabla_{F}$ and the fact that $v$ and $v^{*}$ are bounded, we have

$\displaystyle\|\nabla_{2}\phi_{v}(v,x)\|$
$\displaystyle\leq\|\mathrm{d}z*(x)\|\left[\left\|\frac{1}{2}(\nabla^{3}_{111}G(z^{*}(x),x)|v)v-\frac{1}{2}(\nabla^{3}_{111}G(z^{*}(x),x)|v^{*}(x))v^{*}(x)\right\|\right.$
$\displaystyle\qquad\left.+\|\nabla_{11}^{2}F(z^{*}(x),x)v-\nabla_{11}^{2}F(z^{*}(x),x)v^{*}(x)\|\right]$
$\displaystyle\qquad+\|\frac{1}{2}(\nabla^{3}_{211}G(z^{*}(x),x)|v)v-\frac{1}{2}(\nabla^{3}_{211}G(z^{*}(x),x)|v^{*}(x))v^{*}(x)\|$
$\displaystyle\qquad+\|\nabla_{21}^{2}F(z^{*}(x),x)v-\nabla_{21}^{2}F(z^{*}(x),x)v^{*}(x)\|$
$\displaystyle\leq
L_{*}\left[\left\|\frac{1}{2}(\nabla^{3}_{111}G(z^{*}(x),x)|v-v^{*}(x))v-\frac{1}{2}(\nabla^{3}_{111}G(z^{*}(x),x)|v^{*}(x))(v-v^{*}(x))\right\|\right.$
$\displaystyle\qquad\left.+L^{F}_{2}\|v-v^{*}(x)\|\right]$
$\displaystyle\qquad+\left\|\frac{1}{2}(\nabla^{3}_{211}G(z^{*}(x),x)|v-v^{*}(x))v-\frac{1}{2}(\nabla^{3}_{211}G(z^{*}(x),x)|v^{*}(x))(v-v^{*}(x))\right\|$
$\displaystyle\qquad+L^{F}_{2}\|v-v^{*}(x)\|$ $\displaystyle\leq
L_{*}\left[\left\|\frac{1}{2}(\nabla^{3}_{111}G(z^{*}(x),x)|v-v^{*}(x))v\right\|\right.$
$\displaystyle\qquad\left.+\left\|\frac{1}{2}(\nabla^{3}_{111}G(z^{*}(x),x)|v^{*}(x))(v-v^{*}(x))\right\|+L^{F}_{2}\|v-v^{*}(x)\|\right]$
$\displaystyle\qquad+\left\|\frac{1}{2}(\nabla^{3}_{211}G(z^{*}(x),x)|v-v^{*}(x))v\right\|$
$\displaystyle\qquad+\left\|\frac{1}{2}(\nabla^{3}_{211}G(z^{*}(x),x)|v^{*}(x))(v-v^{*}(x))\right\|+L^{F}_{2}\|v-v^{*}(x)\|$
$\displaystyle\leq
L_{*}\left[\frac{L^{G}_{2}}{2}(\|v\|+\|v^{*}(x)\|)\|v-v^{*}(x))\|+L^{F}_{2}\|v-v^{*}(x)\|\right]$
$\displaystyle\qquad+\frac{L^{G}_{2}}{2}(\|v\|+\|v^{*}(x)\|)\|v-v^{*}(x))\|+L^{F}_{2}\|v-v^{*}(x)\|$
$\displaystyle\leq(L_{*}+1)\left[L^{G}_{2}R+L^{F}_{2}\right]\|v-v^{*}(x)\|\enspace.$

On the other hand, we have by strong convexity

$\displaystyle\|v-v^{*}(x)\|^{2}\leq\frac{2}{\mu_{G}}\phi_{v}(v,x)\enspace.$

As a consequence, we have

$\displaystyle\|\nabla_{2}\phi_{v}(v^{t,k},x^{t,k})\|^{2}\leq
L^{\prime\prime}\phi_{v}(v^{t,k},x^{t,k})$ (42)

with
$L^{\prime\prime}=\frac{2(L_{*}+1)^{2}\left[L^{G}_{2}R+L^{F}_{2}\right]^{2}}{\mu_{G}}$.

Plugging Inequalities (41) and (42) into (40) yields

$\displaystyle\mathbb{E}_{t,k}[\phi_{v}(v^{t,k+1},x^{t,k+1})]$
$\displaystyle\leq\left(1-\frac{\mu_{G}}{8}\rho+cL^{\prime\prime}\gamma\right)\phi_{v}(v^{t,k},x^{t,k})-\frac{\mu_{G}^{2}}{32L^{G}_{1}}\rho\mathbb{E}_{t,k}[\|\mathcal{G}^{t,k}_{v}\|^{2}]$
$\displaystyle\qquad+6\rho\left[\frac{L^{G}_{1}}{\mu_{G}^{2}}+\frac{\mu_{G}^{2}}{16L^{G}_{1}}\right]\mathbb{E}_{t,k}[\|D^{t,k}_{v}-D_{v}(\mathbf{u}^{t,k})\|^{2}]$
$\displaystyle\qquad+6\rho\left[\frac{L^{G}_{1}}{\mu_{G}^{2}}+\frac{\mu_{G}^{2}}{16L^{G}_{1}}\right]\mathbb{E}_{t,k}[\|D_{v}(\mathbf{u}^{t,k})-D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k})\|^{2}]$
$\displaystyle\qquad+\frac{\gamma}{c}\|\mathbb{E}_{t,k}[D^{t,k}_{x}]\|^{2}+\frac{\Lambda_{v}}{2}\gamma^{2}\mathbb{E}_{t,k}[\|D_{x}^{t,k}\|^{2}]\enspace.$

The Lipschitz continuity of $\nabla^{2}_{11}G$ and $\nabla_{1}F$ and the
boundedness of $v$ give us

$\displaystyle\|D_{v}(\mathbf{u}^{t,k})-D_{v}(z^{*}(x^{t,k}),v^{t,k},x^{t,k})\|^{2}$
$\displaystyle\leq\left(\|\nabla^{2}_{11}G(z^{t,k},x^{t,k})v^{t,k}-\nabla^{2}_{11}G(z^{*}(x^{t,k}),x^{t,k})v^{t,k}\|\right.$
$\displaystyle\qquad\left.+\|\nabla_{1}F(z^{t,k},x^{t,k})-\nabla_{1}F(z^{*}(x^{t,k}),x^{t,k})\|\right)^{2}$
$\displaystyle\leq(L^{G}_{2}R+L^{F}_{1})^{2}\|z^{t,k}-z^{*}(x^{t,k})\|^{2}$
$\displaystyle\leq\frac{2(L^{G}_{2}R+L^{F}_{1})^{2}}{\mu_{G}}\phi_{z}(z^{t,k},x^{t,k})\enspace.$

As a consequence

$\displaystyle\mathbb{E}_{t,k}[\phi_{v}(v^{t,k+1},x^{t,k+1})]$
$\displaystyle\leq\left(1-\frac{\mu_{G}}{8}\rho+cL^{\prime\prime}\gamma\right)\phi_{v}(v^{t,k},x^{t,k})-\frac{\mu_{G}^{2}}{32L^{G}_{1}}\rho\mathbb{E}_{t,k}[\|\mathcal{G}^{t,k}_{v}\|^{2}]$
(43)
$\displaystyle\qquad+6\rho\left[\frac{L^{G}_{1}}{\mu_{G}^{2}}+\frac{\mu_{G}^{2}}{16L^{G}_{1}}\right]\mathbb{E}_{t,k}[\|D^{t,k}_{v}-D_{v}(\mathbf{u}^{t,k})\|^{2}]$
$\displaystyle\qquad+6\rho\left[\frac{L^{G}_{1}}{\mu_{G}^{2}}+\frac{\mu_{G}^{2}}{16L^{G}_{1}}\right]\frac{2(L^{G}_{2}R+L^{F}_{1})^{2}}{\mu_{G}}\phi_{z}(z^{t,k},x^{t,k})$
$\displaystyle\qquad+\frac{\gamma}{c}\|\mathbb{E}_{t,k}[D^{t,k}_{x}]\|^{2}+\frac{\Lambda_{v}}{2}\gamma^{2}\mathbb{E}_{t,k}[\|D_{x}^{t,k}\|^{2}]\enspace.$

From Lemma A.2, we have

$\displaystyle\mathbb{E}[\|D_{v}^{t,k}-D_{v}(\mathbf{u}^{t,k})\|^{2}]$
$\displaystyle\leq
4\rho^{2}\left((L^{G}_{2}R)^{2}+(L^{F}_{1})^{2}\right)\sum_{r=1}^{k}\mathbb{E}[\|D^{t,r-1}_{z}\|^{2}]+4\rho^{2}(L^{G}_{1})^{2}\sum_{r=1}^{k}\mathbb{E}[\|\mathcal{G}^{t,r-1}_{v}\|^{2}]$
$\displaystyle\qquad+4\gamma^{2}\left((L^{G}_{2}R)^{2}+(L^{F}_{1})^{2}\right)\sum_{r=1}^{k}\mathbb{E}[\|D_{x}^{t,r-1}\|^{2}]$

Taking the total expectation and plugging the previous inequality in Equation
43 yields

$\displaystyle\phi_{v}^{t,k+1}$
$\displaystyle\leq\left(1-\frac{\mu_{G}}{8}\rho+cL^{\prime\prime}\gamma\right)\phi_{v}^{t,k}-\frac{\mu_{G}^{2}}{32L^{G}_{1}}\rho\mathbb{E}_{t,k}[\|\mathcal{G}^{t,k}_{v}\|^{2}]$
$\displaystyle\qquad+24\rho^{3}\left((L^{G}_{2}R)^{2}+(L^{F}_{1})^{2}\right)\left(\frac{L^{G}_{1}}{\mu_{G}^{2}}+\frac{\mu_{G}^{2}}{16L^{G}_{1}}\right)\sum_{r=1}^{k}\mathbb{E}[\|D^{t,r-1}_{z}\|^{2}]$
$\displaystyle\qquad+24\rho^{3}(L^{G}_{1})^{2}\left(\frac{L^{G}_{1}}{\mu_{G}^{2}}+\frac{\mu_{G}^{2}}{16L^{G}_{1}}\right)\sum_{r=1}^{k}\mathbb{E}[\|\mathcal{G}^{t,r-1}_{v}\|^{2}]$
$\displaystyle\qquad+24\rho\gamma^{2}\left((L^{G}_{2}R)^{2}+(L^{F}_{1})^{2}\right)\left(\frac{L^{G}_{1}}{\mu_{G}^{2}}+\frac{\mu_{G}^{2}}{16L^{G}_{1}}\right)\sum_{r=1}^{k}\mathbb{E}[\|D_{x}^{t,r-1}\|^{2}]$
$\displaystyle\qquad+\left[\frac{L^{G}_{1}}{\mu_{G}^{2}}+\frac{\mu_{G}^{2}}{16L^{G}_{1}}\right]\frac{12(L^{G}_{2}R+L^{F}_{1})^{2}}{\mu_{G}}\rho\phi_{z}^{t,k}$
$\displaystyle\qquad+\frac{\gamma}{c}\mathbb{E}[\|[\mathbb{E}_{t,k}D_{x}^{t,k}]\|^{2}]+\frac{\Lambda_{v}}{2}\gamma^{2}\mathbb{E}[\|D_{x}^{t,k}\|^{2}]\enspace.$

Taking $c=\frac{\mu_{G}\rho}{16L^{\prime\prime}\gamma}$ yields

$\displaystyle\phi_{v}^{t,k+1}$
$\displaystyle\leq\left(1-\frac{\mu_{G}}{16}\rho\right)\phi_{v}^{t,k}-\frac{\mu_{G}^{2}}{32L^{G}_{1}}\rho\mathbb{E}_{t,k}[\|\mathcal{G}^{t,k}_{v}\|^{2}]$
$\displaystyle\qquad+24\rho^{3}\left((L^{G}_{2}R)^{2}+(L^{F}_{1})^{2}\right)\left(\frac{L^{G}_{1}}{\mu_{G}^{2}}+\frac{\mu_{G}^{2}}{16L^{G}_{1}}\right)\sum_{r=1}^{k}\mathbb{E}[\|D^{t,r-1}_{z}\|^{2}]$
$\displaystyle\qquad+24\rho^{3}(L^{G}_{1})^{2}\left(\frac{L^{G}_{1}}{\mu_{G}^{2}}+\frac{\mu_{G}^{2}}{16L^{G}_{1}}\right)\sum_{r=1}^{k}\mathbb{E}[\|\mathcal{G}^{t,r-1}_{v}\|^{2}]$
$\displaystyle\qquad+24\rho\gamma^{2}\left((L^{G}_{2}R)^{2}+(L^{F}_{1})^{2}\right)\left(\frac{L^{G}_{1}}{\mu_{G}^{2}}+\frac{\mu_{G}^{2}}{16L^{G}_{1}}\right)\sum_{r=1}^{k}\mathbb{E}[\|D_{x}^{t,r-1}\|^{2}]$
$\displaystyle\qquad+\left[\frac{L^{G}_{1}}{\mu_{G}^{2}}+\frac{\mu_{G}^{2}}{16L^{G}_{1}}\right]\frac{12(L^{G}_{2}R+L^{F}_{1})^{2}}{\mu_{G}}\rho\phi_{z}^{t,k}$
$\displaystyle\qquad+\frac{16L^{\prime\prime}}{\mu_{G}}\frac{\gamma^{2}}{\rho}\mathbb{E}[\|[\mathbb{E}_{t,k}D_{x}^{t,k}]\|^{2}]+\frac{\Lambda_{v}}{2}\gamma^{2}\mathbb{E}[\|D_{x}^{t,k}\|^{2}]\enspace.$

Combining Equation 21 and Lemma A.2 yields

$\displaystyle\phi_{v}^{t,k+1}$
$\displaystyle\leq\left(1-\frac{\mu_{G}}{16}\rho\right)\phi_{v}^{t,k}-\frac{\mu_{G}^{2}}{32L^{G}_{1}}\rho\mathbb{E}_{t,k}[\|\mathcal{G}^{t,k}_{v}\|^{2}]$
$\displaystyle\qquad+8\rho\left((L^{G}_{2}R)^{2}+(L^{F}_{1})^{2}\right)\left[3\left(\frac{L^{G}_{1}}{\mu_{G}^{2}}+\frac{\mu_{G}^{2}}{16L^{G}_{1}}\right)\rho^{2}+\frac{8L^{\prime\prime}}{\mu_{G}}\gamma^{2}\right]\sum_{r=1}^{k}\mathbb{E}[\|D^{t,r-1}_{z}\|^{2}]$
$\displaystyle\qquad+8\rho(L^{G}_{1})^{2}\left[3\left(\frac{L^{G}_{1}}{\mu_{G}^{2}}+\frac{\mu_{G}^{2}}{16L^{G}_{1}}\right)\rho^{2}+\frac{8L^{\prime\prime}}{\mu_{G}}\gamma^{2}\right]\sum_{r=1}^{k}\mathbb{E}[\|\mathcal{G}^{t,r-1}_{v}\|^{2}]$
$\displaystyle\qquad+8\gamma^{2}\left((L^{G}_{2}R)^{2}+(L^{F}_{1})^{2}\right)\left[3\left(\frac{L^{G}_{1}}{\mu_{G}^{2}}+\frac{\mu_{G}^{2}}{16L^{G}_{1}}\right)\gamma+\frac{8L^{\prime\prime}}{\mu_{G}}\frac{\gamma^{2}}{\rho}\right]\sum_{r=1}^{k}\mathbb{E}[\|D_{x}^{t,r-1}\|^{2}]$
$\displaystyle\qquad+\left[\frac{L^{G}_{1}}{\mu_{G}^{2}}+\frac{\mu_{G}^{2}}{16L^{G}_{1}}\right]\frac{12(L^{G}_{2}R+L^{F}_{1})^{2}}{\mu_{G}}\rho\phi_{z}^{t,k}$
$\displaystyle\qquad+\frac{16L^{\prime\prime}}{\mu_{G}}\frac{\gamma^{2}}{\rho}\mathbb{E}[\|D_{x}(\mathbf{u}^{tk,})\|^{2}]+\frac{\Lambda_{v}}{2}\gamma^{2}\mathbb{E}[\|D_{x}^{t,k}\|^{2}]\enspace.$

By assumption, $\gamma\leq C_{v}\rho$ with
$C_{v}=\sqrt{\frac{\mu_{G}}{8L^{\prime\prime}}\left(\frac{L^{G}_{1}}{\mu_{G}^{2}}+\frac{\mu_{G}^{2}}{16L^{G}_{1}}\right)}$,
therefore

$\displaystyle\phi_{v}^{t,k+1}$
$\displaystyle\leq\left(1-\frac{\mu_{G}}{16}\rho\right)\phi_{v}^{t,k}-\frac{\mu_{G}^{2}}{32L^{G}_{1}}\rho\mathbb{E}_{t,k}[\|\mathcal{G}^{t,k}_{v}\|^{2}]$
$\displaystyle\qquad+32\rho^{3}\left((L^{G}_{2}R)^{2}+(L^{F}_{1})^{2}\right)\left(\frac{L^{G}_{1}}{\mu_{G}^{2}}+\frac{\mu_{G}^{2}}{16L^{G}_{1}}\right)\sum_{r=1}^{k}\mathbb{E}[\|D^{t,r-1}_{z}\|^{2}]$
$\displaystyle\qquad+32\rho^{3}(L^{G}_{1})^{2}\left(\frac{L^{G}_{1}}{\mu_{G}^{2}}+\frac{\mu_{G}^{2}}{16L^{G}_{1}}\right)\sum_{r=1}^{k}\mathbb{E}[\|\mathcal{G}^{t,r-1}_{v}\|^{2}]$
$\displaystyle\qquad+32\gamma^{2}\rho\left((L^{G}_{2}R)^{2}+(L^{F}_{1})^{2}\right)\left(\frac{L^{G}_{1}}{\mu_{G}^{2}}+\frac{\mu_{G}^{2}}{16L^{G}_{1}}\right)\sum_{r=1}^{k}\mathbb{E}[\|D_{x}^{t,r-1}\|^{2}]$
$\displaystyle\qquad+\left[\frac{L^{G}_{1}}{\mu_{G}^{2}}+\frac{\mu_{G}^{2}}{16L^{G}_{1}}\right]\frac{12(L^{G}_{2}R+L^{F}_{1})^{2}}{\mu_{G}}\rho\phi_{z}^{t,k}$
$\displaystyle\qquad+\frac{16L^{\prime\prime}}{\mu_{G}}\frac{\gamma^{2}}{\rho}\mathbb{E}[\|D_{x}(\mathbf{u}^{tk,})\|^{2}]+\frac{\Lambda_{v}}{2}\gamma^{2}\mathbb{E}[\|D_{x}^{t,k}\|^{2}]\enspace.$

We get finally

$\displaystyle\phi_{v}^{t,k+1}$
$\displaystyle\leq\left(1-\frac{\rho\mu_{G}}{16}\right)\phi_{v}^{t,k}-\tilde{\beta}_{vv}\rho
V_{v}^{t}+\rho^{3}\beta_{vz}\mathcal{V}_{z}^{t,k}+2\rho^{3}\beta_{vv}\mathcal{V}_{v}^{t,k}+\gamma^{2}\rho\beta_{vx}\mathcal{V}_{x}^{t,k}$
$\displaystyle\qquad+\rho\alpha_{vz}\phi_{z}^{t,k}+\frac{\Lambda_{v}}{2}\gamma^{2}\mathbb{E}[\|D_{x}^{t,k}\|^{2}]+\frac{\gamma^{2}}{\rho}\overline{\beta}_{vx}\mathbb{E}[\|D_{x}(\mathbf{u}^{t,k})\|^{2}]$

with
$\beta_{vz}=\beta_{vx}=32\left((L^{G}_{2}R)^{2}+(L^{F}_{1})^{2}\right)\left(\frac{L^{G}_{1}}{\mu_{G}^{2}}+\frac{\mu_{G}^{2}}{16L^{G}_{1}}\right)$,
$\beta_{vv}=(L^{G}_{1})^{2}\left(\frac{L^{G}_{1}}{\mu_{G}^{2}}+\frac{\mu_{G}^{2}}{16L^{G}_{1}}\right)$,
$\overline{\beta}_{vx}~{}=~{}\frac{16L^{\prime\prime}}{\mu_{G}}$,
$\tilde{\beta}_{vv}=\frac{\mu_{G}^{2}}{32L^{G}_{1}}$ and
$\alpha_{vz}=\left[\frac{L^{G}_{1}}{\mu_{G}^{2}}+\frac{\mu_{G}^{2}}{16L^{G}_{1}}\right]\frac{12(L^{G}_{2}R+L^{F}_{1})^{2}}{\mu_{G}}$.
∎

###  A.7 Proof of Lemma 3.5

###### Proof.

The smoothness of $h$ (Proposition 2.3) gives us

$\displaystyle h(x^{t,k+1})$ $\displaystyle\leq h(x^{t,k})-\gamma\langle\nabla
h(x^{t,k}),D_{x}^{t,k}\rangle+\gamma^{2}\frac{L^{h}}{2}\|D^{t,k}_{x}\|^{2}\enspace.$

Then, we use the identity $\langle
a,b\rangle=\frac{1}{2}(\|a\|^{2}+\|b\|^{2}-\|a-b\|^{2})$ to get

$\displaystyle h(x^{t,k+1})$ $\displaystyle\leq
h(x^{t,k})-\frac{\gamma}{2}\|\nabla
h(x^{t,k})\|^{2}-\frac{\gamma}{2}\|D^{t,k}_{x}\|^{2}+\frac{\gamma}{2}\|\nabla
h(x^{t,k})-D^{t,k}_{x}\|^{2}+\gamma^{2}\frac{L^{h}}{2}\|D^{t,k}_{x}\|^{2}$
$\displaystyle\leq h(x^{t,k})-\frac{\gamma}{2}\|\nabla
h(x^{t,k})\|^{2}-\frac{\gamma}{2}\|D^{t,k}_{x}\|^{2}+\gamma\|\nabla
h(x^{t,k})-D_{x}(\mathbf{u}^{t,k})\|^{2}$
$\displaystyle\qquad+\gamma\|D_{x}(\mathbf{u}^{t,k})-D^{t,k}_{x}\|^{2}+\gamma^{2}\frac{L^{h}}{2}\|D^{t,k}_{x}\|^{2}\enspace.$

Then taking the expectation gives and using Proposition 2.5 yields

$\displaystyle h^{t,k+1}$ $\displaystyle\leq
h^{t,k}-\frac{\gamma}{2}g^{t,k}+\gamma\mathbb{E}[\|\nabla
h(x^{t,k})-D_{x}(\mathbf{u}^{t,k})\|^{2}]$
$\displaystyle\qquad+\gamma\mathbb{E}[\|D_{x}(\mathbf{u}^{t,k})-D^{t,k}_{x}\|^{2}]-\frac{\gamma}{2}\left(1-L^{h}\gamma\right)\mathbb{E}[\|D^{t,k}_{x}\|^{2}]$
$\displaystyle\leq h^{t,k}-\frac{\gamma}{2}g^{t,k}+\gamma
L_{x}^{2}(\mathbb{E}[\|z^{t,k}-z^{*}(x^{t,k})\|^{2}]+\mathbb{E}[\|v^{t,k}-v^{*}(x^{t,k})\|^{2}])$
$\displaystyle\qquad+\gamma\mathbb{E}[\|D_{x}(\mathbf{u}^{t,k})-D^{t,k}_{x}\|^{2}]-\frac{\gamma}{2}\left(1-L^{h}\gamma\right)\mathbb{E}[\|D^{t,k}_{x}\|^{2}]\enspace.$

The $\mu_{G}$-strong convexity of $G(\,.\,,x)$ ensures that
$\|z-z^{*}(x)\|^{2}\leq\frac{2}{\mu_{G}}\phi_{z}(z,x)$ and
$\|v~{}-~{}v^{*}(x)\|^{2}~{}\leq~{}\frac{2}{\mu_{G}}\phi_{v}(v,x)$. As a
consequence

$\displaystyle h^{t,k+1}$ $\displaystyle\leq
h^{t,k}-\frac{\gamma}{2}g^{t,k}+\gamma\frac{2L_{x}^{2}}{\mu_{G}}(\phi_{z}^{t,k}+\phi_{v}^{t,k})+\gamma\mathbb{E}[\|D_{x}(z^{t,k},v^{t,k},x^{t,k})-D^{t,k}_{x}\|^{2}]$
$\displaystyle\qquad-\frac{\gamma}{2}\left(1-L^{h}\gamma\right)\mathbb{E}[\|D^{t,k}_{x}\|^{2}]\enspace.$

From Lemma A.2, we have

$\displaystyle\mathbb{E}[\|D_{x}^{t,k}-D_{x}(\mathbf{u}^{t,k})\|^{2}]$
$\displaystyle\leq
4\rho^{2}\left((L^{G}_{2}R)^{2}+(L^{F}_{1})^{2}\right)\sum_{r=1}^{k}\mathbb{E}[\|D^{t,r-1}_{z}\|^{2}]+4\rho^{2}(L^{G}_{1})^{2}\sum_{r=1}^{k}\mathbb{E}[\|\mathcal{G}^{t,r-1}_{v}\|^{2}]$
$\displaystyle\qquad+4\gamma^{2}\left((L^{G}_{2}R)^{2}+(L^{F}_{1})^{2}\right)\sum_{r=1}^{k}\mathbb{E}[\|D_{x}^{t,r-1}\|^{2}]\enspace.$

As a consequence

$\displaystyle h^{t,k+1}$ $\displaystyle\leq
h^{t,k}-\frac{\gamma}{2}g^{t,k}+\gamma\frac{2L_{x}^{2}}{\mu_{G}}(\phi_{z}^{t,k}+\phi_{v}^{t,k})$
$\displaystyle\qquad+4\gamma\rho^{2}\left((L^{G}_{2}R)^{2}+2(L^{F}_{1})^{2}\right)\sum_{r=1}^{k}\mathbb{E}[\|D^{t,r-1}_{z}\|^{2}]+4\gamma\rho^{2}(L_{1}^{G})^{2}\sum_{r=1}^{k}\mathbb{E}[\|\mathcal{G}^{t,r-1}_{v}\|^{2}]$
$\displaystyle\qquad+4\gamma^{3}\left((L^{G}_{2}R)^{2}+2(L^{F}_{1})^{2}\right)\sum_{r=1}^{k}\mathbb{E}[\|D^{t,r-1}_{x}\|^{2}]-\frac{\gamma}{2}\left(1-L^{h}\gamma\right)\mathbb{E}[\|D^{t,k}_{x}\|^{2}]$
$\displaystyle\leq
h^{t,k}-\frac{\gamma}{2}g^{t,k}+\gamma\frac{2L_{x}^{2}}{\mu_{G}}(\phi_{z}^{t,k}+\phi_{v}^{t,k})+\gamma\rho^{2}\beta_{hz}\mathcal{V}_{z}^{t,k}+\gamma\rho^{2}\beta_{hv}\mathcal{V}_{v}^{t,k}$
$\displaystyle\qquad+\gamma^{3}\beta_{hx}\mathcal{V}_{x}^{t,k}-\frac{\gamma}{2}\left(1-L^{h}\gamma\right)\mathbb{E}[\|D^{t,k}_{x}\|^{2}]$

with $\beta_{hz}=4\left((L^{G}_{2}R)^{2}+2(L^{F}_{1})^{2}\right)$,
$\beta_{hv}=4(L_{1}^{G})^{2}$ and
$\beta_{hx}=4\left((L^{G}_{2}R)^{2}+2(L^{F}_{1})^{2}\right)$. ∎

###  A.8 Proof of 1 and Corollary 3.6

The constants involved in 1 are

$\displaystyle\psi_{z}=\frac{1}{16\overline{\beta}_{zx}},\quad\psi_{v}=\min\left[\frac{1}{16\overline{\beta}_{vx}},\frac{\alpha_{zv}\mu_{G}}{12}\psi_{z}\right]$
$\displaystyle\overline{\rho}=\min\left[\sqrt{\frac{\psi_{z}}{12q(\psi_{z}\beta_{zz}+\psi_{v}\beta_{zv})}},\sqrt{\frac{1}{6\Lambda_{z}}},\sqrt{\frac{1}{12q\beta_{vv}}},sqrt{\frac{\tilde{\beta}_{vv}}{3\Lambda_{v}}},B_{v}\right]\enspace,$
$\displaystyle\xi=\min\left[C_{z},C_{v},1,\frac{\psi_{v}\mu_{G}^{2}}{16L_{x}^{2}},\sqrt{\frac{\mu_{G}}{8\overline{\beta}_{vx}}},\frac{\psi_{z}\mu_{G}^{2}}{24L_{x}^{2}},\sqrt{\frac{\mu_{G}}{12\overline{\beta}_{zx}}}\right]\enspace,$
$\displaystyle\overline{\gamma}=\min\left[\sqrt{\frac{1}{12q(\psi_{z}\beta_{zx}+\psi_{v}\beta_{vx})}},\sqrt{\frac{1}{12q\beta_{hx}}},\frac{1}{6(L^{h}+\psi_{z}\Lambda_{z}+\psi_{v}\Lambda_{v})},\sqrt{\frac{\psi_{v}\tilde{\beta}_{vv}}{6q(\beta_{hv}+\psi_{z}\beta_{vz})}},\sqrt{\frac{\psi_{z}}{12q\beta_{hz}}}\right]\enspace.$

###### Proof.

The proof is a classical Lyapunov analysis. Consider the following Lyapunov
function
$\mathcal{L}^{t,k}=h^{t,k}+\psi_{z}\phi^{t,k}_{z}+\psi_{v}\phi^{t,k}_{v}$ for
some positive constants $\psi_{z}$ and $\psi_{v}$. We use use Lemmas 3.3 to
3.5 to upper bound $\mathcal{L}^{t,k}-\mathcal{L}^{t,k+1}$.

We have

$\displaystyle\mathcal{L}^{t,k+1}-\mathcal{L}^{t,k}$
$\displaystyle\leq-\frac{\gamma}{2}g^{t,k}+(\psi_{z}\overline{\beta}_{zx}+\psi_{v}\overline{\beta}_{vx})\frac{\gamma^{2}}{\rho}\mathbb{E}[\|D_{x}(\mathbf{u}^{t,k})\|^{2}]$
(44)
$\displaystyle\qquad+\left(\frac{2L_{x}^{2}}{\mu_{G}}\gamma-\psi_{z}\frac{\mu_{G}}{2}\rho+\psi_{v}\alpha_{zv}\rho\right)\phi_{z}^{t,k}+\left(\frac{2L_{x}^{2}}{\mu_{G}}\gamma-\psi_{v}\frac{\mu_{G}}{16}\rho\right)\phi_{v}^{t,k}$
$\displaystyle\qquad+\left(\psi_{z}\frac{\Lambda_{z}}{2}\rho^{2}-\psi_{z}\frac{1}{2}\rho\right)V_{z}^{t,k}-\psi_{v}\tilde{\beta}_{vv}\rho
V_{v}^{t,k}$
$\displaystyle\qquad+\left(\frac{L^{h}}{2}\gamma^{2}+\psi_{z}\frac{\Lambda_{z}}{2}\gamma^{2}+\psi_{v}\frac{\Lambda_{v}}{2}\gamma^{2}-\frac{\gamma}{2}\right)V_{x}^{t,k}$
$\displaystyle\qquad+\left(\beta_{hz}\rho\gamma^{2}+\psi_{z}\beta_{zz}\rho^{3}+\psi_{v}\beta_{zv}\rho^{3}\right)\mathcal{V}_{z}^{t,k}$
$\displaystyle\qquad+\left(\beta_{hv}\rho\gamma^{2}+\psi_{z}\beta_{vz}\gamma^{2}\rho+\psi_{v}\beta_{vv}\rho^{3}\right)\mathcal{V}_{v}^{t,k}$
$\displaystyle\qquad+\left(\beta_{hx}\gamma^{3}+\psi_{z}\beta_{zx}\gamma^{2}\rho+\psi_{v}\beta_{vx}\rho^{3}\right)\mathcal{V}_{x}^{t,k}\enspace.$

We bound $\mathbb{E}[\|D_{x}(\mathbf{u}^{t,k})\|^{2}]$ crudely by using
Proposition 2.5

$\displaystyle\mathbb{E}[\|D_{x}(\mathbf{u}^{t,k})\|^{2}]$ $\displaystyle\leq
2\mathbb{E}[\|\nabla
h(x^{t,k})\|^{2}]+2\mathbb{E}[\|D_{x}(\mathbf{u}^{t,k})-\nabla
h(x^{t,k})\|^{2}]$ $\displaystyle\leq
2g^{t,k}+2(\mathbb{E}[\|z^{t,k}-z^{*}(x^{t,k})\|^{2}]+\mathbb{E}[\|v^{t,k}-v^{*}(x^{t,k})\|^{2}])$
$\displaystyle\leq
2g^{t,k}+\frac{4}{\mu_{G}}(\phi_{z}^{t,k}+\phi_{v}^{t,k})\enspace.$

Summing in (44) for $k=0,\dots,q-1$ yields

$\displaystyle\mathcal{L}^{t,q}-\mathcal{L}^{t,0}$
$\displaystyle\leq-\left(\frac{\gamma}{2}-2\psi_{z}\overline{\beta}_{zx}\frac{\gamma^{2}}{\rho}-2\psi_{v}\overline{\beta}_{vx}\frac{\gamma^{2}}{\rho}\right)\sum_{k=0}^{q-1}g^{t,k}$
(45)
$\displaystyle\qquad+\left(\frac{2L_{x}^{2}}{\mu_{G}}\gamma-\psi_{z}\frac{\mu_{G}}{2}\rho+\psi_{v}\alpha_{zv}\rho+\psi_{z}\overline{\beta}_{zx}\frac{\gamma^{2}}{\rho}\right)\sum_{k=0}^{q-1}\phi_{z}^{t,k}$
$\displaystyle\qquad+\left(\frac{2L_{x}^{2}}{\mu_{G}}\gamma-\psi_{v}\frac{\mu_{G}}{16}\rho+\psi_{v}\overline{\beta}_{vx}\frac{\gamma^{2}}{\rho}\right)\sum_{k=0}^{q-1}\phi_{v}^{t,k}-\psi_{z}\tilde{\beta}_{vv}\rho\sum_{k=0}^{q-1}V_{z}^{t,k}$
$\displaystyle\qquad+\left(\psi_{v}\frac{\Lambda_{v}}{2}\rho^{2}-\psi_{v}\frac{1}{2}\rho\right)\sum_{k=0}^{q-1}V_{v}^{t,k}+\left(\frac{L^{h}}{2}\gamma^{2}+\psi_{z}\frac{\Lambda_{z}}{2}\gamma^{2}+\psi_{v}\frac{\Lambda_{v}}{2}\gamma^{2}-\frac{\gamma}{2}\right)\sum_{k=0}^{q-1}V_{x}^{t,k}$
$\displaystyle\qquad+\left(\beta_{hz}\rho\gamma^{2}+\psi_{z}\beta_{zz}\rho^{3}+\psi_{v}\beta_{zv}\rho^{3}\right)\sum_{k=0}^{q-1}\mathcal{V}_{z}^{t,k}$
$\displaystyle\qquad+\left(\beta_{hv}\rho\gamma^{2}+\psi_{z}\beta_{vz}\gamma^{2}\rho+\psi_{v}\beta_{vv}\rho^{3}\right)\sum_{k=0}^{q-1}\mathcal{V}_{v}^{t,k}$
$\displaystyle\qquad+\left(\beta_{hx}\gamma^{3}+\psi_{z}\beta_{zx}\gamma^{2}\rho+\psi_{v}\beta_{vx}\rho^{3}\right)\sum_{k=0}^{q-1}\mathcal{V}_{x}^{t,k}\enspace.$

Since we have

$\displaystyle\sum_{k=0}^{q-1}\mathcal{V}_{\bullet}^{t,k}$
$\displaystyle=\sum_{k=0}^{q-1}\sum_{r=1}^{k}\mathbb{E}[\|D_{\bullet}^{t,r-1}\|^{2}]=\sum_{r=1}^{q-1}\sum_{k=r}^{q-1}\mathbb{E}[\|D_{\bullet}^{t,r-1}\|^{2}]$
$\displaystyle\qquad=\sum_{r=1}^{q-1}(q-r)\mathbb{E}[\|D_{\bullet}^{t,r-1}\|^{2}]\leq
q\sum_{k=1}^{q-1}\mathbb{E}[\|D_{\bullet}^{t,k-1}\|^{2}]\enspace$

we get

$\displaystyle\mathcal{L}^{t,q}-\mathcal{L}^{t,0}$
$\displaystyle\leq-\left(\frac{\gamma}{2}-2\psi_{z}\overline{\beta}_{zx}\frac{\gamma^{2}}{\rho}-2\psi_{v}\overline{\beta}_{vx}\frac{\gamma^{2}}{\rho}\right)\sum_{k=0}^{q-1}g^{t,k}$
(46)
$\displaystyle\qquad+\left(\frac{2L_{x}^{2}}{\mu_{G}}\gamma-\psi_{z}\frac{\mu_{G}}{2}\rho+\psi_{v}\alpha_{zv}\rho+\psi_{z}\overline{\beta}_{zx}\frac{\gamma^{2}}{\rho}\right)\sum_{k=0}^{q-1}\phi_{z}^{t,k}$
$\displaystyle\qquad+\left(\frac{2L_{x}^{2}}{\mu_{G}}\gamma-\psi_{v}\frac{\mu_{G}}{2}\rho+\psi_{v}\overline{\beta}_{vx}\frac{\gamma^{2}}{\rho}\right)\sum_{k=0}^{q-1}\phi_{v}^{t,k}$
$\displaystyle\qquad+\left(\psi_{z}\frac{\Lambda_{z}}{2}\rho^{2}-\psi_{z}\frac{1}{2}\rho+q\left(\beta_{hz}\rho\gamma^{2}+\psi_{z}\beta_{zz}\rho^{3}+\psi_{v}\beta_{zv}\rho^{3}\right)\right)\sum_{k=0}^{q-1}V_{z}^{t,k}$
$\displaystyle\qquad+\left(\psi_{v}\frac{\Lambda_{v}}{2}\rho^{2}-\psi_{v}\tilde{\beta}_{vv}\rho+q\left(\beta_{hv}\rho\gamma^{2}+\psi_{z}\beta_{vz}\gamma^{2}\rho+\psi_{v}\beta_{vv}\rho^{3}\right)\right)\sum_{k=0}^{q-1}V_{v}^{t,k}$
$\displaystyle\qquad+\left(\frac{L^{h}}{2}\gamma^{2}+\psi_{z}\frac{\Lambda_{z}}{2}\gamma^{2}+\psi_{v}\frac{\Lambda_{v}}{2}\gamma^{2}-\frac{\gamma}{2}+q\left(\beta_{hx}\gamma^{3}+\psi_{z}\beta_{zx}\gamma^{2}\rho+\psi_{v}\beta_{vx}\rho\gamma^{2}\right)\right)\sum_{k=0}^{q-1}V_{x}^{t,k}\enspace.$

Since
$\rho\leq\overline{\rho}\leq\min\left[\sqrt{\frac{\psi_{z}}{12q(\psi_{z}\beta_{zz}+\psi_{v}\beta_{zv})}},\sqrt{\frac{1}{6\Lambda_{z}}}\right]$
and $\gamma\leq\overline{\gamma}\leq\sqrt{\frac{\psi_{z}}{12q\beta_{hz}}}$, we
have

$\displaystyle\psi_{z}\frac{\Lambda_{z}}{2}\rho^{2}-\psi_{z}\frac{1}{2}\rho+q\left(\beta_{hz}\rho\gamma^{2}+\psi_{z}\beta_{zz}\rho^{3}+\psi_{v}\beta_{zv}\rho^{3}\right)<0\enspace.$
(47)

Moreover, the conditions
$\rho\leq\overline{\rho}\leq\min\left[\sqrt{\frac{\tilde{\beta}_{vv}}{6q\beta_{vv}}},\sqrt{\frac{\tilde{\beta}_{vv}}{3\Lambda_{v}}}\right]$
and
$\gamma\leq\overline{\gamma}\leq\sqrt{\frac{\psi_{v}\tilde{\beta}_{vv}}{6q(\beta_{hv}+\psi_{z}\beta_{vz})}}$,
ensure that

$\displaystyle\psi_{v}\frac{\Lambda_{v}}{2}\rho^{2}-\psi_{v}\tilde{\beta}_{vv}\rho+q\left(\beta_{hv}\rho\gamma^{2}+\psi_{z}\beta_{vz}\gamma^{2}\rho+\psi_{v}\beta_{vv}\rho^{3}\right)<0\enspace.$
(48)

The conditions
$\rho\leq\overline{\rho}\leq\sqrt{\frac{1}{12q(\psi_{z}\beta_{zx}+\psi_{v}\beta_{vx})}}$
and
$\gamma\leq\overline{\gamma}\leq\min\left[\sqrt{\frac{1}{12q(\psi_{z}\beta_{zx}+\psi_{v}\beta_{vx})}},\sqrt{\frac{1}{12q\beta_{hx}}},\frac{1}{6(L^{h}+\psi_{z}\Lambda_{z}+\psi_{v}\Lambda_{v})}\right]$
yields

$\displaystyle\frac{L^{h}}{2}\gamma^{2}+\psi_{z}\frac{\Lambda_{z}}{2}\gamma^{2}+\psi_{v}\frac{\Lambda_{v}}{2}\gamma^{2}-\frac{\gamma}{2}+q\left(\beta_{hx}\gamma^{3}+\psi_{z}\beta_{zx}\gamma^{2}\rho+\psi_{v}\beta_{vx}\rho\gamma^{2}\right)<0\enspace.$
(49)

The condition
$\gamma\leq\xi\rho\leq\min\left[\frac{\psi_{v}\mu_{G}^{2}}{16L_{x}^{2}},\sqrt{\frac{\mu_{G}}{8\overline{\beta}_{vx}}}\right]\rho$
ensures

$\displaystyle\frac{2L_{x}^{2}}{\mu_{G}}\gamma-\psi_{v}\frac{\mu_{G}}{2}\rho+\psi_{v}\overline{\beta}_{vx}\frac{\gamma^{2}}{\rho}\leq
0$ (50)

By definition, we have $\psi_{v}\leq\frac{\alpha_{zv}\mu_{G}}{12}\psi_{z}$ and
by assumptions
$\gamma\leq\xi\rho\leq\min\left[\frac{\psi_{z}\mu_{G}^{2}}{24L_{x}^{2}},\sqrt{\frac{\mu_{G}}{12\overline{\beta}_{zx}}}\right]\rho$.
As a consequence

$\displaystyle\frac{2L_{x}^{2}}{\mu_{G}}\gamma-\psi_{z}\frac{\mu_{G}}{2}\rho+\psi_{v}\alpha_{zv}\rho+\psi_{z}\overline{\beta}_{zx}\frac{\gamma^{2}}{\rho}<0\enspace.$
(51)

Plugging Inequalities (47) to (51) into Section A.8 gives

$\displaystyle\mathcal{L}^{t,q}-\mathcal{L}^{t,0}$
$\displaystyle\leq-\left(\frac{\gamma}{2}-2\psi_{z}\overline{\beta}_{zx}\frac{\gamma^{2}}{\rho}-2\psi_{v}\overline{\beta}_{vx}\frac{\gamma^{2}}{\rho}\right)\sum_{k=0}^{q-1}g^{t,k}\enspace.$

Since $\psi_{z}=\frac{\rho}{16\overline{\beta}_{zx}}$ and
$\psi_{v}\leq\frac{\rho}{16\overline{\beta}_{vx}}$ and
$\frac{\gamma^{2}}{\rho}\leq\xi\leq 1$, we get

$\displaystyle\underbrace{\mathcal{L}^{t,q}-\mathcal{L}^{t,0}}_{\mathcal{L}^{t+1,0}-\mathcal{L}^{t,0}}$
$\displaystyle\leq-\frac{\gamma}{4}\sum_{k=0}^{q-1}g^{t,k}\enspace.$

Summing, telescoping and dividing by $Tq$ gives

$\frac{1}{Tq}\sum_{t=0}^{T-1}\sum_{k=0}^{q-1}g^{t,k}\leq\frac{4}{Tq\gamma}\underbrace{\left(h^{0,0}-h^{*}+\psi_{z}\phi^{0,0}+\psi_{v}\phi^{0,0}]\right)}_{\Gamma^{0}}\enspace.$

∎

From 1 we deduce Corollary 3.6.

###### Proof.

Let us take $\rho=\overline{\rho}(n+m)^{-\frac{1}{2}}$,
$\gamma=\min(\xi\rho,\overline{\gamma})$ and $q=n+m$. Then 1 holds:

$\frac{1}{Tq}\sum_{t=0}^{T-1}\sum_{k=0}^{q-1}g^{t,k}\leq\frac{4}{Tq\gamma}\Gamma^{0}\enspace.$

with $\Gamma_{0}=\mathcal{O}(1)$. To get an $\varepsilon$-stationary solution,
we set $T\geq\frac{4}{q\gamma}\Gamma^{0}\varepsilon^{-1}\vee
1=\mathcal{O}\left(\frac{1}{q\gamma\varepsilon}\vee 1\right)$. One iteration
has $\Theta(q)=\Theta(n+m)$ oracle complexity. As a consequence, the sample
complexity to get an $\varepsilon$-stationary point is
$\mathcal{O}\left((n+m)^{\frac{1}{2}}\epsilon^{-1}\vee(n+m)\right)$. ∎

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

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

###  B.1 Preliminary results

In what follows, we provide the building blocks of our worst-case instance.
The proof uses the following quadratic function presented by [Nesterov, 2018].

######  Definition B.1.

Let $d\in\mathbb{N}_{>0}$, $\xi\in[0\,,\,+\infty)$ and $\zeta\leq 1$. We
define $\mathbf{Q}(.;\xi,d):\mathbb{R}^{d}\to\mathbb{R}$ by

$\mathbf{Q}(x;\xi,d)=\frac{\xi}{2}(x_{1}-1)^{2}+\frac{1}{2}\sum_{k=1}^{d-1}(x_{k+1}-x_{k})^{2}\enspace.$

Proposition B.2 proposition comes directly from [Zhou and Gu, 2019,
Proposition 3.5]. The first part of the proposition gives us the regularity of
$\mathbf{Q}$. In the second part shows that a function defined as
$\mathbf{Q}(U\times\,\cdot\,;\xi,d)+\sum_{p=1}^{q}g(\langle
u_{p},\,\cdot\,\rangle)$ verifies the so-called ”zero-chain property” Carmon
et al. [2017b]: if $Ux~{}\in~{}\mathrm{Span}(u_{1},\dots,u_{k})$, then we gain
a non zero coordinate by calling the gradient
$\nabla~{}[\mathbf{Q}~{}(U~{}\times~{}\,\cdot\,;\xi,d)+\sum_{p=1}^{q}g(\langle
u_{p},\,\cdot\,\rangle)](x)$. In other words, that makes us progress in the
problem resolution.

######  Proposition B.2.

For $d\in\mathbb{N}_{>0}$, $\xi\in[0\,,\,+\infty)$ and $\zeta\leq 1$. The
following holds:

  1. 1.

$\mathbf{Q}(\,\cdot\,;\xi,d)$ is convex and 4-smooth.

  2. 2.

Let $q\in\mathbb{N}_{>0}$, $U=[u_{1},\dots,u_{d}]^{\top}\in\mathbb{R}^{d\times
q}$ such that $UU^{\top}=I$ and for $k\leq d$, $U^{(\leq
k)}~{}=~{}[u_{1},\dots,u_{k},0,\dots,0]^{\top}\in\mathbb{R}^{d\times q}$. Let
$g:\mathbb{R}\to\mathbb{R}$ differentiable such that $g^{\prime}(0)=0$. Then
for any $x\in\mathbb{R}^{q}$ such that $Ux=U^{(\leq k)}x$, then

$\nabla\left[\mathbf{Q}(U\times\,\cdot\,;\xi,d)+\sum_{p=1}^{d}g(\langle
u_{p},\,\cdot\,\rangle)\right](x)\in\mathrm{Span}(u_{1},\dots,u_{k},u_{k+1})\enspace.$

###### Proof.

Let $x\in\mathbb{R}^{q}$ such that $Ux=U^{(\leq k)}x$. For $0\geq k\leq d$, we
denote

$\mathbb{R}^{k,d}=\\{v\in\mathbb{R}^{d},v_{k+1}=\cdots=v_{d}=0\\}\enspace.$

Let us write $\mathbf{Q}(x;\xi,d)=\frac{1}{2}x^{\top}Ax+b^{\top}x+c$ with

$A=\begin{bmatrix}1+\xi&-1&0&\cdots&0\\\ -1&2&-1&\ddots&\vdots\\\
0&\ddots&\ddots&\ddots&0\\\ \vdots&\ddots&-1&2&-1\\\
0&\cdots&0&-1&1\end{bmatrix}\in\mathbb{R}^{d\times d}\enspace,$

$b=\xi(1,0,\dots,0)^{\top}$ and $c=\frac{\xi}{2}(1,0,\dots,0)^{\top}$.

On the one hand it is known from [Nesterov, 2018, Lemma 2.5.1] that if
$v\in\mathbb{R}^{k,d}$,

$\displaystyle\nabla\mathbf{Q}(v;\xi,d)\in\mathbb{R}^{k+1,d}$

As a consequence,

$\nabla\mathbf{Q}(Ux;\xi,d)=\nabla\mathbf{Q}(\underbrace{U^{(\leq
k)}x}_{\in\mathbb{R}^{k,d}};\xi,d)\in\mathbb{R}^{k+1,d}$

and

$\nabla[\mathbf{Q}(U\times\,\cdot\,;\xi,d)](x)=U^{\top}\nabla\mathbf{Q}(Ux;\xi,d)\in\mathrm{Span}(u_{1},\dots,u_{k+1})\enspace.$

On the other hand,

$\nabla\left[\sum_{p=1}^{d}g(\langle
u_{p},\,\cdot\,\rangle)\right](x)=\sum_{p=1}^{d}g^{\prime}(\langle
u_{p},x\rangle)u_{p}=\sum_{p=1}^{k}g^{\prime}(\langle
u_{p},x\rangle)u_{p}\in\mathrm{Span}(u_{1},\dots,u_{k+1})\enspace.$

Thus

$\nabla\left[\mathbf{Q}(U\times\,\cdot\,;\xi,d)+\sum_{p=1}^{d}g(\langle
u_{p},\,\cdot\,\rangle)\right](x)\in\mathrm{Span}(u_{1},\dots,u_{k},u_{k+1})\enspace.$

∎

However, the function $\mathbf{Q}$ is convex. That is why we also use the
function $\Gamma$ introduced in Carmon et al. [2017a]. As explained in Carmon
et al. [2017a], this function is essential to lower bound the gradient of our
worst wase instance.

######  Definition B.3.

Let $d\in\mathbb{N}_{>0}$. We define
$\Gamma(\,\cdot\,;d):\mathbb{R}^{d+1}\to\mathbb{R}$ by

$\Gamma(x;d)=120\sum_{k=1}^{d}\int_{1}^{x_{k}}\frac{t^{2}(t-1)}{1+t^{2}}\mathrm{d}t\enspace.$

An important property of $\Gamma$ shown in [Carmon et al., 2017a, Lemma 2] is
the smoothness of the function $\Gamma$.

######  Proposition B.4.

There exists a constant $c>0$ such that $\Gamma(\,\cdot\,;d)$ is $c$-smooth.

Now we introduce the function $f_{\mathrm{nc}}$ we use to build our worst-case
instance. This function comes from [Zhou and Gu, 2019, Definition 3.5]. It is
the sum of the quadratic function defined by B.1 and the non-convex component
given by Definition B.3.

######  Definition B.5.

For $\alpha>0$ and $d\in\mathbb{N}_{>0}$,
$f_{\mathrm{nc}}(\,\cdot\,;\alpha,d):\mathbb{R}^{d+1}\to\mathbb{R}$ is defined
a

$f_{\mathrm{nc}}(x;\alpha,d)=\mathbf{Q}(x;\sqrt{\alpha},d+1)+\alpha\Gamma(x)\enspace.$

The essential properties of $f_{\mathrm{nc}}$ come from [Carmon et al., 2017a,
Lemmas 2, 3, 4]. The first part provides the regularity properties of
$f_{\mathrm{nc}}$. The second part bounds the distance between
$f_{\mathrm{nc}}(\,\cdot\,;\alpha,d)$ and the optimal value of the function.
The third part will be key to the overall proof. In words, it states that as
long $x\in\mathbb{R}^{d+1}$ has its two last components equal to zero, the
norm of the gradient of $f_{\mathrm{nc}}(\,\cdot\,;\alpha,d)$ is higher than a
constant controlled by $\alpha$. As a consequence, if $\alpha$ is properly
chosen, as soon as $x_{d}=x_{d+1}=0$, we are ensured that $\|\nabla
f_{\mathrm{nc}}(x;\alpha,d)\|\geq\epsilon$.

######  Proposition B.6.

For $\alpha\in[0\,,\,1]$, it holds

  1. 1.

$-\alpha c\preceq\nabla^{2}f_{\mathrm{nc}}\preceq 4+\alpha c$.

  2. 2.

$f_{\mathrm{nc}}(0;\alpha,d)-\inf_{x}f_{\mathrm{nc}}(x;\alpha,d)\leq\frac{\sqrt{\alpha}}{2}+10\alpha
d$.

  3. 3.

For $x\in\mathbb{R}^{d+1}$ such that $x_{d}=x_{d+1}=0$, $\|\nabla
f_{\mathrm{nc}}(x;\alpha,d)\|\geq\frac{\alpha^{\frac{3}{4}}}{4}$.

From now we denote

$\mathcal{O}(a,b)=\\{U\in\mathbb{R}^{a\times b},UU^{\top}=I_{a}\\}\enspace.$

The following Lemma adapted from Zhou and Gu [2019] is fundamental for our
lower bound proof.

######  Lemma B.7.

Let $d,m\in\mathbb{N}_{>0}$ and $U\in\mathcal{O}((d+1)m,(d+1)m)$. We denote