# Stochastic approximation with decision-dependent distributions: asymptotic
normality and optimality

Joshua Cutler Department of Mathematics, U. Washington, Seattle, WA 98195;
<EMAIL_ADDRESS>Mateo Díaz Department of Computing and Mathematical Sciences,
Caltech, Pasadena, CA 91125, USA;
http://users.cms.caltech.edu/$\scriptstyle\sim$mateodd/. Dmitriy Drusvyatskiy
Department of Mathematics, U. Washington, Seattle, WA 98195;
www.math.washington.edu/$\scriptstyle\sim$ddrusv. Research of Drusvyatskiy was
supported by the NSF DMS 1651851 and CCF 1740551 awards.

###### Abstract

We analyze a stochastic approximation algorithm for decision-dependent
problems, wherein the data distribution used by the algorithm evolves along
the iterate sequence. The primary examples of such problems appear in
performative prediction and its multiplayer extensions. We show that under
mild assumptions, the deviation between the average iterate of the algorithm
and the solution is asymptotically normal, with a covariance that nicely
decouples the effects of the gradient noise and the distributional shift.
Moreover, building on the work of Hájek and Le Cam, we show that the
asymptotic performance of the algorithm is locally minimax optimal.

##  1 Introduction

The primary role of stochastic optimization in data science is to find a
learning rule (e.g., a classifier) from a limited data sample which enables
accurate prediction on unseen data. Classical theory crucially relies on the
assumption that both the observed data and the unseen data are generated from
the same distribution. Recent literature on strategic classification [12] and
performative prediction [23], however, has highlighted a variety of
contemporary settings where this assumption is grossly violated. One common
reason is that the data seen by a learning system may depend on or react to a
deployed learning rule. For example, members of the population may alter their
features in response to a deployed classifier in order to increase their
likelihood of being positively labeled—a phenomenon called gaming. Even when
the population is agnostic to the learning rule, the decisions made by the
learning system (e.g., loan approval) may inadvertently alter the profile of
the population (e.g., credit score). The goal of the learning system therefore
is to find a classifier that generalizes well under the response distribution.
The situation may be further compounded by a population that reacts to
multiple competing learners simultaneously [22, 33, 24].

In this work, we model decision-dependent problems using variational
inequalities. Namely, let $G(x,z)$ be a map that depends on the decision $x$
and data $z$, and let $\mathcal{X}$ be a convex set of feasible decisions. A
variety of classical learning problems can be posed as solving the variational
inequality

$0\in\mathop{\operatorname{\mathbb{E}}}_{z\sim\mathcal{P}}[G(x,z)]+N_{\mathcal{X}}(x),$
VI($\mathcal{P}$)

where $\mathcal{P}$ is some fixed distribution and $N_{\mathcal{X}}$ is the
normal cone to $\mathcal{X}$. Two examples are worth keeping in mind: $(i)$
standard problems of supervised learning amount to $G(x,z)=\nabla\ell(x,z)$
being the gradient of some loss function, while $(ii)$ stochastic games
correspond to $G(x,z)$ being a stacked gradient of the players’ individual
losses.

Following the recent literature on performative prediction [12, 23, 22], we
will be interested in settings where the distribution $\mathcal{P}$ is not
fixed but rather varies with $x$. With this in mind, let $\mathcal{D}(x)$ be a
family of distributions indexed by $x\in\mathcal{X}$. The interpretation is
that $\mathcal{D}(x)$ is the response of the population to a newly deployed
learning rule $x$. We posit that the goal of a learning system is to find a
point $x^{\star}$ so that $x=x^{\star}$ solves the variational inequality
VI($\mathcal{D}(x^{\star})$), or equivalently:

$0\in\mathop{\operatorname{\mathbb{E}}}_{z\sim\mathcal{D}(x^{\star})}G(x^{\star},z)+N_{\mathcal{X}}(x^{\star}).$

We will say that such points $x^{\star}$ are at equilibrium. In words, a
learning system that deploys an equilibrium point $x^{\star}$ has no incentive
to deviate from $x^{\star}$ based only on the solution of the variational
inequality VI($\mathcal{D}(x^{\star})$) induced by the response distribution
$\mathcal{D}(x^{\star})$. The setting of performative prediction [23]
corresponds to the choice $G(x,z)=\nabla\ell(x,z)$ for some loss function
$\ell$.111In the notation of [23], equilibria coincide with performatively
stable points. More generally, decision-dependent games, proposed by [22, 33,
24], correspond to the choice
$G(x,z)=\big{(}\nabla_{1}\ell_{1}(x,z),\ldots,\nabla_{k}\ell_{k}(x,z)\big{)}$
where $\nabla_{i}\ell_{i}(x,z)$ is the gradient of the $i$’th player’s loss
with respect to their decision $x_{i}$ and
$\mathcal{D}(x)=\mathcal{D}_{1}(x)\times\cdots\times\mathcal{D}_{k}(x)$ is a
product distribution. The specifics of these two examples will not affect our
results, and therefore we work with general maps $G(x,z)$.

Following the prevalent viewpoint in machine learning, we suppose that the
only access to the data distributions $\mathcal{D}(x)$ is by drawing samples
$z\sim\mathcal{D}(x)$. With this in mind, a natural algorithm for finding an
equilibrium point $x^{\star}$ is the stochastic forward-backward algorithm:

$\displaystyle\textrm{Sample}~{}z_{t}\sim\mathcal{D}(x_{t})$ (1)
$\displaystyle\textrm{Set}~{}x_{t+1}={\rm
proj}_{\mathcal{X}}\big{(}x_{t}-\eta_{t}G(x_{t},z_{t})\big{)},$

where ${\rm proj}_{\mathcal{X}}$ is the nearest-point projection onto
$\mathcal{X}$. Specializing to performative prediction [19] and its
multiplayer extension [22], this algorithm reduces to a basic projected
stochastic gradient iteration. The contribution of our paper can be informally
summarized as follows.

We show that algorithm (1) is asymptotically optimal for finding equilibrium
points.

In particular, our results imply asymptotic optimality of the basic stochastic
gradient methods for both single player and multiplayer performative
prediction.

###  1.1 Summary of main results

Arguing optimality of an algorithm is a two step process: $(i)$ estimate the
performance of the specific algorithm and $(ii)$ derive a matching lower-bound
that is valid among all possible estimation procedures. Beginning with the
former, we build on the seminal work of Polyak and Juditsky [25]. Letting
$\bar{x}_{t}=\frac{1}{t}\sum_{i=0}^{t-1}x_{i}$ denote the average iterate, we
show that the deviation $\sqrt{t}(\bar{x}_{t}-x^{\star})$ is asymptotically
normal with an appealingly simple covariance. See Figure 1 for an
illustration.222The code can be found at
https://github.com/mateodd25/Asymptotic-normality-in-performative-prediction.

######  Theorem 1.1 (Informal).

Suppose that $G(\cdot,z)$ is $\alpha$-strongly monotone and Lipschitz,
$G(x,\cdot)$ is $\beta$-Lipschitz, and the distribution $\mathcal{D}(\cdot)$
is $\gamma$-Lipschitz in the Wasserstein-1 distance. Suppose moreover that
$x^{\star}$ lies in the interior of $\mathcal{X}$ and we set
$\eta_{t}=\eta_{0}t^{-\nu}$ with $\eta_{0}>0$ and
$\nu\in\big{(}\frac{1}{2},1\big{)}$. Then in the regime
$\frac{\gamma\beta}{\alpha}<1$, the average iterate $\bar{x}_{t}$ converges to
$x^{\star}$ almost surely and satisfies

$\sqrt{t}(\bar{x}_{t}-x^{\star})\xlongrightarrow{d}\mathsf{N}\big{(}0,W^{-1}\cdot\Sigma\cdot
W^{-\top}\big{)},$

where

$\displaystyle\Sigma=\underset{z\sim\mathcal{D}(x^{\star})}{\operatorname{\mathbb{E}}}\big{[}G(x^{\star},z)G(x^{\star},z)^{\top}\big{]}\qquad\textrm{and}\qquad
W=\underbrace{\underset{z\sim\mathcal{D}(x^{\star})}{\operatorname{\mathbb{E}}}\big{[}\nabla
G(x^{\star},z)\big{]}}_{\rm
static}+\underbrace{\frac{d}{dy}\underset{z\sim\mathcal{D}(y)}{\operatorname{\mathbb{E}}}\big{[}G(x^{\star},z)\big{]}\Big{|}_{y=x^{\star}}}_{\rm
dynamic}.$

A few comments are in order. First, the regime $\frac{\gamma\beta}{\alpha}<1$
is, in essence, optimal because otherwise, equilibrium points may even fail to
exist. Second, the effect of the distributional shift on the asymptotic
covariance is entirely captured by the second “dynamic” term in $W$. Indeed,
when this term is absent, the product $W^{-1}\cdot\Sigma\cdot W^{-\top}$ is
precisely the asymptotic covariance of the stochastic forward-backward
algorithm applied to the static problem VI($\mathcal{D}(x^{\star})$) at
equilibrium.333Of course, this analogy is entirely conceptual, since
$\mathcal{D}(x^{\star})$ is unknown a priori. The proof of Theorem 1.1 follows
by interpreting (1) as a stochastic approximation algorithm for finding a zero
of the nonlinear map
$R(x)=\operatorname{\mathbb{E}}_{z\sim\mathcal{D}(x)}G(x,z)$, and then
applying the results of Polyak and Juditsky [25]. The main technical work is
in verifying the necessary assumptions needed in [25] for this particular
nonlinear map and the corresponding noise sequence.

(a) $\rho=0.25$

(b) $\rho=0.5$

(c) $\rho=0.9$

Figure 1: Consider the problem corresponding to $G(x,z)=\nabla\ell(x,z)$ with
$\ell(x,z)=\frac{1}{2}\|x-z\|^{2}$, and
$\mathcal{D}(x_{1},x_{2})=\mathsf{N}(\rho(x_{2},x_{1}),I_{2})$. A simple
computation shows $\Sigma=I_{2}$ and $W=[1,-\rho;-\rho,1].$ As $\rho$
approaches one, $W^{-1}$ becomes ill conditioned. We run Algorithm (1) $400$
times using $\eta_{t}=t^{-3/4}$ for $10^{6}$ iterations. The first row depicts
the resulting average iterates laid over the confidence regions (plotted in
logarithmic scale) corresponding to the asymptotic normal distribution. The
next two rows depict kernel density estimates from the asymptotic normal
distribution (top) and the deviation $\sqrt{t}(\bar{x}_{t}-x^{\star})$
(bottom).

A reasonable question to ask is whether there exists an algorithm with better
asymptotic guarantees than those of the stochastic forward-backward algorithm.
We will show that in a strong sense, the answer is no; that is, Algorithm (1)
is asymptotically optimal. We will lower-bound the performance of any
estimation procedure for finding an equilibrium point on an adversarially-
chosen sequence of small perturbations of the target problem. That is, we will
define a shrinking neighborhood of distribution maps $B_{1/k}$ around
$\mathcal{D}$. Every map $\mathcal{D}^{\prime}\in B_{1/k}$ defines a new
problem with equilibrium point $x^{\star}_{\mathcal{D}^{\prime}}$ near
$x^{\star}$. Roughly speaking, we will show that for every estimator
$\widehat{x}_{k}$ (i.e. a deterministic function of the observed samples), the
local minimax bound holds:

$\liminf_{k\rightarrow\infty}\sup_{\mathcal{D}^{\prime}\in
B_{1/k}}{\operatorname{\mathbb{E}}}\big{[}\mathcal{L}\big{(}\sqrt{k}(\widehat{x}_{k}-x^{\star}_{\mathcal{D}^{\prime}})\big{)}\big{]}\geq{\operatorname{\mathbb{E}}}\big{[}\mathcal{L}(Z)\big{]},$
(2)

where $\mathcal{L}$ is any nonnegative, symmetric, quasiconvex, lower
semicontinuous function. The end result is the following theorem.

######  Theorem 1.2 (Informal).

Suppose the same setting as in Theorem 1.1 and let $\mathcal{L}$ be any
nonnegative, symmetric, quasiconvex, lower semicontinuous function. Fix any
procedure for finding equilibrium points that outputs $\widehat{x}_{k}$ based
on $k$ samples. As $k\to\infty$, there is a family of perturbed distributions
$\mathcal{D}_{k}(x)$ converging to $\mathcal{D}(x)$, along with corresponding
equilibrium points $x_{k}^{\star}$ converging to $x^{\star}$, such that the
expected error
$\mathbb{E}[\mathcal{L}(\sqrt{k}(\widehat{x}_{k}-x_{k}^{\star}))]$ of the
estimator $\widehat{x}_{k}$ on the perturbed problem is asymptotically lower
bounded by $\mathbb{E}[\mathcal{L}(Z)]$ where
$Z\sim\mathsf{N}(0,W^{-1}\cdot\Sigma\cdot W^{-\top})$.

The formal statement of the theorem and its proof follow closely the classical
work of Hájek and Le Cam [17, 29] on statistical lower bounds and the more
recent work of Duchi and Ruan [8] on asymptotic optimality of the stochastic
gradient method. In particular, the fundamental role of tilt-stability and the
inverse function theorem highlighted in [8] is replaced by the implicit
function theorem paradigm.

###  1.2 Related literature

Our work builds on existing literature in machine learning and stochastic
optimization.

Learning with decision-dependent distributions. The basic setup for decision-
dependent problems that we use is inspired by the performative prediction
framework introduced by [23] and its multiplayer extension developed
independently in [22, 33, 24]. The stochastic gradient method for performative
prediction was first introduced and analyzed in [18], while the stochastic
forward-backward method for games was analyzed in [22]. A related work [7]
showed that a variety of popular gradient-based algorithms for performative
prediction can be understood as the analogous algorithms applied to a certain
static problem corrupted by a vanishing bias. In general, performatively
stable points (equilibria) are not “performatively optimal” in the sense of
[23]. Seeking to develop algorithms for finding performatively optimal points,
the paper [21] provides sufficient conditions for the prediction problem to be
convex; extensions of such conditions to games appear in [33, 22]. Algorithms
for finding performatively optimal points under a variety of different
assumptions and oracle models appear in [15, 21, 33, 22, 14]. The performative
prediction framework is largely motivated by the problem of strategic
classification [12], which has been studied extensively from the perspective
of causal inference [2, 20] and convex optimization [6]. Other lines of work
[3, 26, 4, 34] in performative prediction have focused on the setting in which
the environment evolves dynamically in time or experiences time drift.

Stochastic approximation. There is extensive literature on stochastic
approximation. The most relevant results for us are those of [25] that
quantify the limiting distribution of the average iterate of stochastic
approximation algorithms. Stochastic optimization problems with decision-
dependent uncertainties have appeared in the classical stochastic programming
literature, such as [31, 16, 28, 1, 10]. We refer the reader to the recent
paper [13], which discusses taxonomy and various models of decision-dependent
uncertainties. An important theme of these works is to utilize structural
assumptions on how the decision variables impact the distributions. In
contrast, much of the work on performative prediction [23, 22, 33, 24, 7, 18]
and our current paper are “model-free”.

Local minimax lower bounds in estimation. There is a rich literature on
minimax lower bounds in statistical estimation problems; we refer the reader
to [32, Chapter 15] for a detailed treatment. Typical results of this type
lower-bound the performance of any statistical procedure on a worst-case
instance of that procedure. Minimax lower bounds can be quite loose as they do
not consider the complexity of the particular problem that one is trying to
solve but rather that of an entire problem class to which it belongs. More
precise local minimax lower bounds, as developed by Hájek and Le Cam [17, 29],
provide much finer problem-specific guarantees. Building on this framework,
Duchi and Ruan [8] showed that the stochastic gradient method for standard
single-stage stochastic optimization problems is, in an appropriate sense,
locally minimax optimal. Our paper builds heavily on this line of work.

###  1.3 Outline

The outline of the paper is as follows. Section 2 records some basic notation
that we will use. Section 3 formally introduces/reviews the decision-dependent
framework. In Section 4, we show that the running average of the forward-
backward algorithm is asymptotically normal (Theorem 1.1), and identify its
asymptotic covariance. Finally, Section 5 proves the local minimax lower-bound
(Theorem 1.2). We defer many of the technical proofs to the Appendix.

##  2 Basic notation

Throughout, we let $\mathbf{R}^{d}$ denote a $d$-dimensional Euclidean space,
with inner product $\langle\cdot,\cdot\rangle$ and the induced Euclidean norm
$\|x\|=\sqrt{\langle x,x\rangle}$. For any set
$\mathcal{X}\subset\mathbf{R}^{d}$, the symbol ${\rm proj}_{\mathcal{X}}(x)$
will denote the set of nearest points of $\mathcal{X}$ to
$x\in\mathbf{R}^{d}$. We say that a function
$\mathcal{L}\colon{\mathbf{R}}^{d}\rightarrow{\mathbf{R}}$ is symmetric if it
satisfies $\mathcal{L}(x)=\mathcal{L}(-x)$ for all $x\in\mathbf{R}^{d}$, and
we say that $\mathcal{L}$ is quasiconvex if its sublevel set
$\\{x\mid\mathcal{L}(x)\leq c\\}$ is convex for any $c\in{\mathbf{R}}.$ For
any matrix $A\in{\mathbf{R}}^{m\times n}$, the symbols $\|A\|_{\rm op}$ and
$A^{\dagger}$ stand for the operator norm and Moore-Penrose pseudoinverse of
$A$, respectively. If $A$ is a symmetric matrix, then the notation $A\succeq
0$ will mean that $A$ is positive semidefinite (PSD). Given two symmetric
matrices $A$ and $B$, we write $A\succeq B$ as a shorthand for $A-B\succeq 0$.

A map $F\colon\mathcal{X}\to\mathbf{R}^{d}$ is called $\alpha$-strongly
monotone on $\mathcal{X}\subset\mathbf{R}^{d}$ if it satisfies

$\langle
F(x)-F(x^{\prime}),x-x^{\prime}\rangle\geq\alpha\|x-x^{\prime}\|^{2}\qquad\textrm{for
all }x,x^{\prime}\in\mathcal{X}.$

A map $F\colon\mathcal{X}\to\mathbf{R}^{k}$ is smooth on a set
$\mathcal{X}\subset\mathbf{R}^{d}$ if $F$ extends to a differentiable map on
an open neighborhood of $\mathcal{X}$; further, we say that $F$ is
$\beta$-smooth if the Jacobian of $F$ satisfies the Lipschitz condition

$\|\nabla F(x)-\nabla F(x^{\prime})\|_{\rm
op}\leq\beta\|x-x^{\prime}\|\qquad\textrm{for all
}x,x^{\prime}\in\mathcal{X}.$

Given a metric space $\mathcal{Z}$, we equip $\mathcal{Z}$ with its Borel
$\sigma$-algebra $\mathcal{B}(\mathcal{Z})$ and let $P_{1}(\mathcal{Z})$
denote the set of Radon probability measures on $\mathcal{Z}$ with a finite
first moment. We will measure the deviation between two measures $\mu,\nu\in
P_{1}(\mathcal{Z})$ with the Wasserstein-1 distance:

$W_{1}(\mu,\nu)=\sup_{h\in\textrm{Lip}_{1}(\mathcal{Z})}\left\\{\underset{X\sim\mu}{\operatorname{\mathbb{E}}}\big{[}h(X)\big{]}-\underset{Y\sim\nu}{\operatorname{\mathbb{E}}}\big{[}h(Y)\big{]}\right\\}.$

Here, $\textrm{Lip}_{1}(\mathcal{Z})$ denotes the set of $1$-Lipschitz
continuous functions $h\colon\mathcal{Z}\rightarrow\mathbf{R}$. We will refer
to a Borel measurable map between metric spaces simply as _measurable_. Given
a map $H\colon\mathcal{X}\times\mathcal{Z}\rightarrow\mathbf{R}^{k}$ with
$\mathcal{X}\subset{\mathbf{R}}^{d}$, the symbol $\nabla H(x,z)$ will always
denote the Jacobian of $H(\cdot,z)$ evaluated at $x$.

Let $P$ and $Q$ be probability distributions over $\mathcal{Z}$ such that $P$
is absolutely continuous with respect to $Q$ (denoted $P\ll Q$). Then for any
convex function $f\colon(0,\infty)\rightarrow{\mathbf{R}}$ with $f(1)=0$, the
$f$-divergence between $P$ and $Q$ is given by

$\Delta_{f}\big{(}P\parallel
Q\big{)}=\int_{\mathcal{Z}}f\bigg{(}\frac{dP}{dQ}\bigg{)}\,dQ,$ (3)

where $\frac{dP}{dQ}$ denotes the Radon-Nikodym derivative of $P$ with respect
to $Q$ and $f(0):=\lim_{t\downarrow 0}f(t)$. Abusing notation slightly, if $P$
is not absolutely continuous with respect to $Q$, we set
$\Delta_{f}\big{(}P\parallel Q\big{)}=\infty$.

##  3 Background on learning with decision-dependent distributions

In this section, we formally specify the class of problems that we consider
along with relevant assumptions. In order to model decision-dependence, we fix
a nonempty, closed, convex set $\mathcal{X}\subset\mathbf{R}^{d}$, a nonempty
metric space $(\mathcal{Z},d_{\mathcal{Z}})$, and a map
$\mathcal{D}\colon\mathcal{X}\to P_{1}(\mathcal{Z})$. Thus,
$\\{\mathcal{D}(x)\\}$ is a family of probability distributions indexed by
points $x\in\mathcal{X}$. The variational behavior of the map
$\mathcal{D}(\cdot)$ will play a central role in our work. In particular,
following [23], we will assume that $\mathcal{D}(\cdot)$ is Lipschitz
continuous.

######  Assumption 1 (Lipschitz distribution map).

There is a constant $\gamma>0$ such that

$W_{1}\big{(}\mathcal{D}(x),\mathcal{D}(x^{\prime})\big{)}\leq\gamma\|x-x^{\prime}\|\qquad\textrm{for
all }x,x^{\prime}\in\mathcal{X}.$

