# Stochastic Optimization Schemes for Performative Prediction with Nonconvex
Loss

Qiang Li, Hoi-To Wai Q. Li and H.-T. Wai are with the Department of Systems
Engineering and Engineering Management, The Chinese University of Hong Kong,
Hong Kong SAR of China. Emails: {liqiang<EMAIL_ADDRESS>
###### Abstract

This paper studies a risk minimization problem with decision dependent data
distribution. The problem pertains to the performative prediction setting
where a trained model can affect the outcome that the model estimates. Such
dependency creates a feedback loop that influences the stability of
optimization algorithms such as stochastic gradient descent (SGD). We present
the first study on performative prediction with smooth but possibly non-convex
loss. We analyze a greedy deployment scheme with SGD (SGD-GD). Note that in
the literature, SGD-GD is often studied with strongly convex loss. We first
propose the definition of stationary performative stable (SPS) solutions
through relaxing the popular performative stable condition. We then prove that
SGD-GD converges to a biased SPS solution in expectation. We consider two
conditions of sensitivity on the distribution shifts: (i) the sensitivity is
characterized by Wasserstein-1 distance and the loss is Lipschitz w.r.t. data
samples, or (ii) the sensitivity is characterized by $\chi^{2}$-divergence and
the loss is bounded. In both conditions, the bias levels are proportional to
stochastic gradient’s variance and sensitivity level. Our analysis is extended
to a lazy deployment scheme where models are deployed once per several SGD
updates, and we show that it converges to a bias-free SPS solution. Numerical
experiments corroborate our theories.

##  1 Introduction

When trained models are deployed in social contexts, the outcomes these models
aim to predict can be influenced by the models themselves. Taking email spam
detection as an example. On one hand, email service providers design filters
to protect their users by identifying spam emails. On the other hand, spammers
aim to circumvent these filters to distribute malware and advertisements. Each
time a new classifier is deployed, spammers who are interspersed within the
general population may alter the characteristics of their messages to evade
detection. The above example pertains to the strategic classification problem
(Dalvi et al., 2004; Cai et al., 2015; Hardt et al., 2016; Björkegren et al.,
2020) and can be modelled by dataset shifts (Quiñonero-Candela et al., 2022).

The scenarios described can be captured by the recently proposed performative
prediction problem, which called the above dataset shift phenomena as the
‘performative’ effect. Perdomo et al. (2020) proposed to study the risk
minimization problem with a _decision-dependent_ data distribution:

$\displaystyle\textstyle\min_{{\bm{\theta}}\in\mathbb{R}^{d}}~{}V({\bm{\theta}})\mathrel{\mathop{:}}=\mathbb{E}_{Z\sim{\cal
D}({\bm{\theta}})}[\ell({\bm{\theta}};Z)],$ (1)

where $\ell({\bm{\theta}};z)$ is a loss function that is continuously
differentiable with respect to (w.r.t.) ${\bm{\theta}}$ for any given data
sample $z\in{\sf Z}$, and ${\sf Z}\subseteq\mathbb{R}^{p}$ is the sample
space. The dependence on ${\bm{\theta}}$ in ${\cal D}({\bm{\theta}})$
explicitly captures the distribution shift effect of prediction models on data
samples. The objective function $V({\bm{\theta}})$ is also known as the
performative risk.

The decision variable ${\bm{\theta}}$ in (1) affects simultaneously the
distribution and the loss function. As such, optimizing $V({\bm{\theta}})$
directly is often difficult. Mendler-Dünner et al. (2020) considered the
following stochastic gradient (SGD) recursion: for any $t\geq 0$ and let
$\gamma_{t+1}>0$ be a stepsize,

$\displaystyle{\bm{\theta}}_{t+1}={\bm{\theta}}_{t}-\gamma_{t+1}{\nabla}\ell({\bm{\theta}}_{t};Z_{t+1}),\text{
where }Z_{t+1}\sim{\cal D}({\bm{\theta}}_{t}).$ (2)

The above is known as the _greedy deployment_ scheme with SGD (SGD-GD), where
the learner deploys the current trained model ${\bm{\theta}}_{t}$ before
drawing samples from ${\cal D}({\bm{\theta}}_{t})$. The SGD-GD scheme
describes a training procedure when the learner is _unaware of the
performative phenomena_ with the data distribution ${\cal D}(\cdot)$, which is
plausible in many applications. Relevant studies to (2) include lazy
deployment where the learner deploys a new model only once every few
iterations, or repeated risk minimization; see (Mendler-Dünner et al., 2020;
Perdomo et al., 2020; Zrnic et al., 2021).

Existing convergence analysis of (2) are limited to the case when
$\ell({\bm{\theta}};z)$ is _strongly convex_ 111Note that $V({\bm{\theta}})$
is still non-convex. w.r.t. ${\bm{\theta}}$. Perdomo et al. (2020) introduced
the concept of _performative stable_ (PS) solution as the unique minimizer of
(1) with fixed distribution. The PS solution, while being different from an
optimal or stationary solution to (1), is shown to be the unique limit point
of the recursion (2) provided that the sensitivity of the distribution map
${\cal D}(\cdot)$, measured w.r.t. the Wasserstein-1 ($W_{1}$) distance, is
upper bounded by a factor proportional to the strong convexity modulus of
$\ell({\bm{\theta}};z)$ (Mendler-Dünner et al., 2020). Furthermore, such
convergence condition is proven to be tight (Perdomo et al., 2020) and the
analysis has been extended to proximal algorithm (Drusvyatskiy and Xiao,
2023), online optimization (Cutler et al., 2023), saddle point seeking (Wood
and Dall’Anese, 2023), multi-agent consensus learning (Li et al., 2022), non-
cooperative learning (Wang et al., 2023; Narang et al., 2023; Piliouras and
Yu, 2023), state-dependent learning (Brown et al., 2022; Li and Wai, 2022),
etc.

This paper provides the _first analysis_ of SGD-GD and related stochastic
optimization schemes in performative prediction when $\ell({\bm{\theta}};z)$
is _smooth but possibly non-convex_. This is a more common scenario in machine
learning than the strongly convex loss considered in the prior works, e.g., it
covers the case of training neural network (NN) models. We notice that
existing works are limited to imposing structure on the loss function
$\ell({\bm{\theta}};z)$ and the distribution ${\cal D}({\bm{\theta}})$,
utilizing advanced algorithms that demand extra knowledge on ${\cal
D}({\bm{\theta}})$, etc., as we overview below.

Literature | Ncvx-$\ell$ | Ncvx-$V$ | Sensitivity | Algorithm | Rate | ${\bm{\theta}}_{\infty}$
---|---|---|---|---|---|---
(Izzo et al., 2021) | ✗ | ✓ | Loc.† | 2-Phase | ${\cal O}(T^{-\frac{1}{2}})$ | ${\nabla}V(\cdot)={\bm{0}}$
(Miller et al., 2021) | ✗ | ✗ | $W_{1}$ | 2-Phase | ${\cal O}(T^{-1})$ | $\min V({\bm{\theta}})$
(Mendler-Dünner et al., 2020) | ✗ | ✓ | $W_{1}$ | SGD-GD | ${\cal O}(T^{-1})$ | PS
(Mofakhami et al., 2023) | ✗‡ | ✓ | $\chi^{2}$ | RRM‡ | Linear‡ | PS‡
This Work | ✓ | ✓ |  $\chi^{2}$ or $W_{1}$ | SGD-GD | ${\cal O}(T^{-\frac{1}{2}})$ |  ${\cal O}(\epsilon)$-SPS

Table 1: Comparison of Results in Existing Works. ‘Sensitivity’ indicates the
distance metric imposed on ${\cal D}({\bm{\theta}})$ when the latter is
subject to perturbation, given in the form $d({\cal D}({\bm{\theta}}),{\cal
D}({\bm{\theta}}^{\prime}))\leq\epsilon\|{\bm{\theta}}-{\bm{\theta}}^{\prime}\|$
such that $d(\cdot,\cdot)$ is a distance metric between distributions.
‘$\theta_{\infty}$’ indicates the type of convergent points: ‘PS’ refers to
performative stable solution [cf. (4)], ‘SPS’ refers to Def. 1.
†Izzo et al. (2021) assumed that ${\cal D}({\bm{\theta}})$ belongs to the
location family, i.e., ${\cal D}({\bm{\theta}})={\cal
N}(f({\bm{\theta}});\sigma^{2})$.
‡Mofakhami et al. (2023) considered
$\ell({\bm{\theta}};z)=\tilde{\ell}(f_{{\bm{\theta}}}(x),y)$ with strongly
convex $\tilde{\ell}(\cdot,y)$. The RRM requires solving a non-convex
optimization at each recursion.

#### Related Works

In the non-convex setting, the most related work to ours is (Mofakhami et al.,
2023) which proved that a variant of PS solution can be found when training NN
in the performative prediction setting, i.e., a special case with non-convex
loss. However, their analysis is restrictive: (i) it requires a loss function
of the form $\ell({\bm{\theta}};z)=\tilde{\ell}(f_{\bm{\theta}}(x);y)$ where
$\tilde{\ell}(\hat{y};y)$ is strongly convex w.r.t. $\hat{y}$, (ii) it only
analyzes the case of training NN using a repeated risk minimization (RRM)
procedure which _exactly_ minimizes a non-convex objective function at each
step. In comparison, we concentrate on stochastic (first order) optimization
schemes and require only smoothness for $\ell(\cdot;z)$.

Other works departed from tackling the PS solution and considered alternative
algorithms to directly minimize $V({\bm{\theta}})$. For example, Roy et al.
(2022) assumed that unbiased estimates of ${\nabla}V({\bm{\theta}})$ is
available and studied the convergence of stochastic conditional gradient
algorithms towards a stationary solution of $V({\bm{\theta}})$, Li and Wai
(2022) assumed bounded biasedness w.r.t. ${\nabla}V({\bm{\theta}})$ in (2) and
show that (2) converges to a biased stationary point of $V({\bm{\theta}})$.
Notice that estimating ${\nabla}V({\bm{\theta}})$ requires knowledge on ${\cal
D}({\bm{\theta}})$ which has to be learnt separately. To circumvent this
difficulty, two phases algorithms are studied in (Miller et al., 2021; Izzo et
al., 2021) that learn ${\cal D}({\bm{\theta}})$ via a large batch of samples
at the first stage, then optimize ${\bm{\theta}}$ later (see (Zhu et al.,
2023) for a two-timescale type online algorithm), derivative free optimization
are studied in (Miller et al., 2021; Ray et al., 2022; Liu et al., 2023), and
confidence bound methods in (Jagadeesan et al., 2022). In addition, Miller et
al. (2021) proposed a mixture dominance assumption that can imply the strong
convexity of $V({\bm{\theta}})$. We remark that (Zhao, 2022) studied
conditions to ensure $V({\bm{\theta}})$ to be weakly convex. Our work does not
require such advanced algorithms and show that stochastic (first order)
optimization converges towards a similar solution as the PS solution. We
display a comparison between these related works in Table 1.

#### Our Contributions

This work provides the _first_ analysis of stochastic gradient-based methods
for performative prediction with smooth but possibly _non-convex_ losses. Our
contributions are:

  * •

We propose the concept of _stationary performative stable_ (SPS) solutions
which is a relaxation of the commonly used performative stable (PS) condition
(Perdomo et al., 2020). The relaxation is necessary for handling non-convex
losses using first-order methods.

  * •

We show that the stochastic gradient method with greedy deployment (SGD-GD)
finds a biased SPS solution. Assume that the distribution shift is
$\epsilon$-sensitive, i.e., it holds $d({\cal D}({\bm{\theta}}),{\cal
D}({\bm{\theta}}^{\prime}))\leq\epsilon\|{\bm{\theta}}-{\bm{\theta}}^{\prime}\|$
for some distance measure $d(\cdot,\cdot)$ between the shifted distributions
${\cal D}({\bm{\theta}}),{\cal D}({\bm{\theta}}^{\prime})$, SGD-GD converges
at a rate of ${\cal O}(T^{-\frac{1}{2}})$ in expectation to an ${\cal
O}(\epsilon)$-SPS solution. The bias level is further improved to ${\cal
O}(\epsilon^{2})$ when the gradient is exact.

  * •

Our analysis relies on constructing a time varying Lyapunov function that may
shed new lights for non-gradient stochastic approximation (Dieuleveut et al.,
2023). We studied two alternative conditions on the distance metric between
distributions. When $d(\cdot,\cdot)$ is the Wasserstein-1 distance, SGD-GD
converges to a biased SPS solution for Lipschitz loss function. When
$d(\cdot,\cdot)$ is the $\chi^{2}$-distance, SGD-GD converges to a biased SPS
solution for bounded loss function.

  * •

We extend the analysis to the lazy deployment scheme with SGD (Mendler-Dünner
et al., 2020). The latter scheme finds a _bias-free SPS solution_ as the epoch
length of lazy deployment grows.

Lastly, we provide numerical examples on synthetic and real data to validate
our theoretical findings. The rest of this paper is organized as follows. §2
introduces the problem setup and assumptions for establishing our convergence
results of SGD-GD. Furthermore, we highlight the challenges in analyzing non-
convex performative prediction. §3 introduces the concept of SPS solutions and
presents the convergence results for SGD-GD under two alternative assumptions
on the distribution shifts. We also outline the use of a time varying Lyapunov
function to handle the dynamic nature of SGD-GD. §4 shows the results for the
lazy deployment scheme. Lastly, §5 provides numerical experiments to
illustrate our results.

#### Notations

Let $\mathbb{R}^{d}$ be the $d$-dimensional Euclidean space equipped with
inner product $\left\langle\cdot\,|\,\cdot\right\rangle$ and induced norm
$\left\|x\right\|=\sqrt{\left\langle x\,|\,x\right\rangle}$. Let ${\cal S}$ be
a (measurable) sample space, and $\mu,\nu$ are two probability measures
defined as ${\cal S}$. $\mathbb{E}[\cdot]$ denotes taking expectation w.r.t
all randomness,
$\mathbb{E}_{t}[\cdot]\mathrel{\mathop{:}}=\mathbb{E}_{t}[\cdot|{\cal F}_{t}]$
means taking conditional expectation on the filtration ${\cal
F}_{t}\mathrel{\mathop{:}}=\sigma(\\{{\bm{\theta}}_{0},{\bm{\theta}}_{1},\cdots,{\bm{\theta}}_{t}\\})$,
where $\sigma(\cdot)$ is the sigma-algebra generated by the random variables
in the operand and $\\{{\bm{\theta}}_{t}\\}$ is the sequence of iterates
generated by the SGD-GD scheme (2).

##  2 Stationary Condition for Performative Stability

This section prepares the analysis of (1) with SGD-GD and related schemes in
the non-convex loss setting. To fix idea, we define the decoupled performative
risk and the decoupled partial gradient:

$\displaystyle J({\bm{\theta}}_{1};{\bm{\theta}}_{2})=\mathbb{E}_{Z\sim{\cal
D}({\bm{\theta}}_{2})}\left[\ell({\bm{\theta}}_{1};Z)\right],\quad{\nabla}J({\bm{\theta}}_{1};{\bm{\theta}}_{2})=\mathbb{E}_{Z\sim{\cal
D}({\bm{\theta}}_{2})}\left[{\nabla}\ell({\bm{\theta}}_{1};Z)\right].$ (3)

Observe that while $V({\bm{\theta}})=J({\bm{\theta}};{\bm{\theta}})$,
${\nabla}J({\bm{\theta}};{\bm{\theta}})\neq{\nabla}V({\bm{\theta}})$ in
general since ${\nabla}J({\bm{\theta}};{\bm{\theta}})$ only represents a
partial gradient of $V({\bm{\theta}})$; see (Izzo et al., 2021). In (2), the
conditional expectation of the stochastic gradient update term satisfies
$\mathbb{E}_{t}[{\nabla}\ell({\bm{\theta}}_{t};Z_{t+1})]={\nabla}J({\bm{\theta}}_{t};{\bm{\theta}}_{t})$.

If the loss $\ell({\bm{\theta}};z)$ is strongly convex w.r.t. ${\bm{\theta}}$,
then the decoupled performative risk $J({\bm{\theta}};\bar{{\bm{\theta}}})$
admits a unique minimizer w.r.t. ${\bm{\theta}}$ for any
$\bar{{\bm{\theta}}}$. It has hence motivated (Perdomo et al., 2020) to study
the _performative stable_ (PS) solution ${\bm{\theta}}_{PS}$ which is defined
as a fixed point to the map

$\displaystyle\textstyle{\sf
T}(\bar{\bm{\theta}}):=\mathop{\mathrm{arg\,min}}_{{\bm{\theta}}\in\mathbb{R}^{d}}J({\bm{\theta}};\bar{\bm{\theta}}),~{}~{}\text{i.e.,}~{}~{}{\bm{\theta}}_{PS}={\sf
T}({\bm{\theta}}_{PS}).$ (4)

In the above, the uniqueness and existence of ${\bm{\theta}}_{PS}$ follows by
observing that ${\sf T}(\cdot)$ is a contraction if and only if the
sensitivity of ${\cal D}({\bm{\theta}})$, i.e., the ‘smoothness’ of ${\cal
D}({\bm{\theta}})$ w.r.t. ${\bm{\theta}}$, is upper bounded by the inverse of
condition number of $\ell({\bm{\theta}};z)$. The convergence of SGD-GD follows
by analyzing (2) as a stochastic approximation (SA) scheme for the repeated
risk minimization (RRM) procedure ${\bm{\theta}}\leftarrow{\sf
T}({\bm{\theta}})$ (Mendler-Dünner et al., 2020).

For the case of _non-convex_ loss in this paper, the analysis becomes more
nuanced since the map ${\sf T}(\cdot)$ is no longer well-defined, e.g., there
may exist more than one minimizers to
$\min_{{\bm{\theta}}}J({\bm{\theta}};\bar{{\bm{\theta}}})$. Our remedy is to
concentrate on the following non-convex counter part to the PS solution:

######  Definition 1.

($\delta$ stationary performative stable solution) Let $\delta\geq 0$, the
vector ${\bm{\theta}}_{\delta-SPS}\in\mathbb{R}^{d}$ is said to be an $\delta$
stationary performative stable ($\delta$-SPS) solution of (1) if:

$\displaystyle\left\|{\nabla}J({\bm{\theta}}_{\delta-
SPS};{\bm{\theta}}_{\delta-SPS})\right\|^{2}=\left\|\mathbb{E}_{Z\sim{\cal
D}({\bm{\theta}}_{\delta-SPS})}\left[{\nabla}\ell({\bm{\theta}}_{\delta-
SPS};Z)\right]\right\|^{2}\leq\delta.$ (5)

We also say that ${\bm{\theta}}_{SPS}$ is an (exact) SPS solution if it
satisfies (5) with $\delta=0$. In other words, $\delta\geq 0$ measures the
stationarity of a solution. Notice that if $\ell({\bm{\theta}};z)$ is strongly
convex w.r.t. ${\bm{\theta}}$, then an SPS solution is also a PS solution
defined in (Perdomo et al., 2020).

Although (5) is similar to the usual definitions of stationary solution in
smooth optimization, there is a subtle but critical difference since
${\nabla}J({\bm{\theta}};{\bm{\theta}})$ may not be the _gradient_ of any
function in ${\bm{\theta}}$. For example, consider
$\ell({\bm{\theta}};z)=(1/2)\|{\bm{\theta}}-z\|^{2}$ and ${\cal
D}({\bm{\theta}})\equiv{\cal N}(A{\bm{\theta}},{\bm{I}})$ for some square but
asymmetric matrix $A$, the map
${\nabla}J({\bm{\theta}};{\bm{\theta}})=(I-A){\bm{\theta}}$ has a Jacobian of
$I-A$ which is not symmetric. Furthermore, we observe that the mean field for
SGD-GD scheme (2) is ${\nabla}J({\bm{\theta}};{\bm{\theta}})$, which is not a
gradient. The SGD-GD scheme is thus a special case of non-gradient SA scheme
(Dieuleveut et al., 2023).

To get further insight, as investigated in (Dieuleveut et al., 2023), a common
analysis framework of non-gradient SA scheme is by identifying a smooth
Lyapunov function linked to the recursion (2). When $\ell(\cdot;z)$ is
strongly convex, we may study the Lyapunov functions as the squared distance
$\|{\bm{\theta}}-{\bm{\theta}}_{PS}\|^{2}$ [cf. (4)]. It can be shown that the
properties required in (Dieuleveut et al., 2023) are satisfied under the
conditions analyzed by (Mendler-Dünner et al., 2020). However, in the case of
non-convex loss, identifying a suitable Lyapunov function is hard. In the next
section, we demonstrate how to address this challenge by identifying a time
varying Lyapunov function for SGD-GD.

##  3 Main Results

This section presents theoretical results on the SGD-GD scheme with non-convex
loss. We first show how to construct a time varying Lyapunov function tailor
made for (2). We then show the convergence of SGD-GD under two different sets
of conditions.

We introduce two basic and natural assumptions on the risk minimization
problem (1):

######  A 1.

For any $z\in{\sf Z}$, there exists a constant $L\geq 0$ such that

$\displaystyle\left\|{\nabla}\ell({\bm{\theta}};z)-{\nabla}\ell({\bm{\theta}}^{\prime};z)\right\|\leq
L\left\|{\bm{\theta}}-{\bm{\theta}}^{\prime}\right\|,~{}\forall~{}{\bm{\theta}},{\bm{\theta}}^{\prime}\in\mathbb{R}^{d},$
(6)

where ${\nabla}\ell({\bm{\theta}};z)$ denotes the gradient of
$\ell({\bm{\theta}};z)$ w.r.t. ${\bm{\theta}}$. Moreover, there exists a
constant $\ell^{\star}>-\infty$ such that
$\ell({\bm{\theta}};z)\geq\ell^{\star}$ for any
${\bm{\theta}}\in\mathbb{R}^{d}$.

######  A 2.

For any fixed ${\bm{\theta}}_{1},{\bm{\theta}}_{2}\in\mathbb{R}^{d}$, the
stochastic gradient ${\nabla}\ell({\bm{\theta}}_{1};Z)$, $Z\sim{\cal
D}({\bm{\theta}}_{2})$ is unbiased such that $\mathbb{E}_{Z\sim{\cal
D}({\bm{\theta}}_{2})}[{\nabla}\ell({\bm{\theta}}_{1};Z)]={\nabla}J({\bm{\theta}}_{1};{\bm{\theta}}_{2})$,
and there exists constants $\sigma_{0},\sigma_{1}\geq 0$ such that

$\displaystyle\mathbb{E}_{Z\sim{\cal
D}({\bm{\theta}}_{2})}\left[\left\|{\nabla}\ell({\bm{\theta}}_{1};Z)-{\nabla}J({\bm{\theta}}_{1};{\bm{\theta}}_{2})\right\|^{2}\right]\leq\sigma_{0}^{2}+\sigma_{1}^{2}\left\|{\nabla}J({\bm{\theta}}_{1};{\bm{\theta}}_{2})\right\|^{2}.$
(7)

Note that A1, 2 are standard assumptions that hold for a wide range of
applications and the respective stochastic optimization based training
schemes. For instance, A1 requires the loss function to be smooth, while A2
assumes the stochastic gradient estimates to have a variance that may grow
with $\left\|{\nabla}J({\bm{\theta}}_{1};{\bm{\theta}}_{2})\right\|^{2}$. For
the _non performative prediction_ setting where the data distribution is not
shifted by ${\bm{\theta}}$, i.e., ${\cal D}({\bm{\theta}})\equiv{\cal D}$,
these assumptions guarantee that SGD algorithm (i.e., SGD-GD with
$Z_{t+1}\sim{\cal D}$) to converge to a stationary solution of (1) with a
suitable step size schedule (Ghadimi and Lan, 2013). In particular, in this
case A1 implies that $\mathbb{E}_{Z\sim{\cal D}}[\ell({\bm{\theta}};Z)]$ is a
smooth function and serves as a Lyapunov function for the SGD algorithm.

In this light, it might be tempting to use the performative risk
$V({\bm{\theta}})$ [cf. (1)] as the Lyapunov function for (2) and directly
adopt the analysis in (Dieuleveut et al., 2023). However, the condition A1 is
insufficient to guarantee that $V({\bm{\theta}})$ is smooth, and the mean
field of (2) may not be aligned with ${\nabla}V({\bm{\theta}})$. Instead, from
A1, 2, we proceed with a descent-like lemma for the iterates of SGD-GD:

######  Lemma 3.1.

Under A1, 2. Suppose that the step size satisfies $\sup_{t\geq
1}\gamma_{t}\leq 1/L(1+\sigma_{1}^{2})$, then for any $t\geq 0$, the sequence
of iterates $\\{{\bm{\theta}}_{t}\\}_{t\geq 0}$ generated by SGD-GD (2)
satisfies

$\displaystyle\frac{\gamma_{t+1}}{2}\|{\nabla}J({\bm{\theta}}_{t};{\bm{\theta}}_{t})\|^{2}$
$\displaystyle\leq
J({\bm{\theta}}_{t};{\bm{\theta}}_{t})-\mathbb{E}_{t}[J({\bm{\theta}}_{t+1};{\bm{\theta}}_{t})]+\frac{L}{2}\sigma_{0}^{2}\gamma_{t+1}^{2}.$
(8)

###### Proof.

Fix any $z\in{\sf Z}$, applying A1 and the recursion (2) lead to

$\displaystyle\ell({\bm{\theta}}_{t+1};z)$
$\displaystyle\leq\ell({\bm{\theta}}_{t};z)-\gamma_{t+1}\left\langle{\nabla}\ell({\bm{\theta}}_{t};z)\,|\,{\nabla}\ell({\bm{\theta}}_{t};Z_{t+1})\right\rangle+\frac{L\gamma_{t+1}^{2}}{2}\left\|{\nabla}\ell({\bm{\theta}}_{t};Z_{t+1})\right\|^{2},$
(9)

for any $t\geq 0$. Note that $z\in{\sf Z}$ can be any fixed sample while
$Z_{t+1}$ is the r.v. drawn from ${\cal D}({\bm{\theta}}_{t})$ in (2). Let
$p_{{\bm{\theta}}_{t}}(z)\geq 0$ denotes the pdf of ${\cal
D}({\bm{\theta}}_{t})$. We then multiply $p_{{\bm{\theta}}_{t}}(z)$ on both
sides of the inequality and integrate w.r.t. $z\in{\sf Z}$, i.e., taking the
operator $\int_{\sf Z}(\cdot)p_{{\bm{\theta}}_{t}}(z){\sf d}z$. This yields

$\displaystyle J({\bm{\theta}}_{t+1};{\bm{\theta}}_{t})$ $\displaystyle\leq
J({\bm{\theta}}_{t};{\bm{\theta}}_{t})-\gamma_{t+1}\left\langle{\nabla}J({\bm{\theta}}_{t};{\bm{\theta}}_{t})\,|\,{\nabla}\ell({\bm{\theta}}_{t};Z_{t+1})\right\rangle+\frac{L\gamma_{t+1}^{2}}{2}\left\|{\nabla}\ell({\bm{\theta}}_{t};Z_{t+1})\right\|^{2},$
(10)

since $\int\ell({\bm{\theta}};z)p_{{\bm{\theta}}_{t}}(z){\sf
d}z=J({\bm{\theta}};{\bm{\theta}}_{t})$ according to definition (3). We next
evaluate the conditional expectation, $\mathbb{E}_{t}[\cdot]$, on both sides
of the above inequality

$\displaystyle\mathbb{E}_{t}[J({\bm{\theta}}_{t+1};{\bm{\theta}}_{t})]$
$\displaystyle\leq
J({\bm{\theta}}_{t};{\bm{\theta}}_{t})-\gamma_{t+1}\|{\nabla}J({\bm{\theta}}_{t};{\bm{\theta}}_{t})\|^{2}+\frac{L\gamma_{t+1}^{2}}{2}\mathbb{E}_{t}[\left\|{\nabla}\ell({\bm{\theta}}_{t};Z_{t+1})\right\|^{2}].$
(11)

Using A2, we note that
$\mathbb{E}_{t}[\left\|{\nabla}\ell({\bm{\theta}}_{t};Z_{t+1})\right\|^{2}]\leq\sigma_{0}^{2}+(1+\sigma_{1}^{2})\|{\nabla}J({\bm{\theta}}_{t};{\bm{\theta}}_{t})\|^{2}$.
Reshuffling terms and using the step size condition yield the desired result
(8); see §A for a detailed proof. ∎

For sufficiently small $\gamma_{t+1}>0$ and when ${\bm{\theta}}_{t}$ is not
SPS, eq. (8) implies the descent relation
$\mathbb{E}_{t}[J({\bm{\theta}}_{t+1};{\bm{\theta}}_{t})]<J({\bm{\theta}}_{t};{\bm{\theta}}_{t})$.
This suggests that at the $t$th iteration, the function
$J_{t}({\bm{\theta}}):=J({\bm{\theta}};{\bm{\theta}}_{t})$ may serve as a
Lyapunov function for the SGD-GD scheme. Meanwhile, $J_{t}({\bm{\theta}})$ is
a _time varying_ Lyapunov function. The said descent relation does not
necessarily imply the convergence towards an SPS solution. Instead, the first
term on the right hand side of (8) can be decomposed as:

$\displaystyle\mathbb{E}[J_{t}({\bm{\theta}}_{t})-J_{t}({\bm{\theta}}_{t+1})]=\mathbb{E}[J_{t}({\bm{\theta}}_{t})-J_{t+1}({\bm{\theta}}_{t+1})]+\underbrace{\mathbb{E}[J_{t+1}({\bm{\theta}}_{t+1})-J_{t}({\bm{\theta}}_{t+1})]}_{\text{residual}}.$
(12)

The first part is a difference-of-sequence which is summable, while the second
part is a residual term. The convergence of SGD-GD with non-convex losses
hinges on bounding the latter residual. Taking a closer look, the residual is
the difference of evaluating ${\bm{\theta}}_{t+1}$ on $J_{t}(\cdot)$ and
$J_{t+1}(\cdot)$, i.e.,

$\displaystyle\mathbb{E}[J_{t+1}({\bm{\theta}}_{t+1})-J_{t}({\bm{\theta}}_{t+1})]=\mathbb{E}\left[\mathbb{E}_{Z\sim{\cal
D}({\bm{\theta}}_{t}),Z^{\prime}\sim{\cal
D}({\bm{\theta}}_{t+1})}[\ell({\bm{\theta}}_{t+1};Z^{\prime})-\ell({\bm{\theta}}_{t+1};Z)]\right].$
(13)

The above further depends on the differences between the distributions ${\cal
D}({\bm{\theta}}_{t}),{\cal D}({\bm{\theta}}_{t+1})$, i.e., the _sensitivity_
of the data distribution w.r.t. perturbation in ${\bm{\theta}}$. Next, we
study sufficient conditions that imply the convergence of SGD-GD through
bounding
$\mathbb{E}[J_{t+1}({\bm{\theta}}_{t+1})-J_{t}({\bm{\theta}}_{t+1})]$.

###  3.1 Sufficient Conditions for Convergence of SGD-GD

From (8), we anticipate the convergence of SGD-GD towards a biased SPS
solution if it holds
$\mathbb{E}[J_{t+1}({\bm{\theta}}_{t+1})-J_{t}({\bm{\theta}}_{t+1})]={\cal
O}(\mathbb{E}[\|{\bm{\theta}}_{t}-{\bm{\theta}}_{t+1}\|])={\cal
O}(\gamma_{t+1}\mathbb{E}[\|{\nabla}\ell({\bm{\theta}}_{t};Z_{t+1})]\|)$. Now,
as seen from (13), establishing such relation would require ${\cal
D}({\bm{\theta}})$ to satisfy a certain sensitivity criterion when subject to
perturbation in ${\bm{\theta}}$. Our subsequent discussions are organized
according to various distributional distance measures on sensitivity.

#### Wasserstein-1 Sensitivity

Our first set of conditions uses the Wasserstein-1 distance for measuring the
sensitivity of data distribution. The measure is commonly used in the studies
of performative prediction, e.g., as pioneered by (Perdomo et al., 2020;
Mendler-Dünner et al., 2020):

######  W 1.

($\epsilon$ sensitivity w.r.t. Wasserstein-1 distance) There exists
$\epsilon\geq 0$ such that

$\displaystyle W_{1}({\cal D}({\bm{\theta}}),{\cal
D}({\bm{\theta}}^{\prime}))\leq\epsilon\|{\bm{\theta}}-{\bm{\theta}}^{\prime}\|,$
(14)

for any ${\bm{\theta}},{\bm{\theta}}^{\prime}\in\mathbb{R}^{d}$. Notice that
the Wasserstein-1 distance is defined as
$W_{1}(\cdot,\cdot):=\inf_{P\in\mathcal{P}(\cdot,\cdot)}\mathbb{E}_{\left(z,z^{\prime}\right)\sim
P}\left[\left\|z-z^{\prime}\right\|_{1}\right]$, where
$\mathcal{P}\left(\mathcal{D}({\bm{\theta}}),\mathcal{D}\left({\bm{\theta}}^{\prime}\right)\right)$
is the set of all joint distributions on ${\sf Z}\times{\sf Z}$ whose marginal
distributions are
$\mathcal{D}({\bm{\theta}}),\mathcal{D}\left({\bm{\theta}}^{\prime}\right)$.

In this case, we require the loss function to be Lipschitz continuous w.r.t.
shifts in the data sample $z$.

######  W 2.

There exists a constant $L_{0}>0$ such that for all $z,z^{\prime}\in{\sf Z},$
and ${\bm{\theta}}\in\mathbb{R}^{d}$,

$\displaystyle|\ell({\bm{\theta}};z)-\ell({\bm{\theta}};z^{\prime})|\leq
L_{0}\left\|z-z^{\prime}\right\|.$ (15)

Our key observation is that the above conditions imply the desired Lipschitz
continuity property on $J({\bm{\theta}};\cdot)$. In fact, we have:

######  Lemma 3.2.

Under W1, 2. For any
${\bm{\theta}}_{1},{\bm{\theta}}_{2},{\bm{\theta}}\in\mathbb{R}^{d}$, it holds

$\displaystyle|J({\bm{\theta}};{\bm{\theta}}_{1})-J({\bm{\theta}};{\bm{\theta}}_{2})|\leq
L_{0}\epsilon\left\|{\bm{\theta}}_{1}-{\bm{\theta}}_{2}\right\|.$ (16)

The proof, which is a variant of (Drusvyatskiy and Xiao, 2023, Lemma 2.1), can
be found in §B.

#### Pearson $\chi^{2}$ Sensitivity

Although W1 holds for a number of applications such as location family
distribution (e.g., (Miller et al., 2021; Perdomo et al., 2020; Narang et al.,
2023)), the Lipschitz loss assumption in W2 can be difficult to verify
especially if we want $L_{0}$ to be small. As an alternative, we consider a
slightly stronger sensitivity condition on ${\cal D}({\bm{\theta}})$ via the
Pearson’s $\chi^{2}$ divergence.

######  C 1.

($\epsilon$ sensitivity w.r.t. Pearson $\chi^{2}$ divergence) For any
${\bm{\theta}}\in\mathbb{R}^{d}$, ${\cal D}({\bm{\theta}})$ is not
degenerated, i.e., $p_{{\bm{\theta}}}(\cdot)\neq 0$ almost everywhere. For any
${\bm{\theta}},{\bm{\theta}}^{\prime}\in\mathbb{R}^{d}$, there exists a
constant $\epsilon\geq 0$ such that

$\displaystyle\chi^{2}({\cal D}({\bm{\theta}}),{\cal
D}({\bm{\theta}}^{\prime}))\mathrel{\mathop{:}}=\int\frac{\left(p_{{\bm{\theta}}}(z)-p_{{\bm{\theta}}^{\prime}}(z)\right)^{2}}{p_{{\bm{\theta}}}(z)}{\sf
d}z\leq\epsilon^{2}\left\|{\bm{\theta}}-{\bm{\theta}}^{\prime}\right\|^{2},$
(17)

where $p_{\bm{\theta}}(\cdot),p_{{\bm{\theta}}^{\prime}}(\cdot)$ are p.d.f.s
of induced distributions ${\cal D}({\bm{\theta}}),{\cal
D}({\bm{\theta}}^{\prime})$, respectively.

This setup is modified from (Mofakhami et al., 2023) which used the condition
$\chi^{2}({\cal D}({\bm{\theta}}),{\cal D}({\bm{\theta}}^{\prime}))={\cal
O}(\|f_{\bm{\theta}}(\cdot)-f_{{\bm{\theta}}^{\prime}}(\cdot)\|^{2})$ such
that $f_{\bm{\theta}}(\cdot)$ denotes the output of a neural network
parameterized by ${\bm{\theta}}$. In fact, for bounded sample space ${\sf Z}$,
it holds

$W_{1}({\cal D}({\bm{\theta}}),{\cal D}({\bm{\theta}}^{\prime}))\leq({\rm
diam}({\sf Z})/\sqrt{2})\sqrt{\chi^{2}({\cal D}({\bm{\theta}}),{\cal
D}({\bm{\theta}}^{\prime}))},$

see (Gibbs and Su, 2002, Sec. 2), where ${\rm diam}({\sf
Z}):=\sup_{z,z^{\prime}\in{\sf Z}}\|z-z^{\prime}\|$ is the diameter of the
sample space. We remark that unlike $W_{1}(\cdot,\cdot)$, the measure
$\chi^{2}(\cdot,\cdot)$ is not a metric distance as $\chi^{2}({\cal
D}({\bm{\theta}}),{\cal D}({\bm{\theta}}^{\prime}))\neq\chi^{2}({\cal
D}({\bm{\theta}}^{\prime}),{\cal D}({\bm{\theta}}))$.

Although C1 is slightly strengthened from W1, it allows us to replace the
Lipschitz continuity assumption W2 on the loss. Particularly, we replace W2
by:

######  C 2.

There exists a constant $\ell_{max}\geq 0$ such that
$\sup_{{\bm{\theta}}\in\mathbb{R}^{d},z\in{\sf
Z}}|\ell({\bm{\theta}};z)|\leq\ell_{max}$.

The above condition requires $\ell(\cdot;\cdot)$ to be uniformly bounded.
Compared to W2, it can be easier to verify and $\ell_{max}$ is typically
small. For example, it holds with $\ell_{max}=1$ for the case of sigmoid loss.
Similar to W1, 2, we observe that the above conditions imply
$J({\bm{\theta}};\cdot)$ is Lipschitz continuous:

######  Lemma 3.3.

Under C1, 2. For any
${\bm{\theta}}_{1},{\bm{\theta}}_{2},{\bm{\theta}}\in\mathbb{R}^{d}$, it holds
that

$\displaystyle|J({\bm{\theta}};{\bm{\theta}}_{1})-J({\bm{\theta}};{\bm{\theta}}_{2})|\leq\ell_{max}\epsilon\left\|{\bm{\theta}}_{1}-{\bm{\theta}}_{2}\right\|.$
(18)

See §C. The only difference with Lemma 3.2 is that (18) has a different
Lipschitz constant.

###  3.2 Convergence of SGD-GD with Non-convex Loss

Equipped with Lemmas 3.2, 3.3, we are ready to present the convergence result
for SGD-GD with smooth but non-convex losses. Observe the following theorem
whose proof can be found in §D:

######  Theorem 1.

Under A1, 2. Let the step sizes satisfy $\sup_{t\geq
1}\gamma_{t}\leq{1}/{L(1+\sigma_{1}^{2})}$. Moreover, let
$\displaystyle\tilde{L}=L_{0}~{}~{}\text{if~{}~{}W\ref{assu:w1},
\ref{assu:lip_loss}
hold},~{}~{}\text{or}~{}~{}\tilde{L}=\ell_{max}~{}~{}\text{if~{}~{}C\ref{assu:chi},
\ref{assu:bd_loss} hold}.$ (19) Then, for any $T\geq 1$, the iterates
$\\{{\bm{\theta}}_{t}\\}_{t\geq 0}$ generated by SGD-GD satisfy
$\displaystyle\sum_{t=0}^{T-1}\frac{\gamma_{t+1}}{4}\mathbb{E}[\left\|{\nabla}J({\bm{\theta}}_{t};{\bm{\theta}}_{t})\right\|^{2}]\leq\Delta_{0}+\tilde{L}\epsilon\left(\sigma_{0}+(1+\sigma_{1}^{2})\tilde{L}\epsilon\right)\sum_{t=0}^{T-1}\gamma_{t+1}+\frac{L}{2}\sigma_{0}^{2}\sum_{t=0}^{T-1}\gamma_{t+1}^{2},$
(20) where
$\Delta_{0}\mathrel{\mathop{:}}=J({\bm{\theta}}_{0};{\bm{\theta}}_{0})-\ell_{\star}$
is an upper bound to the initial optimality gap of performative risk.

Using a fixed step size schedule, we simplify the bound as:

######  Corollary 3.1.

Under A1, 2, the alternative conditions W1, 2, or C1, 2. Let $T\geq 1$ be the
maximum number of iterations and set $\gamma_{t}=1/\sqrt{T}$. Let ${\sf T}$ be
a random variable chosen uniformly and independently from
$\\{0,1,\cdots,T-1\\}$. For any $T\geq L^{2}(1+\sigma_{1}^{2})^{2}$, the
iterates by SGD-GD satisfy

$\displaystyle\mathbb{E}\left[\left\|{\nabla}J({\bm{\theta}}_{\sf
T};{\bm{\theta}}_{\sf T})\right\|^{2}\right]\leq
4\left(\Delta_{0}+\frac{L}{2}\sigma_{0}^{2}\right)\cdot\frac{1}{\sqrt{T}}+\underbrace{4\tilde{L}\epsilon\,(\sigma_{0}+(1+\sigma_{1}^{2})\tilde{L}\epsilon)}_{\sf{\bf
bias}}.$ (21)

As $T\to\infty$, the first term in (21) vanishes as ${\cal O}(1/\sqrt{T})$ and
the above shows that the SGD-GD scheme finds an ${\cal
O}(\sigma_{0}\,\epsilon+(1+\sigma_{1}^{2})\,\epsilon^{2})$-SPS solution. This
yields the first convergence guarantee for performative prediction with non-
convex loss via a stochastic optimization scheme.

Lastly, an interesting observation is that the bias level is controlled at
${\cal O}(\sigma_{0}\,\epsilon+(1+\sigma_{1}^{2})\,\epsilon^{2})$. The latter
estimate highlights the role of the stochastic gradient’s variance. To see
this, let us concentrate on the case when $\epsilon$ is small. When stochastic
gradient is used such that $\sigma_{0}>0$, SGD-GD finds an ${\cal
O}(\epsilon)$-SPS solution; while with deterministic gradient, i.e., when
${\nabla}\ell({\bm{\theta}}_{t};Z_{t+1})={\nabla}J({\bm{\theta}}_{t};{\bm{\theta}}_{t})$
with $\sigma_{0}=\sigma_{1}=0$, SGD-GD finds an ${\cal O}(\epsilon^{2})$-SPS
solution. Such a distinction in the bias levels indicate that a _unique
property_ of non-convex performative prediction where the asymptotic
performance of SGD-GD is sensitive to the noise variance. Furthermore, our
result suggests that adjusting the minibatch size in SGD-GD may have a
significant effect on reducing the bias level since $\sigma_{0},\sigma_{1}$
can be controlled by the latter.

######  Remark 1.

Prior analysis in (Mendler-Dünner et al., 2020; Drusvyatskiy and Xiao, 2023)
showed that with $\mu$ strongly convex loss $\ell(\cdot;z)$, both the
existence/uniqueness of the PS solution and the convergence of SGD-GD to the
PS solution critically depend on the condition $\epsilon<\mu/L$ (in addition
to our A1, 2, W1). When $\epsilon>\mu/L$, it is shown that the SGD-GD scheme
may even diverge. In contrary, Theorem 1 does not exhibit such an explicit
condition on $\epsilon$ for the convergence results (20), (21) to hold. This
happens because our result requires the loss function itself to be Lipschitz
[cf. W2] or bounded [cf. C2], which may not be satisfied by their strongly
convex losses.

##  4 Extension: Lazy Deployment Scheme with SGD

Implementing the SGD-GD scheme (2) requires deploying the latest model every
time when drawing samples from ${\cal D}(\cdot)$. This may be difficult to
realize since deploying a new classifier in real time can be time consuming.
As inspired by (Mendler-Dünner et al., 2020), this section studies an
extension of (2) to the _lazy deployment_ scheme where the new (prediction)
models are deployed only once per several SGD updates.

To describe the extended scheme, let $K\geq 1$ denotes the epoch length of
lazy deployment, we have

$\displaystyle\begin{aligned}
&{\bm{\theta}}_{t,k+1}={\bm{\theta}}_{t,k}-\gamma{\nabla}\ell({\bm{\theta}}_{t,k};Z_{t,k+1}),\text{
where }Z_{t,k+1}\sim{\cal D}({\bm{\theta}}_{t}),~{}k=0,...,K-1,\\\
&{\bm{\theta}}_{t+1}={\bm{\theta}}_{t+1,0}={\bm{\theta}}_{t,K}.\end{aligned}$
(22)

For simplicity, we focus on the case with a constant step size $\gamma>0$.
Observe that the lazy deployment scheme is a double-loop algorithm where the
index $t$ denotes the number of deployments and the index $k$ denotes the SGD
update. To analyze the convergence of (22), we further require the stochastic
gradient to be uniformly bounded:

######  A 3.

There exists a constant $G\geq 0$ such that
$\sup_{{\bm{\theta}}\in\mathbb{R}^{d},z\in{\sf
Z}}\left\|{\nabla}\ell({\bm{\theta}};z)\right\|\leq G$.

Despite being a stronger assumption, the above remains valid for practical
non-convex losses, e.g., sigmoid loss. We observe the following convergence
results whose proof is in §E:

######  Theorem 2.

Under A1, 2, 3, and the alternative conditions W1, 2, or C1, 2. Let $T\geq 1$
be the maximum number of deployments to be run, we set $\gamma=1/(K\sqrt{T})$
and let ${\sf T}$ be a random variable chosen uniformly and independently from
$\\{0,1,\cdots,T-1\\}$. For any $T\geq L^{2}(1+\sigma_{1}^{2})^{2}/K^{2}$, the
iterates generated by the lazy deployment scheme with SGD (22) satisfy:
$\displaystyle\mathbb{E}\left[\left\|{\nabla}J({\bm{\theta}}_{\sf
T};{\bm{\theta}}_{\sf
T})\right\|^{2}\right]\leq\frac{8\Delta_{0}}{\sqrt{T}}+\frac{4L\sigma_{0}^{2}}{K\sqrt{T}}+\frac{2LG^{2}}{3T}+\frac{8\tilde{L}\epsilon}{K}\left(\sqrt{2K}\sigma_{0}+2\sqrt{2(K+\sigma_{1}^{2})}\tilde{L}\epsilon\right),$
(23) where we recall that $\tilde{L}$ was defined in (19) and
$\Delta_{0}=J({\bm{\theta}}_{0};{\bm{\theta}}_{0})-\ell_{\star}$.

In (23), the first three terms decay as ${\cal O}(1/\sqrt{T})$ similar to SGD-
GD, the last term simplifies to ${\cal O}(\tilde{L}\epsilon/\sqrt{K})$ to be
controlled with $1/\sqrt{K}$. The lazy deployment scheme (22) finds a _bias-
free SPS solution_ when $T\to\infty,K\to\infty$, contrasting with SGD-GD which
admits a bias level of ${\cal O}(\tilde{L}\epsilon)$.

We remark that the above result can be anticipated. During the $t$th
deployment, (22) runs an SGD recursion for
$\min_{{\bm{\theta}}}J({\bm{\theta}};{\bm{\theta}}_{t})$ where it will find a
stationary solution for the non-convex optimization as $K\to\infty$. The lazy
deployment scheme resembles RRM and we expect that it may find a bias-free SPS
solution as inspired by (Mofakhami et al., 2023) which studied a similar
algorithm.

##  5 Numerical Experiments

We consider two examples of performative prediction with non-convex loss based
on synthetic data and real data. All simulations are performed with Pytorch on
a server using a Intel Xeon 6318 CPU. Additional results can be found in §F.

Synthetic Data with Linear Model. We first consider a binary classification
problem using linear model. To enhance robustness to outliers, we adopt the
sigmoid loss function (Ertekin et al., 2010):

$\displaystyle\ell({\bm{\theta}};z)\mathrel{\mathop{:}}=\left(1+\exp(c\cdot
y\langle
x\,|\,{\bm{\theta}}\rangle)\right)^{-1}+(\beta/2)\left\|{\bm{\theta}}\right\|^{2}.$
(24)

For small regularization $\beta>0$, $\ell(\cdot;z)$ is smooth but non-convex.
To define the data distribution, we have a set of $m$ unshifted samples ${\cal
D}^{o}\equiv\\{(x_{i},y_{i})\\}_{i=1}^{m}$ with feature
$x_{i}\in\mathbb{R}^{d}$ and label $y_{i}\in\\{\pm 1\\}$. For any
${\bm{\theta}}\in\mathbb{R}^{d}$, ${\cal D}({\bm{\theta}})$ is a uniform
distribution on $m$ _shifted samples_
$\\{(x_{i}-\epsilon_{L}{\bm{\theta}},y_{i})\\}_{i=1}^{m}$, where
$\epsilon_{L}>0$ controls the shift magnitude. Applying SGD-GD to the setup
yields a scheme such that A1, 2, W1 (with $\epsilon=\epsilon_{L}$) are
satisfied, and W2 holds as $\|{\bm{\theta}}^{t}\|$ is bounded in practice. To
generate the training data, the unshifted samples ${\cal D}^{o}$ are generated
first as $x_{i}\sim{\cal U}[-1,1]^{d}$, i.e., the uniform distribution,
$\bar{y}_{i}=\text{sgn}(\langle x_{i}\,|\,{\bm{\theta}}^{o}\rangle)\in\\{\pm
1\\}$ such that ${\bm{\theta}}^{o}\sim{\cal N}(0,{\bm{I}})$, then a randomly
selected 10% of the labels are flipped to generate the final $y_{i}$.
Furthermore, we set
$m=800,d=10,c=0.1,\beta=10^{-3},\epsilon\in\\{0,0.1,0.5,2\\}$. For (2), the
batch size is $b=1$ and the stepsize is $\gamma_{t}=\gamma=1/\sqrt{T}$ with
$T=10^{6}$.

First, we validate the convergence behavior of SGD-GD in Theorem 1. In Fig. 1
(left), we compare
$\left\|{\nabla}J({\bm{\theta}}_{t};{\bm{\theta}}_{t})\right\|^{2}$ against
the number of iteration $t$ for the SGD-GD scheme over $10$ repeated runs. The
shaded region indicate the $95\%$ confidence interval. We observe that after a
rapid transient stage, the SPS stationarity
$\left\|{\nabla}J({\bm{\theta}}_{t};{\bm{\theta}}_{t})\right\|^{2}$ saturates
and stay around a constant level, indicating that the SGD-GD converges to a
biased-SPS solution. Increasing $\epsilon_{L}$ leads to an increased bias,
corroborating with Theorem 1 that the bias level is ${\cal O}(\epsilon)$. Fig.
1 (middle) further evaluate the performance of the trained classifier
${\bm{\theta}}_{t}$ in terms of the performative risk value. Second, we
compare the lazy deployment scheme in §4 with $K\in\\{5,10\\}$ and stepsize
$\gamma=1/(K\sqrt{T})$. For fairness, we test SGD-GD with batch size of
$b\in\\{5,25\\}$ and compare the SPS stationarity against the number of
samples accessed. The results in Fig. 1 (right) verifies Theorem 2 where
increasing $K$ effectively reduces the bias level.

Figure 1: Synthetic Data (_left_) SPS measure
$\|{\nabla}J({\bm{\theta}}_{t};{\bm{\theta}}_{t})\|^{2}$ of SGD-GD against
iteration no. $t$. (_middle_) Loss value
$J({\bm{\theta}}_{t};{\bm{\theta}}_{t})$ of SGD-GD against iteration no. $t$.
(_right_) SPS measure $\|{\nabla}J({\bm{\theta}}_{t};{\bm{\theta}}_{t})\|^{2}$
of greedy (SGD-GD) and lazy deployment against number of sample accessed. We
fix $\epsilon_{L}=2$.

Figure 2: Real Data with Neural Network Benchmarking with SPS measure
$\|{\nabla}J({\bm{\theta}}_{t};{\bm{\theta}}_{t})\|^{2}$. (_left_) Against $t$
for SGD-GD with parameters $\epsilon_{\sf NN}\in\\{0,10,100\\}$. (_middle_ &
_right_) Against no. of samples with greedy (SGD-GD) and lazy deployment when
$\epsilon_{\sf NN}=10$ & $\epsilon_{\sf NN}=10^{4}$, respectively.

Real Data with Neural Network. Our second example deals with the task of
training a neural network (NN) on the spambase Hopkins et al. (1999) dataset
with $m=4601$ samples, each with $d=48$ features. We split the training/test
sets as 8 : 2. Our aim is to study the behavior of SGD-GD when training NN
classifier. To specify (1), we let $z\equiv(x,y)$ where $x\in\mathbb{R}^{d}$
is the feature vector, $y\in\\{0,1\\}$ is label (0 for not spam, 1 for spam).
Consider the regularized binary cross entropy loss:

$\displaystyle\ell({\bm{\theta}};z)\equiv\tilde{\ell}(f_{{\bm{\theta}}}(x);y)=-y\log(f_{{\bm{\theta}}}(x))-(1-y)\log(1-f_{{\bm{\theta}}}(x))+({\beta}/{2})\left\|{\bm{\theta}}\right\|^{2},$
(25)

where $f_{{\bm{\theta}}}(x)$ denotes the NN classifier. The unshifted data is
denoted by ${\cal D}^{o}=\\{(x_{i},y_{i})\\}_{i=1}^{m}$. Sampling from the
shifted data distribution ${\cal D}({\bm{\theta}})$ is achieved through (i)
uniformly draw a sample $\bar{z}\equiv(\bar{x},\bar{y})$ from ${\cal D}^{o}$,
(ii) maximize the following utility function:

$\displaystyle\textstyle
x=\arg\max_{x^{\prime}}U(x^{\prime};\bar{x},{\bm{\theta}})\mathrel{\mathop{:}}=\textstyle-
f_{{\bm{\theta}}}(x^{\prime})-\frac{1}{2\epsilon_{\sf
NN}}\left\|x^{\prime}-\bar{x}\right\|^{2},$ (26)

to get $z\equiv(x,\bar{y})\sim{\cal D}({\bm{\theta}})$. In practice, we take
the approximation $x\approx\bar{x}-\epsilon_{\sf
NN}{\nabla}_{x}f_{{\bm{\theta}}}(\bar{x})$.

In our experiment, we set $\epsilon_{\sf NN}\in\\{0,10,100\\}$, batch size as
$b=8$. For SGD-GD, we use $\gamma_{t}=\gamma=200/\sqrt{T}$ and for lazy
deployment, we use $\gamma=200/(K\sqrt{T})$ with $T=10^{5}$. The NN encoded in
$f_{\bm{\theta}}(x)$ consists of three fully-connected layers with $\tanh$
activation and a sigmoid output layer, i.e.,

$f_{{\bm{\theta}}}(x)=Sigmoid\big{(}{\bm{\theta}}_{(1)}^{\top}\cdot\tanh({\bm{\theta}}_{(2)}^{\top}\cdot\tanh({\bm{\theta}}_{(3)}^{\top}x))\big{)},$

where ${\bm{\theta}}_{(i)}\mathrel{\mathop{:}}=[w_{(i)};b_{(i)}]$ concatenates
the weight and bias for each layer with $d_{1}=10,d_{2}=50,d_{3}=57$ neurons,
making a total of $d=3421$ parameters for ${\bm{\theta}}$. For the training,
we initialize these parameters as ${\cal N}(0,1)$ for weights and constant
values for the biases.

Fig. 2 (left) compares the SPS measure
$\left\|{\nabla}J({\bm{\theta}}_{t};{\bm{\theta}}_{t})\right\|^{2}$ against
the iteration number $t$ using SGD-GD. As observed, SGD-GD converges to a near
SPS solution and the behavior seems to be insensitive to $\epsilon_{\sf NN}$.
We speculate that this is due to the shift model (26) but would relegate its
study to future work. Fig. 2 (middle & right) compare the greedy and lazy
deployment schemes with $\epsilon_{\sf NN}\in\\{10,10^{4}\\}$ against the
number of samples used. Compared to SGD-GD, lazy deployment performs
relatively better as $\epsilon_{\sf NN}\uparrow$, as seen from the no. of
samples needed to reach
$\|{\nabla}J({\bm{\theta}}_{t};{\bm{\theta}}_{t})\|^{2}=10^{-4}$ in the plots.
This agrees with (23) which shows the dominant term as ${\cal
O}(\epsilon/\sqrt{K})$ and $\epsilon$ is related to $\epsilon_{\sf NN}$.

#### Conclusions

This paper provides the first study on the performative prediction problem
with smooth but possibly non-convex loss. We proposed a stationary
performative stable (SPS) condition which is the counterpart of performative
stable condition used with strongly convex loss. Using the SPS solution
concept, we studied the convergence of greedy deployment and lazy deployment
schemes with SGD. We prove that SGD-GD finds a biased, ${\cal
O}(\epsilon)$-SPS solution, while the lazy deployment scheme finds a bias-free
SPS solution when the lazy deployment epoch is large. As an initial work on
this subclass of problems, our findings can lead to more general analysis on
algorithms under the non-convex performative prediction framework.

## References

  * Björkegren et al. [2020] Daniel Björkegren, Joshua E Blumenstock, and Samsun Knight.  Manipulation-proof machine learning.  _arXiv preprint arXiv:2004.03865_ , 2020.
  * Brown et al. [2022] Gavin Brown, Shlomi Hod, and Iden Kalemaj.  Performative prediction in a stateful world.  In _International conference on artificial intelligence and statistics_ , pages 6045–6061. PMLR, 2022.
  * Cai et al. [2015] Yang Cai, Constantinos Daskalakis, and Christos Papadimitriou.  Optimum statistical estimation with strategic data sources.  In _Conference on Learning Theory_ , pages 280–296. PMLR, 2015.
  * Cutler et al. [2023] Joshua Cutler, Dmitriy Drusvyatskiy, and Zaid Harchaoui.  Stochastic optimization under distributional drift.  _Journal of machine learning research_ , 24(147):1–56, 2023.
  * Dalvi et al. [2004] Nilesh Dalvi, Pedro Domingos, Mausam, Sumit Sanghai, and Deepak Verma.  Adversarial classification.  In _Proceedings of the tenth ACM SIGKDD international conference on Knowledge discovery and data mining_ , pages 99–108, 2004.
  * Dieuleveut et al. [2023] Aymeric Dieuleveut, Gersende Fort, Eric Moulines, and Hoi-To Wai.  Stochastic approximation beyond gradient for signal processing and machine learning.  _IEEE Transactions on Signal Processing_ , 2023.
  * Drusvyatskiy and Xiao [2023] Dmitriy Drusvyatskiy and Lin Xiao.  Stochastic optimization with decision-dependent distributions.  _Mathematics of Operations Research_ , 48(2):954–998, 2023.
  * Ertekin et al. [2010] Seyda Ertekin, Leon Bottou, and C Lee Giles.  Nonconvex online support vector machines.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 33(2):368–381, 2010.
  * Ghadimi and Lan [2013] Saeed Ghadimi and Guanghui Lan.  Stochastic first-and zeroth-order methods for nonconvex stochastic programming.  _SIAM journal on optimization_ , 23(4):2341–2368, 2013.
  * Gibbs and Su [2002] Alison L Gibbs and Francis Edward Su.  On choosing and bounding probability metrics.  _International statistical review_ , 70(3):419–435, 2002.
  * Hardt et al. [2016] Moritz Hardt, Nimrod Megiddo, Christos Papadimitriou, and Mary Wootters.  Strategic classification.  In _Proceedings of the 2016 ACM conference on innovations in theoretical computer science_ , pages 111–122, 2016.
  * Hopkins et al. [1999] Mark Hopkins, Erik Reeber, George Forman, and Jaap Suermondt.  Spambase.  UCI Machine Learning Repository, 1999.  DOI: https://doi.org/10.24432/C53G6X.
  * Izzo et al. [2021] Zachary Izzo, Lexing Ying, and James 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] Meena Jagadeesan, Tijana Zrnic, and Celestine Mendler-Dünner.  Regret minimization with performative feedback.  In _International Conference on Machine Learning_ , pages 9760–9785. PMLR, 2022.
  * Li and Wai [2022] Qiang Li and Hoi-To Wai.  State dependent performative prediction with stochastic approximation.  In _International Conference on Artificial Intelligence and Statistics_ , pages 3164–3186. PMLR, 2022.
  * Li et al. [2022] Qiang Li, Chung-Yiu Yau, and Hoi-To Wai.  Multi-agent performative prediction with greedy deployment and consensus seeking agents.  _Advances in Neural Information Processing Systems_ , 35:38449–38460, 2022.
  * Liu et al. [2023] Haitong Liu, Qiang Li, and Hoi-To Wai.  Two-timescale derivative free optimization for performative prediction with markovian data.  _ArXiv preprint arXiv:2310.05792_ , 2023.
  * Mendler-Dünner et al. [2020] Celestine Mendler-Dünner, Juan Perdomo, Tijana Zrnic, and Moritz Hardt.  Stochastic optimization for performative prediction.  _Advances in Neural Information Processing Systems_ , 33:4929–4939, 2020.
  * Miller et al. [2021] John P Miller, Juan C Perdomo, and Tijana Zrnic.  Outside the echo chamber: Optimizing the performative risk.  In _International Conference on Machine Learning_ , pages 7710–7720. PMLR, 2021.
  * Mofakhami et al. [2023] Mehrnaz Mofakhami, Ioannis Mitliagkas, and Gauthier Gidel.  Performative prediction with neural networks.  In _International Conference on Artificial Intelligence and Statistics_ , pages 11079–11093. PMLR, 2023.
  * Narang et al. [2023] Adhyyan Narang, Evan Faulkner, Dmitriy Drusvyatskiy, Maryam Fazel, and Lillian J Ratliff.  Multiplayer performative prediction: Learning in decision-dependent games.  _Journal of Machine Learning Research_ , 24(202):1–56, 2023.
  * Perdomo et al. [2020] Juan Perdomo, Tijana Zrnic, Celestine Mendler-Dünner, and Moritz Hardt.  Performative prediction.  In _International Conference on Machine Learning_ , pages 7599–7609. PMLR, 2020.
  * Piliouras and Yu [2023] Georgios Piliouras and Fang-Yi Yu.  Multi-agent performative prediction: From global stability and optimality to chaos.  In _Proceedings of the 24th ACM Conference on Economics and Computation_ , pages 1047–1074, 2023.
  * Quiñonero-Candela et al. [2022] Joaquin Quiñonero-Candela, Masashi Sugiyama, Anton Schwaighofer, and Neil D Lawrence.  _Dataset shift in machine learning_.  Mit Press, 2022.
  * Ray et al. [2022] Mitas Ray, Lillian J Ratliff, Dmitriy Drusvyatskiy, and Maryam Fazel.  Decision-dependent risk minimization in geometrically decaying dynamic environments.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , pages 8081–8088, 2022.
  * Roy et al. [2022] Abhishek Roy, Krishnakumar Balasubramanian, and Saeed Ghadimi.  Constrained stochastic nonconvex optimization with state-dependent markov data.  In _Advances in Neural Information Processing Systems_ , volume 35, pages 23256–23270, 2022.
  * Wang et al. [2023] Xiaolu Wang, Chung-Yiu Yau, and Hoi-To Wai.  Network effects in performative prediction games.  In _International Conference on Machine Learning_ , pages 36514–36540. PMLR, 2023.
  * Wood and Dall’Anese [2023] Killian Wood and Emiliano Dall’Anese.  Stochastic saddle point problems with decision-dependent distributions.  _SIAM Journal on Optimization_ , 33(3):1943–1967, 2023.
  * Zhao [2022] Yulai Zhao.  Optimizing the performative risk under weak convexity assumptions.  _arXiv preprint arXiv:2209.00771_ , 2022.
  * Zhu et al. [2023] Zihan Zhu, Ethan Fang, and Zhuoran Yang.  Online performative gradient descent for learning nash equilibria in decision-dependent games.  _Advances in Neural Information Processing Systems_ , 36, 2023.
  * Zrnic et al. [2021] Tijana Zrnic, Eric Mazumdar, Shankar Sastry, and Michael Jordan.  Who leads and who follows in strategic classification?  In _Advances in Neural Information Processing Systems_ , volume 34, pages 15257–15269, 2021.

