# A regularized variance-reduced modified extragradient method for stochastic
hierarchical games

Shisheng Cui School of Automation, Beijing Institute of Technology
Beijing 100081, China
<EMAIL_ADDRESS> Uday V. Shanbhag Department of Industrial and Manufacturing
Engineering
Pennsylvania State University, University Park, USA
<EMAIL_ADDRESS> Mathias Staudigl Corresponding Author Department of
Mathematics, Universität Mannheim,
B6, 68159, Mannheim, Germany
<EMAIL_ADDRESS>
###### Abstract

We consider an $N$-player hierarchical game in which the $i$th player’s
objective comprises of an expectation-valued term, parametrized by rival
decisions, and a hierarchical term. Such a framework allows for capturing a
broad range of stochastic hierarchical optimization problems, Stackelberg
equilibrium problems, and leader-follower games. We develop an iteratively
regularized and smoothed variance-reduced modified extragradient framework for
iteratively approaching hierarchical equilibria in a stochastic setting. We
equip our analysis with rate statements, complexity guarantees, and almost-
sure convergence results. We then extend these statements to settings where
the lower-level problem is solved inexactly and provide the corresponding rate
and complexity statements. Our model framework encompasses many game theoretic
equilibrium problems studied in the context of power markets. We present a
realistic application to the virtual power plants, emphasizing the role of
hierarchical decision making and regularization.

##  1 Introduction

In this paper we consider a class of stochastic hierarchical optimization
problems and games, generalizing many learning problems involving sequential
optimization. Consider a collection of $N$-agents, where the $i$-th agent
solves the optimization problem parametrized by rival decisions
$\mathbf{x}_{-i}$:

$\min_{\mathbf{x}_{i},\mathbf{y}_{i}}\,\\{\ell_{i}(\mathbf{x}_{i},\mathbf{x}_{-i},\mathbf{y}_{i})\triangleq
f_{i}(\mathbf{x}_{i},\mathbf{x}_{-i})+g_{i}(\mathbf{x}_{i},\mathbf{y}_{i})\\},\text{
s.t.
}\mathbf{x}_{i}\in\mathcal{X}_{i},\mathbf{y}_{i}\in\operatorname{SOL}(\phi_{i}(\mathbf{x}_{i},\cdot),\mathcal{Y}_{i}).$
(P)

We let $i\in\mathcal{I}\triangleq\\{1,\ldots,N\\}$ represent a set of
_leaders_ , characterized by two loss functions: (i)
$f_{i}(\mathbf{x})\triangleq\mathbb{E}_{\xi}[F_{i}(\mathbf{x},\xi)]$, which
depends on the entire action profile
$\mathbf{x}\triangleq(\mathbf{x}_{i},\mathbf{x}_{-i})=(\mathbf{x}_{1},\ldots,\mathbf{x}_{N})\in\mathcal{X}\triangleq\prod_{i\in\mathcal{I}}\mathcal{X}_{i}$;
(ii) $g_{i}(\mathbf{x}_{i},\mathbf{y}_{i})$ is a deterministic function,
jointly controlled by leader $i$’s decision variable
$\mathbf{x}_{i}\in\mathcal{X}_{i}$ and a follower’s decision variable
$\mathbf{y}_{i}\in\mathcal{Y}_{i}$. Each leader’s optimization problem
exhibits two sets of private constraints, the first given by
$\mathbf{x}_{i}\in\mathcal{X}_{i}\subseteq\mathbb{R}^{n_{i}}$, while the
second are equilibrium constraints represented by the solution set of a
parameterized variational inequality
$\operatorname{VI}(\phi_{i}(\mathbf{x}_{i},\cdot),\mathcal{Y}_{i})$, which
reads as

$\text{Find }\mathbf{y}_{i}(\mathbf{x}_{i})\in\mathcal{Y}_{i}\text{ satisfying
}\langle\phi_{i}(\mathbf{x}_{i},\mathbf{y}_{i}(\mathbf{x}_{i})),\mathbf{y}_{i}-\mathbf{y}_{i}(\mathbf{x}_{i})\rangle\geq
0\quad\forall\mathbf{y}_{i}\in\mathcal{Y}_{i}.$ (1.1)

We denote the set of points $\mathbf{y}_{i}$ satisfying this condition by
$\operatorname{SOL}(\phi_{i}(\mathbf{x}_{i},\cdot),\mathcal{Y}_{i})$. This VI
is defined in terms of a closed convex set
$\mathcal{Y}_{i}\subseteq\mathbb{R}^{m_{i}}$ and an expectation-valued mapping
$\phi_{i}(\mathbf{x},\mathbf{y}_{i})\triangleq\mathbb{E}_{\xi}[\Phi_{i}(\mathbf{x},\mathbf{y}_{i},\xi)]$.
All the problem data are affected by random noise represented by a random
variable $\xi:\Omega\to\Xi$, defined on a probability space
$(\Omega,\mathcal{F},\mathbb{P})$ and taking values in a measurable space
$\Xi$. Such hierarchical optimization problems traditionally play a key role
in operations research and engineering, where they are deeply connected to
bilevel programming [16] and mathematical programs under equilibrium
constraints (MPEC) [49]. In fact, the canonical MPEC formulation is obtained
from (P) when $N=1$. The multi-agent formulation (P) also relates to leader-
follower and Stackelberg games (cf. [63, 57, 15, 43]) which are a traditional
model in economics, and also have received increased attention in machine
learning recently [48, 22, 4, 7]. Economic equilibria in power markets have
been extensively studied using a complementarity framework (cf. [30, 31, 32,
24]). More recently, stochastic generalizations have been examined where
uncertainty in price and cost functions have been addressed [62, 40]. As an
immediate application of our algorithmic framework, we present in Section 5 a
model inspired by Hobbs and Pang [32], but suitably modified to account for
uncertainty in prices and costs, multi-period settings with ramping
constraints, and the incorporation of virtual power plants (VPPs) (see [51,
19] for a review of VPPs and power markets).

###  1.1 Our contributions and related work

##### Hierarchical optimization, games and uncertain generalizations.

To date, hierarchical optimization has been studied under the umbrella of
bilevel programming [16, 17] and mathematical programs with equilibrium
constraints (MPECs) [49, 55]. Algorithmic schemes for resolving MPECs where
the lower-level problem is an optimization problem, or a variational
inequality, have largely emphasized either implicit approaches [55] or
regularization/penalization-based techniques [49, 47, 34]. Yet, there appears
to have been a glaring lacuna in non-asymptotic rate and complexity guarantees
for resolving hierarchical optimization and their stochastic and game-
theoretic variants. This gap has been partially addressed in the recent papers
[14] and [13]. Both papers present variance-reduced solution strategies for
various versions of hierarchical optimization problems and games,
respectively, relying on variance reduction via a sequence of increasing mini-
batches. Finite-time and almost sure convergence to solutions is proved under
convexity/monotonicity assumptions on the problem data. However, the framework
for monotone games in [13] requires exact solutions of lower-level problems,
significantly impacting its efficient implementation in large-scale settings.
We complement this literature by developing a novel regularized smoothed
variance reduction method for the family of hierarchical games (P), building
on a disciplined operator splitting approach. Notably, we provide an inexact
generalization allowing for random error-afflicted lower-level solutions,
addressing a significant shortcoming in [13]. Specifically, in this paper we
improve [13] along two important dimensions: (i) First, we allow for inexact
resolution of lower-level problems to accommodate large-scale stochastic
follower problems; (ii) Second, we provide a novel variance-reduction
framework for addressing this problem. Despite the need for inexactness, our
statements match the state-of-the-art both in terms of rate and oracle
complexity. Algorithmically, these advancements are achieved via a novel
stochastic operator splitting approach that combines ideas from iterative
smoothing and regularization [71], with modern variance reduction approaches
originating in machine learning [29].

##### Zeroth-order optimization, smoothing, and regularization.

Zeroth-order (gradient-free) optimization is being increasingly embraced for
solving machine learning problems where explicit expressions of the gradients
are difficult or infeasible to obtain. In hierarchical optimization problems,
this is particularly relevant when solutions of the lower level problem are
injected into the leader’s upper-level problem. In machine learning, this
problem is known as approximating the hypergradient. Various techniques for
estimating this object have been studied recently, ranging from truncated von
Neuman series [26] and fully first-order methods [44]. Instead of
computationally expensive first-order (or higher) information about the
problem data, we develop an online stochastic approximation approach based
solely on function evaluations in order to approximate the directional
derivative of the coupling function between the upper and the lower level.
Approximating the directional derivative of the thus obtained implicit
function has a long history [64] and has been employed for resolving
stochastic optimization [45, 68, 69] and variational inequality problems [70,
71].

##### Variance Reduction.

Variance reduction is a commonly employed method exploiting the finite-sum
structure of variational problems arising in machine learning and engineering.
The classical stochastic variance reduced gradient (SVRG) [37] is embedded
within a double loop structure and tailored to the prototypical finite-sum
structure in empirical risk minimization. Indeed, in the classical SVRG
formulation full gradients are computed “from time-to-time” in the outer loop
while cheap variance reduced gradients are used in the frequently activated
inner loop subroutine. This construction has been extended to saddle-point
problems and stochastic monotone inclusions in [56]. Extensions to monotone
_mixed_ -variational inequality problems were recently provided in [2, 1] and
[10]. In contrast, we consider stochastic hierarchical games over general
sample spaces, complicated by the presence of nested optimization problems
embodied by the interaction between leaders and followers. While the
assumptions we make in this work allow us to recast the problem as a mixed
variational inequality, several challenges persist. First, subgradients of
$g_{i}(\cdot,\mathbf{y}_{i}(\cdot))$ are unavailable; Second, enlisting
smoothing approaches requires $\mathbf{y}_{i}(\mathbf{x}_{i})$, unavailable in
closed form; Third, we are not restricted to finite-sum regimes and allow for
general sample spaces by employing increasingly large batch-sizes to
approximate the gradient in the outer loop.

##  2 Preliminaries

In this section, we articulate the standing assumptions employed in this paper
and introduce our notation. The decision set of leader $i$ is a subset
$\mathcal{X}_{i}$ in $\mathbb{R}^{n_{i}}$. We let
$\mathcal{X}\triangleq\prod_{i\in\mathcal{I}}\mathcal{X}_{i}$ represent the
set of strategy profiles of the leaders and identify it with a subset of
$\mathbb{R}^{n}$, where $n\triangleq\sum_{i\in\mathcal{I}}n_{i}$. For any
$d\geq 1$, we let
$\mathbb{B}_{d}\triangleq\\{\mathbf{x}\in\mathbb{R}^{d}|\;\norm{\mathbf{x}}\leq
1\\}$ denote the unit ball in $\mathbb{R}^{d}$. We start by introducing a
basic assumption on the follower’s problem. It bears reminding that the VI
representation of the follower problem allows for capturing a range of
problems, ranging from smooth convex optimization problems to more intricate
smooth convex games and equilibrium problems; see [21].111Appendix A.2
explains the terminology related to VIs.

######  Assumption 1.

For $i\,\in\,\mathcal{I}$, $\mathcal{Y}_{i}\subseteq\mathbb{R}^{m_{i}}$ is
closed and convex set, and for all $\mathbf{x}_{i}\in\mathcal{X}_{i}$, the
mapping $\phi_{i}(\mathbf{x}_{i},\cdot):\mathcal{Y}_{i}\to\mathbb{R}^{m_{i}}$
is strongly monotone and Lipschitz continuous.

By Assumption 1, for any $i\,\in\,\mathcal{I}$, the set of solutions to
$\operatorname{VI}(\phi_{i}(\mathbf{x}_{i},\cdot),\mathcal{Y}_{i})$, denoted
by $\operatorname{SOL}(\phi_{i}(\mathbf{x}_{i},\cdot),\mathcal{Y}_{i})$, is
single-valued with unique element $\mathbf{y}_{i}(\mathbf{x}_{i})$. Moreover,
$\phi_{i}(\mathbf{x}_{i},\mathbf{y}_{i})=\mathbb{E}_{\xi}\left[\Phi_{i}(\mathbf{x}_{i},\mathbf{y}_{i},\xi)\right]$
for any
$(\mathbf{x}_{i},\mathbf{y}_{i})\in\mathcal{X}_{i}\times\mathcal{Y}_{i}$.

######  Assumption 2.

The following assumptions hold for each leader $i\in\mathcal{I}$:

  * (i)

The set $\mathcal{X}_{i}\subset\mathbb{R}^{n_{i}}$ is nonempty, compact, and
convex. In particular, there exists $C_{i}>0$ such that
$\sup_{\mathbf{x}_{i},\mathbf{x}^{\prime}_{i}\in\mathcal{X}_{i}}\norm{\mathbf{x}_{i}-\mathbf{x}^{\prime}_{i}}\leq
C_{i}$ for all $i\in\mathcal{I}$.

  * (ii)

For some $\delta_{0}>0$, the mapping $\mathbf{x}_{i}\mapsto
g_{i}(\mathbf{x}_{i},\mathbf{y}_{i}(\mathbf{x}_{i}))$ is $L_{1,{i}}$-Lipschitz
on
$\mathcal{X}_{i,\delta_{0}}\triangleq\mathcal{X}_{i}+\delta_{0}\mathbb{B}_{n_{i}}$.

  * (iii)

$\mathbf{x}_{i}\mapsto F_{i}((\mathbf{x}_{i},\mathbf{x}_{-i}),\xi)$ is convex
and continuously differentiable over an open set containing $\mathcal{X}_{i}$,
uniformly for all $\mathbf{x}_{-i}\in\mathcal{X}_{-i}$ and almost every
$\xi\in\Xi$.

  * (iv)

The mapping $\mathbf{y}_{i}\mapsto g_{i}(\mathbf{x}_{i},\mathbf{y}_{i})$ is
$L_{2,i}$-Lipschitz continuous for all
$\mathbf{x}_{i}\in\mathcal{X}_{i,\delta_{0}}$.

  * (v)

The operator $V:\mathbb{R}^{n}\to\mathbb{R}^{n},$ defined by
$V(\mathbf{x})=(V_{i}(\mathbf{x}))_{i\in\mathcal{I}}$ and
$V_{i}(\mathbf{x})\triangleq\nabla_{\mathbf{x}_{i}}f_{i}(\mathbf{x})$, is
$L_{f}$-Lipschitz continuous and monotone on
$\mathcal{X}_{\delta_{0}}=\prod_{i\in\mathcal{I}}\mathcal{X}_{i,\delta_{0}}$.

Assumption 1 is commonly employed in hierarchical optimization problems.
Indeed, in the special case where the VI captures the optimality conditions of
a parametrized convex optimization problem solved by the follower, then strong
monotonicity of $\phi_{i}(\mathbf{x}_{i},\cdot)$, is equivalent to strong
convexity of the follower’s cost function, an assumption that is known in
bilevel optimization literature as the _lower level uniqueness property_ [48].

Given Assumption 1 we define the implicit loss function
$L_{i}:\mathcal{X}\to\mathbb{R}\cup\\{\infty\\}$ by

$L_{i}(\mathbf{x}_{i},\mathbf{x}_{-i})\triangleq\ell_{i}(\mathbf{x}_{i},\mathbf{x}_{-i},\mathbf{y}_{i}(\mathbf{x}_{i}))=\mathbb{E}_{\xi}[F_{i}((\mathbf{x}_{i},\mathbf{x}_{-i}),\xi)]+h_{i}(\mathbf{x}_{i}),$
(2.1)

where $h_{i}(\mathbf{x}_{i})\triangleq
g_{i}(\mathbf{x}_{i},\mathbf{y}_{i}(\mathbf{x}_{i})).$ In terms of the
implicit loss function (2.1), we convert the hierarchical game (P) into a
_stochastic Nash equilibrium problem_ in which each player solves the loss
minimization problem

$(\forall
i\in\mathcal{I}):\quad\min_{\mathbf{x}_{i}\in\mathcal{X}_{i}}L_{i}(\mathbf{x}_{i},\mathbf{x}_{-i}).$
(2.2)

We refer to (2.2) as the _upper-level problem_ , and summarize it as the tuple
$\mathcal{G}^{\text{upper}}\triangleq\\{L_{i},\mathcal{X}_{i}\\}_{i\in\mathcal{I}}$.
Let $\mathbf{NE}(\mathcal{G}^{\text{upper}})$ denote the set of Nash
equilibria of the game $\mathcal{G}^{\text{upper}}$ and
$\mathcal{Y}=\prod_{i\in\mathcal{I}}\mathcal{Y}_{i}$.

######  Definition 2.1.

A $2N$-tuple
$(\mathbf{x}^{\ast},\mathbf{y}^{\ast})\in\mathcal{X}\times\mathcal{Y}$ is
called a _hierarchical equilibrium_ if
$\mathbf{x}^{\ast}\in\mathbf{NE}(\mathcal{G}^{\text{upper}})$ and, for all
$i\in\mathcal{I}$,
$\mathbf{y}_{i}^{\ast}=\mathbf{y}_{i}(\mathbf{x}_{i}^{\ast})$ the unique
solution of
$\operatorname{VI}(\phi_{i}(\mathbf{x}^{\ast}_{i},\cdot),\mathcal{Y}_{i})$.

Typical online learning approaches in game theory employ stochastic
approximation (SA) for iteratively approaching a Nash equilibrium of the game
$\mathcal{G}^{\text{upper}}$. These iterative methods rely on the availability
of a stochastic oracle revealing (noisy) first-order information about the
operators involved (i.e. samples of pseudo-gradients of the objective for each
individual player). Such direct methods are complicated in hierarchical
optimization since the required subgradient is an element of the
subdifferential of the sum of two Lipschitz continuous functions. In the
current setting, this task is even more complicated since the upper level
objective is defined by a function which is available only in an implicit
form, as it depends on the solution of the lower level problem
$\mathbf{y}_{i}(\mathbf{x}_{i})$, and another function given in terms of an
expected value. Thus, even if a sum-rule for a subdifferential applies [11],
it would read as
$\partial_{\mathbf{x}_{i}}L(\mathbf{x}_{i},\mathbf{x}_{-i})=\partial_{\mathbf{x}_{i}}f_{i}(\mathbf{x})+\partial_{\mathbf{x}_{i}}(g_{i}\circ(\operatorname{Id},\mathbf{y}_{i}(\cdot)))(\mathbf{x}_{i}).$
Hence we would need to invoke a non-smooth chain rule for our chosen version
of a subdifferential, in order to evaluate
$\partial_{\mathbf{x}_{i}}(g_{i}\circ(\operatorname{Id},\mathbf{y}_{i}(\cdot)))(\mathbf{x}_{i})$.
In addition, we would require access to the subdifferential of
$\mathbf{y}_{i}(\bullet)$ at $\mathbf{x}_{i}$. We circumvent this
computationally challenging step by developing a random search procedure based
on a finite difference approximation. To develop such scheme, recall that we
defined $h_{i}(\mathbf{x}_{i})\triangleq
g_{i}(\mathbf{x}_{i},\mathbf{y}_{i}(\mathbf{x}_{i}))$ as the loss function
coupling the leader and the follower. The following fact can be found in
Proposition 1 in [14].

######  Lemma 2.2.

Let Assumptions 1-2 hold. Then
$h_{i}(\mathbf{x}_{i})=g_{i}(\mathbf{x}_{i},\mathbf{y}_{i}(\mathbf{x}_{i}))$
is $L_{h_{i}}$-Lipschitz continuous on $\mathcal{X}_{i}$ and directionally
differentiable.

To proceed, we impose a convexity requirement on
$h_{i}(\cdot,\mathbf{y}_{i}(\cdot))$.

######  Assumption 3.

The implicitly defined function $\mathbf{x}_{i}\mapsto
h_{i}(\mathbf{x}_{i})\triangleq
g_{i}(\mathbf{x}_{i},\mathbf{y}_{i}(\mathbf{x}_{i}))$ is convex on
$\mathcal{X}_{i}$.

######  Remark 2.1.

Several papers in the literature provided conditions under which the implicit
function $h_{i}$ is indeed convex in hierarchical settings. A structural model
framework where convexity provably holds is described in [33], and [15].

######  Remark 2.2.

Under Assumptions 2 and 3, an equilibrium exists using classical results due
to [28]. Indeed, since the lower level problem is assumed to have a unique
solution $\mathbf{y}_{i}(\mathbf{x}_{i})$ and Lemma 2.2 guarantees that the
noncooperative game $\mathcal{G}^{\text{upper}}$ satisfies conditions of [50,
Th. I.4.1], implying $\mathbf{NE}(\mathcal{G}^{\text{upper}})\neq\varnothing$.

Assumption 2(v), Assumption 3 and Lemma 2.2, yield a variational
characterization of elements of $\mathbf{NE}(\mathcal{G}^{\text{upper}})$ in
terms of an expectation-valued _mixed-variational inequality_ (cf. Appendix
A.2).

######  Lemma 2.3.

Let $h:\mathcal{X}\to\mathbb{R}$ be defined by
$h(\mathbf{x})\triangleq\sum_{i\in\mathcal{I}}h_{i}(\mathbf{x}_{i})$. Then,
$\mathbf{x}^{\ast}\in\mathbf{NE}(\mathcal{G}^{\text{upper}})$ if and only if
$\mathbf{x}^{\ast}$ solves the mixed variational inequality
$\operatorname{MVI}(V,h)$:

$\langle
V(\mathbf{x}^{\ast}),\mathbf{x}-\mathbf{x}^{\ast}\rangle+h(\mathbf{x})-h(\mathbf{x}^{\ast})\geq
0\qquad\forall\mathbf{x}\in\mathcal{X}.$ (2.3)

Solutions to mixed VIs with expectation-valued operators have been developed
recently in cases where the random variable takes values in a finite set,
and/or when the VI is derived from a zero-sum game displaying a finite-sum
structure [2, 1]. The standard algorithmic approach to iteratively approximate
a solution to such structured VIs are extragradient type of methods. A direct
application of these methods to the mixed VI (2.3) is complicated because of
the following facts:

(i) $L$-smoothness of $h_{i}$. The assumptions made thus far do not guarantee
the differentiability of $h_{i}$ with a Lipschitz continuous gradient. Hence,
a direct application of gradient, or extragradient methods, is a difficult
task in our setting. To cope with this technical difficulty, we develop a
smoothing approach, yielding a family of approximating models enjoying the
typical Lipschitz smoothness requirements.
(ii) Randomness in the operator $V$: Since the operator $V$ is only available
in terms of an expected value, in general, we cannot tractably evaluate it.
Instead we have to use simulation-based methods to obtain random estimators of
this mathematical expectation. To keep this simulation task within a feasible
computational budget, iterative variance reduction ideas are frequently used
in iterative methods for generating the input data. Again, these standard
variance reduction techniques rely on smoothness of the data. In our case,
non-smoothness is present in terms of the implicit function $h_{i}(\cdot)$. In
principle, one could apply splitting techniques to deal with the non-smooth
function via a proximal smoothing. However, this approach requires
$h_{i}(\cdot)$ to be proximable for which we have no a-priori guarantee since
it is the value function of the leader, derived from the solution of the
follower. With these preparatory remarks in mind, we now explain the design of
our algorithmic solution strategy for the hierarchical game problems (P).

##  3 A variance reduced forward-backward-forward algorithm for hierarchical
games

In this section, we present our algorithm for computing an equilibrium of the
hierarchical game (P). As in the seminal SVRG formulation, our method runs in
two loops. Each loop requires as inputs data that are computed in the outer
loop. The inputs of the inner and outer loops are constructed as follows. For
$\eta>0$, we define the (Tikhonov) regularized vector field
$V^{\eta}:\mathcal{X}\to\mathbb{R}^{n}$ by

$V^{\eta}(\mathbf{x})=(V_{i}^{\eta}(\mathbf{x}))_{i\in\mathcal{I}},\mbox{
where }V^{\eta}_{i}(\mathbf{x})\triangleq
V_{i}(\mathbf{x})+\eta\mathbf{x}_{i}\quad\forall i\in\mathcal{I}.$ (3.1)

Tikhonov regularization is a classical tool to obtain stronger convergence
results in numerical schemes. It has been examined for deterministic [38] and
stochastic equilibrium problems [41, 71].

Our next assumption is concerned with the nature of the stochastic oracle
which generated random estimators on the expectation-valued operator $V$ when
queried at a given point $\mathbf{x}$.

######  Assumption 4.

The operator $V$ has a stochastic oracle $\hat{V}(\cdot,\xi)$ that is

  1. 1.