Next, we fix a measurable map
$G\colon\mathcal{X}\times\mathcal{Z}\to\mathbf{R}^{d}$ such that each section
$G(x,\cdot)$ is Lipschitz continuous, and we define the family of maps444Note
that $G_{x}(y)$ is well defined because $\mathcal{D}(x)$ has a finite first
moment and $G(y,\cdot)$ is Lipschitz continuous.

$G_{x}(y)=\underset{z\sim\mathcal{D}(x)}{\operatorname{\mathbb{E}}}G(y,z).$

We impose standard regularity conditions on $G(x,z)$.

######  Assumption 2 (Loss regularity).

There are constants $\bar{L},\beta,\alpha>0$ satisfying the following.

  1. (i)

(Lipschitz) There is a measurable function
$L\colon\mathcal{Z}\rightarrow[0,\infty)$ such that the Lipschitz bounds

$\displaystyle\|G(x,z)-G(x^{\prime},z)\|$ $\displaystyle\leq
L(z)\cdot\|x-x^{\prime}\|,$ $\displaystyle\|G(x,z)-G(x,z^{\prime})\|$
$\displaystyle\leq\beta\cdot d_{\mathcal{Z}}(z,z^{\prime}),$

hold for all $x,x^{\prime}\in\mathcal{X}$ and $z,z^{\prime}\in\mathcal{Z}$.
Further,
${\operatorname{\mathbb{E}}}_{z\sim\mathcal{D}(x)}\big{[}L(z)^{2}\big{]}\leq\bar{L}^{2}$
for all $x\in\mathcal{X}$.

  2. (ii)

(Monotonicity) For all $x\in\mathcal{X}$, the map $G_{x}(\cdot)$ is
$\alpha$-strongly monotone on $\mathcal{X}$.

  3. (iii)

(Compatibility) The inequality $\gamma\beta<\alpha$ holds.

A few comments are in order. Condition (i) simply asserts that the map
$G(x,z)$ is Lipschitz continuous both in $x$ and in $z$. An immediate
consequence is that $G_{x}(\cdot)$ is $\bar{L}$-Lipschitz continuous.
Condition (ii) is a standard convexity and monotonicity requirement. Condition
(iii) ensures that the Lipschitz constant $\gamma$ of $\mathcal{D}(\cdot)$ is
sufficiently small in comparison with the monotonicty constant $\alpha$,
signifying that the dynamics are “mild”. This condition is widely used in the
existing literature; see, e.g., [23, 22, 33, 24].

Assumptions 1 and 2 imply the following useful Lipschitz estimate on the
deviation $G_{x}(y)-G_{x^{\prime}}(y)$, arising from the shift in distribution
from $\mathcal{D}(x)$ to $\mathcal{D}(x^{\prime})$. We will use this estimate
often in what follows. The proof is identical to that of [22, Lemma 2]; a
short argument appears in Appendix A.

######  Lemma 3.1 (Deviation).

Suppose that Assumptions 1 and 2 hold. Then the estimate

$\displaystyle\|G_{x}(y)-G_{x^{\prime}}(y)\|$
$\displaystyle\leq\gamma\beta\cdot\|x-x^{\prime}\|\qquad\textrm{holds for all
}~{}x,x^{\prime},y\in\mathcal{X}.$

The following definition, originating in [23] for performative prediction and
in [22] for its multiplayer extension, is the key solution concept that we
will use.

######  Definition 3.2 (Equilibrium point).

We say that $x^{\star}$ is an equilibrium point if it satisfies:

$0\in G_{x^{\star}}(x^{\star})+N_{\mathcal{X}}(x^{\star}).$

In words, $x^{\star}$ is an equilibrium point if $x=x^{\star}$ solves the
variational inequality VI($\mathcal{P}$) induced by the distribution
$\mathcal{P}=\mathcal{D}(x^{\star})$. Equivalently, these are exactly the
fixed points of the map

$\mathtt{Sol}(x):=\\{y\mid 0\in G_{x}(y)+N_{\mathcal{X}}(y)\\},$ (4)

which is single-valued on $\mathcal{X}$ by strong monotonicity. Equilibrium
points have a clear intuitive meaning: a learning system that deploys a
learning rule $x^{\star}$ that is at equilibrium has no incentive to deviate
from $x^{\star}$ based only on the data drawn from $\mathcal{D}(x^{\star})$.
The key role of equilibrium points in (multiplayer) performative prediction is
by now well documented; see, e.g., [23, 22, 33, 24, 18, 7]. Most importantly,
equilibrium points exist and are unique under Assumptions 1 and 2. The proof
is identical to that of [22, Theorem 1]; we provide a short argument in
Appendix A for completeness.

######  Theorem 3.3 (Existence).

Suppose that Assumptions 1 and 2 hold. Then the map $\mathtt{Sol}(\cdot)$ is
$\frac{\gamma\beta}{\alpha}$-contractive on $\mathcal{X}$ and therefore the
problem admits a unique equilibrium point $x^{\star}$.

We note in passing that when $\gamma\beta\geq\alpha$, equilibrium points may
easily fail to exist; see, e.g., [23, Proposition 3.6]. Therefore, the regime
$\gamma\beta<\alpha$ is the natural setting to consider when searching for
equilibrium points.

##  4 Convergence and asymptotic normality

A central goal of performative prediction is the search for equilibrium
points, which are simply the fixed points of the map $\mathtt{Sol}(\cdot)$
defined in (4). Though the map $\mathtt{Sol}(\cdot)$ is contractive, it cannot
be evaluated directly since it involves evaluating the expectation
$G_{x}(y)={\operatorname{\mathbb{E}}}_{z\sim\mathcal{D}(x)}G(y,z)$. Employing
the standard assumption that the only access to $\mathcal{D}(x)$ is through
sampling, one may instead in iteration $t$ take a single stochastic forward-
backward step on the problem corresponding to $\mathtt{Sol}(x_{t})$. The
resulting procedure is recorded in Algorithm 1. In the setting of performative
prediction [18] and its multiplayer extension [22], the algorithm reduces to
stochastic projected gradient methods.

Algorithm 1 Stochastic Forward-Backward Method (SFB)

Input: initial $x_{0}\in\mathcal{X}$ and step sizes $(\eta_{t})_{t\geq
0}\subset(0,\infty)$

Step $t\geq 0$:

$\displaystyle\textrm{Sample~{}}z_{t}\sim\mathcal{D}(x_{t})$
$\displaystyle\textrm{Set~{}}x_{t+1}={\rm
proj}_{\mathcal{X}}\big{(}x_{t}-\eta_{t}G(x_{t},z_{t})\big{)}$

We will see that under very mild assumptions, the SFB iterates almost surely
converge to the equilibrium point $x^{\star}$. To this end, let
$(\Omega,\mathcal{F},\mathbb{P})$ denote the probability space where this
stochastic process takes place555Formally, in order to define a single
probability space, we may proceed as follows. Suppose that there is
probability space $(\mathcal{S},\mathcal{H},\mu)$ and a measurable map
$\theta\colon\mathcal{S}\times\mathcal{X}\rightarrow\mathcal{Z}$ such that
$\mathcal{D}(x)$ is the pushforward of $\mu$ by the section $\theta(\cdot,x)$.
This means that for every Borel set $V\subset\mathcal{Z}$, the
$\mathcal{D}(x)$-measure of $V$ is equal to the $\mu$-measure of the set
$\\{s\in\mathcal{S}:\theta(s,x)\in V\\}$. Then, we may simply define
$(\Omega,\mathcal{F},\mathbb{P})$ as the countable product
$(\mathcal{S},\mathcal{H},\mu)^{\infty}$., and define the noise vector

$\xi_{x}(z):=G(x,z)-G_{x}(x).$

We impose the following standard second moment bound on $\xi_{x}(z)$.

######  Assumption 3 (Variance bound).

There is a constant $K\geq 0$ such that

$\underset{z\sim\mathcal{D}(x)}{\operatorname{\mathbb{E}}}\|\xi_{x}(z)\|^{2}\leq
K(1+\|x-x^{\star}\|^{2})\qquad\textrm{for all }x\in\mathcal{X}.$

The following proposition shows that the SFB iterates almost surely converge
to the equilibrium point. The proof, which follows from a simple one-step
improvement bound for SFB and the Robbins-Siegmund theorem [27], appears in
Appendix A.3.

######  Proposition 4.1 (Almost sure convergence).

Suppose that Assumptions 1, 2, and 3 hold, and let the step size sequence
satisfy $\eta_{t}>0$, $\sum_{t=0}^{\infty}\eta_{t}=\infty$, and
$\sum_{t=0}^{\infty}\eta_{t}^{2}<\infty.$ Then $x_{t}$ converges to
$x^{\star}$ almost surely. Moreover, if $\eta_{t}=\eta_{0}t^{-\nu}$ with
$\nu\in\big{(}\frac{1}{2},1\big{)}$, then there exists a constant $C\geq 0$
such that

${\operatorname{\mathbb{E}}}\|x_{t}-x^{\star}\|^{2}\leq
C\eta_{t}\qquad\text{for all }t\geq 0.$

The main result of this section is the asymptotic normality of the average
iterates $\bar{x}_{t}:=\frac{1}{t}\sum_{i=0}^{t-1}x_{i}$, for which we require
the following additional assumption.

######  Assumption 4.

The following three conditions hold.

  1. 1.

(Interiority) The equilibrium point $x^{\star}$ lies in the interior of
$\mathcal{X}$.

  2. 2.

(Joint smoothness) The map $(x,y)\mapsto G_{x}(y)$ is smooth and has a
Lipschitz continuous Jacobian on some neighborhood of $(x^{\star},x^{\star})$.

  3. 3.

(Uniform integrability) There is a neighborhood $\mathcal{V}$ around
$x^{\star}$ in $\mathcal{X}$ such that for all $\delta>0$, there exists a
constant $N_{\delta}$ such that:

    1. (i)

$\underset{z\sim\mathcal{D}(x)}{\operatorname{\mathbb{E}}}\big{[}\|G(x^{\star},z)\|^{2}\mathbf{1}\\{\|G(x^{\star},z)\|\geq
N_{\delta}\\}\big{]}\leq\delta$ for all $x\in\mathcal{V}$;

    2. (ii)

$\underset{z\sim\mathcal{D}(x)}{\operatorname{\mathbb{E}}}\big{[}\|\xi_{x}(z)\|^{2}\mathbf{1}\\{\|\xi_{x}(z)\|\geq
N_{\delta}\\}\big{]}\leq\delta$ for all $x\in\mathcal{V}$.

A few comments are in order. First, the interiority condition 1 is a standard
assumption for asymptotic normality results even in static settings [25]. The
joint smoothness condition 2 is fairly mild. For example, it holds
automatically provided the partial derivatives $\nabla_{y}G_{x}(y)$ and
$\nabla_{x}G_{x}(y)$ exist and vary in a locally Lipschitz way with respect to
$(x,y)$; in turn, this holds if, on a neighborhood of $x^{\star}$, each
distribution $\mathcal{D}(x)$ admits a density
$p(x,z)=\frac{d\mathcal{D}(x)}{d\mu}(z)$ with respect to a base measure
$\mu\gg\mathcal{D}(x)$ such that $G(\cdot,z)$ and $p(\cdot,z)$ are
$C^{1,1}$-smooth666That is, differentiable with locally Lipschitz partial
derivatives. and sufficient integrability conditions hold to apply the
dominated convergence theorem. Condition 3 is a standard uniform integrability
condition. For instance, it holds when the random vectors $G(x^{\star},z)$ and
$\xi_{x}(z)$, with $z\sim\mathcal{D}(x)$, are sub-Gaussian with the same
variance proxy $\sigma^{2}$ for all $x$ near $x^{\star}.$

We are now ready to present our main result.

######  Theorem 4.2 (Asymptotic normality).

Suppose that Assumptions 1, 2, 3, and 4 hold, and let the step size sequence
be set to $\eta_{t}=\eta_{0}t^{-\nu}$ with $\eta_{0}>0$ and
$\nu\in\big{(}\frac{1}{2},1\big{)}$. Then the average iterates
$\bar{x}_{t}=\frac{1}{t}\sum_{i=0}^{t-1}x_{i}$ converge to $x^{\star}$ almost
surely, and

$\sqrt{t}(\bar{x}_{t}-x^{\star})\xlongrightarrow{d}\mathsf{N}\big{(}0,\nabla
R(x^{\star})^{-1}\cdot\Sigma\cdot\nabla R(x^{\star})^{-\top}\big{)},$

where

$R(x)=\underset{z\sim\mathcal{D}(x)}{\operatorname{\mathbb{E}}}\big{[}G(x,z)\big{]}\qquad\text{and}\qquad\Sigma=\underset{z\sim\mathcal{D}(x^{\star})}{\operatorname{\mathbb{E}}}\big{[}G(x^{\star},z)G(x^{\star},z)^{\top}\big{]}.$

Thus Theorem 4.2 asserts that under mild assumptions, the deviations
$\sqrt{t}(\bar{x}_{t}-x^{\star})$ converge in distribution to a Gaussian
random vector with a covariance matrix $\nabla
R(x^{\star})^{-1}\cdot\Sigma\cdot\nabla R(x^{\star})^{-\top}$. Explicitly,
under mild regularity conditions we may write

$\nabla
R(x^{\star})=\underbrace{\underset{z\sim\mathcal{D}(x^{\star})}{\operatorname{\mathbb{E}}}\big{[}\nabla
G(x^{\star},z)\big{]}}_{\rm
static}+\underbrace{\frac{d}{dy}\underset{z\sim\mathcal{D}(y)}{\operatorname{\mathbb{E}}}\big{[}G(x^{\star},z)\big{]}\Big{|}_{y=x^{\star}}}_{\rm
dynamic}.$

It is part of the theorem’s conclusion that the matrix $\nabla R(x^{\star})$
is invertible. It is worthwhile to note that the effect of the distributional
shift on the asymptotic variance is entirely captured by the second “dynamic”
term in $\nabla R(x^{\star})$. When the distributions $\mathcal{D}(x)$ admit a
density $p(x,z)=\frac{d\mathcal{D}(x)}{d\mu}(z)$ as before, the Jacobian
$\nabla R(x^{\star})$ admits the simple description:

$\nabla
R(x^{\star})=\underset{z\sim\mathcal{D}(x^{\star})}{\operatorname{\mathbb{E}}}\big{[}\nabla
G(x^{\star},z)\big{]}+\int G(x^{\star},z)\nabla
p(x^{\star},z)^{\top}\,d\mu(z).$

######  Example 4.3 (Performative prediction with location-scale families).

As an explicit example of Theorem 4.2, let us look at the case when
$G(x,z)=\nabla_{x}\ell(x,z)$ is a gradient of a loss function and
$\mathcal{D}(x)$ are affine perturbations of a fixed base distribution
$\mathcal{D}_{0}$. Such distributions are quite reasonable when modeling
performative effects, as explained in [21]. In this case, we have

$z\sim\mathcal{D}(x)\quad\Longleftrightarrow\quad z-Ax\sim\mathcal{D}_{0}$

for some fixed matrix $A\in\mathbf{R}^{n\times d}$. Then a quick computation
shows that we may write

$\nabla
R(x)=\underset{z\sim\mathcal{D}(x)}{\operatorname{\mathbb{E}}}\big{[}\nabla^{2}_{xx}\ell(x,z)+A^{\top}\nabla^{2}_{zx}\ell(x,z)\big{]}$

under mild integrability conditions. Thus the dynamic part of $\nabla
R(x^{\star})$ is governed by the product of the matrix $A^{\top}$ with the
matrix $\nabla^{2}_{zx}\ell(x^{\star},z)$ of mixed partial derivatives. The
former measures the performative effects of the distributional shift, while
the latter measures the sensitivity of the gradient
$\nabla_{x}\ell(x^{\star},z)$ to changes in the data $z$.

######  Example 4.4 (Multiplayer performative prediction with location-scale
families).

More generally, let us look at the problem of multiplayer performative
prediction [22]. In this case, the map $G$ takes the form

$G(x,z)=\big{(}\nabla_{1}\ell_{1}(x,z_{1}),\ldots,\nabla_{k}\ell_{k}(x,z_{k})\big{)}$

where $\ell_{i}$ is a loss for each player $i$ and $\nabla_{i}\ell_{i}$
denotes the gradient of $\ell_{i}$ with respect to the action $x_{i}$ of
player $i$. The distribution $\mathcal{D}(x)$ takes the product form

$\mathcal{D}(x)=\mathcal{D}_{1}(x)\times\cdots\times\mathcal{D}_{k}(x).$

As highlighted in [22], a natural assumption on the distributions is that
there exist some probability distributions $\mathcal{P}_{i}$ and matrices
$A_{i}$, $A_{-i}$ such that the following holds:

$z_{i}\sim\mathcal{D}_{i}(x)\quad\Longleftrightarrow\quad
z_{i}-A_{i}x_{i}-A_{-i}x_{-i}\sim\mathcal{P}_{i}.$

Here $x_{-i}$ denotes the vector obtained from $x$ by deleting the coordinates
$x_{i}$. Thus the distributional assumption asserts that the distribution used
by player $i$ is a “linear perturbation” of some base distribution
$\mathcal{P}_{i}$. We can interpret the matrices $A_{i}$ and $A_{-i}$ as
quantifying the performative effects of player $i$’s decisions and the rest of
the players’ decisions, respectively, on the distribution $\mathcal{D}_{i}$
governing player $i$’s data. It is straightforward to check the expression

$\nabla_{x}R_{i}(x)=\mathop{\mathbb{E}}_{z_{i}\sim\mathcal{D}_{i}(x)}\left[\nabla^{2}_{xx}\ell_{i}(x,z_{i})+\big{[}A_{i}^{\top}\nabla^{2}_{zx}\ell_{i}(x,z_{i}),A_{-i}^{\top}\nabla^{2}_{zx}\ell_{i}(x,z_{i})\big{]}\right]$

under mild integrability conditions. Thus, the dynamic part of $\nabla
R_{i}(x^{\star})$ is governed by $A^{\top}_{i}$, $A^{\top}_{-i}$ and the
matrix $\nabla^{2}_{zx}\ell_{i}(x^{\star},z_{i})$ of mixed partial
derivatives. The former measure the performative effects of the distributional
shift, while the latter measures the sensitivity of the gradient
$\nabla_{x}\ell_{i}(x^{\star},z_{i})$ to changes in the data $z_{i}$.

###  4.1 Proof of Theorem 4.2

The proof of Theorem 4.2 is based on the stochastic approximation result of
Polyak and Juditsky [25, Theorem 2], which we review in Appendix B. For the
remainder of the section, we impose the assumptions of Theorem 4.2.

Consider the map $R\colon\mathcal{X}\rightarrow\mathbf{R}^{d}$ given by

$R(x)=G_{x}(x).$

In light of the interiority condition
$x^{\star}\in\operatorname{int}\mathcal{X}$ of Assumption 4, the equilibrium
point $x^{\star}$ is the solution to the equation $R(x)=0$. To simplify
notation, let us define the direction of motion and noise vectors:

$v_{t}:=G(x_{t},z_{t})\qquad\text{and}\qquad\xi_{t}:=\xi_{x_{t}}(z_{t}).$ (5)

Observe that $v_{t}$ and $\xi_{t}$ satisfy the relation

$v_{t}=R(x_{t})+\xi_{t},$

and so we may write the iterates of Algorithm 1 as

$x_{t+1}=x_{t}-\eta_{t}\big{(}R(x_{t})+\xi_{t}+\zeta_{t}\big{)},$ (6)

where

$\zeta_{t}=\eta_{t}^{-1}\big{(}x_{t}-\eta_{t}v_{t}-{\rm
proj}_{\mathcal{X}}(x_{t}-\eta_{t}v_{t})\big{)}.$ (7)

Our goal is to apply Theorem B.1 to the process (6). In what follows, we
establish the necessary assumptions for the theorem, namely, Assumptions B.1,
B.2, and B.3.

We note first that Proposition 4.1 establishes Assumption B.2. It also follows
from Assumption 3 and Proposition 4.1 that $\xi_{t}$ is square-integrable;
further, letting $\mathbb{F}=(\mathcal{F}_{t})_{t\geq 0}$ denote the
filtration given by

$\mathcal{F}_{t}:=\sigma(z_{0},\dots,z_{t-1})$ (8)

(with $\mathcal{F}_{0}=\left\\{\emptyset,\Omega\right\\}$), it follows that
$\xi_{t}$ is $\mathcal{F}_{t+1}$-measurable and

${\operatorname{\mathbb{E}}}\big{[}\xi_{t}\mid\mathcal{F}_{t}\big{]}=\underset{z\sim\mathcal{D}(x_{t})}{\operatorname{\mathbb{E}}}\big{[}G(x_{t},z)\big{]}-G_{x_{t}}(x_{t})=0.$

Thus, $\xi_{t}$ constitutes a martingale difference process on
$(\Omega,\mathcal{F},{\mathbb{F}},\mathbb{P})$.

To verify condition (i) of Assumption B.1, we need the following lemma.

######  Lemma 4.5 (Lipschitz continuity).

The map $R$ is $(\bar{L}+\gamma\beta)$-Lipschitz continuous.

###### Proof.

Given $x,y\in\mathcal{X}$, we have

$\|R(x)-R(y)\|\leq\|G_{x}(x)-G_{y}(x)\|+\|G_{y}(x)-G_{y}(y)\|\leq(\bar{L}+\gamma\beta)\|x-y\|$

as a consequence of Lemma 3.1 and the $\bar{L}$-Lipschitz continuity of
$G_{y}(\cdot)$. ∎

In light of Assumption 3 and Lemma 4.5, we have

${\operatorname{\mathbb{E}}}\big{[}\|\xi_{t}\|^{2}\mid\mathcal{F}_{t}\big{]}\leq
K(1+\|x_{t}-x^{\star}\|^{2})\qquad\text{and}\qquad\|R(x_{t})\|^{2}\leq(\bar{L}+\gamma\beta)^{2}\|x_{t}-x^{\star}\|^{2};$

thus, upon setting

$\bar{K}:=2\big{(}K+(\bar{L}+\gamma\beta)^{2}\big{)}(1+\|x^{\star}\|^{2}),$