##  Appendix A Proof of Lemma 3.1

###### Proof.

Under A1, for any fixed $z\in{\sf Z}$, we have that

$\displaystyle\ell({\bm{\theta}}_{t+1};z)$
$\displaystyle\leq\ell({\bm{\theta}}_{t};z)+\left\langle{\nabla}\ell({\bm{\theta}}_{t};z)\,|\,{\bm{\theta}}_{t+1}-{\bm{\theta}}_{t}\right\rangle+\frac{L}{2}\left\|{\bm{\theta}}_{t+1}-{\bm{\theta}}_{t}\right\|^{2}$
$\displaystyle\leq\ell({\bm{\theta}}_{t};z)-\gamma_{t+1}\left\langle{\nabla}\ell({\bm{\theta}}_{t};z)\,|\,{\nabla}\ell({\bm{\theta}}_{t};Z_{t+1})\right\rangle+\frac{L\gamma_{t+1}^{2}}{2}\left\|{\nabla}\ell({\bm{\theta}}_{t};Z_{t+1})\right\|^{2},$

where the second inequality is due to the update rule of (2) as we recall that
$Z_{t+1}\sim{\cal D}({\bm{\theta}}_{t})$. Taking integration on $z$ with
weights given by the p.d.f. of ${\cal D}({\bm{\theta}}_{t})$, i.e.,
$\int(\cdot)p_{{\bm{\theta}}_{t}}(z){\sf d}z$, on both sides of above
inequality leads to