unbiased: $V(\mathbf{x})=\mathbb{E}_{\xi}[\hat{V}(\mathbf{x},\xi)]$ for all
$\mathbf{x}\in\mathcal{X}$;

  2. 2.

$\mathcal{L}_{f}(\xi)$-Lipschitz for almost every $\xi\in\Xi$:
$\norm{\hat{V}(\mathbf{x}^{\prime},\xi)-\hat{V}(\mathbf{x},\xi)}\leq\mathcal{L}_{f}(\xi)\norm{\mathbf{x}^{\prime}-\mathbf{x}}$
for all $\mathbf{x},\mathbf{x}^{\prime}\in\mathcal{X}$. The random variable
$\mathcal{L}_{f}(\xi)$ is positive and integrable with
$\mathbb{E}_{\xi}[\mathcal{L}(\xi)]=L_{f}$.

To retrieve in-play information about the value of the implicit loss function
$h_{i}(\cdot)$, we employ a smoothing-based approach, which necessitates
defining another sampling mechanism. We follow the gradient sampling strategy
of [23], though alternative random estimation strategies are certainly
possible (see, e.g. [20, 6, 42]). Specifically, given $\delta>0$, we denote
the finite difference approximation of the directional derivative of $h_{i}$
in direction $\mathbf{w}_{i}\in\mathbb{R}^{n_{i}}$ as

$\nabla_{(\mathbf{w}_{i},\delta)}h_{i}(\mathbf{x}_{i})\,\triangleq\,\frac{h_{i}(\mathbf{x}_{i}+\delta\mathbf{w}_{i})-h_{i}(\mathbf{x})}{\delta}.$

Let $\mathbf{W}_{i}$ be a random vector uniformly distributed on the unit
sphere
$\mathbb{S}_{i}\triangleq\\{\mathbf{x}_{i}\in\mathbb{R}^{n_{i}}|\norm{\mathbf{x}_{i}}=1\\}$.222See
Appendix A.5 for the explicit construction of such an oracle. We then define
the random vector $H_{i,\mathbf{x}_{i}}^{\delta}(\mathbf{W}_{i})$ as a
randomized and suitably rescaled version of the finite difference approximator
reading as

$H_{i,\mathbf{x}_{i}}^{\delta}(\mathbf{W}_{i})\triangleq
n_{i}\mathbf{W}_{i}\nabla_{(\mathbf{W}_{i},\delta)}h_{i}(\mathbf{x}_{i})\in\mathbb{R}^{n_{i}}.$
(3.2)

From eq. (A.8) in App. A.4, we know that
$H^{\delta_{i}}_{i,\mathbf{x}_{i}}(\cdot)$ is an unbiased estimator of the
gradient of the smoothed function

$h^{\delta}_{i}(\mathbf{x}_{i})\triangleq\frac{1}{{\mathbf{Vol}}_{n}(\delta\mathbb{B}_{n_{i}})}\int_{\delta\mathbb{B}_{n_{i}}}h_{i}(\mathbf{x}_{i}+\bf{u})\differential\bf{u},$

where
$\mathbb{B}_{d}\triangleq\\{\mathbf{x}\in\mathbb{R}^{d}|\norm{\mathbf{x}}\leq
1\\}$ for any dimension $d\geq 1$. Furthermore, we discuss in Appendix A.4
that the function $h^{\delta}_{i}$ is continuously differentiable with
gradient

$\nabla
h^{\delta}_{i}(\mathbf{x}_{i})=\frac{n_{i}}{\delta}\mathbb{E}_{\mathbf{W}_{i}\sim\mathsf{U}(\mathbb{S}_{n_{i}})}[\mathbf{W}_{i}\left(h_{i}(\mathbf{x}_{i}+\delta\mathbf{W}_{i})-h_{i}(\mathbf{x}_{i})\right)]=\mathbb{E}_{\mathbf{W}_{i}\sim\mathsf{U}(\mathbb{S}_{n_{i}})}[H^{\delta}_{i,\mathbf{x}_{i}}(\mathbf{W}_{i})],$

and

$\norm{\nabla h_{i}^{\delta}(\mathbf{x}_{i})-\nabla
h_{i}^{\delta}(\mathbf{y}_{i})}\leq\frac{L_{h_{i}}n_{i}}{\delta}\norm{\mathbf{x}_{i}-\mathbf{y}_{i}},\qquad\forall\mathbf{x}_{i},\mathbf{y}_{i}\in\mathbb{R}^{n_{i}},\delta>0$

where $\mathbf{W}_{i}\sim\mathsf{U}(\mathbb{S}_{n_{i}})$ indicates that
$\mathbf{W}_{i}$ is uniformly distributed on the surface of a unit sphere
$\mathbb{S}_{n_{i}}$.

###  3.1 Iterative Regularization methods

Lemma 2.3 shows that the equilibria of our hierarchical game are entirely
captured by the solution set of problem $\operatorname{MVI}(V,h)$. This is a
rich class of variational problems for which the number of contributions is so
numerous that we just point the reader to the monographs [21] and [61].
Deducing convergence results on the _last iterate_ for standard algorithmic
schemes is an important requirement for game-theoretic learning algorithms,
but typically is a rare commodity: Despite some special classes of games [3,
27], first-order methods give only guarantees on a suitably constructed
ergodic average. To obtain last iterate convergence results, we develop an
iterativeTikhonov regularization approach. This leads us to consider the
regularized problem $\operatorname{MVI}(V^{\eta},h)$, which requires to find
$\mathbf{s}(\eta)\in\mathcal{X}$ satisfying

$\langle
V^{\eta}(\mathbf{s}(\eta)),\mathbf{x}-\mathbf{s}(\eta)\rangle+h(\mathbf{x})-h(\mathbf{s}(\eta))\geq
0\qquad\forall\mathbf{x}\in\mathcal{X}.$ (3.3)

Naturally, we would like to understand the nature of the accumulation points
of the sequence $\\{\mathbf{s}_{t}\\}_{t\in\mathbb{N}}$, where
$\mathbf{s}_{t}\equiv\mathbf{s}(\eta_{t})$ and $\eta_{t}\downarrow 0$. This
sequence can be studied in quite some detail, and we summarize some well-known
facts in Proposition 3.1 below. As those results are rather scattered in the
literature, we provide a self-contained proof in Appendix A.3.

######  Proposition 3.1.

Let Assumptions 1, 2 and 3 hold true. Consider the problem
$\operatorname{MVI}(V,h)$ with nonempty solution set
$\operatorname{SOL}(V,h)$. Then the following apply:

  1. (a)

For all $\eta>0$, the set $\operatorname{SOL}(V^{\eta},h)$ is a singleton with
unique element denoted by $\mathbf{s}(\eta)$;

  2. (b)

$(\forall\eta>0):\norm{\mathbf{s}(\eta)}\leq\inf\\{\norm{\mathbf{x}}:\mathbf{x}\in\operatorname{SOL}(V,h)\\}$;

  3. (c)

Let $\\{\eta_{t}\\}_{t\in\mathbb{N}}$ be a positive sequence with
$\eta_{t}\downarrow 0$. Then, the sequence
$\\{\mathbf{s}(\eta_{t})\\}_{t\in\mathbb{N}}$ converges to the least norm
solution
$\arg\min\\{\norm{\mathbf{x}}:\mathbf{x}\in\operatorname{SOL}(V,h)\\}$;

  4. (d)

For any positive sequence $\\{\eta_{t}\\}_{t\in\mathbb{N}}$ satisfying
$\eta_{t}\downarrow 0$, we have

$\left(\frac{\eta_{t}-\eta_{t-1}}{\eta_{t}}\right)\inf_{\mathbf{x}\in\operatorname{SOL}(V,h)}\norm{\mathbf{x}}\geq\norm{\mathbf{s}(\eta_{t})-\mathbf{s}(\eta_{t-1})}.$
(3.4)

###  3.2 The Algorithm

Our algorithm for the hierarchical game setting consists of a double-loop
structure: The outer loop allows the $N$ players to make multiple independent
queries of the stochastic oracle
$\hat{V}_{1}(\cdot,\xi),\ldots,\hat{V}_{N}(\cdot,\xi),$ and draw multiple
independent samples from the surface of a unit sphere, allowing for the
simulation of the random estimator (3.2). However, since the multiple calls
are a negative entry on the oracle complexity of the method, we impose some
control on the number of mini-batches to be constructed by the agents. Within
the inner-loop subroutine, the agents only receive single samples from their
stochastic oracles, and employ this new information in an extragradient-type
algorithm. We give a precise construction in the following paragraphs.

####  3.2.1 The outer loop

Let $t=0,\ldots,T-1$ be the iteration counter for the outer loop. We denote by
$b_{t}\in\mathbb{N}$ the pre-defined sample rate defining the number of random
variables each player is allowed to generate in round $t$. Specifically, each
player generates an iid sample
$\xi^{1:b_{t}}_{i,t}\triangleq\\{\xi_{i,t}^{(s)};1\leq s\leq b_{t}\\}$ and
constructs the mini-batch estimator

$\bar{V}^{t}_{i}\triangleq\frac{1}{b_{t}}\sum_{s=1}^{b_{t}}\hat{V}_{i}(\mathbf{x}^{t},\xi_{i,t}^{(s)}),$
(3.5)

Let $\bar{V}^{t}=(\bar{V}_{1}^{t},\ldots,\bar{V}_{N}^{t})$. A similar
assumption is made to obtain point-estimators of the gradient of the implicit
function $h_{i}(\cdot)$. Hence,
$\mathbf{W}^{1:b_{t}}_{i,t}\triangleq\\{\mathbf{W}^{(s)}_{i,t};1\leq s\leq
b_{t}\\}$ denotes an i.i.d. sample of $b_{t}$ random vectors drawn uniformly
at random from $\mathbb{S}_{i}$ and define the mini-batch estimator
$H_{i,\mathbf{x}_{i}}^{\delta_{t},b_{t}}$ as

$H_{i,\mathbf{x}_{i}}^{\delta_{t},b_{t}}\triangleq\frac{1}{b_{t}}\sum_{s=1}^{b_{t}}H_{i,\mathbf{x}_{i}}^{\delta_{t}}(\mathbf{W}_{i,t}^{(s)}),\quad
H^{\delta_{t},b_{t}}_{\mathbf{x}^{t}}\triangleq(H^{\delta_{t},b_{t}}_{1,\mathbf{x}^{t}_{1}},\ldots,H^{\delta_{t},b_{t}}_{N,\mathbf{x}^{t}_{N}}),$
(3.6)

where $\delta_{t}$ denotes a positive smoothing parameter. Equipped with these
estimators, each player enters the procedure
$\operatorname{\texttt{SFBF}}(\mathbf{x}^{t},\bar{V}^{t},H^{\delta_{t},b_{t}}_{\mathbf{x}^{t}},\gamma_{t},\eta_{t},\delta_{t},K)$,
that relies on steplength $\gamma_{t}$ and regularization parameter
$\eta_{t}$, whose role is explained in the description of the inner loop.

####  3.2.2 The inner loop

Given the inputs
$(\mathbf{x}^{t},\bar{V}^{t},H^{\delta_{t},b_{t}}_{\mathbf{x}_{t}},\gamma_{t},\eta_{t})$
prepared in the outer loop, the inner loop of our method is based on a
stochastic version of Tseng’s modified extragradient method [65], using one-
shot estimators of the relevant data. To be precise, given the current iterate
$\mathbf{x}^{t}$, each player $i$ produces a trajectory
$\\{\mathbf{z}^{(t)}_{i,k}\\}_{k\in\\{0,1/2,1,\ldots,K\\}}$. These interim
strategy profiles are updated recursively by the procedure
$\operatorname{\texttt{SFBF}}(\mathbf{x}^{t},\bar{V}^{t},H^{\delta_{t},b_{t}}_{\mathbf{x}^{t}},\gamma_{t},\eta_{t},\delta_{t},K)$
described in Algorithm 1. Starting with the strategy profile $\mathbf{x}^{t}$,
we choose the initial conditions $\mathbf{z}^{(t)}_{i,0}=\mathbf{x}_{i}^{t}$
for all $i\in\mathcal{I}$. Then, for each $k\in\\{0,1,\ldots,K-1\\}$ each
player queries the stochastic oracle to obtain the feedback signal

$\hat{V}^{\eta_{t}}_{i,t,k+1/2}(\mathbf{z}^{(t)}_{k+1/2})\triangleq\hat{V}_{i}(\mathbf{z}^{(t)}_{k+1/2},\xi_{i,t,k+1/2})+\eta_{t}\mathbf{z}_{i,k+1/2}^{(t)}.$
(3.7)

Similarly, each player obtains the random information
$H^{\delta_{t}}_{\mathbf{z}^{(t)}_{i,k+1/2}}(\mathbf{W}_{i,t,k+1/2})$ and
$H^{\delta_{t}}_{\mathbf{x}^{t}_{i}}(\mathbf{W}_{i,t,k+1/2})$, as defined in
(3.2). These random variables are used to generate the parallel updates

$\displaystyle\mathbf{z}^{(t)}_{i,k+1/2}=\Pi_{\mathcal{X}_{i}}[\mathbf{z}^{(t)}_{i,k}-\gamma_{t}(\bar{V}^{t}_{i}+\eta_{t}\mathbf{x}_{i}^{t}+H^{\delta_{t},b_{t}}_{\mathbf{x}_{i}^{t}})],\text{
and }$
$\displaystyle\mathbf{z}^{(t)}_{i,k+1}=\mathbf{z}^{(t)}_{i,k+1/2}-\gamma_{t}\left(\hat{V}^{\eta_{t}}_{i,t,k+1/2}(\mathbf{z}^{(t)}_{k+1/2})+H^{\delta_{t}}_{\mathbf{z}^{(t)}_{i,k+1/2}}(\mathbf{W}_{i,t,k+1/2})-\hat{V}^{\eta_{t}}_{i,t,k+1/2}(\mathbf{x}^{t})-H^{\delta_{t}}_{\mathbf{x}^{t}_{i}}(\mathbf{W}_{i,t,k+1/2})\right)$

for all $i\in\mathcal{I}$. These iterations correspond to a stochastic
approximation variant of Tseng’s forward-backward-forward method [8] for
solving the time-varying stochastic variational inequality

$0\in V(\bar{\mathbf{x}})+\nabla
h^{\delta_{t}}(\bar{\mathbf{x}})+\eta_{t}\bar{\mathbf{x}}+\operatorname{\mathsf{NC}}_{\mathcal{X}}(\bar{\mathbf{x}}),$

characterized by Tikhonov regularization and smoothing of the implicit
function $h_{i}$.

Result: Iterate $\mathbf{z}^{K}$

Set $\mathbf{z}^{0}=\bar{\mathbf{x}}$;

for  _$k=0,1,\ldots,K-1$_ do

Update
$\mathbf{z}_{k+1/2}=\Pi_{\mathcal{X}}[\mathbf{z}_{k}-\gamma(\bar{\mathbf{v}}+\eta\bar{\mathbf{x}}+\bar{H})]$;

Obtain $\hat{V}^{\eta}_{k+1/2}(\mathbf{z}_{k+1/2})$ and
$\hat{V}^{\eta}_{k+1/2}(\bar{\mathbf{x}})$ as defined in eq. (3.7);

Draw iid direction vectors
$\mathbf{W}_{k+1/2}=\\{\mathbf{W}_{i,k+1/2}\\}_{i\in\mathcal{I}}$, with each
$\mathbf{W}_{i,k+1/2}\sim\mathsf{U}(\mathbb{S}_{i})$. ;

Obtain $H^{\delta}_{\mathbf{z}_{k+1/2}}(\mathbf{W}_{k+1/2})$ and
$H^{\delta}_{\bar{\mathbf{x}}}(\mathbf{W}_{k+1/2})$;

Update
$\mathbf{z}_{k+1}=\mathbf{z}_{k+1/2}-\gamma\left(\hat{V}^{\eta}_{k+1/2}(\mathbf{z}_{k+1/2})+H^{\delta}_{\mathbf{z}_{k+1/2}}(\mathbf{W}_{k+1/2})-\hat{V}^{\eta}_{k+1/2}(\bar{\mathbf{x}})-H^{\delta}_{\bar{\mathbf{x}}}(\mathbf{W}_{k+1/2})\right).$

end for

Algorithm 1
$\operatorname{\texttt{SFBF}}(\bar{\mathbf{x}},\bar{\mathbf{v}},\bar{H},\gamma,\eta,\delta,K)$

Data: $\mathbf{x},T,\\{\gamma_{t}\\}_{t=0}^{T},\\{b_{t}\\}_{t=0}^{T}$

Set $\mathbf{x}^{0}=\mathbf{x}$.

for  _$t=0,1,\ldots,T-1$_ do

For each $i\in\mathcal{I}$ receive the oracle feedback $\bar{V}^{t}$ defined
by
$\bar{V}^{t}_{i}\triangleq\frac{1}{b_{t}}\sum_{s=1}^{b_{t}}\hat{V}_{i}(\mathbf{x}^{t},\xi_{i,t}^{(s)})$.
;

For each $i\in\mathcal{I}$ construct the estimator
$H^{\delta_{t},b_{t}}_{\mathbf{x}^{t}}$ defined by
$H_{i,\mathbf{x}_{i}}^{\delta_{t},b_{t}}\triangleq\frac{1}{b_{t}}\sum_{s=1}^{b_{t}}H_{i,\mathbf{x}_{i}}^{\delta_{t}}(\mathbf{W}_{i,t}^{(s)})$.
;

Update
$\mathbf{x}^{t+1}=\operatorname{\texttt{SFBF}}(\mathbf{x}^{t},\bar{V}^{t},H^{\delta_{t},b_{t}}_{\mathbf{x}^{t}},\gamma_{t},\eta_{t},\delta_{t},K)$.

end for

Algorithm 2 Variance Reduced Hierarchical Game Solver
($\operatorname{\texttt{VRHGS}}$)

##### Discussion