we have

${\operatorname{\mathbb{E}}}\big{[}\|\xi_{t}\|^{2}\mid\mathcal{F}_{t}\big{]}+\|R(x_{t})\|^{2}\leq\bar{K}(1+\|x_{t}\|^{2}).$

Hence condition (i) of Assumption B.1 holds. Moreover, since $x_{t}\rightarrow
x^{\star}$ almost surely, the second uniform integrability condition of
Assumption 4 implies

${\sup}_{t\geq
0}\,{\operatorname{\mathbb{E}}}\big{[}\|\xi_{t}\|^{2}\mathbf{1}\\{\|\xi_{t}\|>N\\}\mid\mathcal{F}_{t}\big{]}\xlongrightarrow{\textrm{a.s.}}0\qquad\text{as
}N\rightarrow\infty,$

so condition (ii) of Assumption B.1 holds.

The next lemma shows that
${\operatorname{\mathbb{E}}}[\xi_{t}\xi_{t}^{\top}\mid\mathcal{F}_{t}]$
converges to the positive semidefinite matrix

$\Sigma:=\underset{z\sim\mathcal{D}(x^{\star})}{\operatorname{\mathbb{E}}}\big{[}G(x^{\star},z)G(x^{\star},z)^{\top}\big{]}$

almost surely, thereby establishing condition (iii) of Assumption B.1.

######  Lemma 4.6 (Asymptotic covariance).

The following limits hold:

${\operatorname{\mathbb{E}}}\big{[}v_{t}v_{t}^{\top}\mid\mathcal{F}_{t}\big{]}\xlongrightarrow{\textrm{a.s.}}\Sigma\qquad\textrm{and}\qquad{\operatorname{\mathbb{E}}}\big{[}\xi_{t}\xi_{t}^{\top}\mid\mathcal{F}_{t}\big{]}\xlongrightarrow{\textrm{a.s.}}\Sigma.$

###### Proof.

Taking into account almost sure convergence of $x_{t}$ to $x^{\star}$, the
first uniform integrability condition of Assumption 4, and the Lipschitz
condition of Assumption 2, we may apply Lemma E.3 with $g=G$ along any sample
path $x_{t}\to x^{\star}$ to obtain
${\operatorname{\mathbb{E}}}\big{[}v_{t}v_{t}^{\top}\mid\mathcal{F}_{t}\big{]}\rightarrow\Sigma$
almost surely. Then we have

$\displaystyle{\operatorname{\mathbb{E}}}\big{[}\xi_{t}\xi_{t}^{\top}\mid\mathcal{F}_{t}\big{]}={\operatorname{\mathbb{E}}}\big{[}v_{t}v_{t}^{\top}\mid\mathcal{F}_{t}\big{]}-R(x_{t})R(x_{t})^{\top}\xlongrightarrow{\textrm{a.s.}}\Sigma,$

as claimed. ∎

Conditions (i), (ii), and (iii) of Assumption B.1 are now established. Now
consider the stochastic process $(\zeta_{t})_{t\geq 0}$ on
$(\Omega,\mathcal{F},{\mathbb{F}},\mathbb{P})$ given by (7). Since
$x_{t}\rightarrow x^{\star}$ almost surely and
$x^{\star}\in\operatorname{int}\mathcal{X}$, it follows that with probability
$1$, the vector $\zeta_{t}$ vanishes for all but finitely many $t$. Thus,
condition (iv) of Assumption B.1 holds, and the verification of Assumption B.1
is complete.

Next, we establish Assumption B.3 by leveraging the joint smoothness condition
of Assumption 4 to compute the Jacobian of $R$ and show that
${\min}_{j}\,\text{Re}\,\lambda_{j}\big{(}\nabla
R(x^{\star})\big{)}\geq\alpha-\gamma\beta>0$. In the following lemma, the
symbol $\nabla G_{x}(x)$ denotes the derivative of $G_{x}(\cdot)$ evaluated at
$x$.

######  Lemma 4.7 (Positivity of the Jacobian).

On a neighborhood of $x^{\star}$, the map $R$ is smooth with Lipschitz
Jacobian satisfying

$\nabla R(x)=\nabla
G_{x}(x)+\frac{d}{dy}\underset{z\sim\mathcal{D}(y)}{\operatorname{\mathbb{E}}}\big{[}G(x,z)\big{]}\Big{|}_{y=x}$
(9)

and hence

$\langle\nabla R(x)v,v\rangle\geq\alpha-\gamma\beta\qquad\text{for all
}v\in\mathbf{S}^{d-1}.$

In particular, the real part of every eigenvalue of $\nabla R(x)$ is no
smaller than $\alpha-\gamma\beta$.

###### Proof.

By the joint smoothness condition of Assumption 4, there exists a neighborhood
$U$ around $x^{\star}$ in $\mathcal{X}$ such that the map $(x,y)\mapsto
G_{x}(y)$ is smooth on $U\times U$ with Lipschitz Jacobian. Thus, the chain
rule reveals that $R$ is smooth on $U$ with Lipschitz Jacobian (9). Now let
$x\in U$ and $v\in\mathbf{S}^{d-1}$. Strong monotonicity of $G_{x}(\cdot)$
implies $\langle\nabla G_{x}(x)v,v\rangle\geq\alpha$. Moreover, Lemma 3.1
shows that the map $y\mapsto G_{y}(x)$ is $\gamma\beta$-Lipschitz on
$\mathcal{X}$, and therefore

$\bigg{|}\bigg{\langle}\bigg{(}\frac{d}{dy}G_{y}(x)\Big{|}_{y=x}\bigg{)}v,v\bigg{\rangle}\bigg{|}\leq\gamma\beta.$

Continuing with (9), we conclude

$\langle\nabla R(x)v,v\rangle\geq\alpha-\gamma\beta\qquad\text{for all
}v\in\mathbf{S}^{d-1}.$ (10)

Finally, let $w\in\mathbf{C}^{d}$ be a normalized eigenvector of $\nabla R(x)$
with associated eigenvalue $\lambda\in\mathbf{C}$. The inequality (10) implies
$\lambda_{\min}\big{(}\nabla R(x)+\nabla R(x)^{\top}\big{)}\geq
2(\alpha-\gamma\beta).$ Letting $w^{\text{H}}$ denote conjugate transpose of
$w$, we conclude

$2(\alpha-\gamma\beta)\leq w^{\text{H}}\big{(}\nabla R(x)+\nabla
R(x)^{\top}\big{)}w=w^{\text{H}}\nabla R(x)w+(w^{\text{H}}\nabla
R(x)w)^{\text{H}}=\lambda+\bar{\lambda}=2(\text{Re}\,\lambda),$

where the first inequality follows from the Rayleigh-Ritz theorem. ∎

Note that since $\nabla R$ is Lipschitz on a neighborhood of $x^{\star}$ and
$R(x^{\star})=0$, we have

$R(x)-\nabla R(x^{\star})(x-x^{\star})=O(\|x-x^{\star}\|^{2})\qquad\text{as
}x\rightarrow x^{\star}.$

Hence Assumption B.3 holds.

The only item left to prove in order to apply Theorem B.1 is that the step
sizes $\eta_{t}$ satisfy the required conditions. Recall that we have set
$\eta_{t}=\eta_{0}t^{-\nu}$ with $\nu\in\big{(}\tfrac{1}{2},1\big{)}$.
Therefore, we compute

$0\leq\frac{\eta_{t}-\eta_{t+1}}{\eta_{t}^{2}}=t^{\nu}\frac{(t+1)^{\nu}-t^{\nu}}{\eta_{0}(t+1)^{\nu}}\leq\frac{(t+1)^{\nu}-t^{\nu}}{\eta_{0}}.$

Using the fact that
${\lim_{t\rightarrow\infty}}\big{(}(t+1)^{\nu}-t^{\nu}\big{)}=0$ for any
$\nu\in(0,1)$, we obtain $\frac{\eta_{t}-\eta_{t+1}}{\eta^{2}_{t}}\rightarrow
0$. Moreover, because $\nu+1/2>1$, we have

$\sum_{t=1}^{\infty}\eta_{t}\cdot
t^{-1/2}=\eta_{0}\sum_{t=1}^{\infty}t^{-(\nu+1/2)}<\infty.$

Thus, the step sizes satisfy the required conditions, and an application of
Theorem B.1 to the process (6) completes the proof of Theorem 4.2.

##  5 Asymptotic optimality

In this section, we establish the local asymptotic optimality of Algorithm 1.
Our result builds on classical ideas from Hájek and Le Cam [17, 29] on lower
bounds for statistical estimation and the more recent work of Duchi and Ruan
[8] on asymptotic optimality of the stochastic gradient method. Throughout the
section, we fix a distribution map $\mathcal{D}\colon\mathcal{X}\rightarrow
P_{1}(\mathcal{Z})$ and a map
$G\colon\mathcal{X}\times\mathcal{Z}\to\mathbf{R}^{d}$ satisfying Assumptions
1 and 2. We will primarily focus on lower-bounding the performance of
estimation procedures when $\mathcal{D}$ is replaced by a small perturbation
$\mathcal{D}^{\prime}$. To this end, we make the following definition.

######  Definition 5.1 (Admissible distribution map).

A distribution map $\mathcal{D}^{\prime}\colon\mathcal{X}\rightarrow
P_{1}(\mathcal{Z})$ will be called admissible if Assumptions 1 and 2 hold for
$\mathcal{D}^{\prime}$ and $G$. The corresponding equilibrium point will then
be denoted by $x^{\star}_{\mathcal{D}^{\prime}}$.

We start with some intuition before delving into the details. Roughly
speaking, we aim to show that the covariance matrix of the asymptotic
distribution in Theorem 4.2 is smaller — in the positive semidefinite order —
than the asymptotic covariance achieved by any other algorithm for finding
equilibrium points. To capture this ordering, we will probe random vectors
with functions $\mathcal{L}\colon{\mathbf{R}}^{d}\rightarrow[0,\infty)$ that
are symmetric, quasiconvex, and lower semicontinuous, saying that
$z_{1}\sim\mathcal{P}_{1}$ is “better” than $z_{2}\sim\mathcal{P}_{2}$ if
${\operatorname{\mathbb{E}}}[\mathcal{L}(z_{1})]\leq{\operatorname{\mathbb{E}}}[\mathcal{L}(z_{2})]$
for all such $\mathcal{L}.$ This notion indeed implies the ordering of
covariances
$\mathbb{E}[z_{1}z_{1}^{\top}]\preceq\mathbb{E}[z_{2}z_{2}^{\top}]$.777If
$\Sigma_{1}$ and $\Sigma_{2}$ are the covariance matrices of $z_{1}$ and
$z_{2}$, respectively, and $\mathcal{L}_{u}(z)=(u^{\top}z)^{2}$, then
$\Sigma_{1}\preceq\Sigma_{2}$ if and only if
${\operatorname{\mathbb{E}}}[\mathcal{L}_{u}(z_{1})]\leq{\operatorname{\mathbb{E}}}[\mathcal{L}_{u}(z_{2})]$
for all $u$. We consider a local notion of complexity that evaluates an
estimator’s performance on problems close to the one we wish to solve. Since
our problem is defined via the distribution map $\mathcal{D}$, we will define
a set of close problems through a collection of admissible distribution maps
$B_{1/k}$ near $\mathcal{D}$ that shrink as $k$ increases. Every map
$\mathcal{D}^{\prime}\in B_{1/k}$ defines a new problem with equilibrium point
$x^{\star}_{\mathcal{D}^{\prime}}$. The rest of the section is devoted to
showing roughly that for every “estimation procedure” $\widehat{x}_{k}$ we
have

$\liminf_{k\rightarrow\infty}\sup_{\mathcal{D}^{\prime}\in
B_{1/k}}{\operatorname{\mathbb{E}}}\big{[}\mathcal{L}\big{(}\sqrt{k}(\widehat{x}_{k}-x^{\star}_{\mathcal{D}^{\prime}})\big{)}\big{]}\geq{\operatorname{\mathbb{E}}}\big{[}\mathcal{L}(Z)\big{]},$
(11)

where $Z\sim\mathsf{N}\big{(}0,W^{-1}\cdot\Sigma\cdot W^{-\top}\big{)}$ with
$\Sigma$ and $W$ as in Theorem 1.1.

To proceed formally, we need to impose the following assumptions on the base
distribution map $\mathcal{D}$ and the map $G$. Despite being technical, these
conditions are very mild and essentially amount to quantifying the smoothness
of $\mathcal{D}$ and $G$. In particular, we will make use of a certain set of
test functions to be integrated against $\mathcal{D}(x)$.

######  Definition 5.2 (Test functions).

Given a compact metric space $\mathcal{K}$, we let
$\mathcal{T}(\mathcal{K},\mathcal{Z})$ consist of all bounded measurable
functions $\phi\colon\mathcal{K}\times\mathcal{Z}\rightarrow\mathbf{R}$ which
admit a constant $L_{\phi}$ such that each section $\phi(\cdot,z)$ is
$L_{\phi}$-Lipschitz on $\mathcal{K}$. For any
$\phi\in\mathcal{T}(\mathcal{K},\mathcal{Z})$, we set $M_{\phi}:=\sup|\phi|$.

######  Assumption 5.

The following three conditions hold.

  1. 1.

(Compactness) The set $\mathcal{X}$ is compact, and the set $\mathcal{Z}$ is
bounded.

  2. 2.

(Smoothness of the distribution) There exists a constant $\vartheta\geq 0$
such that for every function $\phi\in\mathcal{T}(\mathcal{K},\mathcal{Z})$,
the function

$x\mapsto\mathop{\operatorname{\mathbb{E}}}_{z\sim\mathcal{D}(x)}{\phi(y,z)}$

is $C^{1}$-smooth on $\mathcal{X}$ for each $y\in\mathcal{K}$, and the map

$(x,y)\mapsto\nabla_{x}\Big{(}\mathop{\operatorname{\mathbb{E}}}_{z\sim\mathcal{D}(x)}{\phi(y,z)}\Big{)}$

is $\vartheta({L}_{\phi}+M_{\phi})$-Lipschitz on
$\mathcal{X}\times\mathcal{K}$.888The same conclusion then holds for all
measurable maps
$\phi\colon\mathcal{K}\times\mathcal{Z}\rightarrow\mathbf{R}^{n}$ with
$n\in{\mathbf{N}}$,
$L_{\phi}:=\sup_{z}\text{Lip}\left(\phi(\cdot,z)\right)<\infty$, and
$M_{\phi}:=\sup\|\phi\|<\infty$.

  3. 3.

(Lipschitz Jacobian) There exist a measurable function
$\Lambda\colon\mathcal{Z}\rightarrow[0,\infty)$ and constants
$\bar{\Lambda},\beta^{\prime}>0$ such that for every $z\in\mathcal{Z}$ and
$x\in\mathcal{X}$, the section $G(\cdot,z)$ is $\Lambda(z)$-smooth on
$\mathcal{X}$ with
${\operatorname{\mathbb{E}}}_{z\sim\mathcal{D}(x)}\big{[}\Lambda(z)\big{]}\leq\bar{\Lambda}$,
and the section $\nabla_{x}G(x,\cdot)$ is $\beta^{\prime}$-Lipschitz on
$\mathcal{Z}$.

The first condition is imposed mainly for simplicity. The last two smoothness
conditions are required in our arguments to apply a chain of dominated
convergence and implicit function theorems. To illustrate with a concrete
example, suppose that there exists a Borel probability measure $\mu$ on
$\mathcal{Z}$ such that $\mathcal{D}(x)\ll\mu$ for all $x\in\mathcal{X}$, and
consider the density $p(x,z)=\frac{d\mathcal{D}(x)}{d\mu}(z)$. If there exist
constants $\Lambda_{p},L_{p}\geq 0$ such that each section $p(\cdot,z)$ is
$\Lambda_{p}$-smooth and $\sup_{x,z}\|\nabla_{x}p(x,z)\|\leq L_{p}$, then the
second condition of Assumption 5 holds with
$\vartheta=\max\\{\Lambda_{p},L_{p}\\}$.

Next, we describe the “estimation procedures” to which our desired lower bound
will apply. We assume that the estimator $\widehat{x}_{k}$ in the lower bound
depends on samples $z_{0},\dots,z_{k-1}\in\mathcal{Z}$. Unlike in classical
settings, these samples are not necessarily drawn from a fixed distribution.
Instead, they are drawn from distributions
$\mathcal{D}(\tilde{x}_{0}),\ldots,\mathcal{D}(\tilde{x}_{k-1})$ taken along a
possibly random sequence of points $\tilde{x}_{0},\ldots,\tilde{x}_{k-1}$. In
order to model the dynamic nature of the estimation procedure, we introduce
the following definition.

######  Definition 5.3 (Dynamic estimation procedure).

We say that a sequence of measurable maps
$\mathcal{A}_{k}\colon\mathcal{Z}^{k}\times\mathcal{X}^{k}\to\mathcal{X}$ is a
dynamic estimation procedure if for any admissible distribution map
$\mathcal{D}^{\prime}$ and initial point $\tilde{x}_{0}\in\mathcal{X}$, the
recursively defined sequences

$z_{k}\sim\mathcal{D}^{\prime}(\tilde{x}_{k})\qquad\text{and}\qquad\tilde{x}_{k+1}=\mathcal{A}_{k+1}(z_{0},\dots,z_{k},\tilde{x}_{0},\dots,\tilde{x}_{k})$
(12)

are such that $\tilde{x}_{k}$ converges to $x^{\star}_{\mathcal{D}^{\prime}}$
almost surely.

Thus, the dynamic estimation procedure $\mathcal{A}_{k}$ plays the role of the
decision-maker that selects the sequence of points at which to query the
distribution map $\mathcal{D}^{\prime}$. Importantly, $\mathcal{A}_{k}$ is
assumed to be a deterministic function of its arguments. In particular,
Proposition 4.1 shows that maps $\mathcal{A}_{k}$ corresponding to Algorithm
1, i.e.,

$\mathcal{A}_{k+1}(z_{0},\dots,z_{k},x_{0},\dots,x_{k})={\rm
proj}_{\mathcal{X}}\big{(}{x}_{k}-\eta_{k}G({x}_{k},z_{k})\big{)},$

furnish a dynamic estimation procedure under standard conditions on
$\eta_{k}$. Although this map $\mathcal{A}_{k+1}$ depends only on the last
iterate $x_{k}$ and the last sample $z_{k}$, general dynamic estimation
procedures may depend arbitrarily on the previous samples and iterates.

Next, we define the neighborhoods $B_{1/k}$ of $\mathcal{D}$ in (11). These
neighborhoods will be defined based on an $f$-divergence induced by any
$C^{3}$-smooth convex function $f\colon(0,\infty)\rightarrow{\mathbf{R}}$ with
$f(1)=0$. For any distribution map $\mathcal{D}^{\prime}$, we define the
similarity measure

$\Delta_{f}\big{(}\mathcal{D}^{\prime}\parallel\mathcal{D}\big{)}:=\sup_{x\in\mathcal{X}}\Delta_{f}\big{(}\mathcal{D}^{\prime}(x)\parallel\mathcal{D}(x)\big{)},$

where $\Delta_{f}\big{(}\mathcal{D}^{\prime}(x)\parallel\mathcal{D}(x)\big{)}$
denotes the usual $f$-divergence between the probability measures
$\mathcal{D}^{\prime}(x)$ and $\mathcal{D}(x)$, defined in (3). We then define
the admissible neighborhood

$B_{\varepsilon}:=\\{\mathcal{D}^{\prime}\colon\mathcal{X}\rightarrow
P_{1}(\mathcal{Z})\mid\mathcal{D}^{\prime}\text{ is admissible and
}\Delta_{f}(\mathcal{D}^{\prime}\parallel\mathcal{D})\leq\varepsilon\\}.$

Equipped with these tools, we are ready to state the main result of this
section.

######  Theorem 5.4 (Asymptotic optimality).

Suppose that Assumptions 1, 2, and 5 hold and let $\\{\mathcal{A}_{k}\\}$ be a
dynamic estimation procedure. Fix an initial point
$\tilde{x}_{0}\in\mathcal{X}$, and for any admissible distribution map
$\mathcal{D}^{\prime}$, let
${\operatorname{\mathbb{E}}}_{P_{k}^{{}^{\prime}}}[\cdot]$ denote the
expectation with respect to $z_{0},\ldots,z_{k-1}$ produced by the recursion
(12). Let $\mathcal{L}\colon{\mathbf{R}}^{d}\rightarrow[0,\infty)$ be
symmetric, quasiconvex, and lower semicontinuous, and let
$\widehat{x}_{k}\colon\mathcal{Z}^{k}\rightarrow{\mathbf{R}}^{d}$ be a
sequence of estimators. Then for any $c>0$, the inequality

$\liminf_{k\rightarrow\infty}\sup_{\mathcal{D}^{\prime}\in
B_{c/{k}}}{\operatorname{\mathbb{E}}}_{P_{k}^{\prime}}\big{[}\mathcal{L}\big{(}\sqrt{k}({\widehat{x}}_{k}-x^{\star}_{\mathcal{D}^{\prime}})\big{)}\big{]}\geq{\operatorname{\mathbb{E}}}\big{[}\mathcal{L}(Z)\big{]}$

holds, where $Z\sim\mathsf{N}\big{(}0,W^{-1}\cdot\Sigma\cdot W^{-\top}\big{)}$
with

$\displaystyle\Sigma$
$\displaystyle=\underset{z\sim\mathcal{D}(x^{\star})}{\operatorname{\mathbb{E}}}\big{[}G(x^{\star},z)G(x^{\star},z)^{\top}\big{]}\quad\text{and}\quad
W={\underset{z\sim\mathcal{D}(x^{\star})}{\operatorname{\mathbb{E}}}\big{[}\nabla_{x}G(x^{\star},z)\big{]}}+{\frac{d}{dy}\underset{z\sim\mathcal{D}(y)}{\operatorname{\mathbb{E}}}\big{[}G(x^{\star},z)\big{]}\Big{|}_{y=x^{\star}}}.$