$\displaystyle J({\bm{\theta}}_{t+1};{\bm{\theta}}_{t})\leq
J({\bm{\theta}}_{t};{\bm{\theta}}_{t})-\gamma_{t+1}\left\langle{\nabla}J({\bm{\theta}}_{t};{\bm{\theta}}_{t})\,|\,{\nabla}\ell({\bm{\theta}}_{t};Z_{t+1})\right\rangle+\frac{L}{2}\gamma_{t+1}^{2}\left\|{\nabla}\ell({\bm{\theta}}_{t};Z_{t+1})\right\|^{2}.$

As
$\mathbb{E}_{t}[{\nabla}\ell({\bm{\theta}}_{t};Z_{t+1})]={\nabla}J({\bm{\theta}}_{t};{\bm{\theta}}_{t})$,
taking the conditional expectation $\mathbb{E}_{t}[\cdot]$ on both sides yield

$\displaystyle\mathbb{E}_{t}\left[J({\bm{\theta}}_{t+1};{\bm{\theta}}_{t})\right]$
$\displaystyle\leq
J({\bm{\theta}}_{t};{\bm{\theta}}_{t})-\gamma_{t+1}\left\|{\nabla}J({\bm{\theta}}_{t};{\bm{\theta}}_{t})\right\|^{2}+\frac{L}{2}\gamma_{t+1}^{2}\mathbb{E}_{t}\left[\left\|{\nabla}\ell({\bm{\theta}}_{t};Z_{t+1})\right\|^{2}\right]$
(27)
$\displaystyle\overset{(a)}{=}J({\bm{\theta}}_{t};{\bm{\theta}}_{t})-\gamma_{t+1}\left\|{\nabla}J({\bm{\theta}}_{t};{\bm{\theta}}_{t})\right\|^{2}$
$\displaystyle\quad+\frac{L}{2}\gamma_{t+1}^{2}\left(\mathbb{E}_{t}\left\|{\nabla}\ell({\bm{\theta}}_{t};Z_{t+1})-{\nabla}J({\bm{\theta}}_{t};{\bm{\theta}}_{t})\right\|^{2}+\left\|{\nabla}J({\bm{\theta}}_{t};{\bm{\theta}}_{t})\right\|^{2}\right),$
$\displaystyle\overset{(b)}{\leq}J({\bm{\theta}}_{t};{\bm{\theta}}_{t})-\gamma_{t+1}\left\|{\nabla}J({\bm{\theta}}_{t};{\bm{\theta}}_{t})\right\|^{2}+\frac{L}{2}\gamma_{t+1}^{2}\left(\sigma_{0}^{2}+(1+\sigma_{1}^{2})\left\|{\nabla}J({\bm{\theta}};{\bm{\theta}}_{t})\right\|^{2}\right),$