The key innovation of the scheme $\operatorname{\texttt{VRHGS}}$ lies in the
combination of smoothing (to allow for hierarchy), regularization (to contend
with ill-posedness), and variance-reduction (to mitigate bias) within a
stochastic forward-backward-forward framework. Our double-loop solution
strategy mimics the computational architecture of SVRG, which takes a full
gradient sample of the finite sum problem “once in a while", while performing
frequent single-sample updates in between. Our method, adapted to general
probability spaces, proceeds similarly: The “shadow sequence"
$\mathbf{z}^{(t)}_{i,k+1/2}$ uses costly mini-batch estimators computed in the
outer loop; these are maintained in memory while executing the inner loop
(i.e. only “once in a while" updated). The additional forward steps to obtain
the iterates $\mathbf{z}^{(t)}_{i,k+1}$ make use of fresh one-shot estimators
of the payoff gradient and the finite difference estimator. All steps are
overlaid by a Tikhonov regularization, while smoothing facilitates
accommodation with hierarchical objectives. From a computational perspective,
our scheme performs a single projection onto the leaders’ feasible set
$\mathcal{X}_{i}$. This can save considerably on computational time in cases
where the projection operator is costly to evaluate, and constitutes a major
difference compared to viable alternative algorithmic schemes like the
extragradient or optimistic mirror descent. Hence, our method reduces the
sample-complexity of recent mini-batch variance reduction techniques for
stochastic VIs [8], while concomitantly reducing the computational bottlenecks
of double-call algorithms [35, 39] by lifting one projection step. Finally,
similar rate and complexity statements emerge when allowing for inexact
generalizations that allow for $\epsilon$-approximate solutions of lower-level
problem.

##  4 Main Results

In this section we state the main results on the asymptotic convergence of
scheme $\operatorname{\texttt{VRHGS}}$. All technical and lengthy proofs are
collected in Section 6.

In the inner and outer loops of $\operatorname{\texttt{VRHGS}}$, we have two
sources of randomness at each iteration: (i) the sequence of mini-batches
$\xi^{1:b_{t}}_{t}\triangleq\\{\xi^{1:b_{t}}_{i,t}\\}_{i\in\mathcal{I}}$ and
$\mathbf{W}^{1:b_{t}}_{t}\triangleq\\{\mathbf{W}_{i,t}^{1:b_{t}}\\}_{i\in\mathcal{I}}$,
which are used to perform the opening forward-backward step in Algorithm 1;
(ii) the sequences $\xi_{t,k+1/2}=\\{\xi_{i,t,k+1/2}\\}_{i\in\mathcal{I}}$ and
$\mathbf{W}_{t,k+1/2}=\\{\mathbf{W}_{i,t,k+1/2}\\}_{i\in\mathcal{I}}$ for
$k\in\\{0,1,\ldots,K-1\\}$, which are employed in constructing the iterate
$\mathbf{z}^{(t)}_{k+1}$ in Algorithm 1 during the outer epoch $t$. To keep
track of the information structure of the outer and inner loops, we introduce
the filtrations
$\mathcal{F}_{t}\triangleq\sigma(\mathbf{x}^{0},\ldots,\mathbf{x}^{t})$ for
$0\leq t\leq T$, and
$\mathcal{A}_{t,0}\triangleq\sigma(\mathbf{x}^{t},\xi_{t}^{1:b_{t}},\mathbf{W}_{t}^{1:b_{t}})$
as well as
$\mathcal{A}_{t,k}\triangleq\sigma(\mathcal{A}_{t,0}\cup\sigma(\xi_{t,1/2},\mathbf{W}_{t,1/2},\ldots,\xi_{t,k-1/2},\mathbf{W}_{t,k-1/2}))$
for $k\in\\{1,2,\ldots,K-1\\}$. By construction, the iterates
$\mathbf{z}^{(t)}_{k}$ and $\mathbf{z}^{(t)}_{k+1/2}$ are both
$\mathcal{A}_{t,k}$-measurable.

###  4.1 Error Structure of the estimators

We impose a uniform variance bound on the random vector field $\hat{V}$ over
the set $\mathcal{X}$. Compactness of $\mathcal{X}$ implies that such an
assumption comes without loss of generality, and the proof of the variance
bound in Lemma 4.1 is simple to obtain and thus omitted; See [8].

######  Lemma 4.1.

There exists $M_{V}>0$ such that
$\mathbb{E}_{\xi}[\norm{\hat{V}(\mathbf{x},\xi)-V(\mathbf{x})}^{2}]\leq
M^{2}_{V}$ for all $\mathbf{x}\in\mathcal{X}$. Additionally, let $b_{t}\geq 1$
and
$\xi^{1:b_{t}}_{t}=\\{\\{\xi^{(1)}_{i,t}\\}_{i\in\mathcal{I}},\ldots,\\{\xi^{(b_{t})}_{i,t}\\}_{i\in\mathcal{I}}\\}$
denote an i.i.d sample of the random variable $\xi$. Then, for
$\varepsilon_{V}(\xi^{1:b_{t}}_{t})\triangleq\frac{1}{b_{t}}\sum_{s=1}^{b_{t}}\hat{V}(\mathbf{x},\xi^{(s)}_{t})-V(\mathbf{x})$,
we have
$\sqrt{\mathbb{E}\left[\norm{\varepsilon_{V}(\xi^{1:b_{t}}_{t})}^{2}\right]}\leq\frac{M_{V}}{\sqrt{b_{t}}}.$

Concerning the estimator of the gradient of the smoothed lower level function
$h^{\delta}(\cdot)$, we can report the following bounds, which are derived in
Appendix A.5 and proved in Lemma A.12.

######  Lemma 4.2.

Let Assumptions 1-2 hold. Define
$\mathtt{e}_{\mathbf{x}_{i}}(\mathbf{W}_{i}^{1:b})\triangleq
H^{\delta,b}_{i,\mathbf{x}_{i}}-\nabla h^{\delta}_{i}(\mathbf{x}_{i})$, where
$\\{\mathbf{W}^{(s)}_{i}\\}_{s=1}^{b_{t}}$ is an i.i.d sample drawn uniformly
from the unit sphere $\mathbb{S}_{i}$, i.e.
$\mathbf{W}^{1:b}_{i}\sim\mathsf{U}(\mathbb{S}_{i})^{\otimes b}$. Then for all
$i\in\mathcal{I}$,

  * (a)

$\mathbb{E}_{\mathbf{W}^{1:b}_{i}\sim\mathsf{U}(\mathbb{S}_{i})^{\otimes
b}}[\mathtt{e}_{\mathbf{x}_{i}}(\mathbf{W}^{1:b}_{i})]=0$;

  * (b)

$\norm{H_{i,\mathbf{x}_{i}}^{\delta}(\mathbf{w}_{i})}^{2}\leq
L_{h_{i}}^{2}n^{2}_{i}$ for all $\mathbf{w}_{i}\in\mathbb{S}_{i}$;

  * (c)

$\mathbb{E}_{\mathbf{W}^{1:b}_{i}\sim\mathsf{U}(\mathbb{S}_{i})^{\otimes
b}}[\norm{\mathtt{e}_{\mathbf{x}_{i}}(\mathbf{W}^{1:b}_{i})}^{2}]\leq\frac{n^{2}_{i}L^{2}_{h_{i}}}{b_{i}}.$

###  4.2 Almost sure convergence of the last iterate

Our analysis of $\operatorname{\texttt{VRHGS}}$ relies on the following energy
inequality, proved in Section 6.1.1.

######  Lemma 4.3.

Let Assumptions 1,2,3 and 4 hold true. Let $\\{\mathbf{x}^{t}\\}_{t=0}^{T-1}$
be generated by $\operatorname{\texttt{VRHGS}}$ and denote by
$\\{\mathbf{z}^{(t)}_{k}\\}_{k\in\\{0,1/2,\ldots,K\\}}$ the sequence obtained
by executing
$\operatorname{\texttt{SFBF}}(\mathbf{x}^{t},\bar{V}^{t},H^{\delta_{t},b_{t}}_{\mathbf{x}^{t}},\gamma_{t},\eta_{t},\delta_{t},K)$.
Set $L_{h}\triangleq\sum_{i\in\mathcal{I}}L_{h_{i}}$. Then for all
$\mathbf{x}\in\mathcal{X}$ and $k\in\\{0,1,\ldots,K-1\\}$, we have

$\displaystyle\norm{\mathbf{z}^{(t)}_{k+1}-\mathbf{x}}^{2}$
$\displaystyle\leq(1-\gamma_{t}\eta_{t})\norm{\mathbf{z}^{(t)}_{k}-\mathbf{x}}^{2}-(1-2\gamma_{t}\eta_{t})\norm{\mathbf{z}^{(t)}_{k+1/2}-\mathbf{z}^{(t)}_{k}}^{2}$
$\displaystyle+8\gamma_{t}^{2}\sum_{i\in\mathcal{I}}L_{h_{i}}^{2}n^{2}_{i}+4\gamma^{2}_{t}(\mathcal{L}_{f}(\xi_{t,k+1/2})^{2}+\eta^{2}_{0})\norm{\mathbf{z}^{(t)}_{k+1/2}-\mathbf{x}^{t}}^{2}$
$\displaystyle-2\gamma_{t}\langle\hat{V}_{t,k+1/2}(\mathbf{z}^{(t)}_{k+1/2})+H^{\delta_{t}}_{\mathbf{z}^{(t)}_{k+1/2}}(\mathbf{W}_{t,k+1/2})-V(\mathbf{z}^{(t)}_{k+1/2})-\nabla
h^{\delta_{t}}(\mathbf{z}^{(t)}_{k+1/2}),\mathbf{z}^{(t)}_{k+1/2}-\mathbf{x}\rangle$
$\displaystyle-2\gamma_{t}\langle\left(V(\mathbf{x}^{t})+\nabla
h^{\delta_{t}}(\mathbf{x}^{t})\right)-\left(\hat{V}_{t,k+1/2}(\mathbf{x}^{t})+H^{\delta_{t}}_{\mathbf{x}^{t}}(\mathbf{W}_{t,k+1/2})\right),\mathbf{z}^{(t)}_{k+1/2}-\mathbf{x}\rangle$
$\displaystyle-2\gamma_{t}\langle\varepsilon^{t}_{V}(\xi^{1:b_{t}}_{t})+\varepsilon^{t}_{h}(\mathbf{W}^{1:b_{t}}_{t}),\mathbf{z}^{(t)}_{k+1/2}-\mathbf{x}\rangle$
$\displaystyle-2\gamma_{t}\left(\langle
V^{\eta_{t}}(\mathbf{x}),\mathbf{z}^{(t)}_{k+1/2}-\mathbf{x}\rangle+h(\mathbf{z}^{(t)}_{k+1/2})-h(\mathbf{x})\right)+2\gamma_{t}\delta_{t}L_{h}.$

We next prove a.s. convergence of $\\{\mathbf{x}^{t}\\}_{t=0}^{T}$ to the
least-norm solution of $\operatorname{MVI}(V,h)$ as $T\to\infty$. The proof
rests on a fine comparison between the algorithmic sequence
$\\{\mathbf{x}^{t}\\}$ and the sequence of solutions of the regularized
problems $\operatorname{MVI}(V^{\eta_{t}},h)$, denoted as
$\\{\mathbf{s}_{t}\\}_{t=0}^{T}$.

######  Theorem 4.4.

Let Assumptions 1,2,3, and 4 hold. Suppose we are given sequences
$\\{\gamma_{t}\\}_{t\in\mathbb{N}},\\{\delta_{t}\\}_{t\in\mathbb{N}}$ and
$\\{\eta_{t}\\}_{t\in\mathbb{N}}$, satisfying the following conditions:

  * (a)

$\lim_{t\to\infty}\frac{\gamma_{t}}{\eta_{t}}=\lim_{t\to\infty}\frac{\delta_{t}}{\eta_{t}}=0$,
and
$\sum_{t=0}^{\infty}\gamma^{2}_{t}<\infty,\sum_{t=0}^{\infty}\gamma_{t}\eta_{t}=\infty$;

  * (b)

$\gamma_{t}\eta_{t}\in(0,1/2)$ and $\lim_{t\to\infty}\eta_{t}=0$;

  * (c)

$\sum_{t=0}^{\infty}\left(\frac{\eta_{t}-\eta_{t-1}}{\eta_{t}}\right)^{2}(1+\frac{1}{\gamma_{t}\eta_{t}})<\infty$
and
$\lim_{t\to\infty}\left(\frac{\eta_{t}-\eta_{t-1}}{\eta_{t}}\right)^{2}\left(\frac{1+\frac{1}{\gamma_{t}\eta_{t}}}{\gamma_{t}\eta_{t}}\right)=0$.

Then
$\mathbb{P}(\lim_{t\to\infty}\norm{\mathbf{x}^{t}-\mathbf{x}^{\ast}}=0)=1$,
where $\mathbf{x}^{\ast}$ denotes the unique least norm solution of
$\operatorname{MVI}(V,h)$; i.e.
$\\{(\mathbf{x}^{t},\mathbf{y}(\mathbf{x}^{t}))\\}_{t\in\mathbb{N}}$ converges
almost surely to a hierarchical equilibrium of the game (P).

The proof of this Theorem can be found in Section 6.1.2.

###  4.3 Finite-time complexity

The convergence measure usually employed for MVI$(V,h)$ is the gap function

$\Gamma(\mathbf{x})\triangleq\sup_{\mathbf{z}\in\mathcal{X}}\left(\langle
V(\mathbf{z}),\mathbf{x}-\mathbf{z}\rangle+h(\mathbf{x})-h(\mathbf{z})\right),$
(4.1)

Since we work in probabilistic setting, naturally our convergence measure will
be based on $\mathbb{E}[\Gamma(\mathbf{x})]$. Our main finite-time iteration
complexity result in terms of this performance measure is the next Theorem,
whose proof is detailed in Section 6.1.3.

######  Theorem 4.5.

Let Assumptions 1,2,3, and 4 hold and fix $T\in\mathbb{N}$. Consider Algorithm
$\operatorname{\texttt{VRHGS}}$ with the inputs
$\gamma_{t}=\eta_{t}=\delta_{t}=1/T$, as well as $K=T$ and batch size
$b_{t}\geq T^{2}$. Then,
$\mathbb{E}[\Gamma(\bar{\mathbf{z}}^{T})]=\mathcal{O}\left(\frac{C\sigma}{T}\right),$
where
$\sigma\triangleq\sqrt{2M_{V}^{2}+2\sum_{i\in\mathcal{I}}L^{2}_{h_{i}}n^{2}_{i}},C=\max_{i\in\mathcal{I}}C_{i}$
(cf. Assumption (2.i)), and
$\bar{\mathbf{z}}^{T}\triangleq\frac{\sum_{t=0}^{T-1}\gamma_{t}\bar{\mathbf{z}}^{t}}{\sum_{t=0}^{T-1}\gamma_{t}}$
for
$\bar{\mathbf{z}}^{t}\triangleq\frac{1}{K}\sum_{k=0}^{K-1}\mathbf{z}^{(t)}_{k+1/2}$.

We next evaluate the oracle complexity of $\operatorname{\texttt{VRHGS}}$. To
be precise, let $\operatorname{\mathsf{OC}}(T,K,\\{b_{t}\\}_{t=0}^{T-1})$ the
number of random variables method $\operatorname{\texttt{VRHGS}}$ generates in
the inner and outer loop until we achieve a solution that pushes the expected
gap below a target value $\varepsilon$.

######  Remark 4.1.

We point out that this measure of oracle complexity ignores the computational
effort arising from solving the lower level problem attached with player
$i\in\mathcal{I}$. This is consistent as we assume that the solution map is
provided to us in terms of an oracle. A full-fledged complexity analysis can
be done, and will appear in a future publication.

######  Proposition 4.6.

Let $\varepsilon>0$ be given, and set $T=\lceil 1/\varepsilon\rceil$. If we
choose the same sequences as in Theorem 4.5, we have
$\operatorname{\mathsf{OC}}(T,K,\\{b_{t}\\}_{t=0}^{T-1})=\mathcal{O}(2N/\varepsilon^{3}).$

###### Proof.

The number of random variables generated in each inner loop iteration is
$2K\times N$. In each round of the outer loop we sample $2b_{t}\times N$
random variables. Hence, the total oracle complexity is
$\operatorname{\mathsf{OC}}(T,K,\\{b_{t}\\}_{t=0}^{T})=2KTN+2N\sum_{t=0}^{T-1}b_{t}.$
For the specific values of $T,K,\gamma_{t},\delta_{t},\eta_{t}$ defined in
Theorem 4.5 and $b_{t}=T^{2}$, it clearly follows
$\operatorname{\mathsf{OC}}(T,K,\\{b_{t}\\}_{t=0}^{T-1})=\mathcal{O}(2N/\varepsilon^{3}).$
$\blacksquare$

###  4.4 Inexact generalization

A key shortcoming in the implementation of $\operatorname{\texttt{VRHGS}}$ is
the need for exact solutions of the lower-level problem. Naturally, when the
solution map $\mathbf{y}_{i}(\cdot)$ corresponds to the solution of a large-
scale stochastic optimization/VI problem, this claim is hard to justify. In
this section, we allow for an inexact solution
$\mathbf{y}_{i}^{\varepsilon}(\mathbf{x})$ associated with an error level
$\varepsilon$, defined as

$\mathbb{E}[\norm{\mathbf{y}^{\varepsilon}(\mathbf{x})-\mathbf{y}(\mathbf{x})}|\mathbf{x}]\leq\varepsilon\quad\text{
a.s.}$ (4.2)

Under the inexact lower level solution $\mathbf{y}_{i}^{\varepsilon}$, we let
$h_{i}^{\varepsilon}(\mathbf{x}_{i})\triangleq
g_{i}(\mathbf{x}_{i},\mathbf{y}_{i}^{\varepsilon}(\mathbf{x}_{i}))$.

######  Remark 4.2.

We can obtain the inexact solution
$\mathbf{y}_{i}^{\varepsilon}(\mathbf{x}_{i})$ with rather efficient numerical
methods. First, we can parallelize the computation since the problems
$\operatorname{VI}(\phi_{i}(\mathbf{x}_{i},\cdot),\mathcal{Y}_{i})$ are
uncoupled. Second, the mapping $\phi_{i}(\mathbf{x}_{i},\cdot)$ is assumed to
be strongly monotone. Hence, we can solve the VI to $\varepsilon$-accuracy
with exponential rate using for instance the method in [12].

As in the exact regime, we assume that player $i$ has access to an oracle with
which she can construct a spherical approximation of the gradient of the
implicit function $h_{i}^{\varepsilon}$. Hence, for given $\delta>0$, we let
$h^{\varepsilon,\delta}_{i}(\mathbf{x}_{i})\triangleq\int_{\mathbb{B}_{n_{i}}}h^{\varepsilon}_{i}(\mathbf{x}_{i}+\delta\mathbf{w})\frac{\differential\mathbf{w}}{{\mathbf{Vol}}_{n}(\mathbb{B}_{n_{i}})}.$
We denote the resulting estimators by
$H^{\delta,\varepsilon}_{i,\mathbf{x}_{i}}(\mathbf{W}_{i})=n_{i}\mathbf{W}_{i}\nabla_{(\mathbf{W}_{i},\delta)}h^{\varepsilon}_{i}(i,\mathbf{x}_{i})$,
and the mini-batch versions
$H^{\delta,\varepsilon,b}_{i,\mathbf{x}_{i}}\triangleq\frac{1}{b}\sum_{s=1}^{b}H^{\delta,\varepsilon}_{i,\mathbf{x}_{i}}(\mathbf{W}_{i}^{(s)}).$
With these concepts in hand, we can adapt $\operatorname{\texttt{VRHGS}}$ to
run exactly the same way as described in Algorithm 1 and Algorithm 2,
replacing the appearance of quantities involving $h_{i}$ with its inexact
version $h^{\varepsilon}_{i}$; see Section 6.2 for a precise formulation of
the method.

######  Theorem 4.7.

Let Assumptions 2 hold and fix $T\in\mathbb{N}$. Consider Algorithm
I-$\operatorname{\texttt{VRHGS}}$, defined in Section 6.2, with the sequence
$\gamma_{t}=\eta_{t}=\delta_{t}=1/T$, as well as $K=T$, the batch size
$b_{t}\geq T^{2}$ and inexactness regime $\varepsilon_{t}=1/T^{2}$. Then
$\mathbb{E}[\Gamma(\bar{\mathbf{z}}^{T})]=\mathcal{O}\left(\frac{C\sigma}{T}\right),$
where
$\sigma\triangleq\sqrt{2M_{V}^{2}+2\sum_{i\in\mathcal{I}}L^{2}_{h_{i}}n^{2}_{i}},C=\max_{i\in\mathcal{I}}C_{i}$
(cf. Assumption (2.i)).

Notably, tractable resolution of the proposed stochastic hierarchical game is
possible in inexact regimes and such practically motivated schemes are not
adversely affected in terms of either the rate or complexity guarantees.

##  5 Hierarchical games in power markets

In this section we present a model inspired by Hobbs and Pang [32], but
suitably modified to account for uncertainty in prices and costs, multi-period
settings with ramping constraints, and the incorporation of virtual power
plants (VPPs) (see [51, 19] for a review of VPPs and power markets). The model
we present below is at this stage an academic example that demonstrates the
modelling power of our hierarchical games approach. In future studies we aim
for numerical implementations of this model.

Consider a set of nodes ${\cal N}$ of a network and a set of time periods
${\cal T}\,\triangleq\,\left\\{\,1,2,\cdots,T\,\right\\}$. A generation firm
is indexed by $f$, where $f$ belongs to the finite set ${\cal F}$ and each
firm is assumed to have an associated VPP. At a node $i$ in the network, a
firm $f$ may generate $g_{f,i,t}$ units via conventional generation in period
$t$ and sell $s_{f,i,t}$ units during the same period. In addition, at time
period $t$, firm $f$ may generate $P^{\rm pv,S}_{f,t}+P^{\rm pv,L}_{f,t}$
units of power via PV capacity, of which $P^{\rm pv,S}_{f,t}$ is sold and
$P^{\rm pv,L}_{f,t}$ is employed for meeting load. The total amount of power
sold at node $i$ during period $t$ by all generating firms is represented by
$S_{i,t}$, i.e. $S_{i,t}={\displaystyle\sum_{f\in{\cal F}}}s_{f,i,t}$. If the
nodal power price at the $i$th node during period $t$ is a random function
given by $p_{i,t}(\bullet,\xi)$, where $p_{i,t}(\bullet,\xi)$ is a decreasing
function of aggregate nodal sales $S_{i,t}$ for any $\xi\in\Xi$. It follows
that firm $f$’s revenue from non-PV power sales at node $i$ during period $t$
under realization $\xi$ is $p_{i,t}(S_{i,t},\xi)s_{f,i}$. Sales of PV output
by firm $f$ at time $t$ is priced using a function $p^{\rm R}_{f,t}$, earning
a revenue given by $p_{f,t}^{\rm R}\left(P_{f,t}^{\rm pv,S}+P_{f,t}^{\rm
pv,S,V}\right)P^{\rm pv,S}$, where $P_{f,t}^{\rm pv,S,V}$ denotes the sales of
firm $f$’s associated VPP (whose problem is described later in this section).
We observe that renewable power is priced using this price function, distinct
from conventional sources, and is designed to provide incentives for renewable
expansion [60]. The costs incurred by firm $f$ at node $i$ during period $t$
are given by the sum of the cost of generating $g_{f,i,t}$ and the cost of
transmitting the excess $(s_{f,i,t}-g_{f,i,t})$. Let the random cost function
of generation associated with firm $f$ at node $i$ be given by
$c_{fi}(\bullet,\zeta)$ while the cost of transmitting power from an arbitrary
node (referred to as the hub) to node $i$ is given by $w_{i}$. The constraint
set incorporates a balance between aggregate sales, aggregate generation, and
power injection into the VPP at all nodes for every time period $t$. In
addition, we impose nonnegativity bound on sales and generation at any time
period $t$, enforce a capacity limit on generation levels, and introduce
ramping constraints on the change in generation levels. The resulting problem
faced by generating firm $f$, denoted by (Firmf), requires minimizing
generation cost less revenue from conventional and PV sales by optimizing
sales $s_{f,i,t}$ and generation $g_{f,i,t}$ at every node $i$ and every time
period $t$ as well as load-directed PV output $P^{\rm pv,L}_{f,t}$ and PV
sales $P^{\rm pv,S}_{f,t}$ at time $t$. If $P_{f,t}^{\rm
pv,S,V,\epsilon}(\cdot)$ denotes a component of the single-valued solution map
of the $\epsilon$-regularized problem of the VPP associated with firm $f$,
denoted by (VPP${}_{f}(P_{f}^{\rm pv,S})$), then firm $f$’s problem is defined
as follows, where $\widehat{{\cal T}}=\\{1,\cdots,T-1\\}$.

$\begin{array}[]{ll}\displaystyle{{\operatornamewithlimits{\mbox{\bf
maximize}}_{s_{f,i,t},\,g_{f,i,t},\,P^{\rm pv,S}_{f,t}\,P^{\rm
pv,L}_{f,t}}}}&\mathbb{E}\left[\displaystyle{\sum_{t=1}^{T}\sum_{i\in{\cal
N}}}\,\left(\,p_{i,t}(S_{i,t},\xi)s_{f,i,t}-c_{fi}(g_{fi,t},\zeta)-(s_{fi,t}-g_{fi,t})w_{i,t}\,\right)\right]\\\
&+{\displaystyle\sum_{t\in{\cal T}}}\left(p_{f,t}^{\rm R}\left(P_{f,t}^{\rm
pv,S}+P_{f,t}^{\rm pv,S,V,\epsilon}(P_{f,t}^{\rm pv,S})\right)P^{\rm
pv,S}\right)\\\\[14.45377pt] \mbox{\bf subject
to}&\left\\{\begin{array}[]{lll}0&\leq&g_{f,i,t}\,\leq\,{\rm
cap}_{fi}\\\\[5.0pt] 0&\leq&s_{f,i,t}\end{array}\right\\},\hskip
65.04256pt\forall\,t\in{\cal T},\quad\forall i\in{\cal N}\hskip
50.58878pt\left(\mbox{\bf Firm}_{f}(s_{-f},g_{-f})\right)\\\\[14.45377pt]
&-\mbox{RR}^{\rm
down}_{f,i}\,\leq\,g_{f,i,t}-g_{f,i,t-1}\,\leq\,\mbox{RR}^{\rm
up}_{fi},\quad\forall\,t\in\widehat{{\cal T}},\quad\forall i\in{\cal
N}\\\\[14.45377pt] &P^{\rm pv,L}_{f,t}+P^{\rm pv,S}_{f,t}\,\leq\,{\rm
cap}^{\rm pv}_{f,t},\ \qquad\qquad\qquad\forall t\in{\cal T}\\\
&0\,\leq\,P^{\rm pv,L}_{f,t},P^{\rm
pv,S}_{f,t},\qquad\qquad\qquad\qquad\quad\forall t\in{\cal T}\\\\[5.0pt]
\mbox{and}&\displaystyle{\sum_{i\in{\cal
N}}}\,(\,s_{f,i,t}-g_{f,i,t}\,)-P^{\rm pv,L}_{f,t}\,=\,0.\ \
\quad\qquad\forall t\in\ {\cal T}\end{array}$

It bears reminding that the last set of constraints specified in
$\left(\mbox{\bf Firm}_{f}(s_{-f},g_{-f})\right)$ are parametrized by rival
decisions and can be relaxed with Lagrange multiplier $\lambda_{f,t}$, leading
to the following relaxed problem $\left(\mbox{\bf Firm}^{\rm
rel}_{f}(s_{-f},g_{-f})\right)$, defined as follows.

$\begin{array}[]{ll}\displaystyle{{\operatornamewithlimits{\mbox{\bf
maximize}}_{s_{f,i,t},\,g_{f,i,t},\,P^{\rm pv,S}_{f,t}\,P^{\rm
pv,L}_{f,t}}}}&\mathbb{E}\left[\displaystyle{\sum_{t=1}^{T}\sum_{i\in{\cal
N}}}\,\left(\,p_{i,t}(S_{i,t},\xi)s_{f,i,t}-c_{f,i}(g_{f,i,t},\zeta)-(s_{f,i,t}-g_{f,i,t})w_{i,t}\,\right)\right]\\\
&+{\displaystyle\sum_{t\in{\cal T}}}\left(p_{f,t}^{\rm R}\left(P_{f,t}^{\rm
pv,S}+P_{f,t}^{\rm pv,S,V}\right)P_{f,t}^{\rm
pv,S}\right)-{\displaystyle\sum_{t=1}^{T}}\lambda_{f,t}^{\top}\left({\displaystyle\sum_{i\in{\cal
N}}}\left((\,s_{f,i,t}-g_{f,i,t}\,)-P_{f,t}^{\rm
pv,L}\right)\right)\\\\[14.45377pt] \mbox{\bf subject
to}&\left\\{\begin{array}[]{lll}0&\leq&g_{fi,t}\,\leq\,{\rm
cap}_{fi}\\\\[5.0pt] 0&\leq&s_{fi,t}\end{array}\right\\},\hskip
65.04256pt\forall\,t\in{\cal T},\quad\forall i\in{\cal N}\hskip
50.58878pt\left(\mbox{\bf Firm}_{f}(s_{-f},g_{-f})\right)\\\\[14.45377pt]
&-\mbox{RR}^{\rm
down}_{f,i}\,\leq\,g_{f,i,t}-g_{f,i,t-1}\,\leq\,\mbox{RR}^{\rm
up}_{f,i},\quad\forall\,t\in\widehat{{\cal T}},\quad\forall i\in{\cal
N}\\\\[14.45377pt] &P^{\rm pv,L}_{f,t}+P^{\rm pv,S}_{f,t}\,\leq\,{\rm
cap}^{\rm pv}_{f,t},\ \qquad\qquad\qquad\forall t\in{\cal T}\\\
&0\,\leq\,P^{\rm pv,L}_{f,t},P^{\rm
pv,S}_{f,t},\qquad\qquad\qquad\qquad\quad\forall t\in{\cal T}\end{array}$

In addition, we introduce a pricing player $\left(\mbox{\bf
Price}\left(s_{f},g_{f},P_{f,t}^{\rm pv,L}\right)\right)$ corresponding to the
determination of $\lambda_{f,t}$ for $f\in{\cal F}$ and $t\,\in\,{\cal T}$,
defined as follows.

$\begin{array}[]{ll}\displaystyle{{\operatornamewithlimits{\mbox{\bf
minimize}}_{\lambda}}}&{\displaystyle\sum_{f\in{\cal F}}\sum_{t\in{\cal
T}}}\lambda_{f,t}^{\top}\left({\displaystyle\sum_{i\in{\cal
N}}}\,(\,s_{f,i,t}-g_{f,i,t}\,)+P_{f,t}^{\rm pv,L}\right).\hskip
72.26999pt\left(\mbox{\bf Price}\left(s,g,P^{\rm
pv,L}\right)\right)\end{array}$

Note that, the generating firm sees the transmission fee $w_{i,t}$ and the
rival firms’ sales $s_{-fi,t}\equiv\\{s_{hi,t}~{}:~{}h\neq f\\}$ as exogenous
parameters to its optimization problem even though they are endogenous to the
overall equilibrium model as we will see shortly. The ISO sees the
transmission fees $w=(w_{i,t})_{i\in{\cal N},t\in{\cal T}}$ as exogenous and
prescribes flows $y=(y_{i,t})_{i\in{\cal N},t\in{\cal T}}$ as per a solution
of the following linear program

$\begin{array}[]{ll}\displaystyle{{\operatornamewithlimits{\mbox{\bf
maximize}}_{y}}}&\displaystyle{\sum_{i\in{\cal N}}\sum_{t\in{\cal
T}}}\,y_{i,t}w_{i,t}\\\\[14.45377pt] \mbox{\bf subject
to}&\displaystyle{\sum_{i\in{\cal N}}}\,{\rm
PDF}_{ij}y_{i,t}\,\leq\,\hat{T}_{j},\qquad\forall\,j\in{\cal K},\
\forall\,t\in{\cal T},\hskip 93.95122pt\left(\mbox{\bf
ISO}(w)\right)\end{array}$

where ${\cal K}$ is the set of all arcs or links in the network with node set
${\cal N}$, $\hat{T}_{j}$ denotes the transmission capacity of link $j$,
$y_{i,t}$ represents the transfer of power (in MW) by the system operator from
a hub node to node node $i$ and PDFij denotes the power transfer distribution
factor, which specifies the MW flow through link $j$ as a consequence of unit
MW injection at an arbitrary hub node and a unit withdrawal at node $i$.
Finally, to clear the market, the transmission flows $y_{i}$ must must balance
the net sales at each node, as specified next.

$\displaystyle y_{i,t}\,=\,\displaystyle{\sum_{f\in{\cal
F}}}\,\left(\,s_{f,i,t}-g_{f,i,t}\,\right),\qquad\forall\,i\in{\cal
N},\quad\forall\ t\in{\cal T}.$ (5.1)

In fact, this constraint can be recast as a collection of pricing players,
denoted by (Flow${}^{\rm price}(g,s,y)$).

$\begin{array}[]{ll}\displaystyle{{\operatornamewithlimits{\mbox{\bf
minimize}}_{\beta}}}&{\displaystyle\sum_{t\in{\cal T}}\sum_{i\in{\cal
N}}}\beta_{i,t}\left(\,y_{i,t}\,-\,\displaystyle{\sum_{f\in{\cal
F}}}\,\left(\,s_{f,i,t}-g_{f,i,t}\,\right)\,\right)\hskip 144.54pt(\mbox{\bf
Flow}^{\rm price}(g,s,y))\end{array}$

We now extend the scope of the framework of power markets by incorporating
virtual power plants. A virtual power plant (VPP) represents a collections of
distributed energy resources (DERs) (e.g., batteries, smart thermostats,
controllable water heaters, and rooftop solar) that can be coordinated to
enhance the reliability and sustainability of the electric grid. To satisfy
the short-term goals for clean energy technology (CET) deployment, it has been
estimated that U.S. VPP capacity must triple by 2030, leading to potential
savings of $10 billion in annual grid costs [19]. Without loss of generality,
we assume that any firm $f\,\in\,{\cal F}$ has a collection of components,
which collectively provide “virtual power” in addition to conventional
generation. Before proceeding, we model three components in such a VPP, akin
to approaches employed in [36, 25, 9].

(a) Battery storage. Suppose the storage unit associated with firm $f$ has an
associated state of charge (SOC) level at time $t$ by SOCf,t.

$\displaystyle\mbox{SOC}_{f,t+1}\,=\,\mbox{SOC}_{f,t}\,+\,\tfrac{\eta^{\rm
b,ch}_{f}\Delta t}{Q^{\mathrm{b}}_{f}}P^{\rm b,ch}_{f,t}\,-\,\tfrac{\Delta
t}{\eta^{\rm b,ds}_{f}Q^{\mathrm{b}}_{f}}P^{\rm b,ds}_{f,t},\
\forall\,t\in{\cal T}$ (5.2)

where $P^{\rm b,ch}_{f,t}$ and $P^{\rm b,ds}_{f,t}$ represent charging and
discharging power-levels at time $t$, $\eta^{\rm b,ch}_{f}$ and $\eta^{\rm
b,ds}_{f}$ represent charging and discharging efficiencies at time $t$, while
$Q^{\rm b}_{f}$ and $\Delta t$ denote the battery capacity and time interval,
respectively. In addition, SOCf,t is bounded between a minimum value
$\mbox{SOC}^{\rm\min}_{f}$ and maximum value $\mbox{SOC}^{\rm\max}_{f}$ while
at any time $t$, charging and discharging rates cannot exceed $P^{\rm
b,ch,mx}_{f}$ and $P^{\rm b,ds,mx}_{f}$, respectively, as captured by the
following bounds.

$\displaystyle\mbox{SOC}^{\rm\min}_{f}\,$
$\displaystyle\leq\,\mbox{SOC}_{f,t}\,\leq\,\mbox{SOC}^{\rm\max}_{f},\qquad\forall
t\in{\cal T}$ (5.3) $\displaystyle 0\,$ $\displaystyle\leq\,\tfrac{P^{\rm
b,ch}_{f,t}}{P^{\rm b,ch,mx}_{f}}\,\leq\,1,\qquad\qquad\quad\forall t\in{\cal
T}$ (5.4) $\displaystyle 0\,$ $\displaystyle\leq\,\tfrac{P^{\rm
b,ch}_{f,t}}{P^{\rm b,ds,mx}_{f}}\,\,\leq\,1.\qquad\qquad\quad\forall
t\in{\cal T}$ (5.5)

(b) Intermittent resources. We now model intermittency by considering a
photovoltaic (PV) array associated with firm $f$, where at time $t$, $P^{\rm
pv,L,V}_{f,t}$ and $P^{\rm pv,S,V}_{f,t}$ denote the PV output employed for
meeting load and for deriving sales revenue, respectively. Further, $P^{\rm
pv,max}_{f}$ represents maximum PV power at time $t$. Consequently, PV output
is modeled as

$\displaystyle P^{\rm pv,L,V}_{f,t}+P^{\rm pv,S,V}_{f,t}\,$
$\displaystyle=\,(1-U^{\rm pv}_{f})P^{\rm
pv,max}_{f,t}E_{f,t}\qquad\qquad\forall t\in{\cal T}$ (5.6) $\displaystyle
0\,$ $\displaystyle\leq\,U_{f,t}^{\rm
pv}\,\leq\,1,\qquad\qquad\qquad\qquad\,\forall t\in{\cal T}$ (5.7)
$\displaystyle 0\,$ $\displaystyle\leq\,P^{\rm pv,L,V},P^{\rm
pv,S,V}\qquad\,\forall t\in{\cal T}$ (5.8)

where $U^{\rm pv}_{f,t}$ denotes the PV curtailment employed by firm $f$ at
time $t$ while $P^{\rm pv,max}_{f,t}$ scales with the solar irradiance at time
$t$ as seen by firm $f$, denoted by $E_{f,t}$. We observe that $P^{\rm
pv,S,V}_{f,t},P^{\rm pv,L,V}_{f,t}\,\geq\,0$ for any $f\,\in\,{\cal F}$ and
any $t\,\in\,{\cal T}.$

(c) Thermal onsite generation. Often VPPs may incorporate onsite thermal
generation that can be employed. For any $f\,\in\,{\cal F}$, suppose the
generation capacity is denoted by $\mbox{Cap}^{\rm onsite}_{f}$ while the
upward and downward ramping rates are given by $\mbox{RR}^{\rm up}_{f}$ and
$\mbox{RR}^{\rm down}_{f}$, respectively. Consequently, if the generation
output at time $t$ is denoted by $P^{\rm onsite}_{f,t}$, then for any
$t\,\in\,{\cal T}$, we have

$\displaystyle 0\,\leq\,P^{\rm onsite}_{f,t}\,\leq\,\mbox{Cap}^{\rm
onsite}_{f},\qquad\forall\,t\,\in\,{\cal T}$ (5.9)

Furthermore, changes in generation level are bounded by ramping rates, as
captured by the following set of two-sided constraints.

$\displaystyle-\mbox{RR}^{\rm down}_{f}\,\leq\,P^{\rm onsite}_{f,t+1}-P^{\rm
onsite}_{f,t+1}\,\leq\,\mbox{RR}^{\rm up}_{f},\quad\forall\,t\,\in\,{\cal T}$
(5.10)

VPPs are characterized by an idiosyncratic load profile that cannot be
controlled; specifically, $P^{\rm L,V}_{f,t}$ denotes the load associated with
VPP $f$ at time $t$. In more comprehensive models, we may incorporate HVAC and
water heater components that allow for more fine-grained control of such loads
but for purposes of simplicity, we omit such a discussion here. In the current
setting, the effective load emerging from managing the VPP associated with
firm $f$ and time $t$ is given by the sum of the uncontrollable load and the
battery load (charging less discharging level) less the sum of onsite
generation and load-directed PV output is required to be nonpositive, as
specified next.

$\displaystyle P^{\rm L,V}_{f,t}+\left(\,P^{\rm b,ch}_{f,t}-P^{\rm
b,ds}_{f,t}\,\right)-P_{f,t}^{\rm onsite}-P^{\rm
pv,L,V}_{f,t}\,\leq\,0,\qquad\forall\,t\,\in\,{\cal T}$ (5.11)

Note that the satisfaction of this constraint relies on appropriate sizing of
the battery capacity $Q_{f}^{b}$ and the onsite generation capacity
$\mbox{Cap}_{f}^{\rm onsite}$. Suppose the decision vector of firm $f$’s VPP
is denoted by $y^{\rm vpp}_{f}$, defined as

$y_{f}^{\rm vpp}\,=\,\left(\,\mbox{SOC}_{f};P^{\rm b,ch}_{f};P^{\rm
b,ds}_{f};P^{\rm pv,L,V}_{f};P^{\rm pv,S,V}_{f};U_{f}^{\rm pv};P^{\rm
onsite}_{f}\,\right).$

The profit function associated with firm $f$’s VPP is the revenue obtained by
sales revenue derived from PV sales less the VPP’s operational cost (given by
the sum of the costs of onsite generation and the (converted) cost of PV
curtailment), defined as

$\displaystyle\mbox{r}^{\rm vpp}_{f}(y^{\rm vpp}_{f};P^{\rm
pv,S}_{f})\,\triangleq\,\displaystyle{\sum_{t\in{\cal T}}}$
$\displaystyle\,\left(\,\underbrace{p_{f,t}^{\rm R}\left(P^{\rm
pv,S,V}_{f,t}+P^{\rm pv,S}_{f,t}\right)P^{\rm pv,S,V}_{f,t}}_{\tiny\mbox{VPP
revenue from PV sales}}-\underbrace{c_{f}^{\rm onsite}(P^{\rm
onsite}_{f,t})}_{\tiny\mbox{Cost of onsite gen.}}-\underbrace{\beta U^{\rm
pv}_{f,t}P^{\rm max}_{f,t}}_{\tiny\mbox{Env. cost of PV
curtailment}}\,\right),$

where $p_{f,t}^{\rm R}(\cdot)$ denotes the price function of renewables seen
at firm $f$ at time $t$, while the revenue obtained is given by $p_{f,t}^{\rm
R}\left(P^{\rm pv,S,V}_{f,t}+P^{\rm pv,S}_{f,t}\right)P^{\rm pv,S,V}_{f,t}$.
We may then formally define the optimization problem faced by the VPP
associated with firm $f$, where the polyhedral constraints are captured by
$\left\\{\,y^{\rm vpp}_{f}\,\mid\,A_{f}y^{\rm vpp}_{f}\leq d_{f}\,\right\\}$
where $A_{f}\in\mathbb{R}^{m\times n}$ and $d_{f}\in\mathbb{R}^{m}$.

$\begin{array}[]{ll}\displaystyle{{\operatornamewithlimits{\mbox{\bf
maximize}}_{y^{\rm vpp}_{f}}}}&r^{\rm vpp}_{f}(y^{\rm vpp}_{f};P_{f}^{\rm
pv,S})\hskip 166.2212pt\left(\mbox{\bf VPP}_{f}(P_{f}^{\rm
pv,S})\right)\\\\[14.45377pt] \mbox{\bf subject
to}&\eqref{cons1}\,-\,\eqref{cons9}\equiv\left\\{\,y^{\rm
vpp}_{f}\,\mid\,A_{f}y^{\rm vpp}_{f}\leq d_{f}\,\right\\}.\end{array}$

We now observe that the resulting equilibrium problem comprises of a
collection of firms, each of which has a single follower as captured by a VPP,
in addition to the ISO and a set of players that determine prices. To
facilitate analysis of the necessary and sufficient equilibrium conditions of
this hierarchical game, we approximate $\left(\mbox{\bf VPP}_{f}(P_{f}^{\rm
pv,S})\right)$ by employing a smooth (exact) penalized approximation; this
latter formulation is of particular relevance in deriving the concavity of the
function $p_{f,t}^{\rm R}\left(P_{f,t}^{\rm pv,S}+P_{f,t}^{\rm
pv,S,V,\epsilon}(P_{f,t}^{\rm pv,S})\right)P_{f,t}^{\rm pv,S}$ in
$P_{f,t}^{\rm pv,S}$, where $P_{f,t}^{\rm pv,S,V,\epsilon}(P_{f,t}^{\rm
pv,S})$ is a component of the single-valued solution map $y^{\rm
vpp,\epsilon}_{f}(P_{f,t}^{\rm pv,S})$, a solution of the
$\epsilon$-regularized and the $\epsilon$-smoothed (exact) penalized
approximation of $\left(\mbox{\bf VPP}_{f}(P_{f}^{\rm pv,S})\right)$. To this
end, we define the exact penalty function $\varphi$ and its smoothed
counterpart $\varphi_{\epsilon}$ as

$\displaystyle\varphi(A_{f}y^{\rm vpp}_{f}-d_{f})\,$
$\displaystyle\triangleq\,\sum_{i=1}^{m}\max\\{a_{fi}^{\top}y^{\rm
vpp}_{f}-d_{fi},0\\},\quad\varphi_{\epsilon}(A_{f}y^{\rm
vpp}_{f}-d_{f})\,\triangleq\,\sum_{i=1}^{m}\psi_{\epsilon}(a_{fi}^{\top}y^{\rm
vpp}_{f}-d_{fi}),$ (5.12) $\displaystyle\mbox{ and }\psi_{\epsilon}(t)$
$\displaystyle\triangleq\begin{cases}0,&\mbox{ if }t\,\leq\,0\\\ \tfrac{\mu
t^{2}}{2\epsilon},&\mbox{ if }0\,\leq\,t\,\leq\,\epsilon\\\
\mu(t-\tfrac{\epsilon}{2}).&\mbox{ if }t\,\geq\,\epsilon\end{cases}$ (5.13)

This penalty function and its smoothed counterpart are employed in formally
defining the exact penalty reformulation of (VPPf) and its smoothed
counterpart.

$\displaystyle\ \left\\{\begin{aligned} \max_{y^{\rm vpp}_{f}}&\ r^{\rm
vpp}(y^{\rm vpp}_{f},P_{f}^{\rm pv,S})\\\ \mbox{subject to}&\ Ay^{\rm
vpp}_{f}\,\leq\,d_{f}\end{aligned}\right\\}\ \equiv\,\left\\{\begin{aligned}
\max_{y^{\rm vpp}_{f}}&\ \left(r^{\rm vpp}(y^{\rm vpp}_{f},P_{f}^{\rm
pv,S})+\mu\varphi(Ay^{\rm vpp}_{f}-d_{f})\right)\end{aligned}\right\\}$
$\displaystyle\hskip 14.45377pt\ \approx\,\left\\{\begin{aligned} \max_{y^{\rm
vpp}_{f}}&\left(r^{\rm vpp}(y^{\rm vpp}_{f},P_{f}^{\rm
pv,S})+\mu\varphi_{\epsilon}(Ay^{\rm
vpp}_{f}-d_{f})+\tfrac{\epsilon}{2}\|y^{\rm
vpp}_{f}\|^{2}\right).\end{aligned}\right\\}$ (VPP${}_{f}^{\epsilon}(P^{\rm
pv,S}_{f})$)

We observe that the resulting game can be viewed as a noncooperative
hierarchical game, defined by upper-level player problems given by firm
players $\left(\mbox{\bf Firm}^{\rm
rel}_{f}(s_{-f},g_{-f},\lambda_{f})\right)_{f\,\in\,{\cal F}}$, pricing
players $\left(\mbox{\bf Price}\left(s,g,P^{\rm pv,firm-load}\right)\right)$,
the ISO as denoted by (ISO), and the transmission pricing player $(\mbox{\bf
Flow}^{\rm price}(g,s,y))$. In addition, the set of regularized lower-level
VPP problems is given by $\left(\mbox{\bf VPP${}_{f}^{\epsilon}(P^{\rm
pv,S}_{f})$}\right)_{f\in{\cal F}}$. We succinctly represent this
noncooperative game as an $N+2$ player game, in which the first $N$ players
correspond to firm $f$’s problem for $f\in{\cal F}$ while the last two
correspond to pricing players.

$\displaystyle\left\\{\begin{aligned}
\min_{\mathbf{z}^{1}\in\mathbf{Z}^{1}}\,\left(h_{1}(\mathbf{z}^{1},\mathbf{z}^{-1})+g_{1}(\mathbf{z}^{1},\mathbf{u}^{1}(\mathbf{z}^{1}))\right)\end{aligned}\right\\}\quad\cdots\quad\left\\{\begin{aligned}
\min_{\mathbf{z}^{N}\in\mathbf{Z}^{N}}\,\left(h_{N}(\mathbf{z}^{N},\mathbf{z}^{-N})+g_{N}(\mathbf{z}^{N},\mathbf{u}^{N}(\mathbf{z}^{N}))\right)\end{aligned}\right\\},$
$\displaystyle\left\\{\begin{aligned}
\min_{\mathbf{z}^{N+1}\in\mathbf{Z}^{N+1}}\,\left(h_{N+1}(\mathbf{z}^{N+1},\mathbf{z}^{-(N+1)})\right)\end{aligned}\right\\}\quad\left\\{\begin{aligned}
\min_{\mathbf{z}^{N+2}\in\mathbf{Z}^{N+2}}\,\left(h_{N+2}(\mathbf{z}^{N+2},\mathbf{z}^{-(N+2)})\right)\end{aligned}\right\\}.$
(Hier-gamevpp)

Note that the first $N$ players’ objectives are characterized by hierarchical
terms; specifically, the hierarchical terms
$g_{j}(\mathbf{z}^{j},\mathbf{u}^{j}(\mathbf{z}^{j}))$ for any
$j\,\in\\{1,\ldots,N\\}$ correspond to the hierarchical terms in firm $f$’s
problem given by ${\displaystyle\sum_{t\in{\cal T}}}\left(p_{f,t}^{\rm
R}\left(P_{f,t}^{\rm pv,S}+P_{f,t}^{\rm pv,S,V}\right)P_{f,t}^{\rm
pv,S}\right)$ for any $f\,\in\,{\cal F}$.

While convexity of player problems follows in a straightforward fashion from
the definition of firm problems and suitable convexity requirements on the
cost functions as well as affineness requirements on the price functions.
Single-valuedness of the solution map $\mathbf{u}^{j}(\bullet)$ follows from
the observation that the regularized VPP profit function is strongly concave.
Additionally, convexity of $g_{j}(\bullet,\mathbf{u}^{j}(\bullet))$ is a
consequence of analogous analysis for Stackelberg leadership (cf. [66, 15]).
Finally, monotonicity of $F$ can be derived in a fashion similar to that
considered in [30]. Existence of an equilibrium can then be derived in a
fashion similar to that employed in [15]. A comprehensive analysis of this
model is left to future work.

##  6 Proof of the Main Theorem

###  6.1 Analysis of the exact scheme

The proof on the finite time-complexity estimate starts by a technical
derivation of an energy-type inequality that gives us an upper bound on the
change of the energy function
$\frac{1}{2}\norm{\mathbf{z}^{(t)}_{k+1}-\mathbf{x}}^{2}$, computed within an
arbitrary inner loop evaluation, and for an arbitrary anchor point
$\mathbf{x}\in\mathcal{X}$. Via a sequence of tedious, but otherwise
straightforward, manipulations we arrive out our first main result, Lemma 4.3.
From there, we proceed as in the standard analysis of stochastic approximation
schemes [52], and derive an upper bound on the gap function of the mixed
variational inequality.

####  6.1.1 Proof of Lemma 4.3

To simplify notation we omit the dependence on the outer iteration loop $t$,
and thus simply write $\mathbf{z}_{k}$ for $\mathbf{z}^{(t)}_{k}$. The same
notational simplification will be used in all variables that are computed
within the inner loop executed in the $t$-round of the outer loop procedure.
With the hope that the reader agrees that this reduces notational complexity a
bit, we proceed to derive the the postulated energy inequality. To start, we
observe that for each $\mathbf{x}\in\mathcal{X}$ we have

$\displaystyle\norm{\mathbf{z}_{k+1}-\mathbf{x}}^{2}$
$\displaystyle=\norm{\mathbf{z}_{k+1}-\mathbf{z}_{k+1/2}+\mathbf{z}_{k+1/2}-\mathbf{z}_{k}+\mathbf{z}_{k}-\mathbf{x}}^{2}$
$\displaystyle=\norm{\mathbf{z}_{k+1}-\mathbf{z}_{k+1/2}}^{2}-\norm{\mathbf{z}_{k+1/2}-\mathbf{z}_{k}}^{2}+\norm{\mathbf{z}_{k}-\mathbf{x}}^{2}+2\langle\mathbf{z}_{k+1}-\mathbf{z}_{k+1/2},\mathbf{z}_{k+1/2}-\mathbf{x}\rangle$
$\displaystyle+2\langle\mathbf{z}_{k+1/2}-\mathbf{z}_{k},\mathbf{z}_{k+1/2}-\mathbf{x}\rangle$
$\displaystyle=\norm{\mathbf{z}_{k}-\mathbf{x}}^{2}-\norm{\mathbf{z}_{k+1/2}-\mathbf{z}_{k}}^{2}$
$\displaystyle+\norm{\gamma_{t}\left(\hat{V}^{\eta_{t}}_{k+1/2}(\mathbf{z}_{k+1/2})-\hat{V}^{\eta_{t}}_{k+1/2}(\mathbf{x}^{t})+H^{\delta_{t}}_{\mathbf{z}_{k+1/2}}(\mathbf{W}_{k+1/2})-H^{\delta_{t}}_{\mathbf{x}^{t}}(\mathbf{W}_{k+1/2})\right)}^{2}$
$\displaystyle-2\gamma_{t}\langle\hat{V}^{\eta_{t}}_{k+1/2}(\mathbf{z}_{k+1/2})-\hat{V}^{\eta_{t}}_{k+1/2}(\mathbf{x}^{t})+H^{\delta_{t}}_{\mathbf{z}_{k+1/2}}(\mathbf{W}_{k+1/2})-H^{\delta_{t}}_{\mathbf{x}^{t}}(\mathbf{W}_{k+1/2}),\mathbf{z}_{k+1/2}-\mathbf{x}\rangle$
$\displaystyle+2\langle\mathbf{z}_{k+1/2}-\left(\mathbf{z}_{k}-\gamma_{t}(\bar{V}^{t}+H^{\delta_{t},b_{t}}_{\mathbf{x}^{t}}+\eta_{t}\mathbf{x}^{t})\right),\mathbf{z}_{k+1/2}-\mathbf{x}\rangle$
$\displaystyle-2\gamma_{t}\langle\bar{V}^{t}+H^{\delta_{t},b_{t}}_{\mathbf{x}^{t}}+\eta_{t}\mathbf{x}^{t},\mathbf{z}_{k+1/2}-\mathbf{x}\rangle$

Lemma A.1(i) gives
$2\langle\mathbf{z}_{k+1/2}-\left(\mathbf{z}_{k}-\gamma_{t}(\bar{V}^{t}+H^{\delta_{t},b_{t}}_{\mathbf{x}^{t}}+\eta_{t}\mathbf{x}^{t}\right),\mathbf{z}_{k+1/2}-\mathbf{x}\rangle\leq
0.$ Additionally, Assumption 4 and Lemma A.12(b) gives

$\displaystyle\norm{\gamma_{t}\left(\hat{V}^{\eta_{t}}_{k+1/2}(\mathbf{z}_{k+1/2})-\hat{V}^{\eta_{t}}_{k+1/2}(\mathbf{x}^{t})+H^{\delta_{t}}_{\mathbf{z}_{k+1/2}}(\mathbf{W}_{k+1/2})-H^{\delta_{t}}_{\mathbf{x}^{t}}(\mathbf{W}_{k+1/2})\right)}^{2}$
$\displaystyle\leq
2\gamma^{2}_{t}\norm{\hat{V}^{\eta_{t}}_{k+1/2}(\mathbf{z}_{k+1/2})-\hat{V}^{\eta_{t}}_{k+1/2}(\mathbf{x}^{t})}^{2}+2\gamma^{2}_{t}\norm{H^{\delta_{t}}_{\mathbf{z}_{k+1/2}}(\mathbf{W}_{k+1/2})-H^{\delta_{t}}_{\mathbf{x}^{t}}(\mathbf{W}_{k+1/2})}^{2}$
$\displaystyle\leq
4\gamma^{2}_{t}\norm{\hat{V}_{k+1/2}(\mathbf{z}_{k+1/2})-\hat{V}_{k+1/2}(\mathbf{x}^{t})}^{2}+4\eta^{2}_{t}\gamma^{2}_{t}\norm{\mathbf{z}_{k+1/2}-\mathbf{x}^{t}}^{2}$
$\displaystyle+4\gamma^{2}_{t}\left(\norm{H^{\delta_{t}}_{\mathbf{z}_{k+1/2}}(\mathbf{W}_{k+1/2})}^{2}+\norm{H^{\delta_{t}}_{\mathbf{x}^{t}}(\mathbf{W}_{k+1/2})}^{2}\right)$
$\displaystyle\leq
4\gamma^{2}_{t}\left(\mathcal{L}_{f}(\xi_{k+1/2})^{2}+\eta^{2}_{t}\right)\norm{\mathbf{z}_{k+1/2}-\mathbf{x}^{t}}^{2}+8\gamma^{2}_{t}\sum_{i\in\mathcal{I}}L_{h_{i}}^{2}n_{i}^{2}.$

It follows

$\displaystyle\norm{\mathbf{z}_{k+1}-\mathbf{x}}^{2}$
$\displaystyle\leq\norm{\mathbf{z}_{k}-\mathbf{x}}^{2}-\norm{\mathbf{z}_{k+1/2}-\mathbf{z}_{k}}^{2}+8\gamma_{t}^{2}\sum_{i\in\mathcal{I}}L_{h_{i}}^{2}n^{2}_{i}$
$\displaystyle+4\gamma^{2}_{t}(\mathcal{L}_{f}(\xi_{k+1/2})^{2}+\eta^{2}_{0})\norm{\mathbf{z}_{k+1/2}-\mathbf{x}^{t}}^{2}$
$\displaystyle-2\gamma_{t}\langle\hat{V}^{\eta_{t}}_{k+1/2}(\mathbf{z}_{k+1/2})+H^{\delta_{t}}_{\mathbf{z}_{k+1/2}}(\mathbf{W}_{k+1/2})-\hat{V}^{\eta_{t}}_{k+1/2}(\mathbf{x}^{t})-H^{\delta_{t}}_{\mathbf{x}^{t}}(\mathbf{W}_{k+1/2}),z_{k+1/2}-\mathbf{x}\rangle$
$\displaystyle-2\gamma_{t}\langle\bar{V}^{t}+\eta_{t}\mathbf{x}^{t}+H^{\delta_{t},b_{t}}_{\mathbf{x}^{t}},\mathbf{z}_{k+1/2}-\mathbf{x}\rangle$
$\displaystyle=\norm{\mathbf{z}_{k}-\mathbf{x}}^{2}-\norm{\mathbf{z}_{k+1/2}-\mathbf{z}_{k}}^{2}+8\gamma_{t}^{2}\sum_{i\in\mathcal{I}}L_{h_{i}}^{2}n^{2}_{i}+4\gamma^{2}_{t}(\mathcal{L}_{f}(\xi_{k+1/2})^{2}+\eta^{2}_{0})\norm{\mathbf{z}_{k+1/2}-\mathbf{x}^{t}}^{2}$
$\displaystyle-2\gamma_{t}\langle\hat{V}_{k+1/2}(\mathbf{z}_{k+1/2})+H^{\delta_{t}}_{\mathbf{z}_{k+1/2}}(\mathbf{W}_{k+1/2})-V(\mathbf{z}_{k+1/2})-\nabla
h^{\delta_{t}}(\mathbf{z}_{k+1/2}),\mathbf{z}_{k+1/2}-\mathbf{x}\rangle$
$\displaystyle+2\gamma_{t}\langle\hat{V}_{k+1/2}(\mathbf{x}^{t})+H^{\delta_{t}}_{\mathbf{x}^{t}}(\mathbf{W}_{k+1/2})-\bar{V}^{t}-H^{\delta_{t},b_{t}}_{\mathbf{x}^{t}},\mathbf{z}_{k+1/2}-\mathbf{x}\rangle$
$\displaystyle-2\gamma_{t}\langle V(\mathbf{z}_{k+1/2})+\nabla
h^{\delta_{t}}(\mathbf{z}_{k+1/2})+\eta_{t}\mathbf{z}_{k+1/2},\mathbf{z}_{k+1/2}-\mathbf{x}\rangle$
$\displaystyle=\norm{\mathbf{z}_{k}-\mathbf{x}}^{2}-\norm{\mathbf{z}_{k+1/2}-\mathbf{z}_{k}}^{2}+8\gamma_{t}^{2}\sum_{i\in\mathcal{I}}L_{h_{i}}^{2}n^{2}_{i}+4\gamma^{2}_{t}(\mathcal{L}_{f}(\xi_{k+1/2})^{2}+\eta^{2}_{0})\norm{\mathbf{z}_{k+1/2}-\mathbf{x}^{t}}^{2}$
$\displaystyle-2\gamma_{t}\langle\hat{V}_{k+1/2}(\mathbf{z}_{k+1/2})+H^{\delta_{t}}_{\mathbf{z}_{k+1/2}}(\mathbf{W}_{k+1/2})-V(\mathbf{z}_{k+1/2})-\nabla
h^{\delta_{t}}(\mathbf{z}_{k+1/2}),\mathbf{z}_{k+1/2}-\mathbf{x}\rangle$
$\displaystyle+2\gamma_{t}\langle\hat{V}_{k+1/2}(\mathbf{x}^{t})+H^{\delta_{t}}_{\mathbf{x}^{t}}(\mathbf{W}_{k+1/2})-\bar{V}^{t}-H^{\delta_{t},b_{t}}_{\mathbf{x}^{t}},\mathbf{z}_{k+1/2}-\mathbf{x}\rangle$
$\displaystyle-2\gamma_{t}\langle V^{\eta_{t}}(\mathbf{z}_{k+1/2})+\nabla
h^{\delta_{t}}(\mathbf{z}_{k+1/2})-V^{\eta_{t}}(\mathbf{x})-\nabla
h^{\delta_{t}}(\mathbf{x}),\mathbf{z}_{k+1/2}-\mathbf{x}\rangle$
$\displaystyle-2\gamma_{t}\langle V^{\eta_{t}}(\mathbf{x})+\nabla
h^{\delta_{t}}(\mathbf{x}),\mathbf{z}_{k+1/2}-\mathbf{x}\rangle$

Since the operator $\mathbf{x}\mapsto V^{\eta_{t}}(\mathbf{x})+\nabla
h^{\delta_{t}}(\mathbf{x})$ is $\eta_{t}$-strongly monotone, we can further
bound the expression above as

$\displaystyle\norm{\mathbf{z}_{k+1}-\mathbf{x}}^{2}$
$\displaystyle\leq\norm{\mathbf{z}_{k}-\mathbf{x}}^{2}-\norm{\mathbf{z}_{k+1/2}-\mathbf{z}_{k}}^{2}+8\gamma_{t}^{2}\sum_{i\in\mathcal{I}}L_{h_{i}}^{2}n^{2}_{i}+4\gamma^{2}_{t}(\mathcal{L}_{f}(\xi_{k+1/2})^{2}+\eta^{2}_{0})\norm{\mathbf{z}_{k+1/2}-\mathbf{x}^{t}}^{2}$
$\displaystyle-2\gamma_{t}\langle\hat{V}_{k+1/2}(\mathbf{z}_{k+1/2})+H^{\delta_{t}}_{\mathbf{z}_{k+1/2}}(\mathbf{W}_{k+1/2})-V(\mathbf{z}_{k+1/2})-\nabla
h^{\delta_{t}}(\mathbf{z}_{k+1/2}),\mathbf{z}_{k+1/2}-\mathbf{x}\rangle$
$\displaystyle+2\gamma_{t}\langle\hat{V}_{k+1/2}(\mathbf{x}^{t})+H^{\delta_{t}}_{\mathbf{x}^{t}}(\mathbf{W}_{k+1/2})-\bar{V}^{t}-H^{\delta_{t},b_{t}}_{\mathbf{x}^{t}},\mathbf{z}_{k+1/2}-\mathbf{x}\rangle$
$\displaystyle-2\gamma_{t}\eta_{t}\norm{\mathbf{z}_{k+1/2}-\mathbf{x}}^{2}$
$\displaystyle-2\gamma_{t}\langle V^{\eta_{t}}(\mathbf{x})+\nabla
h^{\delta_{t}}(\mathbf{x}),\mathbf{z}_{k+1/2}-\mathbf{x}\rangle$

Next, we split the mini-batch estimator $\bar{V}^{t}$ into its mean component
and its error component so that

$\bar{V}^{t}=V(\mathbf{x}^{t})+\varepsilon^{t}_{V}(\xi^{1:b_{t}}_{t}).$ (6.1)

Similarly, we write

$\varepsilon^{t}_{h}(\mathbf{W}^{1:b_{t}}_{t})\triangleq
H^{\delta_{t},b_{t}}_{\mathbf{x}^{t}}-\nabla h^{\delta_{t}}(\mathbf{x}^{t}).$
(6.2)

Using these error terms, we may further bound the right hand side of the
penultimate display as

$\displaystyle\norm{\mathbf{z}_{k+1}-\mathbf{x}}^{2}$
$\displaystyle\leq\norm{\mathbf{z}_{k}-\mathbf{x}}^{2}-\norm{\mathbf{z}_{k+1/2}-\mathbf{z}_{k}}^{2}+8\gamma_{t}^{2}\sum_{i\in\mathcal{I}}L_{h_{i}}^{2}n^{2}_{i}+4\gamma^{2}_{t}(\mathcal{L}_{f}(\xi_{k+1/2})^{2}+\eta^{2}_{0})\norm{\mathbf{z}_{k+1/2}-\mathbf{x}^{t}}^{2}$
$\displaystyle-2\gamma_{t}\langle\hat{V}_{k+1/2}(\mathbf{z}_{k+1/2})+H^{\delta_{t}}_{\mathbf{z}_{k+1/2}}(\mathbf{W}_{k+1/2})-V(\mathbf{z}_{k+1/2})-\nabla
h^{\delta_{t}}(\mathbf{z}_{k+1/2}),\mathbf{z}_{k+1/2}-\mathbf{x}\rangle$
$\displaystyle-2\gamma_{t}\langle\left(V(\mathbf{x}^{t})+\nabla
h^{\delta_{t}}(\mathbf{x}^{t})\right)-\left(\hat{V}_{k+1/2}(\mathbf{x}^{t})+H^{\delta_{t}}_{\mathbf{x}^{t}}(\mathbf{W}_{k+1/2})\right),\mathbf{z}_{k+1/2}-\mathbf{x}\rangle$
$\displaystyle-2\gamma_{t}\eta_{t}\norm{\mathbf{z}_{k+1/2}-\mathbf{x}}^{2}-2\gamma_{t}\langle\varepsilon^{t}_{V}(\xi^{1:b_{t}}_{t})+\varepsilon^{t}_{h}(\mathbf{W}^{1:b_{t}}_{t}),\mathbf{z}_{k+1/2}-\mathbf{x}\rangle$
$\displaystyle-2\gamma_{t}\langle V^{\eta_{t}}(\mathbf{x})+\nabla
h^{\delta_{t}}(\mathbf{x}),\mathbf{z}_{k+1/2}-\mathbf{x}\rangle$

A simple application of the triangle inequality shows

$-2\gamma_{t}\eta_{t}\norm{\mathbf{z}_{k+1/2}-\mathbf{x}}^{2}\leq
2\gamma_{t}\eta_{t}\norm{\mathbf{z}_{k+1/2}-\mathbf{z}_{k}}^{2}-\gamma_{t}\eta_{t}\norm{\mathbf{z}_{k}-\mathbf{x}}^{2}.$

Using this bound, we continue with the derivations above to arrive at

$\displaystyle\norm{\mathbf{z}_{k+1}-\mathbf{x}}^{2}$
$\displaystyle\leq(1-\gamma_{t}\eta_{t})\norm{\mathbf{z}_{k}-\mathbf{x}}^{2}-(1-2\gamma_{t}\eta_{t})\norm{\mathbf{z}_{k+1/2}-\mathbf{z}_{k}}^{2}$
$\displaystyle+8\gamma_{t}^{2}\sum_{i\in\mathcal{I}}L_{h_{i}}^{2}n^{2}_{i}+4\gamma^{2}_{t}(\mathcal{L}_{f}(\xi_{k+1/2})^{2}+\eta^{2}_{0})\norm{\mathbf{z}_{k+1/2}-\mathbf{x}^{t}}^{2}$
$\displaystyle-2\gamma_{t}\langle\hat{V}_{k+1/2}(\mathbf{z}_{k+1/2})+H^{\delta_{t}}_{\mathbf{z}_{k+1/2}}(\mathbf{W}_{k+1/2})-V(\mathbf{z}_{k+1/2})-\nabla
h^{\delta_{t}}(\mathbf{z}_{k+1/2}),\mathbf{z}_{k+1/2}-\mathbf{x}\rangle$
$\displaystyle-2\gamma_{t}\langle\left(V(\mathbf{x}^{t})+\nabla
h^{\delta_{t}}(\mathbf{x}^{t})\right)-\left(\hat{V}_{k+1/2}(\mathbf{x}^{t})+H^{\delta_{t}}_{\mathbf{x}^{t}}(\mathbf{W}_{k+1/2})\right),\mathbf{z}_{k+1/2}-\mathbf{x}\rangle$
$\displaystyle-2\gamma_{t}\langle\varepsilon^{t}_{V}(\xi^{1:b_{t}}_{t})+\varepsilon^{t}_{h}(\mathbf{W}^{1:b_{t}}_{t}),\mathbf{z}_{k+1/2}-\mathbf{x}\rangle$
$\displaystyle-2\gamma_{t}\langle V^{\eta_{t}}(\mathbf{x})+\nabla
h^{\delta_{t}}(\mathbf{x}),\mathbf{z}_{k+1/2}-\mathbf{x}\rangle$

By convexity of the application $\mathbf{x}\mapsto
h^{\delta_{t}}(\mathbf{x})$, we have

$h^{\delta_{t}}(\mathbf{z}_{k+1/2})\geq
h^{\delta_{t}}(\mathbf{x})+\langle\nabla
h^{\delta_{t}}(\mathbf{x}),\mathbf{z}_{k+1/2}-\mathbf{x}\rangle.$

Hence, the penultimate display turns into

$\displaystyle\norm{\mathbf{z}_{k+1}-\mathbf{x}}^{2}$
$\displaystyle\leq(1-\gamma_{t}\eta_{t})\norm{\mathbf{z}_{k}-\mathbf{x}}^{2}-(1-2\gamma_{t}\eta_{t})\norm{\mathbf{z}_{k+1/2}-\mathbf{z}_{k}}^{2}$
$\displaystyle+8\gamma_{t}^{2}\sum_{i\in\mathcal{I}}L_{h_{i}}^{2}n^{2}_{i}+4\gamma^{2}_{t}(\mathcal{L}_{f}(\xi_{k+1/2})^{2}+\eta^{2}_{0})\norm{\mathbf{z}_{k+1/2}-\mathbf{x}^{t}}^{2}$
$\displaystyle-2\gamma_{t}\langle\hat{V}_{k+1/2}(\mathbf{z}_{k+1/2})+H^{\delta_{t}}_{\mathbf{z}_{k+1/2}}(\mathbf{W}_{k+1/2})-V(\mathbf{z}_{k+1/2})-\nabla
h^{\delta_{t}}(\mathbf{z}_{k+1/2}),\mathbf{z}_{k+1/2}-\mathbf{x}\rangle$
$\displaystyle-2\gamma_{t}\langle\left(V(\mathbf{x}^{t})+\nabla
h^{\delta_{t}}(\mathbf{x}^{t})\right)-\left(\hat{V}_{k+1/2}(\mathbf{x}^{t})+H^{\delta_{t}}_{\mathbf{x}^{t}}(\mathbf{W}_{k+1/2})\right),\mathbf{z}_{k+1/2}-\mathbf{x}\rangle$
$\displaystyle-2\gamma_{t}\langle\varepsilon^{t}_{V}(\xi^{1:b_{t}}_{t})+\varepsilon^{t}_{h}(\mathbf{W}^{1:b_{t}}_{t}),\mathbf{z}_{k+1/2}-\mathbf{x}\rangle$
$\displaystyle-2\gamma_{t}\left(\langle
V^{\eta_{t}}(\mathbf{x}),\mathbf{z}_{k+1/2}-\mathbf{x}\rangle+h^{\delta_{t}}(\mathbf{z}_{k+1/2})-h^{\delta_{t}}(\mathbf{x})\right).$

From Lemma 2 in [68], we know for
$L_{h}\triangleq\sum_{i\in\mathcal{I}}L_{h_{i}}$ that

$h^{\delta_{t}}(\mathbf{z}_{k+1/2})-h^{\delta_{t}}(\mathbf{x})\geq
h(\mathbf{z}_{k+1/2})-h(\mathbf{x})-\delta_{t}L_{h}.$

Therefore,

$\displaystyle\norm{\mathbf{z}_{k+1}-\mathbf{x}}^{2}$
$\displaystyle\leq(1-\gamma_{t}\eta_{t})\norm{\mathbf{z}_{k}-\mathbf{x}}^{2}-(1-2\gamma_{t}\eta_{t})\norm{\mathbf{z}_{k+1/2}-\mathbf{z}_{k}}^{2}$
$\displaystyle+8\gamma_{t}^{2}\sum_{i\in\mathcal{I}}L_{h_{i}}^{2}n^{2}_{i}+4\gamma^{2}_{t}(\mathcal{L}_{f}(\xi_{k+1/2})^{2}+\eta^{2}_{0})\norm{\mathbf{z}_{k+1/2}-\mathbf{x}^{t}}^{2}$
$\displaystyle-2\gamma_{t}\langle\hat{V}_{k+1/2}(\mathbf{z}_{k+1/2})+H^{\delta_{t}}_{\mathbf{z}_{k+1/2}}(\mathbf{W}_{k+1/2})-V(\mathbf{z}_{k+1/2})-\nabla
h^{\delta_{t}}(\mathbf{z}_{k+1/2}),\mathbf{z}_{k+1/2}-\mathbf{x}\rangle$
$\displaystyle-2\gamma_{t}\langle\left(V(\mathbf{x}^{t})+\nabla
h^{\delta_{t}}(\mathbf{x}^{t})\right)-\left(\hat{V}_{k+1/2}(\mathbf{x}^{t})+H^{\delta_{t}}_{\mathbf{x}^{t}}(\mathbf{W}_{k+1/2})\right),\mathbf{z}_{k+1/2}-\mathbf{x}\rangle$
$\displaystyle-2\gamma_{t}\langle\varepsilon^{t}_{V}(\xi^{1:b_{t}}_{t})+\varepsilon^{t}_{h}(\mathbf{W}^{1:b_{t}}_{t}),\mathbf{z}_{k+1/2}-\mathbf{x}\rangle$
$\displaystyle-2\gamma_{t}\left(\langle
V^{\eta_{t}}(\mathbf{x}),\mathbf{z}_{k+1/2}-\mathbf{x}\rangle+h(\mathbf{z}_{k+1/2})-h(\mathbf{x})\right)+2\gamma_{t}\delta_{t}L_{h}$

which is what has been claimed.

####  6.1.2 Proof of Theorem 4.4

Let $\\{\eta_{t}\\}_{t\in\mathbb{N}}$ be a positive sequence with
$\eta_{t}\downarrow 0$. Let $\mathbf{s}_{t}\equiv\mathbf{s}(\eta_{t})$ denote
the corresponding sequence of solutions to
$\operatorname{MVI}(V^{\eta_{t}},h)$. Set $q_{t}\equiv
1-\gamma_{t}\eta_{t}\in(0,1/2)$. Then, iterating the energy inequality
established in Lemma 4.3, we have for $\mathbf{x}=\mathbf{s}_{t}$:

$\displaystyle\norm{\mathbf{z}^{(t)}_{K}-\mathbf{s}_{t}}^{2}\leq
q_{t}^{K}\norm{\mathbf{z}^{(t)}_{0}-\mathbf{s}_{t}}^{2}$
$\displaystyle+\sum_{k=0}^{K-1}q_{t}^{K-k+1}\gamma_{t}^{2}\left(8\left(\sum_{i\in\mathcal{I}}L_{h_{i}}^{2}n^{2}_{i}\right)+4(\mathcal{L}_{f}(\xi_{t,k+1/2})^{2}+\eta^{2}_{0})\norm{\mathbf{z}^{(t)}_{k+1/2}-\mathbf{x}^{t}}^{2}\right)$
$\displaystyle-2\gamma_{t}\sum_{k=0}^{K-1}q_{t}^{K-k+1}\langle\hat{V}_{t,k+1/2}(\mathbf{z}^{(t)}_{k+1/2})+H^{\delta_{t}}_{\mathbf{z}^{(t)}_{k+1/2}}(\mathbf{W}_{t,k+1/2})-V(\mathbf{z}^{(t)}_{k+1/2})-\nabla
h^{\delta_{t}}(\mathbf{z}^{(t)}_{k+1/2}),\mathbf{z}^{(t)}_{k+1/2}-\mathbf{s}_{t}\rangle$
$\displaystyle-2\gamma_{t}\sum_{k=0}^{K-1}q_{t}^{K-k+1}\langle\left(V(\mathbf{x}^{t})+\nabla
h^{\delta_{t}}(\mathbf{x}^{t})\right)-\left(\hat{V}_{t,k+1/2}(\mathbf{x}^{t})+H^{\delta_{t}}_{\mathbf{x}^{t}}(\mathbf{W}_{t,k+1/2})\right),\mathbf{z}^{(t)}_{k+1/2}-\mathbf{s}_{t}\rangle$
$\displaystyle-2\gamma_{t}\sum_{k=0}^{K-1}q_{t}^{K-k+1}\langle\varepsilon^{t}_{V}(\xi^{1:b_{t}}_{t})+\varepsilon^{t}_{h}(\mathbf{W}^{1:b_{t}}_{t}),\mathbf{z}^{(t)}_{k+1/2}-\mathbf{s}_{t}\rangle$
$\displaystyle-2\gamma_{t}\sum_{k=0}^{K-1}q^{K-k+1}_{t}\left(\langle
V^{\eta_{t}}(\mathbf{s}_{t}),\mathbf{z}^{(t)}_{k+1/2}-\mathbf{s}_{t}\rangle+h(\mathbf{z}^{(t)}_{k+1/2})-h(\mathbf{s}_{t})-2\delta_{t}L_{h}\right).$

By definition of the point $\mathbf{s}_{t}$, we have $\langle
V^{\eta_{t}}(\mathbf{s}_{t}),\mathbf{z}^{(t)}_{k+1/2}-\mathbf{s}_{t}\rangle+h(\mathbf{z}^{(t)}_{k+1/2})-h(\mathbf{s}_{t})\geq
0$. Furthermore, the estimators involved are unbiased, which means

$\displaystyle\mathbb{E}[\langle\hat{V}_{t,k+1/2}(\mathbf{z}^{(t)}_{k+1/2})+H^{\delta_{t}}_{\mathbf{z}^{(t)}_{k+1/2}}(\mathbf{W}_{t,k+1/2})-V(\mathbf{z}^{(t)}_{k+1/2})-\nabla
h^{\delta_{t}}(\mathbf{z}^{(t)}_{k+1/2}),\mathbf{z}^{(t)}_{k+1/2}-\mathbf{s}_{t}\rangle|\mathcal{F}_{t}]=0,$
$\displaystyle\mathbb{E}[\left(\langle V(\mathbf{x}^{t})+\nabla
h^{\delta_{t}}(\mathbf{x}^{t})\right)-\left(\hat{V}_{t,k+1/2}(\mathbf{x}^{t})+H^{\delta_{t}}_{\mathbf{x}^{t}}(\mathbf{W}_{t,k+1/2})\right),\mathbf{z}^{(t)}_{k+1/2}-\mathbf{s}_{t}\rangle|\mathcal{F}_{t}]=0,$
$\displaystyle\mathbb{E}[\langle\varepsilon^{t}_{V}(\xi^{1:b_{t}}_{t})+\varepsilon^{t}_{h}(\mathbf{W}^{1:b_{t}}_{t}),\mathbf{z}^{(t)}_{k+1/2}-\mathbf{s}_{t}\rangle|\mathcal{F}_{t}]=0.$

To wit, let us focus on the first line of the above display. Using the law of
iterated expectations, we have

$\displaystyle\mathbb{E}[\langle\hat{V}_{t,k+1/2}(\mathbf{z}^{(t)}_{k+1/2})+H^{\delta_{t}}_{\mathbf{z}^{(t)}_{k+1/2}}(\mathbf{W}_{t,k+1/2})-V(\mathbf{z}^{(t)}_{k+1/2})-\nabla
h^{\delta_{t}}(\mathbf{z}^{(t)}_{k+1/2}),\mathbf{z}^{(t)}_{k+1/2}-\mathbf{s}_{t}\rangle|\mathcal{F}_{t}]=$
$\displaystyle\mathbb{E}\left[\mathbb{E}\left(\langle\hat{V}_{t,k+1/2}(\mathbf{z}^{(t)}_{k+1/2})+H^{\delta_{t}}_{\mathbf{z}^{(t)}_{k+1/2}}(\mathbf{W}_{t,k+1/2})-V(\mathbf{z}^{(t)}_{k+1/2})-\nabla
h^{\delta_{t}}(\mathbf{z}^{(t)}_{k+1/2}),\mathbf{z}^{(t)}_{k+1/2}-\mathbf{s}_{t}\rangle|\mathcal{A}_{t,k}\right)|\mathcal{F}_{t}\right]$
$\displaystyle=0$

This being true because $\mathbf{z}_{k+1/2}$ is
$\mathcal{A}_{t,k}$-measurable. The remaining two equalities can be
demonstrated in the same way. Since $\mathbf{z}_{K}=\mathbf{x}^{t+1}$ and
$\mathbf{z}_{0}=\mathbf{x}^{t}$, and using the results above, we are left with
the estimate

$\mathbb{E}[\norm{\mathbf{x}^{t+1}-\mathbf{s}_{t}}^{2}|\mathcal{F}_{t}]\leq
q_{t}^{K}\norm{\mathbf{x}^{t}-\mathbf{s}_{t}}^{2}+K\gamma_{t}^{2}(L_{f}^{2}+\eta^{2}_{0})C^{2}+8\gamma_{t}^{2}\left(\sum_{i\in\mathcal{I}}L^{2}_{h_{i}}n^{2}_{i}\right)+2K\gamma_{t}\delta_{t}L_{h}.$

From Proposition 3.1(d), we obtain the estimate

$\displaystyle\norm{\mathbf{x}^{t}-\mathbf{s}_{t}}^{2}$
$\displaystyle\leq(1+\gamma_{t}\eta_{t})\norm{\mathbf{x}^{t}-\mathbf{s}_{t-1}}^{2}+(1+\frac{1}{\gamma_{t}\eta_{t}})\norm{\mathbf{s}_{t}-\mathbf{s}_{t-1}}^{2}$
$\displaystyle\leq(1+\gamma_{t}\eta_{t})\norm{\mathbf{x}^{t}-\mathbf{s}_{t-1}}^{2}+(1+\gamma_{t}\eta_{t})\left(\frac{\eta_{t}-\eta_{t-1}}{\eta_{t}}\right)^{2}{\mathbf{A}}^{2}_{x}$

where ${\mathbf{A}}_{x}$ is a constant upper bound of
$\inf_{\mathbf{x}\in\operatorname{SOL}(V,h)}\norm{\mathbf{x}}$ (cf.
Proposition A.8). Moreover,

$\displaystyle(1+\gamma_{t}\eta_{t})q_{t}^{K}=(1+\gamma_{t}\eta_{t})(1-\gamma_{t}\eta_{t})q_{t}^{K-1}=(1-\gamma^{2}_{t}\eta_{t}^{2})q_{t}^{K-1}<q_{t}.$

This allows us to conclude

$\displaystyle\mathbb{E}[\norm{\mathbf{x}^{t+1}-\mathbf{s}_{t}}^{2}|\mathcal{F}_{t}]$
$\displaystyle\leq
q_{t}\norm{\mathbf{x}^{t}-\mathbf{s}_{t-1}}^{2}+q_{t}^{K}(1+\frac{1}{\gamma_{t}\eta_{t}}){\mathbf{A}}^{2}_{x}\left(\frac{\eta_{t}-\eta_{t-1}}{\eta_{t}}\right)^{2}$
$\displaystyle+K\gamma^{2}_{t}(L^{2}_{f}+\eta^{2}_{0})C^{2}+8\sum_{i\in\mathcal{I}}L^{2}_{h_{i}}n^{2}_{i}+2KL_{h}\gamma_{t}\delta_{t}.$

Define

$\mathtt{a}_{t}\triangleq
q_{t}^{K}(1+\frac{1}{\gamma_{t}\eta_{t}}){\mathbf{A}}^{2}_{x}\left(\frac{\eta_{t}-\eta_{t-1}}{\eta_{t}}\right)^{2}+K\gamma^{2}_{t}(L^{2}_{f}+\eta^{2}_{0})C^{2}+8\gamma^{2}_{t}\sum_{i\in\mathcal{I}}L^{2}_{h_{i}}n^{2}_{i}+2KL_{h}\gamma_{t}\delta_{t},$

and
$\psi_{t}\triangleq\mathbb{E}[\norm{\mathbf{x}^{t}-\mathbf{s}_{t-1}}^{2}]$, so
that we obtain the recursion

$\psi_{t+1}\leq q_{t}\psi_{t}+\mathtt{a}_{t}.$

Under the assumptions stated in Theorem 4.4, we have
$\sum_{t=0}^{\infty}\gamma_{t}\eta_{t}=\infty$ and
$\lim_{t\to\infty}\frac{\mathtt{a}_{t}}{\gamma_{t}\eta_{t}}=0$. Using Lemma 3
in [58], it follows $\lim_{t\to\infty}\psi_{t}=0$, and therefore
$\lim_{t\to\infty}\norm{\mathbf{x}^{t}-\mathbf{s}_{t}}=0$ almost surely. Now,
let
$\mathbf{w}\triangleq\inf_{\mathbf{x}\in\operatorname{SOL}(V,h)}\norm{\mathbf{x}}$.
Then, using the triangle inequality we conclude

$\norm{\mathbf{x}^{t+1}-\mathbf{w}}\leq\norm{\mathbf{x}^{t+1}-\mathbf{s}_{t}}+\norm{\mathbf{s}_{t}-\mathbf{w}}\to
0\quad\text{ as }t\to\infty,\text{ a.s.}$

Since $\mathbf{y}(\cdot)$ is Lipschitz continuous (cf. Fact A.4, Appendix
A.2), the claim follows.

####  6.1.3 Proof of Theorem 4.5

We use the energy estimate formulated in Lemma 4.3 to deduce a bound on the
gap function (4.1) relative to a suitably constructed ergodic average.

######  Lemma 6.1.

For any $t\in\\{0,1,\ldots,T-1\\}$, define
$\bar{\mathbf{z}}^{t}\triangleq\frac{1}{K}\sum_{k=0}^{K-1}\mathbf{z}^{(t)}_{k+1/2}$.
Let $\\{\gamma_{t}\\}_{t},\\{\eta_{t}\\}_{t}$ be positive sequences satisfying
$0<\gamma_{t}\eta_{t}<1/2$. For $k\in\\{0,1,\ldots,K-1\\}$ define

$\displaystyle\mathbf{Y}^{1}_{t,k}$
$\displaystyle\triangleq\hat{V}_{t,k+1/2}(\mathbf{z}^{(t)}_{k+1/2})+H^{\delta_{t}}_{\mathbf{z}^{(t)}_{k+1/2}}(\mathbf{W}_{t,k+1/2})-V(\mathbf{z}^{(t)}_{k+1/2})-\nabla
h^{\delta_{t}}(\mathbf{z}^{(t)}_{k+1/2}),\text{ and }$
$\displaystyle\mathbf{Y}^{2}_{t,k}$ $\displaystyle\triangleq
V(\mathbf{x}^{t})+\nabla
h^{\delta_{t}}(\mathbf{x}^{t})-\hat{V}_{t,k+1/2}(\mathbf{x}^{t})-H^{\delta_{t}}_{\mathbf{x}^{t}}(\mathbf{W}_{t,k+1/2}).$

Under the same Assumptions as in Lemma 4.3, we have for all
$\mathbf{x}\in\mathcal{X}$:

$\begin{split}\gamma_{t}&\left(\langle
V(\mathbf{x}),\bar{\mathbf{z}}^{t}-\mathbf{x}\rangle+h(\bar{\mathbf{z}}^{t})-h(\mathbf{x})\right)\leq\frac{1}{2K}\left(\norm{\mathbf{x}^{t}-\mathbf{x}}^{2}-\norm{\mathbf{x}^{t+1}-\mathbf{x}}^{2}\right)\\\
&+2\frac{\gamma^{2}_{t}}{K}(\mathcal{L}_{f}(\xi_{t,k+1/2})^{2}+\eta^{2}_{0})\sum_{k=0}^{K-1}\norm{\mathbf{z}^{(t)}_{k+1/2}-\mathbf{x}}^{2}+4\gamma^{2}_{t}\sum_{i\in\mathcal{I}}L^{2}_{h_{i}}n^{2}_{i}\\\
&-\gamma_{t}\langle\varepsilon^{t}_{V}(\xi^{1:b_{t}}_{t})+\varepsilon^{t}_{h}(\mathbf{W}^{1:b_{t}}_{t}),\bar{\mathbf{z}}^{t}-\mathbf{x}\rangle+\gamma_{t}\delta_{t}L_{h}-\gamma_{t}\eta_{t}\langle\mathbf{x},\bar{\mathbf{z}}^{t}-\mathbf{x}\rangle\\\
&-\frac{\gamma_{t}}{K}\sum_{k=0}^{K-1}\langle\mathbf{Y}^{1}_{t,k},\mathbf{z}^{(t)}_{k+1/2}-\mathbf{x}\rangle-\frac{\gamma_{t}}{K}\sum_{k=0}^{K-1}\langle\mathbf{Y}^{2}_{t,k},\mathbf{z}^{(t)}_{k+1/2}-\mathbf{x}\rangle.\end{split}$
(6.3)

###### Proof.

We depart from the energy bound in Lemma 4.3. Rearranging this inequality and
using $\gamma_{t}\eta_{t}\in(0,1/2)$, it follows

$\displaystyle 2$ $\displaystyle\gamma_{t}\left(\langle
V(\mathbf{x}),\mathbf{z}^{(t)}_{k+1/2}-\mathbf{x}\rangle+h(\mathbf{z}^{(t)}_{k+1/2})-h(\mathbf{x})\right)\leq\norm{\mathbf{z}^{(t)}_{k}-\mathbf{x}}^{2}-\norm{\mathbf{z}^{(t)}_{k+1}-\mathbf{x}}^{2}$
$\displaystyle+8\gamma_{t}^{2}\sum_{i\in\mathcal{I}}L_{h_{i}}^{2}n^{2}_{i}+4\gamma^{2}_{t}(\mathcal{L}_{f}(\xi_{t,k+1/2})^{2}+\eta^{2}_{0})\norm{\mathbf{z}^{(t)}_{k+1/2}-\mathbf{x}^{t}}^{2}$
$\displaystyle-2\gamma_{t}\langle\hat{V}_{t,k+1/2}(\mathbf{z}^{(t)}_{k+1/2})+H^{\delta_{t}}_{\mathbf{z}^{(t)}_{k+1/2}}(\mathbf{W}_{k+1/2})-V(\mathbf{z}^{(t)}_{k+1/2})-\nabla
h^{\delta_{t}}(\mathbf{z}^{(t)}_{k+1/2}),\mathbf{z}^{(t)}_{k+1/2}-\mathbf{x}\rangle$
$\displaystyle-2\gamma_{t}\langle\left(V(\mathbf{x}^{t})+\nabla
h^{\delta_{t}}(\mathbf{x}^{t})\right)-\left(\hat{V}_{t,k+1/2}(\mathbf{x}^{t})+H^{\delta_{t}}_{\mathbf{x}^{t}}(\mathbf{W}_{t,k+1/2})\right),\mathbf{z}^{(t)}_{k+1/2}-\mathbf{x}\rangle$
$\displaystyle-2\gamma_{t}\langle\varepsilon^{t}_{V}(\xi^{1:b_{t}}_{t})+\varepsilon^{t}_{h}(\mathbf{W}^{1:b_{t}}_{t}),\mathbf{z}^{(t)}_{k+1/2}-\mathbf{x}\rangle+2\gamma_{t}\delta_{t}L_{h}-2\gamma_{t}\eta_{t}\langle\mathbf{x},\mathbf{z}^{(t)}_{k+1/2}-\mathbf{x}\rangle.$

Summing from $k=0,\ldots,K-1$ and calling
$\bar{\mathbf{z}}^{t}\triangleq\frac{1}{K}\sum_{k=0}^{K-1}\mathbf{z}^{(t)}_{k+1/2}$,
we get first from Jensen’s inequality

$2\frac{\gamma_{t}}{K}\sum_{k=0}^{K-1}\left(\langle
V(\mathbf{x}),\mathbf{z}^{(t)}_{k+1/2}-\mathbf{x}\rangle+h(\mathbf{z}^{(t)}_{k+1/2})-h(\mathbf{x})\right)\geq
2\gamma_{t}\left(\langle
V(\mathbf{x}),\bar{\mathbf{z}}^{t}-\mathbf{x}\rangle+h(\bar{\mathbf{z}}^{t})-h(\mathbf{x})\right).$

Second, telescoping the expression in the penultimate display and using the
definitions of the process $\\{\mathbf{Y}^{\nu}_{t,k}\\}_{k=0}^{K-1}$, we
deduce the bound

$\displaystyle\gamma_{t}\left(\langle
V(\mathbf{x}),\bar{\mathbf{z}}^{t}-\mathbf{x}\rangle+h(\bar{\mathbf{z}}^{t})-h(\mathbf{x})\right)$
$\displaystyle\leq\frac{1}{2K}\left(\norm{\mathbf{x}^{t}-\mathbf{x}}^{2}-\norm{\mathbf{x}^{t+1}-\mathbf{x}}^{2}\right)$
$\displaystyle+2\frac{\gamma^{2}_{t}}{K}\sum_{k=0}^{K-1}(\mathcal{L}_{f}(\xi_{t,k+1/2})^{2}+\eta^{2}_{0})\norm{\mathbf{z}_{k+1/2}-\mathbf{x}}^{2}+4\gamma^{2}_{t}\sum_{i\in\mathcal{I}}L^{2}_{h_{i}}n^{2}_{i}$
$\displaystyle-\gamma_{t}\langle\varepsilon^{t}_{V}(\xi^{1:b_{t}}_{t})+\varepsilon^{t}_{h}(\mathbf{W}^{1:b_{t}}_{t}),\bar{\mathbf{z}}^{t}-\mathbf{x}\rangle+\gamma_{t}\delta_{t}L_{h}-\gamma_{t}\eta_{t}\langle\mathbf{x},\bar{\mathbf{z}}^{t}-\mathbf{x}\rangle$
$\displaystyle-\frac{\gamma_{t}}{K}\sum_{k=0}^{K-1}\langle\mathbf{Y}^{1}_{t,k},\mathbf{z}^{(t)}_{k+1/2}-\mathbf{x}\rangle-\frac{\gamma_{t}}{K}\sum_{k=0}^{K-1}\langle\mathbf{Y}^{2}_{t,k},\mathbf{z}^{(t)}_{k+1/2}-\mathbf{x}\rangle.$

$\blacksquare$

We can now give the proof of Theorem 4.5. Let us introduce the auxiliary
processes $\\{\mathbf{u}^{\nu}_{t,k}\\}_{k=0}^{K-1}$ for $\nu=1,2$ defined
recursively as

$\mathbf{u}^{\nu}_{t,k+1}=\Pi_{\mathcal{X}}(\mathbf{u}^{\nu}_{t,k}-\gamma_{t}\mathbf{Y}^{\nu}_{t,k}),\quad\mathbf{u}^{\nu}_{t,0}=\mathbf{x}^{t}.$
(6.4)

The definition of the auxiliary sequence gives for $\nu=1,2$ (see e.g. [52])

$\displaystyle\norm{\mathbf{u}^{\nu}_{t,k+1}-\mathbf{x}}^{2}$
$\displaystyle\leq\norm{\mathbf{u}^{\nu}_{t,k}-\mathbf{x}}^{2}-2\gamma_{t}\langle\mathbf{Y}^{\nu}_{t,k},\mathbf{u}^{\nu}_{t,k}-\mathbf{x}\rangle+\gamma^{2}_{t}\norm{\mathbf{Y}^{\nu}_{t,k}}^{2}$
$\displaystyle=\norm{\mathbf{u}^{\nu}_{t,k}-\mathbf{x}}^{2}+2\gamma_{t}\langle\mathbf{Y}^{\nu}_{t,k},\mathbf{z}^{(t)}_{k+1/2}-\mathbf{x}\rangle-2\gamma_{t}\langle\mathbf{Y}^{\nu}_{t,k},\mathbf{u}^{\nu}_{t,k}-\mathbf{z}^{(t)}_{k+1/2}\rangle+\gamma^{2}_{t}\norm{\mathbf{Y}^{\nu}_{t,k}}^{2}.$

Rearranging and telescoping shows

$\displaystyle-2\gamma_{t}\sum_{k=0}^{K-1}\langle\mathbf{Y}^{\nu}_{t,k},\mathbf{z}^{(t)}_{k+1/2}-\mathbf{x}\rangle$
$\displaystyle\leq\norm{\mathbf{u}^{\nu}_{t,0}-\mathbf{x}}^{2}-\norm{\mathbf{u}^{\nu}_{t,K}-\mathbf{x}}^{2}+\gamma^{2}_{t}\sum_{k=0}^{K-1}\norm{\mathbf{Y}^{\nu}_{t,k}}^{2}$
$\displaystyle-2\gamma_{t}\sum_{k=0}^{K-1}\langle\mathbf{Y}^{\nu}_{t,k},\mathbf{z}^{(t)}_{k+1/2}-\mathbf{u}^{\nu}_{t,k}\rangle.$

Plugging this into eq. (6.3), we get

$\displaystyle\gamma_{t}$ $\displaystyle\left(\langle
V(\mathbf{x}),\bar{\mathbf{z}}^{t}-\mathbf{x}\rangle+h(\bar{\mathbf{z}}^{t})-h(\mathbf{x})\right)\leq\frac{1}{2K}\left(\norm{\mathbf{x}^{t}-\mathbf{x}}^{2}-\norm{\mathbf{x}^{t+1}-\mathbf{x}}^{2}\right)$
$\displaystyle+2\frac{\gamma^{2}_{t}}{K}\sum_{k=0}^{K-1}(\mathcal{L}_{f}(\xi_{t,k+1/2})^{2}+\eta^{2}_{0})\norm{\mathbf{z}^{(t)}_{k+1/2}-\mathbf{x}}^{2}+4\gamma^{2}_{t}\sum_{i\in\mathcal{I}}L^{2}_{h_{i}}n^{2}_{i}$
$\displaystyle-\gamma_{t}\langle\varepsilon^{t}_{V}(\xi^{1:b_{t}}_{t})+\varepsilon^{t}_{h}(\mathbf{W}^{1:b_{t}}_{t}),\bar{\mathbf{z}}^{t}-\mathbf{x}\rangle+\gamma_{t}\delta_{t}L_{h}-\gamma_{t}\eta_{t}\langle\mathbf{x},\bar{\mathbf{z}}^{t}-\mathbf{x}\rangle$
$\displaystyle+\frac{\gamma_{t}}{K}\left(\norm{\mathbf{u}^{1}_{t,0}-\mathbf{x}}^{2}-\norm{\mathbf{u}^{1}_{t,K}-\mathbf{x}}^{2}+\gamma^{2}_{t}\sum_{k=0}^{K-1}\norm{\mathbf{Y}^{1}_{t,k}}^{2}-2\gamma_{t}\sum_{k=0}^{K-1}\langle\mathbf{Y}^{1}_{t,k},\mathbf{z}^{(t)}_{k+1/2}-\mathbf{u}^{1}_{t,k}\rangle\right)$
$\displaystyle+\frac{\gamma_{t}}{K}\left(\norm{\mathbf{u}^{2}_{t,0}-\mathbf{x}}^{2}-\norm{\mathbf{u}^{2}_{t,K}-\mathbf{x}}^{2}+\gamma^{2}_{t}\sum_{k=0}^{K-1}\norm{\mathbf{Y}^{2}_{t,k}}^{2}-2\gamma_{t}\sum_{k=0}^{K-1}\langle\mathbf{Y}^{2}_{t,k},\mathbf{z}^{(t)}_{k+1/2}-\mathbf{u}^{2}_{t,k}\rangle\right)$

Summing this expression over the outer-iteration loop and introduce the
averaged iterate

$\bar{\mathbf{z}}^{T}\triangleq\frac{\sum_{t=0}^{T-1}\gamma_{t}\bar{\mathbf{z}}^{t}}{\sum_{t=0}^{T-1}\gamma_{t}}.$

Jensen’s inequality readily implies

$\sum_{t=0}^{T-1}\left(\langle
V(\mathbf{x}),\bar{\mathbf{z}}^{t}-\mathbf{x}\rangle+h(\bar{\mathbf{z}}^{t})-h(\mathbf{x})\right)\geq\left(\sum_{t=0}^{T-1}\gamma_{t}\right)\left(\langle
V(\mathbf{x}),\bar{\mathbf{z}}^{T}-\mathbf{x}\rangle+h(\bar{\mathbf{z}}^{T})-h(\mathbf{x})\right).$

Recall that $C=\max_{i\in\mathcal{I}}C_{i}$ is the upper bound on the diameter
of the set $\mathcal{X}_{i}$ (cf. Assumption 2.(i)). This assumed compactness
of the set $\mathcal{X}$, we derive the a-priori bounds

$\displaystyle\norm{\mathbf{x}^{t}-\mathbf{x}}\leq C\qquad\forall
t=0,1,\ldots,T-1,$
$\displaystyle\norm{\mathbf{z}^{(t)}_{k+1/2}-\mathbf{x}}\leq C\qquad\forall
k=0,1,\ldots,K-1,\text{ and }\norm{\bar{\mathbf{z}}^{t}-\mathbf{x}}\leq
C\qquad\forall t=0,1,\ldots,T-1.$

Using these bounds, we conclude

$\displaystyle\langle
V(\mathbf{x}),\bar{\mathbf{z}}^{T}-\mathbf{x}\rangle+h(\bar{\mathbf{z}}^{T})-h(\mathbf{x})$
$\displaystyle\leq\frac{1}{2K\sum_{t=0}^{T-1}\gamma_{t}}\left(\norm{\mathbf{x}^{0}-\mathbf{x}}^{2}-\norm{\mathbf{x}^{T}-\mathbf{x}}^{2}\right)$
$\displaystyle+\frac{2\sum_{t=0}^{T-1}\gamma^{2}_{t}\sum_{k=0}^{K-1}(\mathcal{L}_{f}(\xi_{t,k+1/2})^{2}+\eta^{2}_{0})\norm{\mathbf{z}^{(t)}_{k+1/2}-\mathbf{x}}^{2}}{K\sum_{t=0}^{T-1}\gamma_{t}}$
$\displaystyle+\frac{\sum_{t=0}^{T-1}4\gamma^{2}_{t}\sum_{i\in\mathcal{I}}L^{2}_{h_{i}}n^{2}_{i}+L_{h}\sum_{t=0}^{T-1}\delta_{t}\gamma_{t}}{\sum_{t=0}^{T-1}\gamma_{t}}$
$\displaystyle+\frac{1}{K\sum_{s=0}^{T-1}\gamma_{s}}\sum_{t=0}^{T-1}\gamma_{t}\left(\norm{\mathbf{u}^{1}_{t,0}-\mathbf{x}}^{2}+\norm{\mathbf{u}^{2}_{t,0}-\mathbf{x}}^{2}\right)$
$\displaystyle+\frac{\sum_{t=0}^{T-1}\gamma^{3}_{t}}{K\sum_{s=0}^{T-1}\gamma_{s}}\sum_{k=0}^{K-1}\left(\norm{\mathbf{Y}^{1}_{t,k}}^{2}+\norm{\mathbf{Y}^{2}_{t,k}}^{2}\right)$
$\displaystyle+\frac{2\sum_{t=0}^{T-1}\gamma^{2}_{t}}{K\sum_{s=0}^{T-1}\gamma_{s}}\sum_{k=0}^{K-1}\left(\langle\mathbf{Y}^{1}_{t,k},\mathbf{u}^{1}_{t,k}-\mathbf{z}^{(t)}_{k+1/2}\rangle+\langle\mathbf{Y}^{2}_{t,k},\mathbf{u}_{t,k}^{2}-\mathbf{z}^{(t)}_{k+1/2}\rangle\right)$
$\displaystyle+\frac{\sum_{t=0}^{T-1}\gamma_{t}\langle\varepsilon^{t}_{V}(\xi^{1:b_{t}}_{t})+\varepsilon^{t}_{h}(\mathbf{W}^{1:b_{t}}_{t}),\mathbf{x}-\bar{\mathbf{z}}^{t}\rangle}{\sum_{t=0}^{T-1}\gamma_{t}}$
$\displaystyle+\frac{\sum_{t=0}^{T}\gamma_{t}\eta_{t}\langle\mathbf{x},\mathbf{x}-\bar{\mathbf{z}}^{t}\rangle}{\sum_{t=0}^{T-1}\gamma_{t}}$

We bound each of the terms above individually as follows:

  1. 1.

$\frac{1}{2K\sum_{t=0}^{T-1}\gamma_{t}}\left(\norm{\mathbf{x}^{0}-\mathbf{x}}^{2}-\norm{\mathbf{x}^{T}-\mathbf{x}}^{2}\right)\leq\frac{C^{2}}{2K\sum_{t=0}^{T-1}\gamma_{t}},$

  2. 2.

$\frac{2\sum_{t=0}^{T-1}\gamma^{2}_{t}\sum_{k=0}^{K-1}(\mathcal{L}_{f}(\xi_{t,k+1/2})^{2}+\eta^{2}_{0})\norm{\mathbf{z}^{(t)}_{k+1/2}-\mathbf{x}}^{2}}{K\sum_{t=0}^{T-1}\gamma_{t}}\leq
2C^{2}\frac{\sum_{t=0}^{T-1}\gamma^{2}_{t}\sum_{k=0}^{K-1}(\mathcal{L}_{f}(\xi_{t,k+1/2})^{2}+\eta^{2}_{0})}{K\sum_{s=0}^{T-1}\gamma_{s}},$

  3. 3.

$\frac{1}{K}\left(\norm{\mathbf{u}^{1}_{t,0}-\mathbf{x}}^{2}+\norm{\mathbf{u}^{2}_{t,0}-\mathbf{x}}^{2}\right)=\frac{2}{K}\norm{\mathbf{x}^{t}-\mathbf{x}}^{2}\leq\frac{2C^{2}}{K},$

  4. 4.

$\frac{\sum_{t=0}^{T}\gamma_{t}\eta_{t}\langle\mathbf{x},\mathbf{x}-\bar{\mathbf{z}}^{t}\rangle}{\sum_{t=0}^{T-1}\gamma_{t}}\leq\frac{3C^{2}\sum_{t=0}^{T-1}\gamma_{t}\eta_{t}}{2\sum_{t=0}^{T-1}\gamma_{t}},$

  5. 5.

$\frac{\sum_{t=0}^{T-1}\gamma_{t}\langle\varepsilon^{t}_{V}(\xi^{1:b_{t}}_{t})+\varepsilon^{t}_{h}(\mathbf{W}^{1:b_{t}}_{t}),\mathbf{x}-\bar{\mathbf{z}}^{t}\rangle}{\sum_{t=0}^{T-1}\gamma_{t}}\leq\frac{\sum_{t=0}^{T-1}\gamma_{t}(C\norm{\varepsilon^{t}_{V}(\xi^{1:b_{t}}_{t})}+C\norm{\varepsilon^{t}_{h}(\mathbf{W}^{1:b_{t}}_{t})})}{\sum_{t=0}^{T-1}\gamma_{t}}$

Plugging all these bounds into the penultimate display gives

$\displaystyle\langle
V(\mathbf{x}),\bar{\mathbf{z}}^{T}-\mathbf{x}\rangle+h(\bar{\mathbf{z}}^{T})-h(\mathbf{x})$
$\displaystyle\leq\frac{C^{2}}{2K\sum_{t=0}^{T-1}\gamma_{t}}+2C^{2}\frac{\sum_{t=0}^{T-1}\gamma^{2}_{t}\sum_{k=0}^{K-1}(\mathcal{L}_{f}(\xi_{t,k+1/2})^{2}+\eta^{2}_{0})}{K\sum_{s=0}^{T-1}\gamma_{s}}$
$\displaystyle+\frac{\sum_{t=0}^{T-1}4\gamma^{2}_{t}\sum_{i\in\mathcal{I}}L^{2}_{h_{i}}n^{2}_{i}+L_{h}\sum_{t=0}^{T-1}\delta_{t}\gamma_{t}}{\sum_{t=0}^{T-1}\gamma_{t}}$
$\displaystyle+\frac{2C^{2}}{K}+\frac{\sum_{t=0}^{T-1}\gamma_{t}(C\norm{\varepsilon^{t}_{V}(\xi^{1:b_{t}}_{t})}+C\norm{\varepsilon^{t}_{h}(\mathbf{W}^{1:b_{t}}_{t})})}{\sum_{t=0}^{T-1}\gamma_{t}}$
$\displaystyle+\frac{\sum_{t=0}^{T-1}\gamma^{3}_{t}}{K\sum_{t=0}^{T-1}\gamma_{t}}\sum_{k=0}^{K-1}\left(\norm{\mathbf{Y}^{1}_{t,k}}^{2}+\norm{\mathbf{Y}^{2}_{t,k}}^{2}\right)$
$\displaystyle+\frac{2\sum_{t=0}^{T-1}\gamma^{2}_{t}}{K\sum_{t=0}^{T-1}\gamma_{t}}\sum_{k=0}^{K-1}\left(\langle\mathbf{Y}^{1}_{t,k},\mathbf{u}^{1}_{t,k}-\mathbf{z}^{(t)}_{k+1/2}\rangle+\langle\mathbf{Y}^{2}_{t,k},\mathbf{u}_{t,k}^{2}-\mathbf{z}_{k+1/2}\rangle\right)$
$\displaystyle+\frac{3C^{2}\sum_{t=0}^{T-1}\gamma_{t}\eta_{t}}{2\sum_{t=0}^{T-1}\gamma_{t}}.$

Let $L_{f}^{2}\triangleq\mathbb{E}_{\xi}[\mathcal{L}_{f}(\xi)^{2}]$, and
define

$\begin{split}\mathcal{D}_{K,T}&\triangleq\frac{C^{2}}{2K\sum_{t=0}^{T-1}\gamma_{t}}+2C^{2}(L_{f}^{2}+\eta^{2}_{0})\frac{\sum_{t=0}^{T-1}\gamma^{2}_{t}}{\sum_{t=0}^{T-1}\gamma_{t}}+\frac{3C^{2}\sum_{t=0}^{T-1}\gamma_{t}\eta_{t}}{2\sum_{t=0}^{T-1}\gamma_{t}}\\\
&+\frac{\sum_{t=0}^{T-1}4\gamma^{2}_{t}\sum_{i\in\mathcal{I}}L^{2}_{h_{i}}n^{2}_{i}+L_{h}\sum_{t=0}^{T-1}\delta_{t}\gamma_{t}}{\sum_{t=0}^{T-1}\gamma_{t}}+\frac{2C^{2}}{K}.\end{split}$
(6.5)

Hence, using the definition of the gap function (4.1), we see

$\displaystyle\mathbb{E}[\Gamma(\bar{\mathbf{z}}^{T})]$
$\displaystyle\leq\mathcal{D}_{K,T}+\mathbb{E}\left[\frac{\sum_{t=0}^{T-1}\gamma^{3}_{t}}{K\sum_{t=0}^{T-1}\gamma_{t}}\sum_{k=0}^{K-1}\left(\norm{\mathbf{Y}^{1}_{t,k}}^{2}+\norm{\mathbf{Y}^{2}_{t,k}}^{2}\right)\right]$
$\displaystyle+\mathbb{E}\left[\frac{2\sum_{t=0}^{T-1}\gamma^{2}_{t}}{K\sum_{t=0}^{T-1}\gamma_{t}}\sum_{k=0}^{K-1}\left(\langle\mathbf{Y}^{1}_{t,k},\mathbf{u}^{1}_{t,k}-\mathbf{z}^{(t)}_{k+1/2}\rangle+\langle\mathbf{Y}^{2}_{t,k},\mathbf{u}_{t,k}^{2}-\mathbf{z}^{(t)}_{k+1/2}\rangle\right)\right]$
$\displaystyle+\mathbb{E}\left[\frac{\sum_{t=0}^{T-1}\gamma_{t}(C\norm{\varepsilon^{t}_{V}(\xi^{1:b_{t}}_{t})}+C\norm{\varepsilon^{t}_{h}(\mathbf{W}^{1:b_{t}}_{t})})}{\sum_{t=0}^{T-1}\gamma_{t}}\right].$

Next, observe that

$\displaystyle\norm{\mathbf{Y}^{1}_{t,k}}\leq\norm{\hat{V}_{t,k+1/2}(\mathbf{z}^{(t)}_{k+1/2})-V(\mathbf{z}_{k+1/2})}+\norm{H^{\delta_{t}}_{\mathbf{z}^{(t)}_{k+1/2}}(\mathbf{W}_{t,k+1/2})-\nabla
h^{\delta_{t}}(\mathbf{z}^{(t)}_{k+1/2})},\text{ and }$
$\displaystyle\norm{\mathbf{Y}^{2}_{t,k}}\leq\norm{\hat{V}_{t,k+1/2}(\mathbf{x}^{t})-V(\mathbf{x}^{t})}+\norm{H^{\delta_{t}}_{\mathbf{x}^{t}}(\mathbf{W}_{t,k+1/2})-\nabla
h^{\delta_{t}}(\mathbf{x}^{t})}.$

Moreover, using compactness of $\mathcal{X}$,

$\langle\mathbf{Y}^{\nu}_{t,k},\mathbf{u}^{\nu}_{t,k}-\mathbf{z}^{(t)}_{k+1/2}\rangle\leq
C\norm{\mathbf{Y}^{\nu}_{t,k}}\quad\forall\nu=1,2.$ (6.6)

Lemma A.12(c) (for $b=1$) gives

$\mathbb{E}_{\mathbf{W}\sim\mathsf{U}(\mathbb{S}_{n})}\left[\norm{H^{\delta_{t}}_{\mathbf{x}^{t}}(\mathbf{W}_{t,k+1/2})-\nabla
h^{\delta_{t}}(\mathbf{x}^{t})}^{2}|\mathcal{A}_{t,k}\right]\leq\sum_{i\in\mathcal{I}}n^{2}_{i}L^{2}_{h_{i}}.$

Lemma 4.1 in turn implies

$\mathbb{E}_{\xi}\left[\norm{\hat{V}_{t,k+1/2}(\mathbf{z}^{(t)}_{k+1/2})-V(\mathbf{z}^{(t)}_{k+1/2})}^{2}|\mathcal{A}_{t,k}\right]\leq
M^{2}_{V}.$

By Jensen’s inequality in tandem with Lemma 4.1 and Lemma 4.2, we conclude
from (6.1) and (6.2)

$\displaystyle\mathbb{E}\left[\norm{\varepsilon_{V}^{t}(\xi^{1:b_{t}}_{t})}|\mathcal{F}_{t}\right]$
$\displaystyle\leq\sqrt{\mathbb{E}\left[\norm{\varepsilon_{V}^{t}(\xi^{1:b_{t}}_{t})}^{2}|\mathcal{F}_{t}\right]}\leq\frac{M_{V}}{\sqrt{b_{t}}},\text{
and }$
$\displaystyle\mathbb{E}\left[\norm{\varepsilon_{h}^{t}(\mathbf{W}^{1:b_{t}}_{t})}|\mathcal{F}_{t}\right]$
$\displaystyle\leq\sqrt{\mathbb{E}\left[\norm{\varepsilon_{h}^{t}(\mathbf{W}^{1:b_{t}}_{t})}^{2}|\mathcal{F}_{t}\right]}\leq\frac{\left(\sum_{i\in\mathcal{I}}L^{2}_{h_{i}}n^{2}_{i}\right)^{1/2}}{\sqrt{b_{t}}}.$

This implies
$\mathbb{E}[\norm{\mathbf{Y}^{\nu}_{t,k}}^{2}|\mathcal{A}_{t,k}]\leq
2M^{2}_{V}+2\sum_{i\in\mathcal{I}}L^{2}_{h_{i}}n^{2}_{i}\equiv\sigma^{2}$ for
all $\nu=1,2$, as well as

$\mathbb{E}[\norm{\mathbf{Y}^{\nu}_{t,k}}|\mathcal{A}_{t,k}]\leq\sqrt{\mathbb{E}[\norm{\mathbf{Y}^{\nu}_{t,k}}^{2}|\mathcal{A}_{t,k}]}\leq\sigma.$

We conclude, via a repeated application of the law of iterated expectations,
that

$\displaystyle\mathbb{E}[\Gamma(\bar{\mathbf{z}}^{T})]$
$\displaystyle\leq\mathcal{D}_{K,T}+\frac{2\sigma^{2}\sum_{t=0}^{T-1}\gamma^{3}_{t}}{\sum_{t=0}^{T-1}\gamma_{t}}+\frac{4C\sigma\sum_{t=0}^{T-1}\gamma^{2}_{t}}{\sum_{t=0}^{T-1}\gamma_{t}}$
$\displaystyle+\frac{\sum_{t=0}^{T-1}\frac{\gamma_{t}}{\sqrt{b_{t}}}\left(CM_{V}+C(\sum_{i\in\mathcal{I}}L_{h_{i}}^{2}n_{i}^{2})^{1/2}\right)}{\sum_{t=0}^{T-1}\gamma_{t}}.$

Using the specification $\gamma_{t}=1/T=\eta_{t}=\delta_{t}$, as well as $K=T$
and $b_{t}\geq T^{2}$ gives

$\displaystyle\mathcal{D}_{K,T}\leq\frac{5C^{2}}{T}+\frac{2C^{2}(L^{2}_{f}+1/T^{2})}{T}+\frac{3C^{2}}{2T}+\frac{4(\sum_{i\in\mathcal{I}}L^{2}_{h_{i}}n^{2}_{i})+L_{h}}{T}\equiv\mathtt{c}_{T}$

and consequently,

$\mathbb{E}[\Gamma(\bar{\mathbf{z}}^{T})]\leq\mathtt{c}_{T}+\frac{2\sigma^{2}}{T^{2}}+\frac{C\sigma}{T}+\frac{C(M_{V}+(\sum_{i\in\mathcal{I}}L^{2}_{h_{i}}n^{2}_{i})^{1/2})}{T}=\mathcal{O}(C\sigma/T)$

###  6.2 Analysis of the inexact scheme

The inexact version of our method $\operatorname{\texttt{VRHGS}}$ is obtained
by replacing the estimates for the implicit function using the inexact
solution map $\mathbf{y}_{i}^{\varepsilon}$. The precise implementation is
summarized in Alg. 3 and Alg. 4.

Result: Iterate $\mathbf{z}^{K}$

Set $\mathbf{z}^{0}=\bar{\mathbf{x}}$;

for  _$k=0,1,\ldots,K-1$_ do

Update
$\mathbf{z}_{k+1/2}=\Pi_{\mathcal{X}}[\mathbf{z}_{k}-\gamma(\bar{\mathbf{v}}+\eta\bar{\mathbf{x}}+\bar{H})]$,;

Obtain $\hat{V}^{\eta}_{k+1/2}(\mathbf{z}_{k+1/2})$ and
$\hat{V}^{\eta}_{k+1/2}(\bar{\mathbf{x}})$ as defined in eq. (3.7);

Draw iid direction vectors
$\mathbf{W}_{k+1/2}=\\{\mathbf{W}_{i,k+1/2}\\}_{i\in\mathcal{I}}$, with each
$\mathbf{W}_{i,k+1/2}\sim\mathsf{U}(\mathbb{S}_{i})$. ;

Obtain $H^{\delta,\varepsilon}_{\mathbf{z}_{k+1/2}}(\mathbf{W}_{k+1/2})$ and
$H^{\delta,\varepsilon}_{\bar{\mathbf{x}}}(\mathbf{W}_{k+1/2})$;

Update
$\mathbf{z}_{k+1}=\mathbf{z}_{k+1/2}-\gamma\left(\hat{V}^{\eta}_{k+1/2}(\mathbf{z}_{k+1/2})+H^{\delta,\varepsilon}_{\mathbf{z}_{k+1/2}}(\mathbf{W}_{k+1/2})-\hat{V}^{\eta}_{k+1/2}(\bar{\mathbf{x}})-H^{\delta,\varepsilon}_{\bar{\mathbf{x}}}(\mathbf{W}_{k+1/2})\right).$
;

end for

Algorithm 3
$\operatorname{\texttt{ISFBF}}(\bar{\mathbf{x}},\bar{\mathbf{v}},\bar{H},\gamma,\eta,\delta,\varepsilon,K)$

Data:
$\mathbf{x},T,\\{\gamma_{t}\\}_{t=0}^{T},\\{b_{t}\\}_{t=0}^{T},\\{\eta_{t}\\}_{t=0}^{T},\\{\varepsilon_{t}\\}_{t=0}^{T}$

Set $\mathbf{x}^{0}=\mathbf{x}$.

for  _$t=0,1,\ldots,T-1$_ do

For each $i\in\mathcal{I}$ receive the oracle feedback $\bar{V}^{t}$ defined
by
$\bar{V}^{t}_{i}\triangleq\frac{1}{b_{t}}\sum_{s=1}^{b_{t}}\hat{V}_{i}(\mathbf{x}^{t},\xi_{i,t}^{(s)})$.
;

For each $i\in\mathcal{I}$ construct the estimator
$H^{\delta_{t},b_{t}}_{\mathbf{x}^{t}}$ defined by
$H_{i,\mathbf{x}_{i}}^{\delta_{t},\varepsilon_{t},b_{t}}\triangleq\frac{1}{b_{t}}\sum_{s=1}^{b_{t}}H_{i,\mathbf{x}_{i}}^{\delta_{t},\varepsilon_{t}}(\mathbf{W}_{i,t}^{(s)})$.
;

Update
$\mathbf{x}^{t+1}=\operatorname{\texttt{ISFBF}}(\mathbf{x}^{t},\bar{V}^{t},H^{\delta_{t},b_{t},\varepsilon_{t}}_{\mathbf{x}^{t}},\gamma_{t},\eta_{t},\delta_{t},\varepsilon_{t},b_{t},K)$

end for

Algorithm 4 Inexact Variance Reduced Hierarchical Game Solver
(I-$\operatorname{\texttt{VRHGS}}$)

The proof of Theorem 4.7 is analogous to the one of Theorem 4.5, with the
simple modification due to inexact feedback from the follower’s problem. We
state the main changes here, leaving the straightforward derivations to the
reader. We begin with the modified energy inequality, similar to Lemma 4.3.
Here, we also follow the same notational simplification by suppressing the
outer iteration counter $t$ from the variables.

######  Lemma 6.2.

Let Assumptions 1,2,3 and 4 hold true. Then, for all
$t\in\\{0,1,\ldots,T-1\\}$ and all anchor points $\mathbf{x}\in\mathcal{X}$,
we have

$\displaystyle\norm{\mathbf{z}_{k+1}-\mathbf{x}}^{2}$
$\displaystyle\leq(1-\gamma_{t}\eta_{t})\norm{\mathbf{z}_{k}-\mathbf{x}}^{2}-(1-2\gamma_{t}\eta_{t})\norm{\mathbf{z}_{k+1/2}-\mathbf{z}_{k}}^{2}$
$\displaystyle+8\gamma_{t}^{2}\left(\norm{H^{\delta_{t},\varepsilon_{t}}_{\mathbf{z}_{k+1/2}}(\mathbf{W}_{k+1/2})}^{2}+\norm{H^{\delta_{t},\varepsilon_{t}}_{\mathbf{x}^{t}}(\mathbf{W}_{k+1/2})}^{2}\right)$
$\displaystyle+4\gamma^{2}_{t}(\mathcal{L}_{f}(\xi_{k+1/2})^{2}+\eta^{2}_{0})\norm{\mathbf{z}_{k+1/2}-\mathbf{x}^{t}}^{2}$
$\displaystyle-2\gamma_{t}\langle\hat{V}_{k+1/2}(\mathbf{z}_{k+1/2})+H^{\delta_{t}}_{\mathbf{z}_{k+1/2}}(\mathbf{W}_{k+1/2})-V(\mathbf{z}_{k+1/2})-\nabla
h^{\delta_{t}}(\mathbf{z}_{k+1/2}),\mathbf{z}_{k+1/2}-\mathbf{x}\rangle$
$\displaystyle-2\gamma_{t}\langle\left(V(\mathbf{x}^{t})+\nabla
h^{\delta_{t}}(\mathbf{x}^{t})\right)-\left(\hat{V}_{k+1/2}(\mathbf{x}^{t})+H^{\delta_{t}}_{\mathbf{x}^{t}}(\mathbf{W}_{k+1/2})\right),\mathbf{z}_{k+1/2}-\mathbf{x}\rangle$
$\displaystyle-2\gamma_{t}\langle\varepsilon^{t}_{V}(\xi^{1:b_{t}})+\varepsilon^{t}_{h}(\mathbf{W}^{1:b_{t}}),\mathbf{z}_{k+1/2}-\mathbf{x}\rangle$
$\displaystyle-2\gamma_{t}\langle
H^{\delta_{t},b_{t}}_{\mathbf{x}^{t}}-H^{\delta_{t},b_{t},\varepsilon_{t}}_{\mathbf{x}^{t}}+H^{\delta_{t},\varepsilon_{t}}_{\mathbf{x}^{t}}(\mathbf{W}_{k+1/2})-H^{\delta_{t}}_{\mathbf{x}^{t}}(\mathbf{W}_{k+1/2}),\mathbf{z}_{k+1/2}-\mathbf{x}\rangle$
$\displaystyle-2\gamma_{t}\langle
H^{\delta_{t},\varepsilon_{t}}_{\mathbf{z}_{k+1/2}}(\mathbf{W}_{k+1/2})-H^{\delta_{t}}_{\mathbf{z}_{k+1/2}}(\mathbf{W}_{k+1/2}),\mathbf{z}_{k+1/2}-\mathbf{x}\rangle$
$\displaystyle-2\gamma_{t}\left(\langle
V^{\eta_{t}}(\mathbf{x}),\mathbf{z}_{k+1/2}-\mathbf{x}\rangle+h(\mathbf{z}_{k+1/2})-h(\mathbf{x})\right)+2\gamma_{t}\delta_{t}L_{h},$

where $L_{h}\triangleq\sum_{i\in\mathcal{I}}L_{h_{i}}$.

The inexact version of Lemma 6.1 reads then as follows.

######  Lemma 6.3.

For any $t\in\\{0,1,\ldots,T-1\\}$, define
$\bar{\mathbf{z}}^{t}\triangleq\frac{1}{K}\sum_{k=0}^{K-1}\mathbf{z}^{(t)}_{k+1/2}$.
Let $\\{\gamma_{t}\\}_{t},\\{\eta_{t}\\}_{t}$ be positive sequences satisfying
$0<\gamma_{t}\eta_{t}<1/2$. For $k\in\\{0,1,\ldots,K-1\\}$ define the process
$\\{\mathbf{Y}^{\nu}_{t,k}\\}_{k=0}^{K-1},\nu\in\\{1,2\\}$ as in Lemma 6.1.
Then, we have for all $\mathbf{x}\in\mathcal{X}$:

$\displaystyle\gamma_{t}$ $\displaystyle\left(\langle
V(\mathbf{x}),\bar{\mathbf{z}}^{t}-\mathbf{x}\rangle+h(\bar{\mathbf{z}}^{t})-h(\mathbf{x})\right)\leq\frac{1}{2K}\left(\norm{\mathbf{x}^{t}-\mathbf{x}}^{2}-\norm{\mathbf{x}^{t+1}-\mathbf{x}}^{2}\right)$
$\displaystyle+\frac{2\gamma^{2}_{t}}{K}\sum_{k=0}^{K-1}(\mathcal{L}_{f}(\xi_{t,k+1/2})^{2}+\eta^{2}_{0})\norm{\mathbf{z}^{(t)}_{k+1/2}-\mathbf{x}}^{2}-\gamma_{t}\eta_{t}\langle\mathbf{x},\bar{\mathbf{z}}^{t}-\mathbf{x}\rangle$
$\displaystyle-\gamma_{t}\langle\varepsilon^{t}_{V}(\xi^{1:b_{t}}_{t})+\varepsilon^{t}_{h}(\mathbf{W}^{1:b_{t}}_{t}),\bar{\mathbf{z}}^{t}-\mathbf{x}\rangle+\gamma_{t}\delta_{t}L_{h}$
$\displaystyle-\frac{\gamma_{t}}{K}\sum_{k=0}^{K-1}\langle\mathbf{Y}^{1}_{t,k},\mathbf{z}^{(t)}_{k+1/2}-\mathbf{x}\rangle-\frac{\gamma_{t}}{K}\sum_{k=0}^{K-1}\langle\mathbf{Y}^{2}_{t,k},\mathbf{z}^{(t)}_{k+1/2}-\mathbf{x}\rangle$
$\displaystyle+\frac{2\gamma^{2}_{t}}{K}\sum_{k=0}^{K-1}\left(\norm{H^{\delta_{t},\varepsilon_{t}}_{\mathbf{z}^{(t)}_{k+1/2}}(\mathbf{W}_{t,k+1/2})}^{2}+\norm{H^{\delta_{t},\varepsilon_{t}}_{\mathbf{x}^{t}}(\mathbf{W}_{t,k+1/2})}^{2}\right)$
$\displaystyle-\frac{\gamma_{t}}{K}\sum_{k=0}^{K-1}\langle
H^{\delta_{t},b_{t}}_{\mathbf{x}^{t}}-H^{\delta_{t},b_{t},\varepsilon_{t}}_{\mathbf{x}^{t}}+H^{\delta_{t},\varepsilon_{t}}_{\mathbf{x}^{t}}(\mathbf{W}_{t,k+1/2})-H^{\delta_{t}}_{\mathbf{x}^{t}}(\mathbf{W}_{t,k+1/2}),\mathbf{z}^{(t)}_{k+1/2}-\mathbf{x}\rangle$
$\displaystyle-\frac{\gamma_{t}}{K}\sum_{k=0}^{K-1}\langle
H^{\delta_{t},\varepsilon_{t}}_{\mathbf{z}^{(t)}_{k+1/2}}(\mathbf{W}_{t,k+1/2})-H^{\delta_{t}}_{\mathbf{z}^{(t)}_{k+1/2}}(\mathbf{W}_{t,k+1/2}),\mathbf{z}^{(t)}_{k+1/2}-\mathbf{x}\rangle.$

Using this bound, we conclude in the same way as in the analysis of the exact
scheme that

$\displaystyle\langle
V(\mathbf{x}),\bar{\mathbf{z}}^{T}-\mathbf{x}\rangle+h(\bar{\mathbf{z}}^{T})-h(\mathbf{x})\leq\frac{1}{2K\sum_{t=0}^{T-1}\gamma_{t}}\norm{\mathbf{x}^{0}-\mathbf{x}}^{2}+\frac{\sum_{t=0}^{T}\gamma_{t}\eta_{t}\langle\mathbf{x},\mathbf{x}-\bar{\mathbf{z}}^{t}\rangle}{\sum_{t=0}^{T-1}\gamma_{t}}+\frac{L_{h}\sum_{t=0}^{T-1}\delta_{t}\gamma_{t}}{\sum_{t=0}^{T-1}\gamma_{t}}$
$\displaystyle+\frac{2\sum_{t=0}^{T-1}\gamma^{2}_{t}\sum_{k=0}^{K-1}(\mathcal{L}_{f}(\xi_{t,k+1/2})^{2}+\eta^{2}_{0})\norm{\mathbf{z}^{(t)}_{k+1/2}-\mathbf{x}}^{2}}{K\sum_{s=0}^{T-1}\gamma_{s}}$
$\displaystyle+\frac{1}{K\sum_{s=0}^{T-1}\gamma_{s}}\sum_{t=0}^{T-1}\gamma_{t}\left(\norm{\mathbf{u}^{1}_{t,0}-\mathbf{x}}^{2}+\norm{\mathbf{u}^{2}_{t,0}-\mathbf{x}}^{2}\right)$
$\displaystyle+\frac{\sum_{t=0}^{T-1}\gamma^{3}_{t}}{K\sum_{t=0}^{T-1}\gamma_{t}}\sum_{k=0}^{K-1}\left(\norm{\mathbf{Y}^{1}_{t,k}}^{2}+\norm{\mathbf{Y}^{2}_{t,k}}^{2}\right)$
$\displaystyle+\frac{2\sum_{t=0}^{T-1}\gamma^{2}_{t}}{K\sum_{t=0}^{T-1}\gamma_{t}}\sum_{k=0}^{K-1}\left(\langle\mathbf{Y}^{1}_{t,k},\mathbf{u}^{1}_{t,k}-\mathbf{z}^{(t)}_{k+1/2}\rangle+\langle\mathbf{Y}^{2}_{t,k},\mathbf{u}_{t,k}^{2}-\mathbf{z}^{(t)}_{k+1/2}\rangle\right)$
$\displaystyle+\frac{\sum_{t=0}^{T-1}\gamma_{t}\langle\varepsilon^{t}_{V}(\xi^{1:b_{t}}_{t})+\varepsilon^{t}_{h}(\mathbf{W}^{1:b_{t}}_{t}),\mathbf{x}-\bar{\mathbf{z}}^{t}\rangle}{\sum_{t=0}^{T-1}\gamma_{t}}$
$\displaystyle+\sum_{t=0}^{T-1}\frac{2\gamma^{2}_{t}}{K\sum_{s=0}^{T-1}\gamma_{s}}\sum_{k=0}^{K-1}\left(\norm{H^{\delta_{t},\varepsilon_{t}}_{\mathbf{z}^{(t)}_{k+1/2}}(\mathbf{W}_{t,k+1/2})}^{2}+\norm{H^{\delta_{t},\varepsilon_{t}}_{\mathbf{x}^{t}}(\mathbf{W}_{t,k+1/2})}^{2}\right)$
$\displaystyle-\sum_{t=0}^{T-1}\frac{\gamma_{t}}{K\sum_{s=0}^{T-1}\gamma_{s}}\sum_{k=0}^{K-1}\langle
H^{\delta_{t},b_{t}}_{\mathbf{x}^{t}}-H^{\delta_{t},b_{t},\varepsilon_{t}}_{\mathbf{x}^{t}}+H^{\delta_{t},\varepsilon_{t}}_{\mathbf{x}^{t}}(\mathbf{W}_{t,k+1/2})-H^{\delta_{t}}_{\mathbf{x}^{t}}(\mathbf{W}_{t,k+1/2}),\mathbf{z}^{(t)}_{k+1/2}-\mathbf{x}\rangle$
$\displaystyle-\sum_{t=0}^{T-1}\frac{\gamma_{t}}{K\sum_{s=0}^{T-1}\gamma_{s}}\sum_{k=0}^{K-1}\langle
H^{\delta_{t},\varepsilon_{t}}_{\mathbf{z}^{(t)}_{k+1/2}}(\mathbf{W}_{t,k+1/2})-H^{\delta_{t}}_{\mathbf{z}^{(t)}_{k+1/2}}(\mathbf{W}_{t,k+1/2}),\mathbf{z}^{(t)}_{k+1/2}-\mathbf{x}\rangle.$

Performing the same bounding steps as done in the exact case, we readily
arrive at the expression

$\displaystyle\langle
V(\mathbf{x}),\bar{\mathbf{z}}^{T}-\mathbf{x}\rangle+h(\bar{\mathbf{z}}^{T})-h(\mathbf{x})\leq\frac{C^{2}}{2K\sum_{t=0}^{T-1}\gamma_{t}}+2C^{2}\frac{\sum_{t=0}^{T-1}\gamma^{2}_{t}\sum_{k=0}^{K-1}(\mathcal{L}_{f}(\xi_{t,k+1/2})^{2}+\eta^{2}_{0})}{K\sum_{t=0}^{T-1}\gamma_{t}}$
$\displaystyle+\frac{3C^{2}\sum_{t=0}^{T-1}\gamma_{t}\eta_{t}}{2\sum_{t=0}^{T-1}\gamma_{t}}+\frac{2C^{2}}{K}+\frac{L_{h}\sum_{t=0}^{T-1}\delta_{t}\gamma_{t}}{\sum_{t=0}^{T-1}\gamma_{t}}+\frac{\sum_{t=0}^{T-1}\gamma_{t}\left(C\norm{\varepsilon^{t}_{V}(\xi^{1:b_{t}}_{t})}+C\norm{\varepsilon^{t}_{h}(\mathbf{W}^{1:b_{t}}_{t})}\right)}{\sum_{t=0}^{T-1}\gamma_{t}}$
$\displaystyle+\frac{\sum_{t=0}^{T-1}\gamma^{3}_{t}}{K\sum_{t=0}^{T-1}\gamma_{t}}\sum_{k=0}^{K-1}\left(\norm{\mathbf{Y}^{1}_{t,k}}^{2}+\norm{\mathbf{Y}^{2}_{t,k}}^{2}\right)+\frac{2\sum_{t=0}^{T-1}\gamma^{2}_{t}}{K\sum_{t=0}^{T-1}\gamma_{t}}\sum_{k=0}^{K-1}\left(\langle\mathbf{Y}^{1}_{t,k},\mathbf{u}^{1}_{t,k}-\mathbf{z}^{(t)}_{k+1/2}\rangle+\langle\mathbf{Y}^{2}_{t,k},\mathbf{u}_{t,k}^{2}-\mathbf{z}^{(t)}_{k+1/2}\rangle\right)$
$\displaystyle+\sum_{t=0}^{T-1}\frac{2\gamma^{2}_{t}}{K\sum_{s=0}^{T-1}\gamma_{s}}\sum_{k=0}^{K-1}\left(\norm{H^{\delta_{t},\varepsilon_{t}}_{\mathbf{z}^{(t)}_{k+1/2}}(\mathbf{W}_{t,k+1/2})}^{2}+\norm{H^{\delta_{t},\varepsilon_{t}}_{\mathbf{x}^{t}}(\mathbf{W}_{t,k+1/2})}^{2}\right)$
$\displaystyle+\sum_{t=0}^{T-1}\frac{C\gamma_{t}}{K\sum_{s=0}^{T-1}\gamma_{s}}\sum_{k=0}^{K-1}\norm{H^{\delta_{t},b_{t}}_{\mathbf{x}^{t}}-H^{\delta_{t},b_{t},\varepsilon_{t}}_{\mathbf{x}^{t}}+H^{\delta_{t},\varepsilon_{t}}_{\mathbf{x}^{t}}(\mathbf{W}_{t,k+1/2})-H^{\delta_{t}}_{\mathbf{x}^{t}}(\mathbf{W}_{t,k+1/2})}$
$\displaystyle+\sum_{t=0}^{T-1}\frac{C\gamma_{t}}{K\sum_{s=0}^{T-1}\gamma_{s}}\sum_{k=0}^{K-1}\norm{H^{\delta_{t},\varepsilon_{t}}_{\mathbf{z}_{k+1/2}}(\mathbf{W}_{t,k+1/2})-H^{\delta_{t}}_{\mathbf{z}^{(t)}_{k+1/2}}(\mathbf{W}_{t,k+1/2})}$

We next estimate the error terms appearing because of the inexact feedback map
in the coupling function. Lemma A.14 yields

$\mathbb{E}\left[\norm{H^{\delta_{t},\varepsilon_{t}}_{\mathbf{z}_{k+1/2}}(\mathbf{W}_{t,k+1/2})}^{2}+\norm{H^{\delta_{t},\varepsilon_{t}}_{\mathbf{x}^{t}}(\mathbf{W}_{t,k+1/2})}^{2}|\mathcal{A}_{t,k}\right]\leq
6\sum_{i\in\mathcal{I}}\left(\frac{n_{i}}{\delta_{t}}\right)^{2}(2L^{2}_{2,i}\varepsilon^{2}_{t}+L_{1,i}^{2}\delta^{2}_{t})\equiv\alpha^{(1)}_{t}.$

Furthermore, the triangle inequality, Jensen’s inequality, and (A.11) gives

$\displaystyle\mathbb{E}\left[\norm{H^{\delta_{t},b_{t}}_{\mathbf{x}^{t}}-H^{\delta_{t},b_{t},\varepsilon_{t}}_{\mathbf{x}^{t}}+H^{\delta_{t},\varepsilon_{t}}_{\mathbf{x}^{t}}(\mathbf{W}_{k+1/2})-H^{\delta_{t}}_{\mathbf{x}^{t}}(\mathbf{W}_{k+1/2})}|\mathcal{A}_{t,k}\right]\leq\mathbb{E}\left[\norm{H^{\delta_{t},b_{t}}_{\mathbf{x}^{t}}-H^{\delta_{t},b_{t},\varepsilon_{t}}_{\mathbf{x}^{t}}}|\mathcal{A}_{t,k}\right]$
$\displaystyle+\mathbb{E}\left[\norm{H^{\delta_{t},\varepsilon_{t}}_{\mathbf{x}^{t}}(\mathbf{W}_{k+1/2})-H^{\delta_{t}}_{\mathbf{x}^{t}}(\mathbf{W}_{k+1/2})}|\mathcal{A}_{t,k}\right]$
$\displaystyle\leq\sqrt{\mathbb{E}\left[\norm{H^{\delta_{t},b_{t}}_{\mathbf{x}^{t}}-H^{\delta_{t},b_{t},\varepsilon_{t}}_{\mathbf{x}^{t}}}^{2}|\mathcal{A}_{t,k}\right]}+\sqrt{\mathbb{E}\left[\norm{H^{\delta_{t},\varepsilon_{t}}_{\mathbf{x}^{t}}(\mathbf{W}_{k+1/2})-H^{\delta_{t}}_{\mathbf{x}^{t}}(\mathbf{W}_{k+1/2})}^{2}|\mathcal{A}_{t,k}\right]}$
$\displaystyle\leq\frac{2\varepsilon_{t}}{\delta_{t}}\left(\frac{1}{\sqrt{b_{t}}}+1\right)\sqrt{\sum_{i\in\mathcal{I}}L_{2,i}^{2}n_{i}^{2}}\equiv\alpha^{(2)}_{t}.$

Lastly, we bound

$\mathbb{E}\left[\norm{H^{\delta_{t},\varepsilon_{t}}_{\mathbf{z}_{k+1/2}}(\mathbf{W}_{k+1/2})-H^{\delta_{t}}_{\mathbf{z}_{k+1/2}}(\mathbf{W}_{k+1/2})}|\mathcal{A}_{t,k}\right]\leq\frac{2\varepsilon_{t}}{\delta_{t}}\sqrt{\sum_{i\in\mathcal{I}}L^{2}_{2,i}n^{2}_{i}}\equiv\alpha^{(3)}_{t}.$

We now set

$\displaystyle\mathcal{D}_{K,T}$
$\displaystyle\triangleq\frac{C^{2}}{2K\sum_{t=0}^{T-1}\gamma_{t}}+2C^{2}(L_{f}^{2}+\eta^{2}_{0})\frac{\sum_{t=0}^{T-1}\gamma^{2}_{t}}{\sum_{t=0}^{T-1}\gamma_{t}}$
$\displaystyle+\frac{3C^{2}\sum_{t=0}^{T-1}\gamma_{t}\eta_{t}}{2\sum_{t=0}^{T-1}\gamma_{t}}+\frac{2C^{2}}{K}+\frac{L_{h}\sum_{t=0}^{T-1}\delta_{t}\gamma_{t}}{\sum_{t=0}^{T-1}\gamma_{t}}.$

Using the definition of the gap function (4.1), we deduce

$\displaystyle\mathbb{E}[\Gamma(\bar{\mathbf{z}}^{T})]$
$\displaystyle\leq\mathcal{D}_{K,T}+\mathbb{E}\left[\frac{\sum_{t=0}^{T-1}\gamma^{3}_{t}}{K\sum_{t=0}^{T-1}\gamma_{t}}\sum_{k=0}^{K-1}\left(\norm{\mathbf{Y}^{1}_{t,k}}^{2}+\norm{\mathbf{Y}^{2}_{t,k}}^{2}\right)\right]$
$\displaystyle+\mathbb{E}\left[\frac{2\sum_{t=0}^{T-1}\gamma^{2}_{t}}{K\sum_{t=0}^{T-1}\gamma_{t}}\sum_{k=0}^{K-1}\left(\langle\mathbf{Y}^{1}_{t,k},\mathbf{u}^{1}_{t,k}-\mathbf{z}^{(t)}_{k+1/2}\rangle+\langle\mathbf{Y}^{2}_{t,k},\mathbf{u}_{t,k}^{2}-\mathbf{z}^{(t)}_{k+1/2}\rangle\right)\right]$
$\displaystyle+\mathbb{E}\left[\frac{\sum_{t=0}^{T-1}\gamma_{t}(C\norm{\varepsilon^{t}_{V}(\xi^{1:b_{t}}_{t})}+C\norm{\varepsilon^{t}_{h}(\mathbf{W}^{1:b_{t}}_{t})})}{\sum_{t=0}^{T-1}\gamma_{t}}\right]$
$\displaystyle+\mathbb{E}\left[\sum_{t=0}^{T-1}\frac{2\gamma^{2}_{t}}{K\sum_{s=0}^{T-1}\gamma_{s}}\sum_{k=0}^{K-1}\left(\norm{H^{\delta_{t},\varepsilon_{t}}_{\mathbf{z}_{k+1/2}}(\mathbf{W}_{t,k+1/2})}^{2}+\norm{H^{\delta_{t},\varepsilon_{t}}_{\mathbf{x}^{t}}(\mathbf{W}_{t,k+1/2})}^{2}\right)\right]$
$\displaystyle+\mathbb{E}\left[\sum_{t=0}^{T-1}\frac{C\gamma_{t}}{K\sum_{s=0}^{T-1}\gamma_{s}}\sum_{k=0}^{K-1}\norm{H^{\delta_{t},b_{t}}_{\mathbf{x}^{t}}-H^{\delta_{t},b_{t},\varepsilon_{t}}_{\mathbf{x}^{t}}+H^{\delta_{t},\varepsilon_{t}}_{\mathbf{x}^{t}}(\mathbf{W}_{t,k+1/2})-H^{\delta_{t}}_{\mathbf{x}^{t}}(\mathbf{W}_{t,k+1/2})}\right]$
$\displaystyle+\mathbb{E}\left[\sum_{t=0}^{T-1}\frac{C\gamma_{t}}{K\sum_{s=0}^{T-1}\gamma_{s}}\sum_{k=0}^{K-1}\norm{H^{\delta_{t},\varepsilon_{t}}_{\mathbf{z}^{(t)}_{k+1/2}}(\mathbf{W}_{t,k+1/2})-H^{\delta_{t}}_{\mathbf{z}^{(t)}_{k+1/2}}(\mathbf{W}_{t,k+1/2})}\right]$
$\displaystyle\leq\mathcal{D}_{K,T}+\frac{2\sigma^{2}\sum_{t=0}^{T-1}\gamma^{3}_{t}}{\sum_{t=0}^{T-1}\gamma_{t}}+\frac{4C\sigma\sum_{t=0}^{T-1}\gamma^{2}_{t}}{\sum_{t=0}^{T-1}\gamma_{t}}+\frac{\sum_{t=0}^{T-1}\frac{\gamma_{t}}{\sqrt{b_{t}}}C(M_{V}+(\sum_{i\in\mathcal{I}}L^{2}_{h_{i}}n_{i}^{2})^{1/2}}{\sum_{t=0}^{T-1}\gamma_{t}}$
$\displaystyle+\sum_{t=0}^{T-1}\frac{2\gamma^{2}_{t}\alpha^{(1)}_{t}}{\sum_{s=0}^{T-1}\gamma_{s}}+\sum_{t=0}^{T-1}\frac{C\alpha^{(2)}_{t}\gamma_{t}}{\sum_{s=0}^{T-1}\gamma_{s}}+\sum_{t=0}^{T-1}\frac{C\gamma_{t}\alpha^{(3)}_{t}}{\sum_{s=0}^{T-1}\gamma_{s}}$

Making the choice $\gamma_{t}=\eta_{t}=\delta_{t}=1/T$ as well as
$K=T,b_{t}\geq T^{2}$ and $\varepsilon_{t}=1/T^{2}$, we see that
$\alpha^{(1)}_{t}=\mathcal{O}(1/T^{2}),\alpha^{(2)}_{t}=\mathcal{O}(1/T)$ and
$\alpha^{(3)}_{t}=\mathcal{O}(1/T)$. It follows
$\mathbb{E}[\Gamma(\bar{\mathbf{z}}^{T})]=\mathcal{O}(C\sigma/T)$, which
completes the proof of Theorem 4.7.

##  7 Conclusion

In this work, we proposed a new solution approach to solve a fairly large
class of stochastic hierarchical games. Using a combination of smoothing,
zeroth-order gradient approximation, and iterative regularization, we develop
a novel variance reduction method for stochastic VIs affected by general
stochastic noise. We demonstrate consistency of the method by proving that
solution trajectory converges almost surely to a particular equilibrium of the
game and derive a $\mathcal{O}(1/T)$ convergence rate in terms of the expected
gap function, using a suitably defined averaged trajectory. This rate result
is robust to inexact solutions of the lower level problem of the follower and
aligns with state-of-the-art variance reduction methods tailored to finite-sum
problems. Furthermore, our approach is based on Tseng’s splitting technique,
which shares the same number of function calls as the popular extragradient
method, but saves on one projection step. This implies that our scheme reduces
the oracle complexity relative to vanilla mini-batch approaches and at the
same time reduces the computational bottlenecks in every single iteration.
This leaves open the door for many future investigations, involving bias and
non-convexities, that we leave for future research.

#### acknowledgements

The research of MST benefited from the support of the FMJH Program Gaspard
Monge for optimization and operations research and their interactions with
data science. CS acknowledges supported in part by NSFC under Grant 62373050.

##  Appendix A Auxiliary Facts

###  A.1 Generalities

Given a closed convex set $\mathcal{X}\subset\mathbb{R}^{n}$, we denote by
$\Pi_{\mathcal{X}}:\mathbb{R}^{n}\to\mathcal{X}$ the orthogonal projector
defined as

$\Pi_{\mathcal{X}}(\mathbf{w}):=\operatorname*{argmin}_{\mathbf{x}\in\mathcal{X}}\frac{1}{2}\norm{\mathbf{x}-\mathbf{w}}^{2}.$

This is the solution map of a strongly convex optimization problem with the
following well-known properties.

######  Lemma A.1.

Let $\mathcal{X}\subset\mathbb{R}^{n}$ be a nonempty closed convex set. Then:

  * (i)

$\Pi_{\mathcal{X}}(\mathbf{w})$ is the unique point satisfying
$\langle\mathbf{w}-\Pi_{\mathcal{X}}(\mathbf{w}),\mathbf{x}-\Pi_{\mathcal{X}}(\mathbf{w})\rangle\leq
0$ for all $\mathbf{x}\in\mathcal{X}$;

  * (ii)

For all $\mathbf{w}\in\mathbb{R}^{n}$ and $\mathbf{x}\in\mathcal{X}$, we have
$\norm{\Pi_{\mathcal{X}}(\mathbf{w})-\mathbf{x}}^{2}+\norm{\Pi_{\mathcal{X}}(\mathbf{w})-\mathbf{w}}^{2}\leq\norm{\mathbf{w}-\mathbf{x}}^{2}$;

  * (iii)

For all $\mathbf{w},\mathbf{v}\in\mathbb{R}^{n}$,
$\norm{\Pi_{\mathcal{X}}(\mathbf{w})-\Pi_{\mathcal{X}}(\mathbf{v})}^{2}\leq\norm{\mathbf{w}-\mathbf{v}}^{2}$;

Let $(\Omega,\mathcal{F},\mathbb{P})$ be a given probability space carrying a