Most importantly, observe that the distribution of $Z$ in Theorem 5.4
coincides with the asymptotic distribution of
$\sqrt{t}(\bar{x}_{t}-x^{\star})$ in Theorem 4.2, thereby justifying the
asymptotic optimality of the stochastic forward backward algorithm.

###  5.1 Proof of Theorem 5.4

Our strategy for proving Theorem 5.4 follows a similar general pattern as in
[8] and crucially relies on the Hájek-Le Cam minimax theorem. Throughout the
section, we impose the assumptions and notation of Theorem 5.4. In order to
simplify notation, we will use the shorthand $\mathcal{D}_{x}=\mathcal{D}(x)$.

####  5.1.1 Tilted distributions

The proof proceeds by constructing “tilt perturbations” of $\mathcal{D}$, in
the language of [8], which are contained in $B_{1/k}$ and realize the claimed
lower bound of Theorem 5.4. These distributions encode the difficult instances
near the target problem. The tilted distributions will be parametrized by a
function $h\colon\mathbf{R}\to[-1,1]$, a map
$g\colon\mathcal{X}\times\mathcal{Z}\to\mathbf{R}^{d}$, and a vector
$u\in\mathbf{R}^{d}$. We start by defining the relevant class of maps
$g\in\mathcal{G}$.

######  Definition 5.5.

Let $\mathcal{G}$ consist of all measurable maps
$g\colon\mathcal{X}\times\mathcal{Z}\rightarrow{\mathbf{R}}^{d}$ satisfying
the following.

  1. 1.

(Unbiased)
${\operatorname{\mathbb{E}}}_{z\sim\mathcal{D}(x)}\big{[}g(x,z)\big{]}=0$ for
all $x\in\mathcal{X}$.

  2. 2.

(Lipschitz) There exists a constant $\beta_{g}>0$ such that for every
$x\in\mathcal{X}$, the section $g(x,\cdot)$ is $\beta_{g}$-Lipschitz on
$\mathcal{Z}$.

  3. 3.

(Smoothness) There exist a measurable function
$\Lambda_{g}\colon\mathcal{Z}\rightarrow[0,\infty)$ and constants
${\bar{\Lambda}_{g},\beta^{\prime}_{g}>0}$ such that for every
$z\in\mathcal{Z}$ and $x\in\mathcal{X}$, the section $g(\cdot,z)$ is
$\Lambda_{g}(z)$-smooth on $\mathcal{X}$ with
${\operatorname{\mathbb{E}}}_{z\sim\mathcal{D}(x)}\big{[}\Lambda_{g}(z)\big{]}\leq\bar{\Lambda}_{g}$,
and the section $\nabla_{x}g(x,\cdot)$ is $\beta^{\prime}_{g}$-Lipschitz on
$\mathcal{Z}$.

For our purposes, the most important map in $\mathcal{G}$ is the noise

$\xi(x,z):=G(x,z)-G_{x}(x),$

which belongs to $\mathcal{G}$ as a consequence of Assumptions 2 and 5 and
Lemma C.4. For any $g\in\mathcal{G}$, we will employ the notation
$g_{x}(z):=g(x,z)$.

Next, we fix $g\in\mathcal{G}$ and an arbitrary $C^{3}$-smooth function
$h\colon{\mathbf{R}}\rightarrow[-1,1]$ such that its first three derivatives
are bounded and $h(t)=t$ for all $t\in[-1/2,1/2]$. Given
$u\in{\mathbf{R}}^{d}$, we define a new map
$\mathcal{D}^{u}\colon\mathcal{X}\rightarrow P_{1}(\mathcal{Z})$ that maps $x$
to the probability distribution $\mathcal{D}_{x}^{u}$ given by

$d\mathcal{D}^{u}_{x}(z):=\frac{1+h\big{(}u^{\top}g_{x}(z)\big{)}}{C^{u}_{x}}\,d\mathcal{D}_{x}(z),$
(13)

where $C_{x}^{u}$ is the normalizing constant $C_{x}^{u}:=1+\int
h\big{(}u^{\top}g_{x}(z)\big{)}\,d\mathcal{D}_{x}(z)$. Reassuringly, for all
sufficiently small $u$, the perturbed distribution map $\mathcal{D}^{u}$ is
admissible. This is the content of the following lemma, whose proof appears in
Appendix C.1.

######  Lemma 5.6 (Perturbed distributions are admissible).

There exists a convex neighborhood $\mathcal{U}$ of zero such that for all
$u\in\mathcal{U}$, the map $\mathcal{D}^{u}$ is admissible.999We actually show
that constants $\gamma^{u}$, $\bar{L}^{u}$, and $\alpha^{u}$ for
$\mathcal{D}^{u}$ from Assumptions 1 and 2 deviate from $\gamma$, $\bar{L}$,
and $\alpha$ by at most $O(\|u\|)$ as $u\rightarrow 0$.

With this lemma in mind, we set $x^{\star}_{u}:=x^{\star}_{\mathcal{D}^{u}}$
for all $u\in\mathcal{U}$. For each $u\in\mathbf{R}^{d}$, we define the
sequence of perturbed distributions $P_{k,u}$ on $\mathcal{Z}^{k}$ through the
recursion

$z_{i}\sim\mathcal{D}_{\tilde{x}_{i}}^{u}\qquad\text{and}\qquad\tilde{x}_{i+1}=\mathcal{A}(z_{0},\dots,z_{i},\tilde{x}_{0},\dots,\tilde{x}_{i})\qquad\text{for
all }i\in\\{0,\ldots,k-1\\}.$

Upon taking $g=\xi$, this collection of tilted distributions achieves the
desired lower bound.

######  Theorem 5.7.

Let $\mathcal{L}\colon{\mathbf{R}}^{d}\rightarrow[0,\infty)$ be symmetric,
quasiconvex, and lower semicontinuous, let
$\widehat{x}_{k}\colon\mathcal{Z}^{k}\rightarrow{\mathbf{R}}^{d}$ be a
sequence of estimators, and suppose $g=\xi$. Then for any $c>0$, the
inequality

$\liminf_{k\rightarrow\infty}\sup_{\|u\|\leq\sqrt{c/k}}{\operatorname{\mathbb{E}}}_{P_{k,u}}\big{[}\mathcal{L}\big{(}\sqrt{k}(\widehat{x}_{k}-x^{\star}_{u})\big{)}\big{]}\geq{\operatorname{\mathbb{E}}}\big{[}\mathcal{L}(Z)\big{]}$
(14)

holds, where $Z\sim\mathsf{N}\big{(}0,W^{-1}\cdot\Sigma\cdot
W^{-\top}\big{)}$.

Theorem 5.4 follows directly from the above result. To see this, it suffices
to show that for any given $c>0$, there exists $\bar{c}>0$ such that the
condition $\|u\|\leq\sqrt{\bar{c}/k}$ implies $\mathcal{D}^{u}\in B_{c/k}$ for
all sufficiently large $k$. To this end, we have to check two conditions:
first, that $\mathcal{D}^{u}$ is admissible, and second, that
$\Delta_{f}(\mathcal{D}^{u}\parallel\mathcal{D})\leq c/k.$ It follows from
Lemma 5.6 that for any $\bar{c}>0$, the condition $\|u\|\leq\sqrt{\bar{c}/k}$
implies that $\mathcal{D}^{u}$ is admissible for all sufficiently large $k$.
Let us now consider the second condition.

Without loss of generality, we may assume that $f^{\prime}(1)=0,$ since we can
take $\tilde{f}(t)=f(t)-f^{\prime}(1)(t-1)$ and
$\Delta_{f}=\Delta_{\tilde{f}}.$ Then for all $x\in\mathcal{X}$ and
sufficiently small $u\in\mathbf{R}^{d}$, we have

$\displaystyle\Delta_{f}\big{(}\mathcal{D}_{x}^{u}\parallel\mathcal{D}_{x}\big{)}$
$\displaystyle=\int
f\bigg{(}\frac{1+h\big{(}u^{\top}g_{x}(z)\big{)}}{C_{x}^{u}}\bigg{)}\,d\mathcal{D}_{x}(z)$
$\displaystyle=\int
f\bigg{(}\frac{1+u^{\top}g_{x}(z)}{C_{x}^{u}}\bigg{)}\,d\mathcal{D}_{x}(z)$
(15)
$\displaystyle=\frac{f^{\prime\prime}(1)}{2}u^{\top}\bigg{(}\mathop{\operatorname{\mathbb{E}}}_{z\sim\mathcal{D}_{x}}g_{x}(z)g_{x}(z)^{\top}\bigg{)}u+r_{x}(u),$
(16)

where $r_{x}\colon{\mathbf{R}}^{d}\rightarrow{\mathbf{R}}$ is such that
$\sup_{x\in\mathcal{X}}|r_{x}(u)|=o(\|u\|^{2})$ as $u\to 0$. The equality (15)
follows since $g_{x}(z)$ is uniformly bounded over
$\mathcal{X}\times\mathcal{Z}$ (see Lemma C.1) and hence for small $u$ we have
$h\big{(}u^{\top}g_{x}(z)\big{)}=u^{\top}g_{x}(z)$ for all $x\in\mathcal{X}$
and $z\in\mathcal{Z}$. The last equality (16) follows by the dominated
convergence theorem; we defer the details to Lemma C.3. Since $g_{x}(z)$ is
uniformly bounded, there is a constant $a>0$ for which
$\sup_{x\in\mathcal{X}}\|{\operatorname{\mathbb{E}}}_{z\sim\mathcal{D}_{x}}[g_{x}(z)g_{x}(z)^{\top}]\|_{\rm
op}\leq a$. Further, given any $b>0$, there is a neighborhood $U$ of zero such
that $\sup_{x\in\mathcal{X},\,u\in U}\|u\|^{-2}|r_{x}(u)|\leq b.$ Therefore,
by setting $\bar{c}=c/(\frac{a}{2}f^{\prime\prime}(1)+b)$, it follows that the
condition $\|u\|^{2}\leq\bar{c}/k$ implies
$\Delta_{f}(\mathcal{D}^{u}\parallel\mathcal{D})\leq c/k$ for all sufficiently
large $k$. Thus, Theorem 5.4 is a corollary of Theorem 5.7. The remainder of
this section focuses on establishing Theorem 5.7.

####  5.1.2 Proof of Theorem 5.7

The proof of this result is based on the classical Hájek-Le Cam minimax
theorem. To state this result, we need two standard definitions from
statistics. In what follows, we say that a sequence of random vectors $X_{k}$
in $\mathbf{R}^{n}$ on probability spaces $(\Omega_{k},\mathcal{F}_{k},P_{k})$
is $o_{P_{k}}(1)$ if $X_{k}$ tends to zero in $P_{k}$ probability, i.e.,
$\lim_{k\rightarrow\infty}P_{k}\\{\|X_{k}\|\geq\varepsilon\\}=0$ for all
$\varepsilon>0$. We also employ the notation

$X_{k}\mathrel{\stackunder[.5pt]{\xrightarrow{\makebox[18.26393pt]{$\scriptstyle
d$}}}{\scriptstyle P_{k}\,}}\mathcal{L}_{X}$

to indicate that $X_{k}\xlongrightarrow{d}X$ where $X\sim\mathcal{L}_{X}$ and
the laws of $X_{k}$ are taken with respect to $P_{k}$.

######  Definition 5.8 (Locally asymptotically normal).

Let $U\subset{\mathbf{R}}^{d}$ be a neighborhood of zero. For each
$k\in\mathbf{N}$ and $u\in U$, let $Q_{k,u}$ be a probability measure on a
measurable space $(\Omega_{k},\mathcal{F}_{k})$. The family $\\{Q_{{k,u}}\\}$
is _locally asymptotically normal (at zero)_ with precision $V\succeq 0$ if
there exist
$Z_{k}\mathrel{\stackunder[.5pt]{\xrightarrow{\makebox[22.0556pt]{$\scriptstyle
d$}}}{\scriptstyle Q_{k,0}\,}}\mathsf{N}(0,V)$ such that for every bounded
sequence $\\{u_{k}\\}$ in $U$, we have

$\log\frac{dQ_{k,u_{k}}}{dQ_{k,0}}=u_{k}^{\top}Z_{k}-\frac{1}{2}u_{k}^{\top}Vu_{k}+o_{Q_{k,0}}(1).$
(17)

######  Definition 5.9 (Regular mapping sequence).

Let $U\subset{\mathbf{R}}^{n}$ be a neighborhood of zero. A sequence of
mappings $\\{\Gamma_{k}\colon U\rightarrow{\mathbf{R}}^{d}\\}$ is _regular_
with derivative $D\in{\mathbf{R}}^{d\times n}$ if

$\lim_{k\rightarrow\infty}\sqrt{k}\big{(}\Gamma_{k}(u)-\Gamma_{k}(0)\big{)}=Du\qquad\text{for
all }u\in U.$

Equipped with these definitions, we are ready to state the minimax theorem;
which appears for example in [17, Theorem 6.6.2], [30, Theorem 3.11.5], and
[8, Lemma 8.2].

######  Theorem 5.10.

Let $\\{\Omega_{k},\mathcal{F}_{k},Q_{{k,u}}\\}_{u\in U}$ be a locally
asymptotically normal family with precision $V\succeq 0$ and let
$\\{\Gamma_{k}\colon U\rightarrow{\mathbf{R}}^{d}\\}$ be regular sequence with
derivative $D$. Let $\mathcal{L}\colon{\mathbf{R}}^{d}\rightarrow[0,\infty)$
be symmetric, quasiconvex, and lower semicontinuous. Then, for any sequence of
estimators $T_{k}\colon\Omega_{k}\rightarrow{\mathbf{R}}^{d}$, we have

$\sup_{U_{0}\subset U,|U_{0}|<\infty}\liminf_{k\rightarrow\infty}\max_{u\in
U_{0}}\,{\operatorname{\mathbb{E}}}_{Q_{k,u}}\big{[}\mathcal{L}\big{(}\sqrt{k}\big{(}T_{k}-\Gamma_{k}(u)\big{)}\big{)}\big{]}\geq{\operatorname{\mathbb{E}}}\big{[}\mathcal{L}(Z)\big{]},$
(18)

where $Z\sim\mathsf{N}\big{(}0,DV^{{-1}}D^{\top}\big{)}$ when $V$ is
invertible; if $V$ is singular, then (18) holds with
$Z\sim\mathsf{N}\big{(}0,D(V+\lambda I)^{{-1}}D^{\top}\big{)}$ for all
$\lambda>0$.

We take $\Omega_{k}=\mathcal{Z}^{k}$ and set
$\mathcal{F}_{k}=\mathcal{B}(\Omega_{k})$ to be its Borel $\sigma$-algebra. To
prove Theorem 5.7 we will apply Theorem 5.10 to the sequence of probability
measures $Q_{k,u}:=P_{k,u/\sqrt{k}}$ and regular mapping sequence
$\Gamma_{k}\colon\mathcal{U}\rightarrow\mathbf{R}^{d}$ given by $u\mapsto
x^{\star}_{u/\sqrt{k}}$, where $\mathcal{U}$ is the neighborhood of zero
furnished by Lemma 5.6. We now state two key lemmas that will allow us to
apply Theorem 5.10; we defer their proofs to Appendix C.2.

The following lemma shows that $\Gamma_{k}$ constitutes a regular mapping
sequence.

######  Lemma 5.11.

The sequence $\\{\Gamma_{k}\colon\mathcal{U}\rightarrow\mathbf{R}^{d}\\}$ is
regular with derivative $-W^{{-1}}\Sigma_{g,G}^{\top}$, where

$\Sigma_{g,G}:=\underset{z\sim\mathcal{D}_{x^{\star}}}{\operatorname{\mathbb{E}}}\big{[}g_{x^{\star}}(z)G(x^{\star},z)^{\top}\big{]}.$

Furthermore, the family of distributions $Q_{k,u}$ is locally asymptotically
normal.

######  Lemma 5.12.

The family $\\{\Omega_{k},\mathcal{F}_{k},Q_{{k,u}}\\}_{u\in\mathbf{R}^{d}}$
is locally asymptotically normal with precision

$\Sigma_{g}:=\underset{z\sim\mathcal{D}_{x^{\star}}}{\mathop{\operatorname{\mathbb{E}}}}\big{[}g_{x^{\star}}(z)g_{x^{\star}}(z)^{\top}\big{]}.$

We now apply Theorem 5.10. Let
$\mathcal{L}\colon{\mathbf{R}}^{d}\rightarrow[0,\infty)$ be symmetric,
quasiconvex, and lower semicontinuous,
$\widehat{x}_{k}\colon\mathcal{Z}^{k}\rightarrow{\mathbf{R}}^{d}$ be a
sequence of estimators, and $c>0$. To connect the inequality (18) to (14),
take $U=\left\\{u\in\mathcal{U}\mid\|u\|\leq\sqrt{c}\right\\}$ and observe
that for any finite subset $U_{0}\subset U$, we have

$\liminf_{k\rightarrow\infty}\sup_{\|u\|\leq\sqrt{c/k}}{\operatorname{\mathbb{E}}}_{P_{k,u}}\big{[}\mathcal{L}\big{(}\sqrt{k}(\widehat{x}_{k}-x^{\star}_{u})\big{)}\big{]}\geq\liminf_{k\rightarrow\infty}\,\max_{u\in
U_{0}}\,{\operatorname{\mathbb{E}}}_{Q_{k,u}}\big{[}\mathcal{L}\big{(}\sqrt{k}\big{(}\widehat{x}_{k}-\Gamma_{k}(u)\big{)}\big{)}\big{]}.$

Taking the supremum over all finite $U_{0}\subset U$ and applying Theorem 5.10
with Lemmas 5.11 and 5.12 with $g=\xi$ (noting
$\Sigma_{\xi,G}=\Sigma_{\xi}=\Sigma$) yields

$\liminf_{k\rightarrow\infty}\sup_{\|u\|\leq\sqrt{c/k}}{\operatorname{\mathbb{E}}}_{P_{k,u}}\big{[}\mathcal{L}\big{(}\sqrt{k}(\widehat{x}_{k}-x^{\star}_{u})\big{)}\big{]}\geq{\operatorname{\mathbb{E}}}\big{[}\mathcal{L}(Z_{\lambda})\big{]}$
(19)

where $Z_{\lambda}\sim\mathsf{N}\big{(}0,W^{{-1}}\Sigma(\Sigma+\lambda
I)^{{-1}}\Sigma W^{-\top}\big{)}$ for any $\lambda>0$. Letting
$\lambda\downarrow 0$ in (19) establishes (14). Indeed, let $\Sigma=AA^{\top}$
be a Cholesky decomposition of $\Sigma$ and observe that the pseudoinverse
identities $A^{\dagger}=\lim_{\lambda\downarrow
0}A^{\top}\big{(}AA^{\top}+\lambda I\big{)}^{-1}$ and $AA^{\dagger}A=A$ imply

$\lim_{\lambda\downarrow 0}\Sigma\big{(}\Sigma+\lambda
I\big{)}^{{-1}}\Sigma=A\Big{(}\lim_{\lambda\downarrow
0}A^{\top}\big{(}AA^{\top}+\lambda
I\big{)}^{{-1}}\Big{)}AA^{\top}=\big{(}AA^{\dagger}A\big{)}A^{\top}=AA^{\top}=\Sigma.$

Thus, upon setting
$\widetilde{\Sigma}_{\lambda}:=W^{{-1}}\Sigma(\Sigma+\lambda I)^{{-1}}\Sigma
W^{-\top}$ and $\widetilde{\Sigma}:=W^{{-1}}\Sigma W^{-\top}$, we have
$\widetilde{\Sigma}_{\lambda}\rightarrow\widetilde{\Sigma}$ as
$\lambda\downarrow 0$. Further, for all $0<\lambda_{1}\leq\lambda_{2}$, we
have $\widetilde{\Sigma}_{\lambda_{1}}\succeq\widetilde{\Sigma}_{\lambda_{2}}$
and hence
${\exp}\big{(}{-}\tfrac{1}{2}u^{\top}\widetilde{\Sigma}_{\lambda_{1}}^{\dagger}u\big{)}\geq{\exp}\big{(}{-}\tfrac{1}{2}u^{\top}\widetilde{\Sigma}_{\lambda_{2}}^{\dagger}u\big{)}$
for all $u\in\mathbf{R}^{d}$. Since the densities corresponding to
$Z_{\lambda}\sim\mathsf{N}\big{(}0,\widetilde{\Sigma}_{\lambda}\big{)}$ and
$Z\sim\mathsf{N}\big{(}0,\widetilde{\Sigma}\big{)}$ with respect to the
Lebesgue measure restricted to $S:=\operatorname{range}\widetilde{\Sigma}$ are
given by

$p_{\lambda}(u):=\frac{{\exp}\big{(}{-}\tfrac{1}{2}u^{\top}\widetilde{\Sigma}_{\lambda}^{\dagger}u\big{)}}{\sqrt{(2\pi)^{r}\,{\det}^{*}\big{(}\widetilde{\Sigma}_{\lambda}\big{)}}}\qquad\text{and}\qquad
p(u):=\frac{{\exp}\big{(}{-}\tfrac{1}{2}u^{\top}\widetilde{\Sigma}^{\dagger}u\big{)}}{\sqrt{(2\pi)^{r}\,{\det}^{*}\big{(}\widetilde{\Sigma}\big{)}}}$

where $r$ is the rank of $\Sigma$, we may therefore apply the monotone
convergence theorem to obtain

$\displaystyle\lim_{\lambda\downarrow
0}{\operatorname{\mathbb{E}}}\big{[}\mathcal{L}(Z_{\lambda})\big{]}$
$\displaystyle=\lim_{\lambda\downarrow
0}\frac{1}{\sqrt{(2\pi)^{r}\,{\det}^{*}\big{(}\widetilde{\Sigma}_{\lambda}\big{)}}}\int_{S}\mathcal{L}(u)\,{\exp}\big{(}{-}\tfrac{1}{2}u^{\top}\widetilde{\Sigma}_{\lambda}^{\dagger}u\big{)}\,du$
$\displaystyle=\frac{1}{\sqrt{(2\pi)^{r}\,{\det}^{*}\big{(}\widetilde{\Sigma}\big{)}}}\int_{S}\mathcal{L}(u)\,{\exp}\big{(}{-}\tfrac{1}{2}u^{\top}\widetilde{\Sigma}^{\dagger}u\big{)}\,du$
$\displaystyle={\operatorname{\mathbb{E}}}\big{[}\mathcal{L}(Z)\big{]}.$