where $(a)$ used A2 and the property:

$\mathbb{E}_{t}\left[\left\|{\nabla}\ell({\bm{\theta}}_{t};Z_{t+1})\right\|^{2}\right]=\mathbb{E}_{t}\left[\left\|{\nabla}\ell({\bm{\theta}}_{t};Z_{t+1})-{\nabla}J({\bm{\theta}}_{t};{\bm{\theta}}_{t})\right\|^{2}\right]+\|{\nabla}J({\bm{\theta}}_{t};{\bm{\theta}}_{t})\|^{2},$

and $(b)$ is due to the variance bound in A2. Rearranging terms in (27) leads
to

$\displaystyle\left(1-\frac{L}{2}(1+\sigma_{1}^{2})\gamma_{t+1}\right)\gamma_{t+1}\left\|{\nabla}J({\bm{\theta}}_{t};{\bm{\theta}}_{t})\right\|^{2}$
$\displaystyle\leq
J({\bm{\theta}}_{t};{\bm{\theta}}_{t})-\mathbb{E}_{t}[J({\bm{\theta}}_{t+1};{\bm{\theta}}_{t})]+\frac{L}{2}\sigma_{0}^{2}\gamma_{t+1}^{2}.$

The step size condition implies
$1-\frac{L}{2}(1+\sigma_{1}^{2})\gamma_{t+1}\geq 1/2$. This concludes the
proof. ∎