Hence (19) entails (14) and the proof of Theorem 5.7 is complete.

## References

  * Ahmed [2000] S. Ahmed.  _Strategic planning under uncertainty: Stochastic integer programming approaches_.  University of Illinois at Urbana-Champaign, 2000.
  * Bechavod et al. [2020] Y. Bechavod, K. Ligett, Z. S. Wu, and J. Ziani.  Causal feature discovery through strategic modification.  _arXiv:2002.07024_ , 2020.
  * Brown et al. [2020] G. Brown, S. Hod, and I. Kalemaj.  Performative prediction in a stateful world.  _arXiv preprint arXiv:2011.03885_ , 2020.
  * Cutler et al. [2021] J. Cutler, D. Drusvyatskiy, and Z. Harchaoui.  Stochastic optimization under time drift: iterate averaging, step-decay schedules, and high probability guarantees.  _Advances in Neural Information Processing Systems_ , 34, 2021.
  * Davis et al. [2021] D. Davis, D. Drusvyatskiy, and L. Jiang.  Subgradient methods near active manifolds: saddle point avoidance, local convergence, and asymptotic normality.  _arXiv preprint arXiv:2108.11832_ , 2021.
  * Dong et al. [2018] J. Dong, A. Roth, Z. Schutzman, B. Waggoner, and Z. S. Wu.  Strategic classification from revealed preferences.  In _Proceedings of the 2018 ACM Conference on Economics and Computation_ , pages 55–70, 2018.
  * Drusvyatskiy and Xiao [2020] D. Drusvyatskiy and L. Xiao.  Stochastic optimization with decision-dependent distributions.  _arXiv:2011.11173_ , 2020.
  * Duchi and Ruan [2021] J. C. Duchi and F. Ruan.  Asymptotic optimality in stochastic optimization.  _The Annals of Statistics_ , 49(1):21–48, 2021\.
  * Duflo [1997] M. Duflo.  _Random iterative models_ , volume 34.  Berlin: Springer, 1997.  ISBN 3-540-57100-0.
  * Dupacová [2006] J. Dupacová.  Optimization under exogenous and endogenous uncertainty.  _University of West Bohemia in Pilsen_ , 2006.
  * Feller [1971] W. Feller.  _An introduction to probability theory and its applications. Vol II. 2nd ed._ Wiley Ser. Probab. Math. Stat. John Wiley & Sons, Hoboken, NJ, 1971.
  * Hardt et al. [2016] M. Hardt, N. Megiddo, C. H. Papadimitriou, and M. Wootters.  Strategic classification.  In _Proceedings of the 2016 ACM Conference on Innovations in Theoretical Computer Science_ , pages 111–122. ACM, 2016.
  * Hellemo et al. [2018] L. Hellemo, P. I. Barton, and A. Tomasgard.  Decision-dependent probabilities in stochastic programs with recourse.  _Computational Management Science_ , 15(3):369–395, 2018.
  * Izzo et al. [2021] Z. Izzo, L. Ying, and J. Zou.  How to learn when data reacts to your model: performative gradient descent.  In _International Conference on Machine Learning_ , pages 4641–4650. PMLR, 2021.
  * Jagadeesan et al. [2022] M. Jagadeesan, T. Zrnic, and C. Mendler-Dünner.  Regret minimization with performative feedback.  _arXiv preprint arXiv:2202.00628_ , 2022.
  * Jonsbråten et al. [1998] T. W. Jonsbråten, R. J. Wets, and D. L. Woodruff.  A class of stochastic programs withdecision dependent random elements.  _Annals of Operations Research_ , 82:83–106, 1998.
  * Le Cam et al. [2000] L. Le Cam, L. M. LeCam, and G. L. Yang.  _Asymptotics in statistics: some basic concepts_.  Springer Science & Business Media, 2000.
  * Mendler-Dünner et al. [2020] C. Mendler-Dünner, J. Perdomo, T. Zrnic, and M. Hardt.  Stochastic optimization for performative prediction.  In _Advances in Neural Information Processing Systems_ , volume 33, pages 4929–4939. Curran Associates, Inc., 2020.
  * Mendler-Dünner et al. [2020] C. Mendler-Dünner, J. C. Perdomo, T. Zrnic, and M. Hardt.  Stochastic optimization for performative prediction.  _arXiv preprint arXiv:2006.06887_ , 2020.
  * Miller et al. [2020] J. Miller, S. Milli, and M. Hardt.  Strategic classification is causal modeling in disguise.  In _International Conference on Machine Learning_ , pages 6917–6926. PMLR, 2020.
  * Miller et al. [2021] J. P. Miller, J. C. Perdomo, and T. Zrnic.  Outside the echo chamber: Optimizing the performative risk.  In _International Conference on Machine Learning_ , pages 7710–7720. PMLR, 2021.
  * Narang et al. [2022] A. Narang, E. Faulkner, D. Drusvyatskiy, M. Fazel, and L. J. Ratliff.  Multiplayer performative prediction: Learning in decision-dependent games.  _arXiv preprint arXiv:2201.03398_ , 2022.
  * Perdomo et al. [2020] J. Perdomo, T. Zrnic, C. Mendler-Dünner, and M. Hardt.  Performative prediction.  In _Proceedings of the 37th International Conference on Machine Learning_ , volume 119 of _Proceedings of Machine Learning Research_ , pages 7599–7609. PMLR, 2020.
  * Piliouras and Yu [2022] G. Piliouras and F.-Y. Yu.  Multi-agent performative prediction: From global stability and optimality to chaos.  _arXiv preprint arXiv:2201.10483_ , 2022.
  * Polyak and Juditsky [1992] B. T. Polyak and A. B. Juditsky.  Acceleration of stochastic approximation by averaging.  _SIAM Journal on Control and Optimization_ , 30(4):838–855, 1992.  ISSN 0363-0129.  doi: 10.1137/0330046.
  * Ray et al. [2022] M. Ray, L. J. Ratliff, D. Drusvyatskiy, and M. Fazel.  Decision-dependent risk minimization in geometrically decaying dynamic environments.  _Proceedings of the AAAI International Conference on Artificial Intelligence (AAAI)_ , 2022.
  * Robbins and Siegmund [1971] H. Robbins and D. Siegmund.  A convergence theorem for non negative almost supermartingales and some applications.  1971\.
  * Rubinstein and Shapiro [1993] R. Y. Rubinstein and A. Shapiro.  Sensitivity analysis and stochastic optimization by the score function method, 1993.
  * Van der Vaart [2000] A. W. Van der Vaart.  _Asymptotic statistics_ , volume 3.  Cambridge university press, 2000.
  * van der Vaart and Wellner [1996] A. W. van der Vaart and J. A. Wellner.  _Convolution and Minimax Theorems_ , pages 412–422.  Springer New York, New York, NY, 1996.  ISBN 978-1-4757-2545-2.  doi: 10.1007/978-1-4757-2545-2_37.  URL https://doi.org/10.1007/978-1-4757-2545-2_37.
  * Varaiya and Wets [1988] P. Varaiya and R.-B. Wets.  Stochastic dynamic optimization approaches and computation.  1988\.
  * Wainwright [2019] M. J. Wainwright.  _High-dimensional statistics: A non-asymptotic viewpoint_ , volume 48.  Cambridge University Press, 2019.
  * Wood and Dall’Anese [2022] K. Wood and E. Dall’Anese.  Stochastic saddle point problems with decision-dependent distributions.  _arXiv preprint arXiv:2201.02313_ , 2022.
  * Wood et al. [2021] K. Wood, G. Bianchin, and E. Dall’Anese.  Online projected gradient descent for stochastic optimization with decision-dependent distributions.  _IEEE Control Systems Letters_ , 6:1646–1651, 2021.

##  Appendix A Proofs from Sections 3 and 4

###  A.1 Proof of Lemma 3.1

We successively estimate

$\displaystyle\|G_{x}(y)-G_{x^{\prime}}(y)\|$
$\displaystyle=\left\|\mathop{\operatorname{\mathbb{E}}}_{z\sim\mathcal{D}(x)}G(y,z)-\mathop{\operatorname{\mathbb{E}}}_{z\sim\mathcal{D}(x^{\prime})}G(y,z)\right\|$
$\displaystyle=\sup_{\|v\|\leq
1}\left\\{\mathop{\operatorname{\mathbb{E}}}_{z\sim\mathcal{D}(x)}\langle
G(y,z),v\rangle-\mathop{\operatorname{\mathbb{E}}}_{z\sim\mathcal{D}(x^{\prime})}\langle
G(y,z),v\rangle\right\\}$ $\displaystyle\leq\beta\cdot
W_{1}\big{(}\mathcal{D}(x),\mathcal{D}(x^{\prime})\big{)}$ (20)
$\displaystyle\leq\beta\gamma\cdot\|x-x^{\prime}\|,$

where inequality (20) uses that the map $z\mapsto\langle G(y,z),v\rangle$ is
$\beta$-Lipschitz continuous and the definition of the $W_{1}$ distance. The
proof is complete.

###  A.2 Proof of Theorem 3.3

Fix two points $x$ and $x^{\prime}$ in $\mathcal{X}$ and set
$y:=\mathtt{Sol}(x)$ and $y^{\prime}:=\mathtt{Sol}(x^{\prime})$. Note that the
definition of the normal cone implies

$\langle G_{x}(y),y-y^{\prime}\rangle\leq 0\qquad\textrm{and}\qquad\langle
G_{x^{\prime}}(y^{\prime}),y^{\prime}-y\rangle\leq 0.$

Strong monotonicity therefore ensures

$\displaystyle\alpha\|y-y^{\prime}\|^{2}$ $\displaystyle\leq\langle
G_{x}(y)-G_{x}(y^{\prime}),y-y^{\prime}\rangle$ $\displaystyle\leq\langle
G_{x^{\prime}}(y^{\prime})-G_{x}(y^{\prime}),y-y^{\prime}\rangle$
$\displaystyle\leq\|G_{x^{\prime}}(y^{\prime})-G_{x}(y^{\prime})\|\cdot\|y-y^{\prime}\|\
$ $\displaystyle\leq\gamma\beta\|x-x^{\prime}\|\cdot\|y-y^{\prime}\|,$

where the last inequality follows from Lemma 3.1. Dividing through by
$\alpha\|y-y^{\prime}\|$ guarantees that $\mathtt{Sol}(\cdot)$ is indeed a
contraction on $\mathcal{X}$ with parameter $\frac{\gamma\beta}{\alpha}$. The
result follows immediately from the Banach fixed point theorem.

###  A.3 Proof of Proposition 4.1

We will use the following classical convergence result [9, Theorem 1.3.12].

######  Lemma A.1 (Robbins-Siegmund).

Let $(A_{t}),(B_{t}),(C_{t}),(D_{t})$ be sequences of nonnegative finite
random variables on a filtered probability space
$(\Omega,\mathcal{F},\mathbb{F},\mathbb{P})$ adapted to the filtration
$\mathbb{F}=(\mathcal{F}_{t})$ and satisfying

${\operatorname{\mathbb{E}}}[A_{t+1}\mid\mathcal{F}_{t}]\leq(1+B_{t})A_{t}+C_{t}-D_{t}$

for all $t$. Then on the event
$\left\\{\sum_{t}B_{t}<\infty,\sum_{t}C_{t}<\infty\right\\}$, there is a
nonnegative finite random variable $A_{\infty}$ such that $A_{t}\rightarrow
A_{\infty}$ and $\sum_{t}D_{t}<\infty$ almost surely.

Set $\Delta_{t}=x_{t}-x^{\star}$ and $\bar{\alpha}=\alpha-\gamma\beta$; note
that $\bar{\alpha}>0$ by Assumption 2. Let $v_{t}$ and $\xi_{t}$ be the
direction of motion and noise vectors given by (5), and let
$(\mathcal{F}_{t})$ be the filtration given by (8). Observe that Lemma 3.1 and
Assumption 3 immediately yield the bias/variance bounds

$\|{\operatorname{\mathbb{E}}}_{t}[v_{t}]-G_{x^{\star}}(x_{t})\|\leq\gamma\beta\|x_{t}-x^{\star}\|\qquad\text{and}\qquad{\operatorname{\mathbb{E}}}_{t}\|\xi_{t}\|^{2}\leq
K(1+\|x_{t}-x^{\star}\|^{2}),$ (21)

where $\operatorname{\mathbb{E}}_{t}[\cdot]$ denotes the conditional
expectation with respect to $\mathcal{F}_{t}$. In light of (21) and the fact
that $\eta_{t}\rightarrow 0$ (since $\sum_{t}\eta_{t}^{2}<\infty$), we see
that for all sufficiently large $t$, the result [22, Theorem 8 (Benign bias)]
directly implies the one-step improvement bound

$\displaystyle{\operatorname{\mathbb{E}}}_{t}\|\Delta_{t+1}\|^{2}$
$\displaystyle\leq\|\Delta_{t}\|^{2}+4K\eta^{2}_{t}-\tfrac{1}{2}\bar{\alpha}\eta_{t}\|\Delta_{t}\|^{2}.$
(22)

We now apply Lemma A.1 with $A_{t}=\|\Delta_{t}\|^{2}$, $B_{t}=0$,
$C_{t}=4K\eta_{t}^{2}$, and
$D_{t}=\tfrac{1}{2}\bar{\mu}\eta_{t}\|\Delta_{t}\|^{2}$. Notice that the event
$\left\\{\sum_{t}B_{t}<\infty,\sum_{t}C_{t}<\infty\right\\}$ holds since
$\sum_{t}\eta_{t}^{2}<\infty$. Thus, there exists a nonnegative finite random
variable $A_{\infty}$ such that $A_{t}\rightarrow A_{\infty}$ and
$\sum_{t}D_{t}<\infty$ almost surely. Hence $\|\Delta_{t}\|^{2}\rightarrow
A_{\infty}$ and $\sum_{t}\eta_{t}\|\Delta_{t}\|^{2}<\infty$ almost surely.
Since $\sum_{t}\eta_{t}$ diverges, we conclude that
$\lim_{t}\|\Delta_{t}\|^{2}=A_{\infty}=0$ almost surely.

Next, to establish the in-expectation rate, we note that by (22) and the tower
rule, we have

${\operatorname{\mathbb{E}}}\|\Delta_{t+1}\|^{2}\leq\big{(}1-\tfrac{1}{2}\bar{\alpha}\eta_{t}\big{)}{\operatorname{\mathbb{E}}}\|\Delta_{t}\|^{2}+4K\eta_{t}^{2}$

for all sufficiently large $t$. Thus, upon supposing
$\eta_{t}=\eta_{0}t^{-\nu}$ with $\nu\in\big{(}\frac{1}{2},1\big{)}$, a
standard inductive argument (see [5, Lemma 3.11.8]) yields a constant $C\geq
0$ such that ${\operatorname{\mathbb{E}}}\|\Delta_{t}\|^{2}\leq C\eta_{t}$ for
all $t\geq 0$. This completes the proof.

##  Appendix B Review of asymptotic normality

In this section, we briefly review the asymptotic normality results of Polyak
and Juditsky [25]. Consider a set $\mathcal{X}\subset\mathbf{R}^{d}$ and a map
$R\colon\mathcal{X}\rightarrow{\mathbf{R}}^{d}$. Suppose that there exists a
solution $x^{\star}\in\operatorname{int}\mathcal{X}$ to the equation $R(x)=0$.
The goal is to approximate $x^{\star}$ while only having access to noisy
evaluations of $R$. Given $x_{0}\in\mathcal{X}$, consider the iterative
process

$x_{t+1}=x_{t}-\eta_{t}\big{(}R(x_{t})+\xi_{t}+\zeta_{t}\big{)},$ (23)

where $\eta_{t}$ is a deterministic positive step size and $\xi_{t}$ and
$\zeta_{t}$ are random vectors in $\mathbf{R}^{d}$ such that
$x_{t}\in\mathcal{X}$ for all $t$. We impose the following three assumptions.

######  Assumption B.1 (Stochastic setting).

The sequences $(\xi_{t})_{t\geq 0}$ and $(\zeta_{t})_{t\geq 0}$ are stochastic
processes defined on a filtered probability space
$(\Omega,\mathcal{F},{\mathbb{F}},\mathbb{P})$ with filtration
${\mathbb{F}}=(\mathcal{F})_{t\geq 0}$ such that $\xi_{t}$ and $\zeta_{t}$ are
$\mathcal{F}_{t+1}$-measurable, and $\xi_{t}$ constitutes a martingale
difference process satisfying
$\operatorname{\mathbb{E}}[\xi_{t}\,|\,\mathcal{F}_{t}]=0$. Further, the
following four conditions hold.

  1. (i)

There is a constant $\bar{K}\geq 0$ such that

${\operatorname{\mathbb{E}}}\big{[}\|\xi_{t}\|^{2}\mid\mathcal{F}_{t}\big{]}+\|R(x_{t})\|^{2}\leq\bar{K}(1+\|x_{t}\|^{2})\qquad\text{for
all }t\geq 0.$

  2. (ii)

The variables $\sup_{t\geq
0}{\operatorname{\mathbb{E}}}\big{[}\|\xi_{t}\|^{2}\mathbf{1}\\{\|\xi_{t}\|>N\\}\mid\mathcal{F}_{t}\big{]}$
converge to zero in probability as $N\rightarrow\infty$.

  3. (iii)

The sequence
${\operatorname{\mathbb{E}}}\big{[}\xi_{t}\xi_{t}^{\top}\mid\mathcal{F}_{t}\big{]}$
converges to a matrix $\Sigma\succeq 0$ in probability as
$t\rightarrow\infty$.

  4. (iv)

With probability $1$, the vector $\zeta_{t}$ vanishes for all but finitely
many $t$.

######  Assumption B.2 (Convergence).

The iterates $x_{t}$ converge to $x^{\star}$ almost surely, and there exists a
constant $C\geq 0$ such that

${\operatorname{\mathbb{E}}}\|x_{t}-x^{\star}\|^{2}\leq
C\eta_{t}\qquad\text{for all }t\geq 0.$

######  Assumption B.3 (Smoothness).

The map $R$ is differentiable at $x^{\star}$, and the Jacobian $\nabla
R(x^{\star})$ satisfies

$R(x)-\nabla R(x^{\star})(x-x^{\star})=O(\|x-x^{\star}\|^{2})\qquad\text{as
}x\rightarrow x^{\star}.$

Further, the real parts of the eigenvalues of $\nabla R(x^{\star})$ are
positive, i.e., ${\min}_{j}\,\text{Re}\,\lambda_{j}\big{(}\nabla
R(x^{\star})\big{)}>0$.

######  Theorem B.1 ([25, Theorem 2]).

Suppose that Assumptions B.1, B.2, and B.3 hold. Further, suppose that the
step sizes satisfy $\frac{\eta_{t}-\eta_{t+1}}{\eta_{t}^{2}}\to 0$ and
$\sum_{t=1}^{\infty}\eta_{t}\cdot t^{-1/2}<\infty$. Then the average iterates
$\bar{x}_{t}=\frac{1}{t}\sum_{i=0}^{t-1}x_{i}$ converge to $x^{\star}$ almost
surely, and

$\sqrt{t}(\bar{x}_{t}-x^{\star})\xlongrightarrow{d}\mathsf{N}\big{(}0,\nabla
R(x^{\star})^{-1}\cdot\Sigma\cdot\nabla R(x^{\star})^{-\top}\big{)}.$

We remark that the assumptions of Theorem B.1 differ slightly from those in
[25]. The proof of the theorem follows by slightly modifying the arguments in
[25]. We provide a quick sketch of the small modifications that are needed.
There are two differences between the assumptions of Theorem B.1 and the one
in [25]: $(a)$ Assumption B.2 replaces a Lyapunov function assumption
(Assumption 3.1 in [25]) and $(b)$ our process (23) includes the additional
term $\zeta_{t}$. The proof of [25, Theorem 2] is divided into four parts:

##### Parts 1 and 2:

use the existence of a Lyapunov function to establish almost sure convergence
of $x_{t}$ to $x^{\star}$ and to derive an in-expectation rate of convergence.
Both of these conclusions follow directly from Assumption B.2.

##### Part 3:

considers a linearized version of the stochastic process (23) — which we dub
$x_{t}^{1}$ — and uses a martingale central limit theorem to show that it has
the desired asymptotic distribution. To do so, the error
$\sqrt{t}(\bar{x}_{t}^{1}-x^{\star})$ is decomposed into three sums, two of
which tend to zero in mean square and the third realizing the stated
asymptotic normal distribution. The addition of $\zeta_{t}$ is
straightforward: it results in extra sums involving $\zeta_{t}$ in the error
decomposition, and these additional sums tend to zero almost surely by virtue
of condition (iv) of Assumption B.1.

##### Part 4:

shows that the linearized stochastic process $x_{t}^{1}$ and (23) are
asymptotically equivalent. This step follows without changes since it is
solely based on the in-expectation rate obtained in Parts 1 and 2 and an
application of Kronecker’s lemma.

##  Appendix C Proofs from Section 5

This section presents all the missing proofs from Section 5. We assume
throughout that the assumptions used in Section 5 are valid; in particular,
$\mathcal{X}$ is compact, $\mathcal{Z}$ is bounded, and $g\in\mathcal{G}$.

######  Lemma C.1.

We have

$\sup_{x\in\mathcal{X},\,z\in\mathcal{Z}}\|g_{x}(z)\|<\infty\qquad\text{and}\qquad\sup_{x\in\mathcal{X},\,z\in\mathcal{Z}}\|\nabla_{x}g_{x}(z)\|_{\textup{op}}<\infty.$

###### Proof.

Fix $x^{\circ}\in\mathcal{X}$ and $z^{\circ}\in\mathcal{Z}$. Since
$\mathcal{X}$ and $\mathcal{Z}$ are bounded, we compute