##  Appendix B Proof of Lemma 3.2

###### Proof.

Our proof is modified from Lemma 2.1 of [Drusvyatskiy and Xiao, 2023]. By W2,
since $\ell({\bm{\theta}};z)$ is $L_{0}$-Lipchitz in $z$, we have

$\displaystyle|J({\bm{\theta}};{\bm{\theta}}_{1})-J({\bm{\theta}};{\bm{\theta}}_{2})|=|\mathbb{E}_{Z\sim{\cal
D}({\bm{\theta}}_{1})}\ell({\bm{\theta}};Z)-\mathbb{E}_{Z^{\prime}\sim{\cal
D}({\bm{\theta}}_{2})}\ell({\bm{\theta}};Z^{\prime})|\leq L_{0}W_{1}({\cal
D}({\bm{\theta}}_{1}),{\cal D}({\bm{\theta}}_{2})).$

Applying W1 gives

$\displaystyle|J({\bm{\theta}};{\bm{\theta}}_{1})-J({\bm{\theta}};{\bm{\theta}}_{2})|\leq
L_{0}\epsilon\left\|{\bm{\theta}}_{1}-{\bm{\theta}}_{2}\right\|,$

which finishes the proof. ∎

##  Appendix C Proof of Lemma 3.3

###### Proof.

Observe the following chain

$\displaystyle|J({\bm{\theta}};{\bm{\theta}}_{1})-J({\bm{\theta}};{\bm{\theta}}_{2})|$
$\displaystyle=\left|\int\ell({\bm{\theta}};z)\left(p_{{\bm{\theta}}_{1}}(z)-p_{{\bm{\theta}}_{2}}(z)\right){\sf
d}z\right|$
$\displaystyle\leq\int\left|\ell({\bm{\theta}};z)\frac{p_{{\bm{\theta}}_{2}}(z)-p_{{\bm{\theta}}_{1}}(z)}{p_{{\bm{\theta}}_{1}}(z)}\right|p_{{\bm{\theta}}_{1}}(z){\sf
d}z$
$\displaystyle\overset{(a)}{\leq}\left(\int\left\|\ell({\bm{\theta}};z)\right\|^{2}p_{{\bm{\theta}}_{1}}(z){\sf
d}z\right)^{\frac{1}{2}}\cdot\left(\int\left(\frac{p_{{\bm{\theta}}_{1}}(z)-p_{{\bm{\theta}}_{2}}(z)}{p_{{\bm{\theta}}_{1}}(z)}\right)^{2}p_{{\bm{\theta}}_{1}}(z){\sf
d}z\right)^{\frac{1}{2}}$
$\displaystyle\overset{(b)}{\leq}\ell_{max}\epsilon\left\|{\bm{\theta}}_{1}-{\bm{\theta}}_{2}\right\|,$

where the inequality (a) is due to the Cauchy-Schwarz inequality, (b) is due
to the stated assumptions C1, 2. ∎

##  Appendix D Proof of Theorem 1

###### Proof.

We recall from Lemma 3.1 the following relation:

$\displaystyle\frac{\gamma_{t+1}}{2}\|{\nabla}J({\bm{\theta}}_{t};{\bm{\theta}}_{t})\|^{2}\leq\mathbb{E}_{t}[J({\bm{\theta}}_{t};{\bm{\theta}}_{t})-J({\bm{\theta}}_{t+1};{\bm{\theta}}_{t})]+\frac{L}{2}\sigma_{0}^{2}\gamma_{t+1}^{2}.$
(28)

We notice that Lemmas 3.2, 3.3 imply

$\displaystyle|J(\bar{{\bm{\theta}}};{\bm{\theta}})-J(\bar{{\bm{\theta}}};{\bm{\theta}}^{\prime})|\leq\tilde{L}\epsilon\,\|{\bm{\theta}}-{\bm{\theta}}^{\prime}\|,$

where $\tilde{L}=L_{0}$ if W1, 2 hold, or $\tilde{L}=\ell_{max}$ if C1, 2
hold. Subsequently, the first term on the right hand side of (28) can be
bounded by

$\begin{split}\mathbb{E}_{t}[J({\bm{\theta}}_{t};{\bm{\theta}}_{t})-J({\bm{\theta}}_{t+1};{\bm{\theta}}_{t})]&\leq\mathbb{E}_{t}[J({\bm{\theta}}_{t};{\bm{\theta}}_{t})-J({\bm{\theta}}_{t+1};{\bm{\theta}}_{t+1})]+\mathbb{E}_{t}[|J({\bm{\theta}}_{t+1};{\bm{\theta}}_{t+1})-J({\bm{\theta}}_{t+1};{\bm{\theta}}_{t})|]\\\
&\leq\mathbb{E}_{t}[J({\bm{\theta}}_{t};{\bm{\theta}}_{t})-J({\bm{\theta}}_{t+1};{\bm{\theta}}_{t+1})]+\tilde{L}\epsilon\,\mathbb{E}_{t}[\|{\bm{\theta}}_{t+1}-{\bm{\theta}}_{t}\|]\\\
&=\mathbb{E}_{t}[J({\bm{\theta}}_{t};{\bm{\theta}}_{t})-J({\bm{\theta}}_{t+1};{\bm{\theta}}_{t+1})]+\gamma_{t+1}\tilde{L}\epsilon\,\mathbb{E}_{t}[\|{\nabla}\ell({\bm{\theta}}_{t};Z_{t+1})\|].\end{split}$

Notice that

$\displaystyle\begin{aligned}
\gamma_{t+1}\tilde{L}\epsilon\,\mathbb{E}_{t}\left[\left\|{\nabla}\ell({\bm{\theta}}_{t};Z_{t+1})\right\|\right]&\overset{(a)}{\leq}\gamma_{t+1}\tilde{L}\epsilon\sqrt{\mathbb{E}_{t}\left[\left\|{\nabla}\ell({\bm{\theta}}_{t};Z_{t+1})\right\|^{2}\right]}\\\
&\overset{(b)}{\leq}\gamma_{t+1}\tilde{L}\epsilon\left(\sigma_{0}+\sqrt{1+\sigma_{1}^{2}}\left\|{\nabla}J({\bm{\theta}}_{t};{\bm{\theta}}_{t})\right\|\right)\\\
&\overset{(c)}{\leq}\gamma_{t+1}\tilde{L}\epsilon\left(\sigma_{0}+{(1+\sigma_{1}^{2})\tilde{L}\epsilon}+\frac{1}{4\tilde{L}\epsilon}\left\|{\nabla}J({\bm{\theta}}_{t};{\bm{\theta}}_{t})\right\|^{2}\right),\end{aligned}$

where $(a)$ is due to the Cauchy-Schwarz inequality
$\mathbb{E}[\|X\|]\leq\sqrt{\mathbb{E}[\|X\|^{2}]}$, $(b)$ is due to the
chain:

$\begin{split}\mathbb{E}_{t}\left[\left\|{\nabla}\ell({\bm{\theta}}_{t};Z_{t+1})\right\|^{2}\right]&=\|{\nabla}J({\bm{\theta}}_{t};{\bm{\theta}}_{t})\|^{2}+\mathbb{E}_{t}\left[\left\|{\nabla}\ell({\bm{\theta}}_{t};Z_{t+1})-{\nabla}J({\bm{\theta}}_{t};{\bm{\theta}}_{t})\right\|^{2}\right]\\\
&\leq\sigma_{0}^{2}+(1+\sigma_{1}^{2})\|{\nabla}J({\bm{\theta}}_{t};{\bm{\theta}}_{t})\|^{2}\leq\left(\sigma_{0}+\sqrt{1+\sigma_{1}^{2}}\|{\nabla}J({\bm{\theta}}_{t};{\bm{\theta}}_{t})\|\right)^{2}\end{split}$

and $(c)$ is due to the Young’s inequality. Substituting back into (28) gives

$\frac{\gamma_{t+1}}{4}\|{\nabla}J({\bm{\theta}}_{t};{\bm{\theta}}_{t})\|^{2}\leq\mathbb{E}_{t}[J({\bm{\theta}}_{t};{\bm{\theta}}_{t})-J({\bm{\theta}}_{t+1};{\bm{\theta}}_{t+1})]+\gamma_{t+1}\tilde{L}\epsilon\left(\sigma_{0}+{(1+\sigma_{1}^{2})\tilde{L}\epsilon}\right)+\frac{L}{2}\sigma_{0}^{2}\gamma_{t+1}^{2}.$

Notice that taking full expectation and summing both sides of the inequality
from $t=0$ to $T-1$ yields the theorem. ∎

##  Appendix E Proof of Theorem 2

###### Proof.

The first steps of our proof resemble that of Lemma 3.1 and is repeated here
for completeness. Under A1, for any fixed $z\in{\sf Z}$, we have that

$\displaystyle\ell({\bm{\theta}}_{t,k+1};z)$
$\displaystyle\leq\ell({\bm{\theta}}_{t,k};z)+\left\langle{\nabla}\ell({\bm{\theta}}_{t,k};z)\,|\,{\bm{\theta}}_{t,k+1}-{\bm{\theta}}_{t,k}\right\rangle+\frac{L}{2}\left\|{\bm{\theta}}_{t,k+1}-{\bm{\theta}}_{t,k}\right\|^{2}$
$\displaystyle\leq\ell({\bm{\theta}}_{t,k};z)-\gamma\left\langle{\nabla}\ell({\bm{\theta}}_{t,k};z)\,|\,{\nabla}\ell({\bm{\theta}}_{t,k};Z_{t,k+1})\right\rangle+\frac{L\gamma^{2}}{2}\left\|{\nabla}\ell({\bm{\theta}}_{t,k};Z_{t,k+1})\right\|^{2},$

where the second inequality is due to the update rule of (2) as we recall that
$Z_{t+1}\sim{\cal D}({\bm{\theta}}_{t})$. Taking integration on $z$ with
weights given by the p.d.f. of ${\cal D}({\bm{\theta}}_{t})$, i.e.,
$\int(\cdot)p_{{\bm{\theta}}_{t}}(z){\sf d}z$, on both sides of above
inequality leads to

$\displaystyle J({\bm{\theta}}_{t,k+1};{\bm{\theta}}_{t,0})\leq
J({\bm{\theta}}_{t,k};{\bm{\theta}}_{t,0})-\gamma\left\langle{\nabla}J({\bm{\theta}}_{t,k};{\bm{\theta}}_{t,0})\,|\,{\nabla}\ell({\bm{\theta}}_{t,k};Z_{t,k+1})\right\rangle+\frac{L}{2}\gamma^{2}\left\|{\nabla}\ell({\bm{\theta}}_{t,k};Z_{t,k+1})\right\|^{2}.$

As $Z_{t,k+1}\sim{\cal D}({\bm{\theta}}_{t,0})$, we have
$\mathbb{E}_{t,k}[{\nabla}\ell({\bm{\theta}}_{t,k};Z_{t,k+1})]={\nabla}J({\bm{\theta}}_{t,k};{\bm{\theta}}_{t,0})$,
where $\mathbb{E}_{t,k}[\cdot]$ denotes the conditional expectation on the
filtration

${\cal
F}_{t,k}=\sigma(\\{{\bm{\theta}}_{0},{\bm{\theta}}_{0,1},\cdots,{\bm{\theta}}_{0,K},{\bm{\theta}}_{1,1},\cdots,{\bm{\theta}}_{t},{\bm{\theta}}_{t,1},\cdots,{\bm{\theta}}_{t,k}\\}).$

Taking the conditional expectation $\mathbb{E}_{t,k}[\cdot]$ on both sides
yield

$\displaystyle\mathbb{E}_{t,k}\left[J({\bm{\theta}}_{t,k+1};{\bm{\theta}}_{t,0})\right]$
$\displaystyle\\!\leq\\!J({\bm{\theta}}_{t,k};{\bm{\theta}}_{t,0})\\!-\\!\gamma\left\|{\nabla}J({\bm{\theta}}_{t,k};{\bm{\theta}}_{t,0})\right\|^{2}\\!+\\!\frac{L\gamma^{2}}{2}\mathbb{E}_{t,k}\left\|{\nabla}\ell({\bm{\theta}}_{t,k};Z_{t,k+1})\right\|^{2}$
(29)
$\displaystyle\overset{(a)}{=}J({\bm{\theta}}_{t,k};{\bm{\theta}}_{t,0})-\gamma\left\|{\nabla}J({\bm{\theta}}_{t,k};{\bm{\theta}}_{t,0})\right\|^{2}$
$\displaystyle\quad+\frac{L}{2}\gamma^{2}\left(\mathbb{E}_{t,k}\left\|{\nabla}\ell({\bm{\theta}}_{t,k};Z_{t,k+1})-{\nabla}J({\bm{\theta}}_{t,k};{\bm{\theta}}_{t,0})\right\|^{2}+\left\|{\nabla}J({\bm{\theta}}_{t,k};{\bm{\theta}}_{t,0})\right\|^{2}\right),$
$\displaystyle\overset{(b)}{\leq}J({\bm{\theta}}_{t,k};{\bm{\theta}}_{t,0})-\gamma\left\|{\nabla}J({\bm{\theta}}_{t,k};{\bm{\theta}}_{t,0})\right\|^{2}+\frac{L}{2}\gamma^{2}\left(\sigma_{0}^{2}+(1+\sigma_{1}^{2})\left\|{\nabla}J({\bm{\theta}}_{t,k};{\bm{\theta}}_{t,0})\right\|^{2}\right),$

where $(a)$ used A2 and the property:

$\mathbb{E}_{t,k}\left[\left\|{\nabla}\ell({\bm{\theta}}_{t,k};Z_{t,k+1})\right\|^{2}\right]=\mathbb{E}_{t,k}\left[\left\|{\nabla}\ell({\bm{\theta}}_{t,k};Z_{t+1})-{\nabla}J({\bm{\theta}}_{t,k};{\bm{\theta}}_{t,0})\right\|^{2}\right]+\|{\nabla}J({\bm{\theta}}_{t,k};{\bm{\theta}}_{t,0})\|^{2},$

and $(b)$ is due to the variance bound in A2. Rearranging terms in (29),

$\displaystyle\left(1-\frac{L}{2}(1+\sigma_{1}^{2})\gamma\right)\gamma\left\|{\nabla}J({\bm{\theta}}_{t,k};{\bm{\theta}}_{t,0})\right\|^{2}$
$\displaystyle\leq
J({\bm{\theta}}_{t,k};{\bm{\theta}}_{t,0})-\mathbb{E}_{t,k}[J({\bm{\theta}}_{t,k+1};{\bm{\theta}}_{t,0})]+\frac{L}{2}\sigma_{0}^{2}\gamma^{2}.$
(30)

The step size condition implies
$1-\frac{L}{2}(1+\sigma_{1}^{2})\gamma\geq\frac{1}{2}$.

$\displaystyle\frac{\gamma}{2}\left\|{\nabla}J({\bm{\theta}}_{t,k};{\bm{\theta}}_{t,0})\right\|^{2}$
$\displaystyle\leq
J({\bm{\theta}}_{t,k};{\bm{\theta}}_{t,0})-\mathbb{E}_{t,k}[J({\bm{\theta}}_{t,k+1};{\bm{\theta}}_{t,0})]+\frac{L}{2}\sigma_{0}^{2}\gamma^{2}.$

Taking summation on $k$ from 0 to $K-1$ leads to

$\displaystyle\frac{\gamma}{2}\sum_{k=0}^{K-1}\left\|{\nabla}J({\bm{\theta}}_{t,k};{\bm{\theta}}_{t,0})\right\|^{2}$
$\displaystyle\leq
J({\bm{\theta}}_{t,0};{\bm{\theta}}_{t,0})-\mathbb{E}_{t,k}[J({\bm{\theta}}_{t+1,0};{\bm{\theta}}_{t,0})]+\frac{LK}{2}\sigma_{0}^{2}\gamma^{2}.$
(31)

Recall that ${\bm{\theta}}_{t+1,0}={\bm{\theta}}_{t+1}={\bm{\theta}}_{t,K}$.
Subsequently, the first term on the right hand side of (31) can be bounded by

$\begin{split}\mathbb{E}_{t}[J({\bm{\theta}}_{t};{\bm{\theta}}_{t})-J({\bm{\theta}}_{t+1};{\bm{\theta}}_{t})]&\leq\mathbb{E}_{t}[J({\bm{\theta}}_{t};{\bm{\theta}}_{t})-J({\bm{\theta}}_{t+1};{\bm{\theta}}_{t+1})]+\mathbb{E}_{t}[|J({\bm{\theta}}_{t+1};{\bm{\theta}}_{t+1})-J({\bm{\theta}}_{t+1};{\bm{\theta}}_{t})|]\\\
&\leq\mathbb{E}_{t}[J({\bm{\theta}}_{t};{\bm{\theta}}_{t})-J({\bm{\theta}}_{t+1};{\bm{\theta}}_{t+1})]+\tilde{L}\epsilon\,\mathbb{E}_{t}[\|{\bm{\theta}}_{t+1}-{\bm{\theta}}_{t}\|]\\\
&=\mathbb{E}_{t}[J({\bm{\theta}}_{t};{\bm{\theta}}_{t})-J({\bm{\theta}}_{t+1};{\bm{\theta}}_{t+1})]+\gamma\tilde{L}\epsilon\,\mathbb{E}_{t}\left[\left\|\sum_{k=0}^{K-1}{\nabla}\ell({\bm{\theta}}_{t,k};Z_{t,k+1})\right\|\right].\end{split}$