$\displaystyle
M^{\prime}_{g}:=\sup_{x\in\mathcal{X},\,z\in\mathcal{Z}}\|\nabla_{x}g_{x}(z)\|_{\text{op}}$
$\displaystyle\leq\begin{aligned}
\|\nabla_{x}g_{x^{\circ}}(z^{\circ})\|_{\text{op}}&+\sup_{x\in\mathcal{X}}\|\nabla_{x}g_{x}(z^{\circ})-\nabla_{x}g_{x^{\circ}}(z^{\circ})\|_{\text{op}}\\\
&+\sup_{x\in\mathcal{X},\,z\in\mathcal{Z}}\|\nabla_{x}g_{x}(z)-\nabla_{x}g_{x}(z^{\circ})\|_{\text{op}}\end{aligned}$
$\displaystyle\leq\|\nabla_{x}g_{x^{\circ}}(z^{\circ})\|_{\text{op}}+\Lambda_{g}(z^{\circ})\operatorname{diam}(\mathcal{X})+\beta^{\prime}_{g}\operatorname{diam}(\mathcal{Z})<\infty.$

Hence every section $g(\cdot,z)$ is $M^{\prime}_{g}$-Lipschitz on
$\mathcal{X}$, and the estimate

$\displaystyle M_{g}:=\sup_{x\in\mathcal{X},\,z\in\mathcal{Z}}\|g_{x}(z)\|$
$\displaystyle\leq\|g_{x^{\circ}}(z^{\circ})\|+\sup_{x\in\mathcal{X}}\|g_{x}(z^{\circ})-g_{x^{\circ}}(z^{\circ})\|+\sup_{x\in\mathcal{X},\,z\in\mathcal{Z}}\|g_{x}(z)-g_{x}(z^{\circ})\|$
$\displaystyle\leq\|g_{x^{\circ}}(z^{\circ})\|+M^{\prime}_{g}\operatorname{diam}(\mathcal{X})+\beta_{g}\operatorname{diam}(\mathcal{Z})$

completes the proof. ∎

######  Lemma C.2.

For any $x\in\mathcal{X}$, the function
$C_{x}(u):=1+\mathop{\operatorname{\mathbb{E}}}_{z\sim\mathcal{D}_{x}}h\big{(}u^{\top}g_{x}(z)\big{)}$
is $C^{2}$-smooth on $\mathbf{R}^{d}$ with $A_{h}M_{g}^{3}$-Lipschitz
continuous Hessian, where $A_{h}:=\sup|h^{\prime\prime\prime}|$ and
$M_{g}:=\sup_{x\in\mathcal{X},\,z\in\mathcal{Z}}\|g_{x}(z)\|$. Moreover, we
have
$\sup_{x\in\mathcal{X}}\big{|}\mathop{\operatorname{\mathbb{E}}}_{z\sim\mathcal{D}_{x}}h\big{(}u^{\top}g_{x}(z)\big{)}\big{|}=o(\|u\|^{2})$
as $u\to 0$.

###### Proof.

Given $x\in\mathcal{X}$, the dominated convergence theorem yields

$\nabla_{u}\Big{(}\mathop{\operatorname{\mathbb{E}}}_{z\sim\mathcal{D}_{x}}h\big{(}u^{\top}g_{x}(z)\big{)}\Big{)}=\mathop{\operatorname{\mathbb{E}}}_{z\sim\mathcal{D}_{x}}h^{\prime}\big{(}u^{\top}g_{x}(z)\big{)}g_{x}(z)$

and

$\nabla_{u}^{2}\Big{(}\mathop{\operatorname{\mathbb{E}}}_{z\sim\mathcal{D}_{x}}h\big{(}u^{\top}g_{x}(z)\big{)}\Big{)}=\mathop{\operatorname{\mathbb{E}}}_{z\sim\mathcal{D}_{x}}h^{\prime\prime}\big{(}u^{\top}g_{x}(z)\big{)}g_{x}(z)g_{x}(z)^{\top}.$

Thus, $C_{x}(\cdot)$ is $C^{2}$-smooth on ${\mathbf{R}}^{d}$, and the claimed
Lipschitz constant of the Hessian follows trivially.

To prove the final claim of the lemma, recall that $h^{\prime}(0)=1$ and
$h^{\prime\prime}(0)=0$. Therefore

$\nabla_{u}\Big{(}\mathop{\operatorname{\mathbb{E}}}_{z\sim\mathcal{D}_{x}}h\big{(}u^{\top}g_{x}(z)\big{)}\Big{)}\Big{|}_{u=0}=\mathop{\operatorname{\mathbb{E}}}_{z\sim\mathcal{D}_{x}}g_{x}(z)=0\qquad\text{and}\qquad\nabla_{u}^{2}\Big{(}\mathop{\operatorname{\mathbb{E}}}_{z\sim\mathcal{D}_{x}}h\big{(}u^{\top}g_{x}(z)\big{)}\Big{)}\Big{|}_{u=0}=0.$

Consequently, the second-order Taylor polynomial of
$\mathop{\operatorname{\mathbb{E}}}_{z\sim\mathcal{D}_{x}}h\big{(}u^{\top}g_{x}(z)\big{)}$
in $u$ based at $0$ is identically zero. Lipschitzness of the Hessian
therefore implies

$\displaystyle\Big{|}\mathop{\operatorname{\mathbb{E}}}_{z\sim\mathcal{D}_{x}}h\big{(}u^{\top}g_{x}(z)\big{)}\Big{|}\leq\frac{A_{h}M_{g}^{3}}{6}\|u\|^{3}$

for all $x\in\mathcal{X}$ and $u\in\mathbf{R}^{d}$. Hence
$\sup_{x\in\mathcal{X}}\big{|}\mathop{\operatorname{\mathbb{E}}}_{z\sim\mathcal{D}_{x}}h\big{(}u^{\top}g_{x}(z)\big{)}\big{|}=o(\|u\|^{2})$
as $u\to 0$. ∎

######  Lemma C.3.

Let $f\colon(0,\infty)\rightarrow{\mathbf{R}}$ be a function that is
$C^{3}$-smooth around $t=1$ and satisfies $f(1)=f^{\prime}(1)=0.$ Then for all
$x\in\mathcal{X}$ and sufficiently small $u\in\mathbf{R}^{d}$, we have

$\int
f\bigg{(}\frac{1+u^{\top}g_{x}(z)}{C_{x}^{u}}\bigg{)}\,d\mathcal{D}_{x}(z)=\frac{f^{\prime\prime}(1)}{2}u^{\top}\bigg{(}\mathop{\operatorname{\mathbb{E}}}_{z\sim\mathcal{D}_{x}}g_{x}(z)g_{x}(z)^{\top}\bigg{)}u+r_{x}(u)$
(24)

with $r_{x}\colon{\mathbf{R}}^{d}\rightarrow{\mathbf{R}}$ satisfying
$\sup_{x\in\mathcal{X}}|r_{x}(u)|=o(\|u\|^{2})$ as $u\rightarrow 0$.

###### Proof.

Fix $x\in\mathcal{X}$ and define
$\varphi_{x}(u):={\operatorname{\mathbb{E}}}_{z\sim\mathcal{D}_{x}}f\Big{(}\frac{1+u^{\top}g_{x}(z)}{C_{x}^{u}}\Big{)}$.
By the dominated convergence theorem,

$\nabla_{u}\varphi_{x}(u)={\mathop{\operatorname{\mathbb{E}}}_{z\sim\mathcal{D}_{x}}}\bigg{[}f^{\prime}\bigg{(}\frac{1+u^{\top}g_{x}(z)}{C_{x}^{u}}\bigg{)}\bigg{(}\frac{g_{x}(z)C_{x}^{u}-\big{(}1+u^{\top}g_{x}(z)\big{)}\nabla_{u}C_{x}^{u}}{(C_{x}^{u})^{2}}\bigg{)}\bigg{]}.$

Similarly, the Hessian $\nabla_{uu}^{2}\varphi_{x}(u)$ of $\varphi_{x}$ at $u$
is equal to

$\displaystyle\frac{1}{(C_{x}^{u})^{4}}\mathop{\operatorname{\mathbb{E}}}_{z\sim\mathcal{D}_{x}}\left[f^{\prime\prime}\bigg{(}\frac{1+u^{\top}g_{x}(z)}{C_{x}^{u}}\bigg{)}\Big{(}g_{x}(z)C_{x}^{u}-\big{(}1+u^{\top}g_{x}(z)\big{)}\nabla_{u}C_{x}^{u}\Big{)}\Big{(}g_{x}(z)C_{x}^{u}-\big{(}1+u^{\top}g_{x}(z)\big{)}\nabla_{u}C_{x}^{u}\Big{)}^{\top}\right.$
$\displaystyle\hskip
62.59596pt+\left.f^{\prime}\bigg{(}\frac{1+u^{\top}g_{x}(z)}{C_{x}^{u}}\bigg{)}\bigg{(}(C_{x}^{u})^{2}\Big{(}g_{x}(z)(\nabla_{u}C_{x}^{u})^{\top}-(\nabla_{u}C_{x}^{u})g_{x}(z)^{\top}-\big{(}1+u^{\top}g_{x}(z)\big{)}\nabla^{2}_{uu}C_{x}^{u}\Big{)}\right.$
$\displaystyle\hskip
170.71652pt-\left.2C_{x}^{u}\Big{(}g_{x}(z)C_{x}^{u}-\big{(}1+u^{\top}g_{x}(z)\big{)}\nabla_{u}C_{x}^{u}\Big{)}(\nabla_{u}C_{x}^{u})^{\top}\bigg{)}\right].$

Thus, $\varphi_{x}$ is $C^{2}$-smooth. Taking a Taylor expansion at $u=0$ and
using the equalities $C_{x}^{0}=1$, $\nabla_{u}C_{x}^{u}|_{u=0}=0$, and
$f(1)=f^{\prime}(1)=0$ yields (24) with second-order Taylor remainder
$r_{x}\colon{\mathbf{R}}^{d}\rightarrow{\mathbf{R}}$. It remains to verify
that $\sup_{x\in\mathcal{X}}|r_{x}(u)|=o(\|u\|^{2})$ as $u\rightarrow 0$.

Lemmas C.1 and C.2 ensure that $C_{x}^{u},\nabla_{u}C_{x}^{u},$ and
$\nabla_{uu}^{2}C_{x}^{u}$ are Lipschitz and bounded on a compact neighborhood
of $u=0$, with Lipschitz constants and bounds independent of $x$. Further,
since $f$ is $C^{3}$-smooth around $t=1$, we have that $f^{\prime}$ and
$f^{\prime\prime}$ are Lipschitz and bounded on a neighborhood of $t=1$. It
follows that $\nabla^{2}_{uu}\varphi_{x}$ is $\tilde{L}$-Lipschitz on a
neighborhood $U$ of $u=0$ with constant $\tilde{L}$ independent of $x$. Thus
we deduce $|r_{x}(u)|\leq\frac{\tilde{L}}{6}\|u\|^{3}$ for all
$(x,u)\in\mathcal{X}\times U$, and the result follows. ∎

###  C.1 Proof of Lemma 5.6

The proof of this lemma is divided into four steps: the first step verifies
Assumption 1 and the next three steps establish Assumption 2. The strategy in
all steps is to prove that various quantities of interest change continuously
with $u$ near zero. The main tool we will use to infer continuity is the
following elementary lemma. Its proof consists of several applications of the
dominated convergence theorem and is deferred to Appendix D.

######  Lemma C.4 (Inferring smoothness).

Suppose that $n\in{\mathbf{N}}$ and
$T\colon\mathcal{X}\times\mathcal{Z}\rightarrow{\mathbf{R}}^{n}$ is a map
satisfying the following two conditions.

  1. 1.

There exists $\beta_{T}>0$ such that for every $x\in\mathcal{X}$, the section
$T(x,\cdot)$ is $\beta_{T}$-Lipschitz on $\mathcal{Z}$.

  2. 2.

There exist a measurable function
$\Lambda_{T}\colon\mathcal{Z}\rightarrow[0,\infty)$ and constants
$\bar{\Lambda}_{T},\beta^{\prime}_{T}>0$ such that for every $z\in\mathcal{Z}$
and $x\in\mathcal{X}$, the section $T(\cdot,z)$ is $\Lambda_{T}(z)$-smooth on
$\mathcal{X}$ with
$\operatorname{\mathbb{E}}_{z\sim\mathcal{D}_{x}}[\Lambda_{T}(z)]\leq\bar{\Lambda}_{T}$,
and the section $\nabla_{x}T(x,\cdot)$ is $\beta^{\prime}_{T}$-Lipschitz on
$\mathcal{Z}$.

Set

$M_{T}:=\sup_{x\in\mathcal{X},\,z\in\mathcal{Z}}\|T(x,z)\|\qquad\text{and}\qquad
M^{\prime}_{T}:=\sup_{x\in\mathcal{X},\,z\in\mathcal{Z}}\|\nabla_{x}T(x,z)\|_{\textup{op}}.$

Then $M_{T}$ and $M^{\prime}_{T}$ are finite. Further, given any compact
subset $\mathcal{W}\subset\mathbf{R}^{d}$, the map
$H\colon\mathcal{X}\times\mathcal{W}\rightarrow\mathbf{R}^{n}$ given by

$H(x,u)=\mathop{\operatorname{\mathbb{E}}}_{z\sim\mathcal{D}_{x}^{u}}T(x,z)$

is $\Lambda_{H}$-smooth, where $\Lambda_{H}$ is a constant depending on $T$
only through $\beta_{T},\bar{\Lambda}_{T},\beta^{\prime}_{T}$, $M_{T}$, and
$M^{\prime}_{T}$.

##### Step 1 (Assumption 1)

First, we show that the perturbed distribution map $\mathcal{D}^{u}$ satisfies
Assumption 1 with Lipschitz constant $\gamma^{u}=\gamma+O(\|u\|)$ as
$u\rightarrow 0$, where $\gamma$ is the Lipschitz constant for $\mathcal{D}$.
To this end, we take $\mathcal{W}$ to be the unit ball in $\mathbf{R}^{d}$ and
apply Lemma C.4 to identify a constant $L_{1}$ such that for every
$1$-Lipschitz function $\phi\colon\mathcal{Z}\rightarrow{\mathbf{R}}$ and
every $u\in\mathcal{W}$, the function

$\rho_{\phi}(x,u):=\mathop{\operatorname{\mathbb{E}}}_{z\sim\mathcal{D}_{x}^{u}}\phi(z)$

is Lipschitz in the $x$-component with constant
$\gamma^{u}:=\gamma+L_{1}\|u\|$. Without loss of generality, we may assume
that $\phi$ is bounded by $\operatorname{diam}(\mathcal{Z})$, since we can
always shift the function $\phi\leftarrow\phi-\inf\phi$ and this will not
change the Lipschitz constant of $\rho_{\phi}(\cdot,u)$. Then Lemma C.4 yields
a constant $L_{1}$ independent of $\phi$ such that $\rho_{\phi}$ is
$L_{1}$-smooth on $\mathcal{X}\times\mathcal{W}$.

To establish the Lipschitz constant of $\rho_{\phi}(\cdot,u)$, we derive a
bound on the gradient of $\rho_{\phi}$. By the triangle inequality, the
following holds for all $(x,u)\in\mathcal{X}\times\mathcal{W}$:

$\displaystyle\|\nabla_{x}\rho_{\phi}(x,u)\|$
$\displaystyle\leq\|\nabla_{x}\rho_{\phi}(x,0)\|+\|\nabla_{x}\rho_{\phi}(x,u)-\nabla_{x}\rho_{\phi}(x,0)\|$
$\displaystyle\leq\gamma+L_{1}\|u\|.$

Therefore $\mathcal{D}^{u}$ satisfies Assumption 1 with Lipschitz constant
$\gamma^{u}=\gamma+O(\|u\|)$ as $u\rightarrow 0$.

##### Step 2 (Lipschitz)

Next, we establish Assumption 2(i) for the problem with the perturbed
distribution map $\mathcal{D}^{u}$. Observe that the Lipschitz bounds in
Assumption 2(i) remain unchanged, and that we only need to identify for all
sufficiently small $u$ a constant $\bar{L}^{u}>0$ such that
$\sup_{x\in\mathcal{X}}{\operatorname{\mathbb{E}}}_{z\sim\mathcal{D}_{x}^{u}}\big{[}L(z)^{2}\big{]}\leq(\bar{L}^{u})^{2}$.
We will show that we can select $(\bar{L}^{u})^{2}=\bar{L}^{2}+O(\|u\|)$ as
$u\rightarrow 0$, where $\bar{L}$ is the constant satisfying
$\sup_{x\in\mathcal{X}}{\operatorname{\mathbb{E}}}_{z\sim\mathcal{D}_{x}}\big{[}L(z)^{2}\big{]}\leq\bar{L}^{2}$
furnished by Assumption 2(i). To this end, observe first that for all
$x\in\mathcal{X}$ and $u\in\mathbf{R}^{d}$, we have

${\mathop{\operatorname{\mathbb{E}}}_{z\sim\mathcal{D}_{x}^{u}}}\big{[}L(z)^{2}\big{]}=\frac{1}{C_{x}^{u}}\,{\mathop{\operatorname{\mathbb{E}}}_{z\sim\mathcal{D}_{x}}}\big{[}L(z)^{2}\big{(}1+h\big{(}u^{\top}g_{x}(z)\big{)}\big{)}\big{]}.$

Using the Taylor series of $(1+t)^{-1}$ for $|t|<1$, it follows from Lemma C.2
that

$\sup_{x\in\mathcal{X}}\frac{1}{C_{x}^{u}}=1+o(\|u\|^{2})\qquad\text{as~{}}u\rightarrow
0.$ (25)

Moreover, upon noting that $h(0)=0$ and the function $u\mapsto
h\big{(}u^{\top}g_{x}(z)\big{)}$ is Lipschitz on $\mathbf{R}^{d}$ with
constant $(\sup|h^{\prime}|)M_{g}$ for each $x\in\mathcal{X}$ and
$z\in\mathcal{Z}$, we obtain the following for all $u\in\mathbf{R}^{d}$:

$\sup_{x\in\mathcal{X},\,z\in\mathcal{Z}}\big{|}h\big{(}u^{\top}g_{x}(z)\big{)}\big{|}\leq(\sup|h^{\prime}|)M_{g}\|u\|$

and hence

$\sup_{x\in\mathcal{X}}{\mathop{\operatorname{\mathbb{E}}}_{z\sim\mathcal{D}_{x}}}\Big{[}L(z)^{2}\big{(}1+h\big{(}u^{\top}g_{x}(z)\big{)}\big{)}\Big{]}\leq\bar{L}^{2}\cdot\big{(}1+(\sup|h^{\prime}|)M_{g}\|u\|\big{)}.$

We deduce

$\sup_{x\in\mathcal{X}}{\mathop{\operatorname{\mathbb{E}}}_{z\sim\mathcal{D}_{x}^{u}}}\big{[}L(z)^{2}\big{]}\leq\big{(}1+o(\|u\|^{2})\big{)}\cdot\bar{L}^{2}\cdot\big{(}1+(\sup|h^{\prime}|)M_{g}\|u\|\big{)}=\bar{L}^{2}+O(\|u\|)\qquad\text{as~{}}u\rightarrow
0.$

##### Step 3 (Monotonicity)

We prove that for all $x\in\mathcal{X}$, the map $G^{u}_{x}(\cdot)$ is
strongly monotone on $\mathcal{X}$ with constant
$\alpha^{u}=\alpha+o(\|u\|^{2})$ as $u\rightarrow 0$, where $\alpha$ is the
strong monotonicity constant of $G_{x}(\cdot)$. Given $x\in\mathcal{X}$ and
$u\in\mathbf{R}^{d}$, we have

$\displaystyle\big{\langle}G_{x}^{u}(y)-G_{x}^{u}(y^{\prime}),y-y^{\prime}\big{\rangle}$
$\displaystyle=\big{\langle}G_{x}(y)-G_{x}(y^{\prime}),y-y^{\prime}\big{\rangle}+\big{\langle}\big{(}G_{x}^{u}(y)-G_{x}(y)\big{)}-\big{(}G_{x}^{u}(y^{\prime})-G_{x}(y^{\prime})\big{)},y-y^{\prime}\big{\rangle}$
$\displaystyle\geq\alpha\|y-y^{\prime}\|^{2}-\big{\|}\big{(}G_{x}^{u}(y)-G_{x}(y)\big{)}-\big{(}G_{x}^{u}(y^{\prime})-G_{x}(y^{\prime})\big{)}\big{\|}\cdot\big{\|}y-y^{\prime}\big{\|}$

for all $y,y^{\prime}\in\mathcal{X}$ by the $\alpha$-strong monotonicity of
$G_{x}(\cdot)$. We claim that for all sufficiently small $u$, there exists
$\ell^{u}=o(\|u\|^{2})$ independent of $x$ such that the map $y\mapsto
G_{x}^{u}(y)-G_{x}(y)$ is $\ell^{u}$-Lipschitz on $\mathcal{X}$ for all
$x\in\mathcal{X}$. Indeed, upon noting
$\sup_{x\in\mathcal{X},\,z\in\mathcal{Z}}\|\nabla_{x}G(x,z)\|_{\textup{op}}<\infty$
(as in Lemma C.4) and applying the dominated convergence theorem together with
(25) and Lemma C.2, we obtain