Notice that through a careful use of A2 and the independence between
stochastic gradients, we have

$\displaystyle\begin{aligned}
&\mathbb{E}_{t}\left[\left\|\sum_{k=0}^{K-1}{\nabla}\ell({\bm{\theta}}_{t,k};Z_{t,k+1})\right\|\right]\leq\sqrt{\mathbb{E}_{t}\left[\left\|\sum_{k=0}^{K-1}{\nabla}\ell({\bm{\theta}}_{t,k};Z_{t,k+1})\right\|^{2}\right]}\\\
&\leq\sqrt{2\mathbb{E}_{t}\left[\left\|\sum_{k=0}^{K-1}({\nabla}\ell({\bm{\theta}}_{t,k};Z_{t,k+1})-{\nabla}J({\bm{\theta}}_{t,k};{\bm{\theta}}_{t}))\right\|^{2}\right]+2\mathbb{E}_{t}\left[\left\|\sum_{k=0}^{K-1}{\nabla}J({\bm{\theta}}_{t,k};{\bm{\theta}}_{t})\right\|^{2}\right]}\\\
&=\sqrt{2\sum_{k=0}^{K-1}\mathbb{E}_{t}\left[\left\|{\nabla}\ell({\bm{\theta}}_{t,k};Z_{t,k+1})-{\nabla}J({\bm{\theta}}_{t,k};{\bm{\theta}}_{t})\right\|^{2}\right]+2\mathbb{E}_{t}\left[\left\|\sum_{k=0}^{K-1}{\nabla}J({\bm{\theta}}_{t,k};{\bm{\theta}}_{t})\right\|^{2}\right]}\\\
&\leq\sqrt{2K\sigma_{0}^{2}+2\sigma_{1}^{2}\sum_{k=0}^{K-1}\mathbb{E}_{t}\left[\left\|{\nabla}J({\bm{\theta}}_{t,k};{\bm{\theta}}_{t})\right\|^{2}\right]+2\mathbb{E}_{t}\left[\left\|\sum_{k=0}^{K-1}{\nabla}J({\bm{\theta}}_{t,k};{\bm{\theta}}_{t})\right\|^{2}\right]}\\\
&\leq\sqrt{2K\sigma_{0}^{2}+2(K+\sigma_{1}^{2})\sum_{k=0}^{K-1}\mathbb{E}_{t}\left[\left\|{\nabla}J({\bm{\theta}}_{t,k};{\bm{\theta}}_{t})\right\|^{2}\right]}.\end{aligned}$
(32)

Using $\sqrt{a^{2}+b^{2}}\leq a+b$ for $a,b\geq 0$, we have

$\displaystyle\begin{aligned}
\mathbb{E}_{t}\left[\left\|\sum_{k=0}^{K-1}{\nabla}\ell({\bm{\theta}}_{t,k};Z_{t,k+1})\right\|\right]&\leq\sqrt{2K}\sigma_{0}+\sqrt{2(K+\sigma_{1}^{2})}\sqrt{\sum_{k=0}^{K-1}\mathbb{E}_{t}\left[\left\|{\nabla}J({\bm{\theta}}_{t,k};{\bm{\theta}}_{t})\right\|^{2}\right]}\\\
&\leq\sqrt{2K}\sigma_{0}+2\sqrt{2(K+\sigma_{1}^{2})}\,\tilde{L}\epsilon+\frac{1}{4\tilde{L}\epsilon}\sum_{k=0}^{K-1}\mathbb{E}_{t}\left[\left\|{\nabla}J({\bm{\theta}}_{t,k};{\bm{\theta}}_{t})\right\|^{2}\right],\end{aligned}$
(33)

where the last inequality used the property
$\sqrt{x}\leq\frac{c}{2}+\frac{x}{2c}$ for any $c>0$. Substituting above
results to (31) and taking full expectation on both sides give us

$\displaystyle\frac{\gamma}{4}\sum_{k=0}^{K-1}\mathbb{E}\left\|{\nabla}J({\bm{\theta}}_{t,k},{\bm{\theta}}_{t,0})\right\|^{2}$
$\displaystyle\leq\mathbb{E}\left[J({\bm{\theta}}_{t},{\bm{\theta}}_{t})-J({\bm{\theta}}_{t+1},{\bm{\theta}}_{t+1})\right]$
(34)
$\displaystyle\quad+\gamma\tilde{L}\epsilon\left(\sqrt{2K}\sigma_{0}+2\sqrt{2(K+\sigma_{1}^{2})}\tilde{L}\epsilon\right)+\frac{LK}{2}\sigma_{0}^{2}\gamma^{2}.$

Next, we lower bound the left hand side by observing:

$\displaystyle\sum_{k=0}^{K-1}\mathbb{E}\left\|{\nabla}J({\bm{\theta}}_{t,k};{\bm{\theta}}_{t,0})\right\|^{2}$
$\displaystyle\overset{(a)}{\geq}\sum_{k=0}^{K-1}\mathbb{E}\left[\frac{1}{2}\left\|{\nabla}J({\bm{\theta}}_{t},{\bm{\theta}}_{t})\right\|^{2}-\left\|{\nabla}J({\bm{\theta}}_{t,k};{\bm{\theta}}_{t})-{\nabla}J({\bm{\theta}}_{t},{\bm{\theta}}_{t})\right\|^{2}\right]$
$\displaystyle\overset{(b)}{\geq}\frac{1}{2}K\mathbb{E}\left\|{\nabla}J({\bm{\theta}}_{t},{\bm{\theta}}_{t})\right\|^{2}-L\sum_{k=0}^{K-1}\mathbb{E}\left\|{\bm{\theta}}_{t,k}-{\bm{\theta}}_{t}\right\|^{2}$
$\displaystyle\overset{(c)}{=}\frac{1}{2}K\mathbb{E}\left\|{\nabla}J({\bm{\theta}}_{t},{\bm{\theta}}_{t})\right\|^{2}-L\sum_{k=0}^{K-1}\mathbb{E}\left\|\sum_{\ell=0}^{k-1}\gamma{\nabla}\ell({\bm{\theta}}_{t,\ell};Z_{t,\ell})\right\|^{2}$
$\displaystyle\geq\frac{1}{2}K\mathbb{E}\left\|{\nabla}J({\bm{\theta}}_{t},{\bm{\theta}}_{t})\right\|^{2}-L\gamma^{2}\sum_{k=0}^{K-1}k\sum_{\ell=0}^{k-1}\mathbb{E}\left\|{\nabla}\ell({\bm{\theta}}_{t,\ell};Z_{t,\ell})\right\|^{2}$
$\displaystyle\overset{(d)}{\geq}\frac{1}{2}K\mathbb{E}\left\|{\nabla}J({\bm{\theta}}_{t},{\bm{\theta}}_{t})\right\|^{2}-L\gamma^{2}\sum_{k=0}^{K-1}k\sum_{\ell=0}^{k-1}G^{2}$
$\displaystyle=\frac{1}{2}K\mathbb{E}\left\|{\nabla}J({\bm{\theta}}_{t},{\bm{\theta}}_{t})\right\|^{2}-LG^{2}\gamma^{2}\cdot\frac{K(K-1)(2K-1)}{6}$
$\displaystyle\geq\frac{1}{2}K\mathbb{E}\left\|{\nabla}J({\bm{\theta}}_{t},{\bm{\theta}}_{t})\right\|^{2}-LG^{2}\gamma^{2}\cdot\frac{K^{3}}{3},$

where $(a)$ is due to the fact that
$\left\|a\right\|^{2}\geq\frac{1}{2}\left\|a+b\right\|^{2}-\left\|b\right\|^{2}$,
for any $a,b\in\mathbb{R}^{n}$, $(b)$ is due to A1, $(c)$ is obtained from the
updating rule (22). In $(d)$, we used the additional assumption A3. The last
chain is due to
$\sum_{k=0}^{K-1}k^{2}=\frac{K(K-1)(2K-1)}{6}\leq\frac{K^{3}}{3}$, when $K\geq
1$. Substituting the above lower bound to (34) and rearrange terms lead to

$\displaystyle\frac{\gamma
K}{8}\mathbb{E}\left\|{\nabla}J({\bm{\theta}}_{t},{\bm{\theta}}_{t})\right\|^{2}$
$\displaystyle\leq\mathbb{E}\left[J({\bm{\theta}}_{t},{\bm{\theta}}_{t})-J({\bm{\theta}}_{t+1},{\bm{\theta}}_{t+1})\right]+\frac{1}{12}\gamma^{3}LG^{2}K^{3}$
$\displaystyle\quad+\gamma\tilde{L}\epsilon\left(\sqrt{2K}\sigma_{0}+2\sqrt{2(K+\sigma_{1}^{2})}\tilde{L}\epsilon\right)+\frac{LK}{2}\sigma_{0}^{2}\gamma^{2}.$

Taking summation from $t=0,1,\cdots,T-1$ gives us

$\displaystyle\frac{\gamma
K}{8}\sum_{t=0}^{T-1}\mathbb{E}\left\|{\nabla}J({\bm{\theta}}_{t},{\bm{\theta}}_{t})\right\|^{2}$
$\displaystyle\leq\mathbb{E}\left[J({\bm{\theta}}_{0},{\bm{\theta}}_{0})-J({\bm{\theta}}_{T},{\bm{\theta}}_{T})\right]+\frac{T}{12}\gamma^{3}LG^{2}K^{3}$
$\displaystyle\quad+\gamma
T\tilde{L}\epsilon\left(\sqrt{2K}\sigma_{0}+2\sqrt{2(K+\sigma_{1}^{2})}\tilde{L}\epsilon\right)+\frac{TLK}{2}\sigma_{0}^{2}\gamma^{2}.$

Dividing $\gamma KT/8$ on both sides, we have

$\displaystyle\frac{1}{T}\sum_{t=0}^{T-1}\mathbb{E}\left\|{\nabla}J({\bm{\theta}}_{t},{\bm{\theta}}_{t})\right\|^{2}$
$\displaystyle\leq\frac{8\Delta_{0}}{\gamma
KT}+4L\sigma_{0}^{2}\gamma+\frac{2}{3}\gamma^{2}LG^{2}K^{2}+\frac{8\tilde{L}\epsilon}{K}\left(\sqrt{2K}\sigma_{0}+2\sqrt{2(K+\sigma_{1}^{2})}\tilde{L}\epsilon\right),$

where we recall
$\Delta_{0}\mathrel{\mathop{:}}=J({\bm{\theta}}_{0};{\bm{\theta}}_{0})-\ell_{max}$.
∎

##  Appendix F Additional Numerical Results

This section provides additional details for the numerical experiments that
were omitted due to space limitation.

#### Synthetic Data with Linear Model

Fig. 3 shows the trajectories of training and testing accuracy with the SGD-GD
scheme under different shift parameters, using the same settings as in Fig. 1
(left & right). Note that the testing dataset with 200 samples is generated
from the same procedure described in the main paper with the same ground truth
${\bm{\theta}}^{o}$, but without the label flipping step. In this case,
although increasing the shift leads to a larger risk value
$J({\bm{\theta}}_{t};{\bm{\theta}}_{t})$ and more biased stationary solution
in terms of $\|{\nabla}J({\bm{\theta}}_{t};{\bm{\theta}}_{t})\|^{2}$ [cf. Fig.
1 (left & middle)], the test/train accuracy remain relatively stable
regardless of the shift parameter. We remark that as observed from [Miller et
al., 2021, Fig. 2], increasing the shift parameter $\epsilon$ does not always
lead to a deteriorated or improved model accuracy. Importantly, the effects
can be unpredictable in general, especially when only biased SPS solutions are
guaranteed.

Figure 3: Synthetic Data (_Left_) Training accuracy under different
sensitivity parameter $\epsilon_{L}$. (_Right_) Testing accuracy under
different $\epsilon_{L}$.

Fig. 4 shows the trajectories of loss values
$J({\bm{\theta}}_{t};{\bm{\theta}}_{t})$, training and testing accuracy with
the greedy and lazy deployment scheme using the same settings as in Fig. 1
(right). We observe similar behaviors as indicated in Fig. 3. Moreover, we
notice that although the lazy deployment scheme converges to a less biased SPS
solution than greedy deployment scheme utilizing the same number of samples,
the initial convergence speed is slower. This can be predicted from Theorem 2
as the lazy deployment scheme is simulated with a larger noise variance
$\sigma_{0}$.

Figure 4: Synthetic Data (_Left_) Loss $V({\bm{\theta}})$ against no of sample
accessed. (_Middle_) Training accuracy under different sensitivity parameter
$\epsilon_{L}$. (_Right_) Testing accuracy under different $\epsilon_{L}$.

#### Real Data with Neural Network Model

Similar to the above paragraph, Fig. 5, 6, 7 show the trajectories of
train/test accuracy, for greedy/lazy deployment scheme when
$\epsilon\in\\{10,10^{4}\\}$ for completeness. The figures demonstrate similar
behavior as described in the main paper. Moreover, we observe that the
sensitivity parameter $\epsilon_{\sf NN}$ has a small effect in the
training/testing acccuracies of the trained models.

Figure 5: Real Data with Neural Network (_Left_) Training accuracy under
different sensitivity parameter $\epsilon$. (_Right_) Testing accuracy under
different $\epsilon$.

Figure 6: Real Data with Neural Network (_left_ & _right_) Training accuracy
under different deployment scheme when $\epsilon_{\sf NN}=10$.

Figure 7: Real Data with Neural Network (_left_ & _right_) Training accuracy
under different deployment scheme when $\epsilon_{\sf NN}=10^{4}$.


[a]Ed Bennett

# Status of reproducibility and open science in hep-lat in 2021

###### Abstract

As a fully computational discipline, Lattice Field Theory has the potential to
give results that anyone with sufficient computational resources can
reproduce, going from input parameters to published numbers and plots correct
to the last byte. After briefly motivating and outlining some of the key steps
in making lattice computations reproducible, this contribution presents the
results of a survey of all 1,229 submissions to the hep-lat arXiv in 2021 of
how explicitly reproducible each is. Areas where LFT has historically been
well ahead of the curve are highlighted, as are areas where there are
opportunities to do more.

##  1 Introduction

There has been increasing pressure from many sides in recent years to embrace
the principle of open science, and to ensure that data analyses are
reproducible. As the software used to analyse data has become increasingly
complex, it has become less and less feasible to unambiguously explain in a
traditional journal publication exactly what steps were carried out, to the
point that another researcher could reproduce them. As enabling this is one of
the key aims of academic publication, this trend poses a clear problem, which
research funders, journals, and others are looking to address.

As a purely computational field, Lattice Field Theory (LFT) stands to be more
severely affected than most fields. However, LFT is also one of the very early
adopters of many principles of Open Science. The latter fact, specifically the
ubiquitous use of freely-available preprints for written publications, has
enabled the work presented here: an analysis of all 1,229 submissions to the
arXiv [1] preprint server in the hep-lat category in 2021, to assess where the
field is strong and where there remain opportunities for improvement.

##  2 Terminology

In this work, _reproducible_ is used as defined by the Turing Way project [2]:
that given the same data, another researcher can repeat the same analysis and
obtain the same results. This definition may seem trivial at first glance, as
it seems obvious that the same analysis on the same data should give the same
result; however, it is a prerequisite of the more challenging concepts of
_replicability_ (being able to apply the same analysis to fresh data and
obtain the same results), and _robustness_ (being able to apply different
analyses to the same data and obtain compatible conclusions).

Since the mapping from human languages to programming ones is not one-to-one,
for any non-trivial software then access to the code is an essential part of
reproducibility. Similarly, for data beyond those that can be presented in a
handful of tables, access to the raw data is also needed. Removing the
possibility for human error is another way to improve reproducibility; every
manual step in a process is a place that could be done inconsistently (either
within a work, or between the original work and attempts to reproduce it). The
primary purpose of computers is to automatically perform repeated tasks in a
consistent way, so encoding manual workflows in software significantly reduces
the possibility for human error, and where errors do creep into code, then
they can be inspected and discovered after the fact.

_Open science_ is the movement that all research—including not only
publications, but also physical samples, data, software, and other
outputs—should be accessible to all by default (barring specific ethical or
legal constraints, for example when working with private personal data).
Theoretical particle physics pioneered the use of open-access preprints for
publications, and remains a leading field—in many disciplines preprints are
still seen by many as a novelty. (For example, in a survey of 3759 researchers
[3], 46% of respondents in medical and health sciences had never viewed or
downloaded a preprint, and 63.6% of the same cohort had never authored one.)

Data and software published openly can (and should) be FAIR: specifically,
_findable_ , _accessible_ , _interoperable_ , and _reusable_. This is a
separate consideration to the data being made available publicly—data that
will never leave a research group still benefits from new researchers to the
group being able to easily find it, and conversely it is entirely possible to
share data in un-FAIR ways. Making data and software FAIR increases the uses
that others can make of it, and hence makes it more valuable.

##  3 Survey results

Figure 1: Breakdown of submissions to the hep-lat arXiv in 2021 indicating
whether they had hep-lat as their primary arXiv or were crosslists, and
whether they presented new numerical results.

The 1,229 submissions to the hep-lat arXiv (including those cross-listed into
the category) were surveyed for some aspects affecting the reproducibility of
the work done. The questions were answered based on a combination of full-text
search for relevant keywords and a brief skim read of the contents. The
complete survey data, and the analysis scripts used to prepare the plots
presented here, are released separately [4]. The analysis was performed using
Python [5], pandas [6], Matplotlib [7], and seaborn [8], via a Jupyter
Notebook [9].

Following the principle of findability, and that a paper should include
sufficient information to be able to reproduce it without hunting elsewhere,
information not included in the preprint was not considered when completing
the survey. This excluded information added in published manuscripts but not
updated on the arXiv, and also excluded information that could be found by
digging through a collaboration’s website.

Figure 1 shows an initial breakdown of the submissions; a small majority had
hep-lat as their primary category, with the remainder being cross-lists. A
substantial majority of both cases—1013 in all—included numerical results
(e.g. a plot that was not merely schematic, or numbers with uncertainties not
quoted from another publication); it is this subset that will be the focus of
the remainder of the analysis.

Figure 2: Breakdown of whether submissions to the hep-lat arXiv in 2021 that
included new numerical results acknowledged an HPC resource (left panel) or
specified any of the software used to generate their results (right panel).

Since one step towards reproducibility is to specify the software used, an
initial question to ask is what proportion of submissions do this. As a
baseline, this is compared with the proportion of submissions that acknowledge
the use of a high-performance computing (HPC) facility. These are of similar
difficulty—adding one or two sentences—and fulfil a similar purpose, in
demonstrating the utility of the software/facility and ensuring that its
development or activities continue to be funded. Figure 2 presents the
breakdown of submissions for these questions. The vast majority of submissions
presenting new numerical results with hep-lat as the primary arXiv (referred
to hereafter as _hep-lat numerical submissions_) acknowledge the use of an HPC
facility. However, fewer than half of hep-lat numerical submissions specify or
acknowledge any of the software used. This is potentially the simplest change
for authors to make; while citing software isn’t sufficient to ensure
reproducibility if custom code has been written, it significantly reduces the
barrier while also recognising software developers’ contributions to the
research.

The survey takes a relatively reductive view towards categorising the work
done in LFT: it is assumed that a typical piece of work may generate field
configurations using a Monte Carlo-like algorithm, may compute observables on
such configurations (generated as part of the same work or otherwise), and
will do some kind of statistical analysis on data resulting from one or both
of the previous two steps. (Alternative approaches including tensor networks
and quantum simulation have their software effort considered in the “analysis”
category; their count is sufficiently small that this is not considered to
introduce significant bias.)

0px

Figure 3: Breakdown of whether submissions to the hep-lat arXiv in 2021 that
generate new field configurations (left panel) or measure observables on field
configurations (right panel) specify the software that was used to do this.

Generation of configurations is the most computationally expensive activity in
LFT, and the associated software suites have had corresponding amounts of
effort put into ensuring they are robust and efficient. Steps taken to ensure
the reproducibiliy of generation include having deterministic seeded random
number generators, and having detailed log files and storing metadata in
configurations recording what version of a piece of software was used and on
what machine. Computation of observables from these ensembles is typically the
next most expensive activity, and frequently uses some of the same software
infrastructure.

Figure 3 shows the proportion of submissions that specify the software used
for each of these activities. Less than one in three of the 262 submissions
generating new configurations specified the software, and a slightly larger
proportion of the 462 submissions measuring observables on configurations
did.111In each case, some of the software specified was toolkits such as Grid
[10] and Chroma [11], where significant amounts of additional code (e.g. XML
or C++) is needed to perform the work being reported; this additional code was
not typically found to be shared or specified. Given that due to differences
in normalisations and other undocumented conventions, running two different
codes with the same parameters can give differing results, this is a barrier
to reproducibility of this effort.

Figure 4: Plot showing the proportion of submissions to the hep-lat arXiv in
2021 making use of existing gauge configurations that acknowledge the
infrastructure used to host and serve them.

Figure 5: Breakdown of the infrastructure acknowledged by those submissions in
Fig. 4 that do acknowledge hosting infrastructure.

Figure 6: Breakdown of submissions to the hep-lat arXiv in 2021 that publish
any of the data they generate.

Figure 7: Breakdown of the repositories used to host data published by the
authors of the submissions indicated in Fig. 6. Unibi is the institutional
repository of the University of Bielefeld.

Where configurations are not generated as part of the work, they must be
obtained from somewhere. This may be a collaborations own storage, or
publicly-available ensembles may be used. In either case storage
infrastructure must be used, and in the latter case a facility for finding and
sharing the configurations. The International Lattice Data Grid (ILDG)
provides schemas and specifications for such services, which have been
deployed in a number of Regional Grids (RGs). This was a very early example of
FAIR data sharing not only in LFT but in all of academia; it in fact pre-dated
the coining of the “FAIR” acronym itself. Being the first in the space has
however meant that tooling development in the surrounding world has moved on,
and left the ILDG tooling unmaintained. Only 14 of 230 submissions using
existing configurations acknowledge the hosting and sharing infrastructure
used (Fig. 4), all of which acknowledge the Japan Lattice Data Grid (Fig. 5),
reflecting the dormant state of the ILDG project and the RGs. Working groups
have now however resumed work on ILDG to modernise it, with the support of new
funding, as discussed in ILDG’s joint contribution [12].

Publications can also share data other than field configurations. Typically
being smaller in size, these require less infrastructure so can be shared with
generally-available data repositories such as CERN Zenodo. An advantage of
using such a repository over using a regular code or web hosting service is
that they will typically provide a persistent identifier (PID) such as a
Digital Object Identifier (DOI), which gives a reference that will not change
over time (e.g. when a researcher changes their GitHub username, or
institution), and a commitment to remain available in the moderately long term
(unlike services such as GitLab, which recently announced older repositories
would be made unavailable after a period of inactivity, with only a few weeks’
notice222This decision was later overturned due to public outcry, but there is
no barrier to the situation occurring again.). Sharing data in this format,
rather than relying on tables in a PDF file, makes collating data from
multiple sources and re-using it in other work significantly easier, and
removes the significant likelihood of errors occuring when transcribing data.
Figure 6 illustrates that this is not yet common practice in LFT, while Fig. 7
shows a significant reliance on code hosting services like GitHub and GitLab
that do not provide persistent identifiers or a guarantee of longevity. The
phrase “data available on request” is commonly used in some areas of science
where sharing of data is required by publisher policy but authors have not
prepared it; in some cases this has led to significant delays in access to
data as requests are ignored [13]. This phrase has yet to catch on in LFT.

Figure 8: Breakdown of submissions to the hep-lat arXiv in 2021 that generate
new numerical results that specify any of the software used for the analysis
of the data presented (left panel), and those that publish the full or partial
software workflow for this analysis (right panel). Figure 9: Breakdown of the
top ten specific pieces of analysis software referred to by submissions to the
hep-lat arXiv in 2021 that acknowledged such software.

The final part of the survey considers the data analysis of LFT data, i.e. the
process that takes in the data produced from configurations on HPC facilities
and outputs the plots and tables shown in the submissions. Many areas of
experimental science have reproducibility efforts, and these have experimental
facilities in place of deterministically-generated field configurations, and
take measurements on these facilities that are not expected to give machine-
precision identical results when repeated; the reproducibility effort focuses
on ensuring that the data from experimental measurements always gives the same
numerical results once analysed. As discussed above, encoding these procedures
into software and sharing that software is the only feasible way to do this
unambiguously.

Figure 8 illustrates that fewer than one in six hep-lat numerical submissions
specified any of the software used. The most popular tools referred to are
shown in Fig. 9; most are programming languages or frameworks rather than
computation-specific tools. Also shown in Fig. 8 is that only two of these
submissions included the workflow allowing another researcher to fully
reproduce the analysis performed. This represents a significant opportunity
for growth. A parallel survey, the initial results of which are reported in a
separate contribution [14] indicates that significantly more work is enabled
by workflows that are at least partially automated; publication of these would
significantly boost the reproducibility of the data analysis phase of LFT
computations.

##  4 Conclusions

While LFT has in many cases been at the forefront of open science, being one
of the first disciplines to embrace open-access preprints and being ahead of
its time with FAIR data for field configurations, there are some areas that
suffer from a “first-mover disadvantage” where work is ongoing to realign with
more modern tooling, and others where there is an opportunity to learn from
the example of other disciplines.

Some low-hanging fruit to improve the reproducibility of publications include
specifying and acknowledging publicly-available software that has enabled the
research presented, and making existing software workflows available (and
citing them in work that uses them). Taking manual processes and automating
them so they can be published is more challenging, in particular doing so in a
way that can be run end-to-end without intervention. There is work to be done
to write tools that will enable this to be done more easily; this work will be
informed by a parallel survey of individual researchers’ practices in this
area, whose initial results are reported in a separate contribution [14].

##  5 Acknowledgements

This work has been funded by the UKRI Science and Technologies Facilities
Council (STFC) Research Software Engineering Fellowship EP/V052489/1. The
author acknowledges the support of the Supercomputing Wales project, which is
part-funded by the European Regional Development Fund (ERDF) via Welsh
Government. The author would like to thank the organisers of the UKLFT Annual
Meeting for the invitation to present on this topic, which inspired this work.
The author also thanks Julian Lenz for reviewing the draft of this
contribution.

## References

  * [1] Cornell University, _arXiv_ , https://arxiv.org.
  * Community [2022] T. T. W. Community, _The Turing Way: A handbook for reproducible, ethical and collaborative research_ (2022), URL https://doi.org/10.5281/zenodo.6909298.
  * Soderberg et al. [2020] C. K. Soderberg, T. M. Errington, and B. A. Nosek, Royal Society Open Science 7, 201520 (2020), https://royalsocietypublishing.org/doi/pdf/10.1098/rsos.201520, URL https://royalsocietypublishing.org/doi/abs/10.1098/rsos.201520.
  * Bennett [2022] E. Bennett, _Survey of reproducibility in hep-lat publications in 2021_ (2022), URL https://doi.org/10.5281/zenodo.7373413.
  * Van Rossum and Drake [2009] G. Van Rossum and F. L. Drake, _Python 3 Reference Manual_ (CreateSpace, Scotts Valley, CA, 2009), ISBN 1441412697.
  * pandas development team [2022] T. pandas development team, _pandas-dev/pandas: Pandas 1.4.2_ (2022), URL https://doi.org/10.5281/zenodo.6408044.
  * Hunter [2007] J. D. Hunter, Computing in Science & Engineering 9, 90 (2007).
  * Waskom [2021] M. L. Waskom, Journal of Open Source Software 6, 3021 (2021), URL https://doi.org/10.21105/joss.03021.
  * Granger and Pérez [2021] B. E. Granger and F. Pérez, Computing in Science & Engineering 23, 7 (2021).
  * Boyle et al. [2016] P. A. Boyle, G. Cossu, A. Yamaguchi, and A. Portelli, in _Proceedings of The 33rd International Symposium on Lattice Field Theory — PoS(LATTICE 2015)_ (2016), vol. 251, p. 023.
  * Edwards and Joo [2005] R. G. Edwards and B. Joo (SciDAC, LHPC, UKQCD), Nucl. Phys. B Proc. Suppl. 140, 832 (2005), hep-lat/0409003.
  * Karsch et al. [2023] F. Karsch, H. Simma, et al., in _To appear in Proceedings of The 39th International Symposium on Lattice Field Theory — PoS(LATTICE 2022)_ (2023).
  * Smart [2018] A. G. Smart, Physics Today (2018), URL https://physicstoday.scitation.org/do/10.1063/PT.6.1.20180822a/full/.
  * Athenodorou et al. [2023] A. Athenodorou, E. Bennett, J. Lenz, and E. Papadopoulou, in _To appear in Proceedings of The 39th International Symposium on Lattice Field Theory — PoS(LATTICE 2022)_ (2023).