$\displaystyle\ell^{u}$
$\displaystyle:=\sup_{x,y\in\mathcal{X}}\big{\|}\nabla_{y}\big{(}G_{x}^{u}(y)-G_{x}(y)\big{)}\big{\|}_{\text{op}}$
$\displaystyle=\sup_{x,y\in\mathcal{X}}\left\|\frac{1}{C_{x}^{u}}\underset{z\sim\mathcal{D}_{x}}{\operatorname{\mathbb{E}}}\big{[}\nabla_{y}G(y,z)\big{(}1+h\big{(}u^{\top}g_{x}(z)\big{)}\big{)}\big{]}-\underset{z\sim\mathcal{D}_{x}}{\operatorname{\mathbb{E}}}\big{[}\nabla_{y}G(y,z)\big{]}\right\|_{\text{op}}$
$\displaystyle\leq\underbrace{\sup_{x,y\in\mathcal{X}}\left\|\frac{1}{C_{x}^{u}}\underset{z\sim\mathcal{D}_{x}}{\operatorname{\mathbb{E}}}\big{[}\nabla_{y}G(y,z)\big{]}-\underset{z\sim\mathcal{D}_{x}}{\operatorname{\mathbb{E}}}\big{[}\nabla_{y}G(y,z)\big{]}\right\|_{\text{op}}}_{o(\|u\|^{2})}+\underbrace{\sup_{x,y\in\mathcal{X}}\left\|\frac{1}{C_{x}^{u}}\underset{z\sim\mathcal{D}_{x}}{\operatorname{\mathbb{E}}}\big{[}\nabla_{y}G(y,z)h\big{(}u^{\top}g_{x}(z)\big{)}\big{]}\right\|_{\text{op}}}_{\left(1+o(\|u\|^{2})\right)\cdot
o(\|u\|^{2})}$ $\displaystyle=o(\|u\|^{2})\qquad\text{as~{}}u\rightarrow 0.$

Setting $\alpha^{u}:=\alpha-\ell^{u}$ for all $u$ in a neighborhood of zero,
we conclude that for all $x,y,y^{\prime}\in\mathcal{X}$,

$\big{\langle}G_{x}^{u}(y)-G_{x}^{u}(y^{\prime}),y-y^{\prime}\big{\rangle}\geq\alpha^{u}\|y-y^{\prime}\|^{2}$

and hence $G^{u}_{x}(\cdot)$ is strongly monotone on $\mathcal{X}$ with
constant $\alpha^{u}=\alpha+o(\|u\|^{2})$ as $u\rightarrow 0$.

##### Step 4 (Compatibility)

Finally, we verify that Assumption 2(iii) holds for the perturbed problem
corresponding to $\mathcal{D}^{u}$. Indeed, the previous steps show
$\gamma^{u}=\gamma+O(\|u\|)$ and $\alpha^{u}=\alpha+o(\|u\|^{2})$ as
$u\rightarrow 0$, so the compatibility inequality $\gamma\beta<\alpha$
corresponding to $\mathcal{\mathcal{D}}$ implies $\gamma^{u}\beta<\alpha^{u}$
for all sufficiently small $u$.

###  C.2 Lemmas used in the proof of Theorem 5.7

In this section, we collect the proofs of the two key lemmas used in the proof
of Theorem 5.7.

####  C.2.1 Proof of Lemma 5.11

Let $F\colon\mathcal{X}\times\mathbf{R}^{d}\rightarrow\mathbf{R}^{d}$ be the
map given by

$F(x,u)={\mathop{\operatorname{\mathbb{E}}}_{z\sim\mathcal{D}_{x}^{u}}}G(x,z)=\frac{1}{C_{x}^{u}}\,{\mathop{\operatorname{\mathbb{E}}}_{z\sim\mathcal{D}_{x}}}\big{[}G(x,z)\big{(}1+h\big{(}u^{\top}g_{x}(z)\big{)}\big{)}\big{]},$

where we recall
$C_{x}^{u}=1+{\mathop{\operatorname{\mathbb{E}}}_{z\sim\mathcal{D}_{x}}}h\big{(}u^{\top}g_{x}(z)\big{)}$.
Lemma C.4 directly implies that $F$ is $C^{1}$-smooth. Consider now the family
of smooth nonlinear equations

$F(x,u)=0$ (26)

parametrized by $u\in\mathcal{\mathbf{R}}^{d}$. Observe that
$F(x^{\star},0)=0$ since $x^{\star}\in\operatorname{int}\mathcal{X}$. More
generally, the equality (26) with $x\in\rm{int}\,\mathcal{X}$ holds precisely
when $x$ is equal to $x_{u}^{\star}$. We will apply the implicit function
theorem to $x_{u}^{\star}$ as a smooth implicit function of $u$. To this end,
we compute the the partial Jacobian

$\nabla_{x}F(x^{\star},0)=\nabla_{x}\Big{(}\mathop{\operatorname{\mathbb{E}}}_{z\sim\mathcal{D}_{x}}G(x,z)\Big{)}\Big{|}_{x=x^{\star}}=W.$

Lemma 4.7 therefore implies that $\nabla_{x}F(x^{\star},0)$ is invertible.
Consequently, the implicit function theorem yields open sets
$V\subset\operatorname{int}\mathcal{X}$ and $U\subset\mathbf{R}^{d}$ such that
$x^{\star}\in V$, $0\in U$, and for each $u\in U$ there exists a unique point
$x_{u}^{\star}\in V$ satisfying $F(x_{u}^{\star},u)=0$. Further, the map
$u\mapsto x_{u}^{\star}$ is $C^{1}$-smooth at $u=0$ with Jacobian
$-W^{-1}\nabla_{u}F(x^{\star},0)$. This yields the first-order Taylor
expansion

$x^{\star}_{u}=x^{\star}-W^{-1}\nabla_{u}F(x^{\star},0)u+o(\|u\|)\qquad\textrm{as
}u\to 0.$ (27)

To compute $\nabla_{u}F(x^{\star},0)$, we apply the quotient rule and the
dominated convergence theorem to obtain

$\nabla_{u}F(x,u)=\frac{1}{(C_{x}^{u})^{2}}\,{\mathop{\operatorname{\mathbb{E}}}_{z\sim\mathcal{D}_{x}}}\big{[}C_{x}^{u}h^{\prime}\big{(}u^{\top}g_{x}(z)\big{)}G(x,z)g_{x}(z)^{\top}-\big{(}1+h\big{(}u^{\top}g_{x}(z)\big{)}\big{)}G(x,z)(\nabla_{u}C_{x}^{u})^{\top}\big{]}$

for all $x\in\mathcal{X}$ and $u\in\mathbf{R}^{d}$. Then the equalities
$C_{x}^{0}=1$, $h^{\prime}(0)=1$, and

$\nabla_{u}C_{x}^{u}\big{|}_{u=0}={\mathop{\operatorname{\mathbb{E}}}_{z\sim\mathcal{D}_{x}}}g_{x}(z)=0$

imply

$\nabla_{u}F(x,0)={\mathop{\operatorname{\mathbb{E}}}_{z\sim\mathcal{D}_{x}}}\big{[}G(x,z)g_{x}(z)^{\top}\big{]}$

for all $x\in\mathcal{X}$. In particular,
$\nabla_{u}F(x^{\star},0)=\Sigma_{g,G}^{\top}$. Thus, the expansion (27)
becomes

$\displaystyle x^{\star}_{u}$
$\displaystyle=x^{\star}-W^{-1}\Sigma_{g,G}^{\top}u+o(\|u\|)\qquad\textrm{as
}u\to 0.$

Consequently, for any $u\in\mathbf{R}^{d}$, we have

$\sqrt{k}\big{(}x^{\star}_{u/\sqrt{k}}-x^{\star}\big{)}=-W^{-1}\Sigma_{g,G}^{\top}u+\sqrt{k}\cdot
o\bigg{(}\frac{\|u\|}{\sqrt{k}}\bigg{)}\rightarrow-W^{-1}\Sigma_{g,G}^{\top}u\qquad\text{as~{}}k\rightarrow\infty.$

The proof is complete.

####  C.2.2 Proof of Lemma 5.12

We start by deriving the expression

$dQ_{k,u}(z_{0},\dots,z_{k-1})=d\mathcal{D}_{\tilde{x}_{k-1}}^{u/\sqrt{k}}(z_{k-1})\cdots
d\mathcal{D}_{\tilde{x}_{0}}^{u/\sqrt{k}}(z_{0}).$ (28)

To see this, let $B\in\mathcal{B}(\mathcal{Z}^{k})$ be a Borel measurable set.
By repeatedly applying the tower rule, we get

$\displaystyle\int\mathbf{1}_{B}\,dQ_{k,u}$
$\displaystyle=\operatorname{\mathbb{E}}_{P_{k,{u/\sqrt{k}}}}[\mathbf{1}_{B}]$
$\displaystyle={\mathop{\operatorname{\mathbb{E}}}_{z_{0}\sim\mathcal{D}^{u/\sqrt{k}}_{\tilde{x}_{0}}}}\big{[}\operatorname{\mathbb{E}}[\mathbf{1}_{B}\mid
z_{0}]\big{]}$ $\displaystyle=\int\operatorname{\mathbb{E}}[\mathbf{1}_{B}\mid
z_{0}]\,d\mathcal{D}^{u/\sqrt{k}}_{\tilde{x}_{0}}(z_{0})$
$\displaystyle=\int{\mathop{\operatorname{\mathbb{E}}}_{z_{1}\sim\mathcal{D}^{u/\sqrt{k}}_{\tilde{x}_{1}}}}\big{[}\operatorname{\mathbb{E}}[\mathbf{1}_{B}\mid
z_{0},z_{1}]\big{]}\,d\mathcal{D}^{u/\sqrt{k}}_{\tilde{x}_{0}}(z_{0})$
$\displaystyle=\int\operatorname{\mathbb{E}}[\mathbf{1}_{B}\mid
z_{0},z_{1}]\,d\mathcal{D}^{u/\sqrt{k}}_{\tilde{x}_{1}}(z_{1})\,d\mathcal{D}^{u/\sqrt{k}}_{\tilde{x}_{0}}(z_{0})$
$\displaystyle\hskip 5.69046pt\vdots$
$\displaystyle=\int\mathbf{1}_{B}\,d\mathcal{D}^{u/\sqrt{k}}_{\tilde{x}_{k-1}}(z_{k-1})\cdots
d\mathcal{D}^{u/\sqrt{k}}_{\tilde{x}_{1}}(z_{1})\,d\mathcal{D}^{u/\sqrt{k}}_{\tilde{x}_{0}}(z_{0})$

using that $\tilde{x}_{j}$ is a deterministic function of the samples
$(z_{0},\ldots,z_{j-1})$.

Now let $\\{u_{k}\\}$ be any bounded sequence in $\mathbf{R}^{d}$ and observe
that (28) together with (13) yield

$\log\frac{dQ_{k,u_{k}}}{dQ_{k,0}}(z_{0},\dots,z_{k-1})=-\sum_{j=0}^{k-1}\log
C_{\tilde{x}_{j}}^{u_{k}/\sqrt{k}}+\sum_{j=0}^{k-1}{\log}\Big{(}1+h\big{(}u_{k}^{\top}g_{\tilde{x}_{j}}(z_{j})/\sqrt{k}\big{)}\Big{)}.$
(29)

Setting
$e_{x}(u):=\mathop{\operatorname{\mathbb{E}}}_{z\sim\mathcal{D}_{x}}h\big{(}u^{\top}g_{x}(z)\big{)}$
gives

$C_{x}^{u}=1+e_{x}(u),$

and Lemma C.2 shows $\sup_{x\in\mathcal{X}}|e_{x}(u)|=o(\|u\|^{2})$ as
$u\rightarrow 0$. Using the fact that $\log(1+t)=t+o(t)$ as $t\rightarrow 0$,
we derive

$\displaystyle\sum_{j=0}^{k-1}\log C_{\tilde{x}_{j}}^{u_{k}/\sqrt{k}}$
$\displaystyle=\sum_{j=0}^{k-1}\left(e_{\tilde{x}_{j}}\big{(}u_{k}/\sqrt{k}\big{)}+o\big{(}e_{\tilde{x}_{j}}\big{(}u_{k}/\sqrt{k}\big{)}\big{)}\right)=k\cdot
o\bigg{(}\frac{\|u_{k}\|^{2}}{k}\bigg{)}\to 0\qquad\text{as~{}}k\to\infty.$
(30)

We can therefore ignore the first term in (29).

Next, note that the boundedness of $\\{u_{k}\\}$ and $g$ (Lemma C.1) implies
that $\sup_{j,k}|u_{k}^{\top}g_{\tilde{x}_{j}}(z_{j})|<\infty$. Therefore
$h\big{(}{u_{k}^{\top}g_{\tilde{x}_{j}}(z_{j})}/{\sqrt{k}}\big{)}={u_{k}^{\top}g_{\tilde{x}_{j}}(z_{j})}/{\sqrt{k}}\in[-1/2,1/2]$
for all sufficiently large $k$. Then, using the second-order mean value
theorem for $\log(1+t)$ around zero, we have

$\sum_{j=0}^{k-1}{\log}\bigg{(}1+\frac{u_{k}^{\top}g_{\tilde{x}_{j}}(z_{j})}{\sqrt{k}}\bigg{)}=u_{k}^{\top}\bigg{(}\frac{1}{\sqrt{k}}\sum_{j=0}^{k-1}g_{\tilde{x}_{j}}(z_{j})\bigg{)}-\frac{1}{2}u_{k}^{\top}\bigg{(}\frac{1}{k}\sum_{j=0}^{k-1}g_{\tilde{x}_{j}}(z_{j})g_{\tilde{x}_{j}}(z_{j})^{\top}\bigg{)}u_{k}+\sum_{j=0}^{k-1}\tau_{j}\frac{\big{(}u_{k}^{\top}g_{\tilde{x}_{j}}(z_{j})\big{)}^{3}}{3k^{3/2}}$
(31)

where the coefficients $\tau_{j}$ satisfy $|\tau_{j}|\leq 8$ for all $j$. The
first term on the right-hand side matches the corresponding term in (17) upon
setting $Z_{k}:=k^{-1/2}\sum_{j=0}^{k-1}{g_{\tilde{x}_{j}}(z_{j})}$. Before
establishing the distributional convergence of $Z_{k}$, let us study the last
two terms of (31).

We start with the middle term on the right-hand side of (31); we claim

$\frac{1}{k}\sum_{j=0}^{k-1}g_{\tilde{x}_{j}}(z_{j})g_{\tilde{x}_{j}}(z_{j})^{\top}=\Sigma_{g}+o_{Q_{k,0}}(1)$
(32)

as a consequence of Theorem E.1. Indeed, let $\mathbb{F}=(\mathcal{F}_{k})$ be
the filtration given by $\mathcal{F}_{k}=\sigma(z_{0},\ldots,z_{k-1})$ where
$(z_{0},\ldots,z_{k-1})\sim Q_{k,0}$, and for each $j\geq 0$, set

$\displaystyle X_{j+1}$
$\displaystyle=g_{\tilde{x}_{j}}(z_{j})g_{\tilde{x}_{j}}(z_{j})^{\top}-{\operatorname{\mathbb{E}}}\big{[}g_{\tilde{x}_{j}}(z_{j})g_{\tilde{x}_{j}}(z_{j})^{\top}\mid\mathcal{F}_{j}\big{]}$
$\displaystyle=g_{\tilde{x}_{j}}(z_{j})g_{\tilde{x}_{j}}(z_{j})^{\top}-{\mathop{\operatorname{\mathbb{E}}}_{z_{j}\sim\mathcal{D}_{\tilde{x}_{j}}}}\big{[}g_{\tilde{x}_{j}}(z_{j})g_{\tilde{x}_{j}}(z_{j})^{\top}\big{]},$

thereby defining a martingale difference sequence $X$ adapted to $\mathbb{F}$;
note that Lemma C.1 implies $\sup_{j}\|X_{j}\|_{\text{F}}^{2}<\infty$, so $X$
is square-integrable and

$\sum_{j=1}^{\infty}j^{-2}{\operatorname{\mathbb{E}}}\|X_{j}\|_{\text{F}}^{2}\leq\tfrac{\pi^{2}}{6}\sup_{j}\|X_{j}\|_{\text{F}}^{2}<\infty.$

Therefore Theorem E.1 implies

$\frac{1}{k}\sum_{j=0}^{k-1}\Big{(}g_{\tilde{x}_{j}}(z_{j})g_{\tilde{x}_{j}}(z_{j})^{\top}-{\mathop{\operatorname{\mathbb{E}}}_{z_{j}\sim\mathcal{D}_{\tilde{x}_{j}}}}\big{[}g_{\tilde{x}_{j}}(z_{j})g_{\tilde{x}_{j}}(z_{j})^{\top}\big{]}\Big{)}=\frac{1}{k}\sum_{j=1}^{k}X_{j}\xlongrightarrow{\textrm{a.s.}}0.$
(33)

On the other hand, we have $\tilde{x}_{j}\rightarrow x^{\star}$ almost surely
by Assumption 5.3, so

${\mathop{\operatorname{\mathbb{E}}}_{z_{j}\sim\mathcal{D}_{\tilde{x}_{j}}}}\big{[}g_{\tilde{x}_{j}}(z_{j})g_{\tilde{x}_{j}}(z_{j})^{\top}\big{]}\xlongrightarrow{\textrm{a.s.}}\underset{z\sim\mathcal{D}_{x^{\star}}}{\mathop{\operatorname{\mathbb{E}}}}\big{[}g_{x^{\star}}(z)g_{x^{\star}}(z)^{\top}\big{]}=\Sigma_{g}$

by Lemma E.3. Thus, $\Sigma_{g}$ is a Cesàro mean of covariances:

$\frac{1}{k}\sum_{j=0}^{k-1}{\mathop{\operatorname{\mathbb{E}}}_{z_{j}\sim\mathcal{D}_{\tilde{x}_{j}}}}\big{[}g_{\tilde{x}_{j}}(z_{j})g_{\tilde{x}_{j}}(z_{j})^{\top}\big{]}\xlongrightarrow{\textrm{a.s.}}\Sigma_{g}.$
(34)

Combining (33) and (34), we deduce

$\frac{1}{k}\sum_{j=0}^{k-1}g_{\tilde{x}_{j}}(z_{j})g_{\tilde{x}_{j}}(z_{j})^{\top}\xlongrightarrow{\textrm{a.s.}}\Sigma_{g},$

thereby establishing (32). As for the last term in (31), we have

$\displaystyle\left|\sum_{j=0}^{k-1}\tau_{j}\frac{(u_{k}^{\top}g_{\tilde{x}_{j}}(z_{j}))^{3}}{3k^{3/2}}\right|\leq\frac{8\big{(}\sup_{i,j}|u_{i}^{\top}g_{\tilde{x}_{j}}(z_{j})|^{3}\big{)}}{3\sqrt{k}}\to
0\qquad\text{as~{}}k\to\infty.$ (35)

Hence, it only remains to demonstrate that
$Z_{k}\mathrel{\stackunder[.5pt]{\xrightarrow{\makebox[22.0556pt]{$\scriptstyle
d$}}}{\scriptstyle Q_{k,0}\,}}\mathsf{N}(0,\Sigma_{g})$. To see this, we apply
the martingale central limit theorem (Theorem E.2). Set $M_{0}=0$ and
$M_{k}=\sum_{j=0}^{k-1}g_{\tilde{x}_{j}}(z_{j})$ for all $k\geq 1$; then $M$
is a square-integrable martingale in $\mathbf{R}^{d}$ adapted to the
filtration $\mathbb{F}$. Indeed, $M_{k}$ is clearly bounded (Lemma C.1),
$\mathcal{F}_{k}$-measurable, and

$\operatorname{\mathbb{E}}\big{[}M_{k+1}\mid\mathcal{F}_{k}\big{]}=M_{k}+{\mathop{\operatorname{\mathbb{E}}}_{z_{k}\sim\mathcal{D}_{\tilde{x}_{k}}}}\big{[}g_{\tilde{x}_{k}}(z_{k})\big{]}=M_{k}$

by the unbiasedness assumption of Definition 5.5, which also facilitates the
inductive calculation of the predictable square variation $\langle M\rangle$
of $M$:

$\langle
M\rangle_{k}=\sum_{j=0}^{k-1}{\mathop{\operatorname{\mathbb{E}}}_{z_{j}\sim\mathcal{D}_{\tilde{x}_{j}}}}\big{[}g_{\tilde{x}_{j}}(z_{j})g_{\tilde{x}_{j}}(z_{j})^{\top}\big{]}.$

Thus, by (34), we have

$k^{-1}\langle
M\rangle_{k}=\frac{1}{k}\sum_{j=0}^{k-1}{\mathop{\operatorname{\mathbb{E}}}_{z_{j}\sim\mathcal{D}_{\tilde{x}_{j}}}}\big{[}g_{\tilde{x}_{j}}(z_{j})g_{\tilde{x}_{j}}(z_{j})^{\top}\big{]}\xlongrightarrow{\textrm{a.s.}}\Sigma_{g}.$

The assumptions of Theorem E.2 are therefore fulfilled with $a_{k}=k$ (note
that Lindeberg’s condition holds trivially since
$\sup_{k}\|M_{k}-M_{k-1}\|<\infty$ by Lemma C.1). Hence

$Z_{k}=k^{-1/2}M_{k}\mathrel{\stackunder[.5pt]{\xrightarrow{\makebox[22.0556pt]{$\scriptstyle
d$}}}{\scriptstyle Q_{k,0}\,}}\mathsf{N}(0,\Sigma_{g}).$

This completes the proof.

##  Appendix D Proof of Lemma C.4

First, we recall that the quantities
$M^{\prime}_{g}:=\sup_{x\in\mathcal{X},\,z\in\mathcal{Z}}\|\nabla_{x}g_{x}(z)\|_{\text{op}}$
and $M_{g}:=\sup_{x\in\mathcal{X},\,z\in\mathcal{Z}}\|g_{x}(z)\|$ are finite
by Lemma C.1. Exactly the same argument shows that $M^{\prime}_{T}$ and
$M_{T}$ are finite.

Next, we turn to establishing that the map
$H\colon\mathcal{X}\times\mathcal{W}\rightarrow\mathbf{R}^{n}$ given by

$H(x,u)=\frac{1}{C_{x}^{u}}\,{\mathop{\operatorname{\mathbb{E}}}_{z\sim\mathcal{D}_{x}}}\big{[}T(x,z)\big{(}1+h\big{(}u^{\top}g_{x}(z)\big{)}\big{)}\big{]}$

is smooth with Lipschitz Jacobian on the compact set
$\mathcal{K}:=\mathcal{X}\times\mathcal{W}$. By Lemma E.5, it is enough to
show that $(x,u)\mapsto C_{x}^{u}$ and

$\hat{H}(x,u):={\mathop{\operatorname{\mathbb{E}}}_{z\sim\mathcal{D}_{x}}}\big{[}T(x,z)\big{(}1+h\big{(}u^{\top}g_{x}(z)\big{)}\big{)}\big{]}$

are smooth with Lipschitz Jacobians on $\mathcal{K}$; in turn, it suffices to
establish this fact for $\hat{H}$ since we can then take $T\equiv 1$ to derive
the result for $C_{x}^{u}$.

We reason this via the chain rule. Namely, consider the map
$\bar{H}\colon\mathcal{X}\times\mathcal{X}\times\mathcal{W}\rightarrow{\mathbf{R}}^{n}$
given by

$\bar{H}(x,y,u)={\mathop{\operatorname{\mathbb{E}}}_{z\sim\mathcal{D}_{x}}}\big{[}T(y,z)\big{(}1+h\big{(}u^{\top}g_{y}(z)\big{)}\big{)}\big{]}.$

Clearly $\hat{H}=\bar{H}\circ J$ with $J(x,u):=(x,x,u)$ and therefore the
chain rule implies $\nabla\hat{H}(x,u)=\nabla\bar{H}(x,x,u)\nabla J(x,u)$
provided $\bar{H}$ is smooth. Thus, it suffices to show that $\bar{H}$ is
smooth with Lipschitz Jacobian. To this end, we demonstrate that the three
partial derivatives of $\bar{H}$ are all Lipschitz with constants depending on
$T$ only through $\beta_{T},\bar{\Lambda}_{T},\beta^{\prime}_{T}$, $M_{T}$,
and $M^{\prime}_{T}$.

We begin with the partial derivative of $\bar{H}$ with respect to $x$.
Consider the function
$\phi\colon\mathcal{K}\times\mathcal{Z}\rightarrow\mathbf{R}^{n}$ given by

$\phi(y,u,z)=T(y,z)\big{(}1+h\big{(}u^{\top}g_{y}(z)\big{)}\big{)}.$

Let us verify that $\phi$ is a test function to which the second item of
Assumption 5 applies. Clearly $\phi$ is measurable and bounded with
$\sup\|\phi\|\leq 2M_{T}$. Further, for each $z\in\mathcal{Z}$, it follows
readily that the section $\phi(\cdot,z)$ is Lipschitz on $\mathcal{K}$ with
constant

$L_{\phi}:=2M^{\prime}_{T}+M_{T}L_{h}\big{(}\operatorname{diam}(\mathcal{W})M^{\prime}_{g}+M_{g}\big{)},$

where $L_{h}:=\sup|h^{\prime}|$. Thus, the second item of Assumption 5 implies
that the map

$x\mapsto\mathop{\operatorname{\mathbb{E}}}_{z\sim\mathcal{D}_{x}}{\phi(y,u,z)}=\bar{H}(x,y,u)$

is smooth on $\mathcal{X}$ for each $(y,u)\in\mathcal{K}$, and that the map

$(x,y,u)\mapsto\nabla_{x}\bar{H}(x,y,u)$

is Lipschitz on $\mathcal{X}\times\mathcal{X}\times\mathcal{W}$ with constant
$\vartheta(L_{\phi}+2M_{T})$, which depends on $T$ only through $M_{T}$ and
$M^{\prime}_{T}$.

Next, we consider the partial derivative of $\bar{H}$ with respect to $y$.
Given $(x,u)\in\mathcal{X}\times\mathcal{W}$, the dominated convergence
theorem ensures that $\bar{H}(x,y,u)$ is smooth in $y$ with

$\nabla_{y}\bar{H}(x,y,u)={\mathop{\operatorname{\mathbb{E}}}_{z\sim\mathcal{D}_{x}}}\nabla_{y}\phi(y,u,z)$
(36)

provided $\|\nabla_{y}\phi(y,u,z)\|_{\text{op}}$ is dominated by a
$\mathcal{D}_{x}$-integrable random variable independent of $y$. Using the
product rule, we have

$\nabla_{y}\phi(y,u,z)=\big{(}\nabla_{y}T(y,z)\big{)}\big{(}1+h\big{(}u^{\top}g_{y}(z)\big{)}\big{)}+h^{\prime}\big{(}u^{\top}g_{y}(z)\big{)}\big{(}T(y,z)u^{\top}\big{)}\nabla_{y}g_{y}(z)$
(37)

and hence

$\displaystyle\|\nabla_{y}\phi(y,u,z)\|_{\text{op}}$ $\displaystyle\leq
2\|\nabla_{y}T(y,z)\|_{\text{op}}+(\sup|h^{\prime}|)\|T(y,z)\|\|u\|\|\nabla_{y}g_{y}(z)\|_{\text{op}}$
$\displaystyle\leq
2M^{\prime}_{T}+\operatorname{diam}(\mathcal{W})L_{h}M_{T}M^{\prime}_{g},$

so $\nabla_{y}\phi$ is in fact uniformly bounded. Therefore $\bar{H}(x,y,u)$
is smooth in $y$ and (36) holds. Moreover, it follows from (37) that the map

$(x,y,u)\mapsto\nabla_{y}\bar{H}(x,y,u)$

is Lipschitz on $\mathcal{X}\times\mathcal{X}\times\mathcal{W}$; we will
verify this by computing Lipschitz constants separately in $x$, $y$, and $u$.
To begin, note that it follows from (37) that $z\mapsto\nabla_{y}\phi(y,u,z)$
is Lipschitz on $\mathcal{Z}$ with constant

$a:=2\beta^{\prime}_{T}+\operatorname{diam}(\mathcal{W})M^{\prime}_{T}L_{h}\beta_{g}+\operatorname{diam}(\mathcal{W})L_{h}(\beta_{T}M^{\prime}_{g}+M_{T}\beta^{\prime}_{g})+\operatorname{diam}(\mathcal{W})^{2}M_{T}M^{\prime}_{g}L_{h^{\prime}}\beta_{g}.$

Hence (36) and Assumption 1 imply that $x\mapsto\nabla_{y}\bar{H}(x,y,u)$ is
Lipschitz on $\mathcal{X}$ with constant $\gamma a$, which depends on $T$ only
through $\beta_{T},\beta^{\prime}_{T},M_{T}$, and $M^{\prime}_{T}$. Likewise,
it follows from (37) that $y\mapsto\nabla_{y}\phi(y,u,z)$ is Lipschitz on
$\mathcal{X}$ with constant

$2\Lambda_{T}(z)+\operatorname{diam}(\mathcal{W})M^{\prime}_{T}L_{h}M^{\prime}_{g}+\operatorname{diam}(\mathcal{W})L_{h}\big{(}M_{T}\Lambda_{g}(z)+M^{\prime}_{T}M^{\prime}_{g}\big{)}+\operatorname{diam}(\mathcal{W})^{2}M_{T}L_{h^{\prime}}(M^{\prime}_{T})^{2},$

where $L_{h^{\prime}}:=\sup|h^{\prime\prime}|$. Hence (36) implies that
$y\mapsto\nabla_{y}\bar{H}(x,y,u)$ is Lipschitz on $\mathcal{X}$ with constant

$2\bar{\Lambda}_{T}+\operatorname{diam}(\mathcal{W})L_{h}\big{(}M_{T}\bar{\Lambda}_{g}+2M^{\prime}_{T}M^{\prime}_{g}\big{)}+\operatorname{diam}(\mathcal{W})^{2}M_{T}L_{h^{\prime}}(M^{\prime}_{T})^{2}.$

Similarly, it follows from (37) that $u\mapsto\nabla_{y}\phi(y,u,z)$ is
Lipschitz on $\mathcal{W}$ with constant

$M^{\prime}_{T}M_{g}L_{h}+M_{T}M^{\prime}_{g}\big{(}L_{h}+\operatorname{diam}(\mathcal{W})M_{g}L_{h^{\prime}}\big{)},$

so (36) implies that $u\mapsto\nabla_{y}\bar{H}(x,y,u)$ is Lipschitz on
$\mathcal{W}$ with the same constant. We conclude therefore that the map
$(x,y,u)\mapsto\nabla_{y}\bar{H}(x,y,u)$ is Lipschitz on
$\mathcal{X}\times\mathcal{X}\times\mathcal{W}$ with constant depending on $T$
only through $\beta_{T},\bar{\Lambda}_{T},\beta^{\prime}_{T}$, $M_{T}$, and
$M^{\prime}_{T}$.

Finally, we consider the partial derivative of $\bar{H}$ with respect to $u$.
Given $(x,y)\in\mathcal{X}\times\mathcal{X}$, the dominated convergence
theorem ensures that $\bar{H}(x,y,u)$ is smooth in $u$ with

$\nabla_{u}\bar{H}(x,y,u)={\mathop{\operatorname{\mathbb{E}}}_{z\sim\mathcal{D}_{x}}}\nabla_{u}\phi(y,u,z)$
(38)

provided $\|\nabla_{u}\phi(y,u,z)\|_{\text{op}}$ is dominated by a
$\mathcal{D}_{x}$-integrable random variable independent of $u$. In this case,
we have

$\nabla_{u}\phi(y,u,z)=h^{\prime}\big{(}u^{\top}g_{y}(z)\big{)}T(y,z)g_{y}(z)^{\top}$
(39)

and hence

$\displaystyle\|\nabla_{u}\phi(y,u,z)\|_{\text{op}}$
$\displaystyle\leq(\sup|h^{\prime}|)\|T(y,z)\|\|g_{y}(z)\|\leq
L_{h}M_{T}M_{g}.$

Therefore $\bar{H}(x,y,u)$ is smooth in $u$ and (38) holds. Moreover, it
follows from (39) that the map

$(x,y,u)\mapsto\nabla_{u}\bar{H}(x,y,u)$

is Lipschitz on $\mathcal{X}\times\mathcal{X}\times\mathcal{W}$; as before, we
will verify this by computing Lipschitz constants separately in $x$, $y$, and
$u$. First, note that it follows from (39) that
$z\mapsto\nabla_{u}\phi(y,u,z)$ is Lipschitz on $\mathcal{Z}$ with constant

$b:=L_{h}(\beta_{T}M_{g}+M_{T}\beta_{g})+\operatorname{diam}(\mathcal{W})M_{T}M_{g}L_{h^{\prime}}\beta_{g}.$

Hence (38) and Assumption 1 imply that $x\mapsto\nabla_{u}\bar{H}(x,y,u)$ is
Lipschitz on $\mathcal{X}$ with constant $\gamma b$, which depends on $T$ only
through $\beta_{T}$ and $M_{T}$. Likewise, it follows from (39) that
$y\mapsto\nabla_{u}\phi(y,u,z)$ is Lipschitz on $\mathcal{X}$ with constant

$L_{h}(M^{\prime}_{T}M_{g}+M_{T}M^{\prime}_{g})+\operatorname{diam}(\mathcal{W})M_{T}M_{g}L_{h^{\prime}}M^{\prime}_{g},$

hence so is $y\mapsto\nabla_{u}\bar{H}(x,y,u)$ by (38). Similarly, it follows
from (39) that $u\mapsto\nabla_{u}\phi(y,u,z)$ is
$L_{h^{\prime}}M_{T}M_{g}^{2}$-Lipschitz on $\mathcal{W}$, hence so is
$u\mapsto\nabla_{u}\bar{H}(x,y,u)$ by (38). We conclude therefore that the map
$(x,y,u)\mapsto\nabla_{u}\bar{H}(x,y,u)$ is Lipschitz on
$\mathcal{X}\times\mathcal{X}\times\mathcal{W}$ with constant depending on $T$
only through $\beta_{T},M_{T}$, and $M^{\prime}_{T}$.

The preceding reveals that $\bar{H}$ and hence $\hat{H}=\bar{H}\circ J$ are
smooth, with Lipschitz Jacobians with constants depending on $T$ only through
$\beta_{T},\bar{\Lambda}_{T},\beta^{\prime}_{T}$, $M_{T}$, and
$M^{\prime}_{T}$. Taking $T\equiv 1$, we conclude that $(x,u)\mapsto
C_{x}^{u}$ is smooth, with Lipschitz Jacobian with constant independent of
$T$. Upon observing in the same way as above that $\bar{H}$ and hence
$\hat{H}$ are Lipschitz with constants depending on $T$ only on
$\beta_{T},M_{T}$, and $M^{\prime}_{T}$, it follows from Lemma E.5 and its
proof that $H$ is smooth, with Lipschitz Jacobian with constant depending on
$T$ only through $\beta_{T},\bar{\Lambda}_{T},\beta^{\prime}_{T}$, $M_{T}$,
and $M^{\prime}_{T}$. This completes the proof.

##  Appendix E Supplementary results

In this section, we provide some supplementary results required in our proofs.
First, we record fairly general versions of the Strong Law of Large Numbers
and the Central Limit Theorem.

######  Theorem E.1 (Martingale Strong Law of Large Numbers [11, Theorem
VII.9.3]).

Let $X$ be a square-integrable martingale difference sequence in
$\mathbf{R}^{n}$, and let $(a_{k})$ be a deterministic sequence of nonzero
real numbers increasing to $\infty$. If

$\sum_{j=1}^{\infty}a_{j}^{-2}{\operatorname{\mathbb{E}}}\|X_{j}\|^{2}<\infty,$

then

$a_{k}^{-1}\sum_{j=1}^{k}X_{j}\xlongrightarrow{\textrm{a.s.}}0.$

######  Theorem E.2 (Martingale Central Limit Theorem [9, Corollary 2.1.10]).

Let $M$ be a square-integrable martingale in $\mathbf{R}^{n}$ adapted to a
filtration $(\mathcal{F}_{k})$ and let $\langle M\rangle$ denote the
predictable square variation of $M$, given by

$\langle M\rangle_{k}-\langle
M\rangle_{k-1}={\operatorname{\mathbb{E}}}\big{[}M_{k}M_{k}^{\top}\mid\mathcal{F}_{k-1}\big{]}-M_{k-1}M_{k-1}^{\top}$

with $\langle M\rangle_{0}=0$. Further, let $(a_{k})$ be a deterministic
sequence of nonzero real numbers increasing to $\infty$, and suppose that the
following two assumptions hold.

  1. 1.

(Square variation) There exists a deterministic positive semidefinite matrix
$\Sigma$ such that

$a_{k}^{-1}\langle M\rangle_{k}\xlongrightarrow{p}\Sigma.$

  2. 2.

(Lindeberg’s condition) For all $\varepsilon>0$,

$a_{k}^{-1}\sum_{j=1}^{k}{\operatorname{\mathbb{E}}}\big{[}\|M_{j}-M_{j-1}\|^{2}\mathbf{1}\\{\|M_{j}-M_{j-1}\|\geq\varepsilon
a_{k}^{1/2}\\}\mid\mathcal{F}_{j-1}\big{]}\xlongrightarrow{p}0.$

Then

$a_{k}^{-1}M_{k}\xlongrightarrow{\textrm{a.s.}}0\quad\textrm{and}\quad
a_{k}^{-1/2}M_{k}\xlongrightarrow{\textrm{d}}\mathsf{N}(0,\Sigma).$

The following lemma is used multiple times in our arguments to control
covariance matrices.

######  Lemma E.3 (Asymptotic covariance).

Let $x_{t}\in\mathcal{X}$ be a sequence in some set
$\mathcal{X}\subset\mathbf{R}^{d}$ that converges to some point
$x^{\star}\in\mathcal{X}$. Let $\mu_{t}\in P_{1}(\mathcal{Z})$ be a sequence
of probability measures on a metric space $\mathcal{Z}$ converging to some
measure $\mu^{\star}\in P_{1}(\mathcal{Z})$ in the Wasserstein-1 distance. Let
$g\colon\mathcal{X}\times\mathcal{Z}\rightarrow{\mathbf{R}}^{n}$ be a map
satisfying the following two conditions.

  1. 1.

For every $\delta>0$, there exists a constant $N_{\delta}$ such that

$\displaystyle\mathop{\operatorname{\mathbb{E}}}_{z\sim\mu^{\star}}\big{[}\|g(x^{\star},z)\|^{2}\mathbf{1}\\{\|g(x^{\star},z)\|^{2}\geq
N_{\delta}\\}\big{]}\leq\delta,$ $\displaystyle\limsup_{t\to\infty}$
$\displaystyle\mathop{\operatorname{\mathbb{E}}}_{z\sim\mu_{t}}\big{[}\|g(x^{\star},z)\|^{2}\mathbf{1}\\{\|g(x^{\star},z)\|^{2}\geq
N_{\delta}\\}\big{]}\leq\delta.$

  2. 2.

There exist a neighborhood $\mathcal{V}$ of $x^{\star}$, a measurable function
$L\colon\mathcal{Z}\rightarrow[0,\infty)$, and constants $\bar{L},\beta>0$
such that for every $z\in\mathcal{Z}$ and $x\in\mathcal{V}$, the section
$g(\cdot,z)$ is $L(z)$-Lipschitz on $\mathcal{V}$ with
$\limsup_{t\to\infty}{\operatorname{\mathbb{E}}}_{z\sim\mu_{t}}\big{[}L(z)^{2}\big{]}\leq\bar{L}^{2}$,
and the section $g(x^{\star},\cdot)$ is $\beta$-Lipschitz on $\mathcal{Z}$.

Then

$\mathop{\operatorname{\mathbb{E}}}_{z\sim\mu_{t}}\big{[}g(x_{t},z)g(x_{t},z)^{\top}\big{]}\xlongrightarrow{t\to\infty}\mathop{\operatorname{\mathbb{E}}}_{z\sim\mu^{\star}}\big{[}g(x^{\star},z)g(x^{\star},z)^{\top}\big{]}.$

###### Proof.

For notational convenience, we set $g_{x}(z)=g(x,z)$ and

$\Sigma={\mathop{\operatorname{\mathbb{E}}}_{z\sim\mu^{\star}}}\big{[}g_{x^{\star}}(z)g_{x^{\star}}(z)^{\top}\big{]}.$

Since $x_{t}\rightarrow x^{\star}$, it follows that all but finitely many
$x_{t}$ lie in $\mathcal{V}$. Suppose $x_{t}\in\mathcal{V}$. Condition 1 of
the lemma implies a uniform upper bound on the second moment of
$g_{x^{\star}}$:

$\displaystyle\begin{split}\underset{z\sim\mu_{t}}{\operatorname{\mathbb{E}}}\big{[}\|g_{x^{\star}}(z)\|^{2}\big{]}&=\underset{z\sim\mu_{t}}{\operatorname{\mathbb{E}}}\big{[}\|g_{x^{\star}}(z)\|^{2}\mathbf{1}\\{\|g_{x^{\star}}(z)\|<N_{1}\\}\big{]}+\underset{z\sim\mu_{t}}{\operatorname{\mathbb{E}}}\big{[}\|g_{x^{\star}}(z)\|^{2}\mathbf{1}\\{\|g_{x^{\star}}(z)\|\geq
N_{1}\\}\big{]},\end{split}$

and therefore

$\limsup_{t\to\infty}\underset{z\sim\mu_{t}}{\operatorname{\mathbb{E}}}\big{[}\|g_{x^{\star}}(z)\|^{2}\big{]}\leq
N_{1}^{2}+1.$ (40)

Now observe that we have the decomposition

$\displaystyle\underset{z\sim\mu_{t}}{\operatorname{\mathbb{E}}}\big{[}g_{x_{t}}(z)g_{x_{t}}(z)^{\top}\big{]}$
$\displaystyle=\underset{z\sim\mu_{t}}{\operatorname{\mathbb{E}}}\big{[}g_{x^{\star}}(z)g_{x^{\star}}(z)^{\top}\big{]}+\underset{z\sim\mu_{t}}{\operatorname{\mathbb{E}}}\big{[}g_{x^{\star}}(z)\big{(}g_{x_{t}}(z)-g_{x^{\star}}(z)\big{)}^{\top}\big{]}$
$\displaystyle~{}~{}~{}+\underset{z\sim\mu_{t}}{\operatorname{\mathbb{E}}}\big{[}\big{(}g_{x_{t}}(z)-g_{x^{\star}}(z)\big{)}g_{x_{t}}(z)^{\top}\big{]}.$

The last two summands may be bounded as follows:

$\displaystyle\left\|\underset{z\sim\mu_{t}}{\operatorname{\mathbb{E}}}\big{[}g_{x^{\star}}(z)\big{(}g_{x_{t}}(z)-g_{x^{\star}}(z)\big{)}^{\top}\big{]}\right\|_{\text{op}}$
$\displaystyle\leq\underset{z\sim\mu_{t}}{\operatorname{\mathbb{E}}}\big{[}\|g_{x^{\star}}(z)\|\cdot\|g_{x_{t}}(z)-g_{x^{\star}}(z)\|\big{]}$
$\displaystyle\leq\|x_{t}-x^{\star}\|\sqrt{\underset{z\sim\mu_{t}}{\operatorname{\mathbb{E}}}\big{[}\|g_{x^{\star}}(z)\|^{2}\big{]}\underset{z\sim\mu_{t}}{\operatorname{\mathbb{E}}}\big{[}L(z)^{2}\big{]}}\to
0$

and

$\displaystyle\left\|\underset{z\sim\mu_{t}}{\operatorname{\mathbb{E}}}\big{[}\big{(}g_{x_{t}}(z)-g_{x^{\star}}(z)\big{)}g_{x_{t}}(z)^{\top}\big{]}\right\|_{\text{op}}$
$\displaystyle\leq\underset{z\sim\mu_{t}}{\operatorname{\mathbb{E}}}\big{[}\|g_{x_{t}}(z)-g_{x^{\star}}(z)\|\cdot\|g_{x_{t}}(z)\|\big{]}$
$\displaystyle\leq\|x_{t}-x^{\star}\|\sqrt{\underset{z\sim\mu_{t}}{\operatorname{\mathbb{E}}}\big{[}L(z)^{2}\big{]}\underset{z\sim\mu_{t}}{\operatorname{\mathbb{E}}}\big{[}\|g_{x_{t}}(z)\|^{2}\big{]}}$