SPRINT: Stochastic Performative Prediction With Variance Reduction

## Abstract

Abstract Performative prediction (PP) is an algorithmic framework for optimizing machine learning (ML) models where the model’s deployment affects the distribution of the data it is trained on. Compared to traditional ML with fixed data, designing algorithms in PP converging to a stable point – known as a stationary performative stable (SPS) solution – is more challenging than the counterpart in conventional ML tasks due to the model-induced distribution shifts. While considerable efforts have been made to find SPS solutions using methods such as repeated gradient descent (RGD) and greedy stochastic gradient descent (SGD-GD), most prior studies assumed a strongly convex loss until a recent work established 𝒪 ​ ( 1 / T ) \mathcal{O}(1/\sqrt{T}) convergence of SGD-GD to SPS solutions under smooth, non-convex losses. However, this latest progress is still based on the restricted bounded variance assumption in stochastic gradient estimates and yields convergence bounds with a non-vanishing error neighborhood that scales with the variance. This limitation motivates us to improve convergence rates and reduce error in stochastic optimization for PP, particularly in non-convex settings. Thus, we propose a new algorithm called s tochastic p erformative p r ediction w i th varia n ce reduc t ion ( 𝖲𝖯𝖱𝖨𝖭𝖳 \mathsf{SPRINT} ) and establish its convergence to an SPS solution at a rate of 𝒪 ​ ( 1 / T ) \mathcal{O}(1/T) . Notably, the resulting error neighborhood is independent of the variance of the stochastic gradients. Experiments on multiple real datasets with non-convex models demonstrate that 𝖲𝖯𝖱𝖨𝖭𝖳 \mathsf{SPRINT} outperforms SGD-GD in both convergence rate and stability.

## 1 Introduction

Machine learning (ML) models are typically trained under the assumption that both training and testing data are drawn independently from a static distribution.
However, in many real-world applications where an ML system is used to make predictions about humans, this assumption is often violated, as the deployment of the model can alter the behavior of the population interacting with it, thereby inducing changes in the data distribution that the model is meant to predict.
This phenomenon, known as “performative effects” or “model-induced distribution shift” (Perdomo et al., 2021; Drusvyatskiy and Xiao, 2023), presents distinctive challenges to conventional learning paradigms. For example, in strategic classification scenarios, individuals can manipulate their features (e.g., by modifying their resumes or financial profiles) to obtain more favorable results (Hardt et al., 2015); In digital platforms, users may change their engagement behavior based on their beliefs about how the model operates (Zhang et al., 2019; Chi et al., 2022); In adversarial settings such as spam detection, spammers continuously adapt their tactics to evade newly deployed filters. Understanding and mitigating these performative effects is essential for building robust ML models.

Thus, Perdomo et al. (2021) introduced performative prediction (PP) as the first optimization framework for model-dependent distribution shifts.
In this framework, the data distribution $\mathcal{D}(\boldsymbol{\theta})$ is explicitly modeled as a mapping from the ML model parameters $\boldsymbol{\theta}$ to the space of distributions, and the goal is to minimize risk for this distribution.
Since the data distribution itself depends on the model being optimized, the objective becomes minimizing the performative risk (PR), i.e.,

$$ $\displaystyle\boldsymbol{\theta}^{\text{PO}}=\underset{\boldsymbol{\theta}}{\operatorname{argmin}}~\mathcal{V}(\boldsymbol{\theta})\stackrel{{\scriptstyle\text{def}}}{{=}}\mathbb{E}_{z\sim\mathcal{D}({\boldsymbol{\theta}})}[\ell(z;\boldsymbol{\theta})],$ (1) $$

where $\ell(\boldsymbol{\theta};z)$ is the loss function and $z=(x,y)$ is sampled from the distribution $\mathcal{D}(\boldsymbol{\theta})$. The coupling nature of PP implies that $\boldsymbol{\theta}$ determines both the model and the data distribution, and the minimizer $\boldsymbol{\theta}^{\text{PO}}$ is referred to as the performative optimal (PO) solution. Since the data distribution is itself a function of the variable $\boldsymbol{\theta}$ being optimized, finding $\boldsymbol{\theta}^{\text{PO}}$ is often intractable (Perdomo et al., 2021; Izzo et al., 2021).
To address this challenge, instead of finding PO solution, many existing works focused on performative stable (PS) solution $\boldsymbol{\theta}^{\text{PS}}$ to minimize the decoupled performative risk $\mathcal{J}(\boldsymbol{\theta};\boldsymbol{\theta}^{\text{PS}})$ in the following form:

$$ $\displaystyle\displaystyle\boldsymbol{\theta}^{\text{PS}}=\underset{\boldsymbol{\theta}}{\operatorname{argmin}}~\mathcal{J}(\boldsymbol{\theta};\boldsymbol{\theta}^{\text{PS}})\stackrel{{\scriptstyle\text{def}}}{{=}}\mathbb{E}_{z\sim\mathcal{D}(\boldsymbol{\theta}^{\text{PS}})}[\ell(z;\boldsymbol{\theta})].$ (2) $$

Different from $\mathcal{V}(\boldsymbol{\theta})$ where the data distribution $\mathcal{D}(\boldsymbol{\theta})$ depends on the variable $\boldsymbol{\theta}$ being optimized, $\mathcal{J}(\boldsymbol{\theta};\boldsymbol{\theta}^{\text{PS}})$ decouples the two: the data distribution is induced by a fixed parameter $\boldsymbol{\theta}^{\text{PS}}$, and the the learner only needs to minimize the loss over this fixed distribution.
This enables us to find $\boldsymbol{\theta}^{\text{PS}}$ iteratively using repeated optimization schemes—by repeatedly updating $\boldsymbol{\theta}_{t+1}$ to minimize risk on a fixed distribution $\mathcal{D}(\boldsymbol{\theta}_{t})$.
If the procedure converges, then the minimizer $\boldsymbol{\theta}_{t+1}$ approaches $\boldsymbol{\theta}^{\text{PS}}$. Examples of such approaches commonly used in the literature include repeated risk minimization (Perdomo et al., 2021), repeated gradient descent (Perdomo et al., 2021), and stochastic gradient descent-greedy deploy (SGD-GD) (Mendler-Dünner et al., 2020).
Since the population size can be large in practice, it is not feasible to repeatedly run full gradient descent at each iteration (Perdomo et al., 2021). In this paper, we consider stochastic optimization in PP with a primary focus on the SGD-GD approach (Li and Wai, 2024) since it is more computationally efficient and requires now knowledge of the performative phenomena.
In general, the SGD-GD approach can be written as in the following form:

$$ $\displaystyle\boldsymbol{\theta}_{t+1}=\boldsymbol{\theta}_{t}-\gamma_{t+1}\nabla\ell(z_{t+1};\boldsymbol{\theta}_{t}),~z_{t+1}\sim\mathcal{D}(\boldsymbol{\theta}_{t}).$ (3) $$

It has been shown that Eqn. [3](https://arxiv.org/html/2509.17304v2#S1.E3) could converge to a unique PS solution (Mendler-Dünner et al., 2020), but under two strong assumptions:
(i) the loss function $\ell$ is strongly convex and smooth in $\boldsymbol{\theta}$;
and (ii) the sensitivity of the distribution map $\mathcal{D}$ measured by the Wasserstein-1 distance is upper bounded.
However, the strong convexity requirement significantly limits the practical applicability of PP.
A very recent work (Li and Wai, 2024) addressed this limitation by establishing the first convergence guarantees to a stationary performative stable (SPS) solution (cf. Def. [1.1](https://arxiv.org/html/2509.17304v2#S1.Thmtheorem1)) under the more general setting where $\ell$ is smooth and non-convex.

###### Definition 1.1 ( $\delta$ -Stationary Performative Stable Solution (Li and Wai, 2024 ) ) .

For a given $\delta\geq 0$, the vector $\boldsymbol{\theta}_{\delta\text{-SPS}}\in\mathbb{R}^{d}$ is said to be a $\delta$-stationary performative stable ($\delta$-SPS) solution of Eqn. [2](https://arxiv.org/html/2509.17304v2#S1.E2) if

$$ $\displaystyle\left\|\nabla\mathcal{J}(\boldsymbol{\theta}_{\delta\text{-SPS}};\boldsymbol{\theta}_{\delta\text{-SPS}})\right\|^{2}=\left\|\mathbb{E}_{z\sim\mathcal{D}(\boldsymbol{\theta}_{\delta\text{-SPS}})}\left[\nabla\ell(z;\boldsymbol{\theta}_{\delta\text{-SPS}})\right]\right\|^{2}\leq\delta.$ $$

The parameter $\delta$ in Def. [1.1](https://arxiv.org/html/2509.17304v2#S1.Thmtheorem1) measures the stability of a solution. When $\delta=0$, $\boldsymbol{\theta}_{\delta\text{-SPS}}$ coincides with PS solution $\boldsymbol{\theta}^{\text{PS}}$. As showed by Li and Wai (2024), SGD-GD can converge to $\boldsymbol{\theta}_{\delta\text{-SPS}}$ at a rate of $\mathcal{O}(1/\sqrt{T})$ for non-convex, smooth loss function $\ell$ under an additional assumption of bounded variance in the stochastic gradient estimates. Moreover, the convergence bound includes a non-vanishing error neighborhood that scales with the variance. This naturally motivates the following fundamental question:

We note that, in the literature,
several variance reduction techniques have been developed for the conventional (i.e., non-PP) stochastic gradient descent (SGD) method in ML (e.g., the stochastic average gradient method (SAG) (Roux et al., 2012), the stochastic average gradient augmented (SAGA) method (Defazio et al., 2014), the stochastic variance reduced gradient method (SVRG) (Reddi et al., 2016), and the stochastic path-integrated differential estimator (SPIDER) (Fang et al., 2018)), which achieve better convergence rates.
However, their applicability and effectiveness in the PP settings remain largely unknown.
In light of the increasing importance of PP and its significant gap between the theory and practice, our goal in this paper is to develop new variance-reduced SGD-GD-based algorithms for PP.

As a starting point in the new variance-reduced PP paradigm, we propose a stochastic performative prediction with variance reduction ($\mathsf{SPRINT}$) method, which is the first variance-reduced performative prediction (PP) framework and is inspired by the SVRG approach(Reddi et al., 2016). Our key idea is to divide the updating iterations of SGD-GD into multiple epochs and store a snapshot of the full gradient $\nabla\mathcal{J}$ at the end of each epoch. At each iteration of repeated optimization, the most recent snapshot can reduce the variance of the current stochastic gradient descent with a slight bias.

We note, however, that our proposed $\mathsf{SPRINT}$ approach is far from a straightforward application of SVRG in the PP paradigm.
Unlike conventional stochastic optimization, the unique nature of PP settings introduces two major challenges in algorithm design and theoretical analysis:
(i) The effects change the data distribution and the resulting decoupled performative loss gradient across iterations, thereby rendering the full gradient snapshot taken in the previous iteration a biased estimate of the current stochastic gradient evaluated at the same $\boldsymbol{\theta}$ on a different data distribution; and (ii) The effects introduce new difficulties, absent in the standard SVRG theoretical analysis, in constructing a suitable Lyapunov function and establishing a negative Lyapunov drift, since the change from $\boldsymbol{\theta}_{k}^{s}$ to $\boldsymbol{\theta}_{k+1}^{s}$ incurs extra distribution shifts unseen in the non-PP settings.

The major contribution of this paper is that we overcome the above technical challenges and establish the convergence of our $\mathsf{SPRINT}$ method to an SPS solution with a variance-independent fast convergence rate.
The main results and key contributions of this work are summarized as follows:

- •
We propose the $\mathsf{SPRINT}$ algorithm to reduce the variance of the SGD-GD-Based approach in non-convex PP settings.
Compared to existing works, $\mathsf{SPRINT}$ converges to an SPS solution at an accelerated rate of $\mathcal{O}(1/T)$, and, more importantly, features an error neighborhood that is independent of the stochastic gradient variance.
To our knowledge, all these results are first in the literature.
- •
To establish the aforementioned $\mathcal{O}(1/T)$ convergence rate of our $\mathsf{SPRINT}$ method, we propose a series of new Lyapunov function construction techniques, which are of general and independent interests in the PP literature.
Moreover, we derive the incremental first-order oracle (IFO) complexity (the number of operations of taking a sample and calculating its gradient to achieve a $\mathcal{O}(\delta)$-SPS (Reddi et al., 2016)) of the algorithm.
- •
We conduct extensive numerical experiments to validate the theoretical results on three real-world datasets.
In addition to the credit data (Kaggle, 2012) and MNIST (Deng, 2012), we present additional experimental results for training MLP/CNN models using SGD-GD-Based and $\mathsf{SPRINT}$ approaches on the full CIFAR-10 (Krizhevsky et al., 2009) dataset.
These broad experiments verify and enrich the practical applicability of optimization-based approaches for PP.

The remainder of this paper is organized as follows: In Sec. [2](https://arxiv.org/html/2509.17304v2#S2), we review the related literature and formally present the problem formulation along with a discussion of prior results in Sec. [3](https://arxiv.org/html/2509.17304v2#S3). We then present our $\mathsf{SPRINT}$ algorithm in Sec. [4](https://arxiv.org/html/2509.17304v2#S4) and prove its convergence and IFO complexity in Sec. [5](https://arxiv.org/html/2509.17304v2#S5).
Lastly, we present the numerical results in Sec. [6](https://arxiv.org/html/2509.17304v2#S6) and conclude the paper in Sec. [7](https://arxiv.org/html/2509.17304v2#S7).

###### Definition 1.1 ( δ \delta -Stationary Performative Stable Solution (Li and Wai, 2024 ) ) .

## 2 Related work

In this section, we provide a brief overview on three research areas that are closely related to our work: 1) finding SPS solutions for PP; 2) other related solution metrics in PP; and 3) variance reduction (VR) techniques for SGD-based optimization methods.

1) Finding PS/SPS solutions for PP: Performative prediction (PP) was first formulated as an optimization framework by Perdomo et al. (2021) to handle endogenous data distribution shifts, based on which
an iterative optimization procedure named repeated risk minimization (RRM) to find a performative stable point and also bound the distance between the PS solution and the PO solution in (Perdomo et al., 2021). Perdomo et al. (2021) also proposed a repeated gradient descent (RGD) method, but a full gradient is required in each iteration. In contrast, Mendler-Dünner et al. (2020) designed the first algorithm to find the PS solution under the online setting. Mendler-Dünner et al. (2022) later established the convergence rate of greedy deployment and lazy deployment after each random update under the assumptions of smoothness and strong convexity.
Later, Mofakhami et al. (2023) investigated training neural network in the PP settings, but they assumed the model as $\hat{y}=f_{\boldsymbol{\theta}}$. The loss function is written as $\ell(\hat{y},y)$ and assumed to be strongly convex to $\hat{y}$. The distribution map $D(\hat{y})$ must also be $\epsilon$-sensitive to $\hat{y}$. Most recently, Li and Wai (2024) provided the first convergence results of the PP settings where the loss function in PP is non-convex and no particular solution structure is assumed.
However, this result is based on the bounded variance assumption limitation and suffers from a non-vanishing error neighborhood that scales linearly with the variance. Our work overcomes this limitation by proposing a variance reduction approach to accelerate the convergence and mitigate the error neighborhood.

2) Related solution metrics for PP:
It is worth noting that, besides finding SPS solutions, there also exist other solution metrics for PP.
For example, Miller et al. (2021) tackled PP problems by directly optimizing performance risk (PR) to find the PO solution for a restricted set of the distribution maps.
In contrast, Ray et al. (2022); Liu et al. (2023) utilized derivative-free optimization to find PO solutions. Izzo et al. (2021) also designed an algorithm based on performative gradient descent for finding PO solutions under a convex PR assumption, which, however, is hard to verify even with a strongly convex loss function. Zhu et al. (2023); Zheng et al. (2024) also studied PP under weakly convex conditions. Brown et al. (2022); Li and Wai (2022) focused on state-dependent PP settings and proposed algorithms that converge to PS solutions.
There also exist other works on PP that focused on fairness(Jin et al., 2024), social welfare (Kim and Perdomo, 2022) and privacy (Li et al., 2024) metrics.

3) Variance reduction techniques for SGD-based optimization: In the literature, variance reduction (VR) techniques have been developed to accelerate the convergence of the SGD method
(Roux et al., 2012; Reddi et al., 2016; Nguyen et al., 2017; Defazio et al., 2014; Fang et al., 2018).
For example, stochastic average gradient (SAG) (Roux et al., 2012) and stochastic average gradient augmented (SAGA) (Defazio et al., 2014) maintain an average of the stochastic gradients of all data samples with infrequent updates.
However, these approaches could incur high memory costs in the large dataset regime.
In contrast, the stochastic variance-reduced gradient (SVRG) (Reddi et al., 2016) eliminates the need for storing per-sample stochastic gradients by periodically computing a full gradient at the beginning of each epoch and using it to construct a variance-reduced update.
SVRG enjoys a linear convergence rate with a modest memory requirement under the strong convexity setting.
Subsequently, the stochastic recursive gradient (SARAH) (Nguyen et al., 2017) and the stochastic path-integrated differential estimator (SPIDER)(Fang et al., 2018) methods further improve SVRG by incorporating fresher recursive iterates to achieve optimal convergence rate in terms of both stationary gap and dataset size dependencies.
We note, however, that all these VR techniques were only designed for non-PP settings.
To date, the development of VR techniques for PP remains largely underexplored.

## 3 Preliminaries and state of the art of performative prediction (PP)

Following the convention in the variance reduction literature (e.g., (Reddi et al., 2016; Roux et al., 2012)), we assume the full population contains $n$ samples $\{z_{1},...,z_{n}\}\subset\mathcal{Z}$ in total (^1^11This finite sum setting has been widely used in variance reduction literature (Roux et al., 2012; Reddi et al., 2016; Nguyen et al., 2017) and is practical in PP settings because typical applications (e.g., credit scoring (Perdomo et al., 2021), college admission (Xie et al., 2024)) often involve a finite population. Moreover, our framework can be extended to infinite sum settings as discussed in App. [G](https://arxiv.org/html/2509.17304v2#A7).) and the ML model is parameterized by $\boldsymbol{\theta}\in\Theta\subseteq\mathbb{R}^{d}$.
Then, the expected gradient of the decoupled performative risk among the full population can be computed as: $\nabla\mathcal{J}(\boldsymbol{\theta};\boldsymbol{\theta})=\mathbb{E}_{z\sim\mathcal{D}(\boldsymbol{\theta})}[\nabla\ell(\boldsymbol{\theta};z)]=\frac{1}{n}\sum_{i=1}^{n}\nabla\ell(\boldsymbol{\theta};z_{i})$, where $\nabla\ell(\boldsymbol{\theta};z)$ denotes the gradient of $\ell(\boldsymbol{\theta};z)$ with respect to $\boldsymbol{\theta}$. Unlike most prior works that studied the convergence of SGD-GD (Eqn. [3](https://arxiv.org/html/2509.17304v2#S1.E3)) in performative prediction (PP) under the assumption that the loss function $\ell$ is strongly convex (Mendler-Dünner et al., 2020; Izzo et al., 2021; Mofakhami et al., 2023), we consider non-convex $\ell$ and focus on the problem setup studied by Li and Wai (2024), a recent work that showed the convergence of SGD-GD to an
SPS solution in non-convex PP settings. To familiarize the readers with necessary background and facilitate the subsequent discussions, we first present the technical assumptions and then the theoretical results of (Li and Wai, 2024).
Throughout the paper, we use $\|\cdot\|$ to denote $\ell_{2}$ norm.

###### Assumption 3.1 (Smoothness and Lower Bound of the Gradients) .

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

$$ $\displaystyle\|\nabla\ell(z;\boldsymbol{\theta})-\nabla\ell(z;\boldsymbol{\theta^{\prime}})\|\leq L\|\boldsymbol{\theta}-\boldsymbol{\theta^{\prime}}\|,\quad\forall\boldsymbol{\theta},\boldsymbol{\theta^{\prime}}\in\mathbb{R}^{d}.$ $$

Moreover, there exists a constant $\ell^{\star}>-\infty$ such that $\ell(z;\boldsymbol{\theta})\geq\ell^{\star}$ for all $\boldsymbol{\theta},z$.

###### Definition 3.2 ( $\epsilon$ -Sensitivity) .

Distribution map $\mathcal{D}(\theta)$ is $\epsilon$-sensitive if there exists an $\epsilon\geq 0$ such that, for any $\boldsymbol{\theta},\boldsymbol{\theta^{\prime}}\in\mathbb{R}^{d}$, we have

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

where $W_{1}(\cdot,\cdot)$ is Wasserstein-1 distance (Li and Wai, 2024).

###### Assumption 3.3 (Lipschitzness) .

There exists a constant $L_{0}>0$ such that,

$$ $\displaystyle|\ell(z;\boldsymbol{\theta})-\ell(z^{\prime};\boldsymbol{\theta})|\leq L_{0}\|z-z^{\prime}\|,~~~\forall z,z^{\prime}\in\mathcal{Z},\boldsymbol{\theta}\in\mathbb{R}^{d}.$ $$

###### Assumption 3.4 (Variance of Stochastic Gradient Estimates) .

The stochastic gradient is unbiased. i.e., $\nabla\mathcal{J}(\boldsymbol{\theta_{1}};\boldsymbol{\theta_{2}})=\mathbb{E}_{Z\sim\mathcal{D}(\boldsymbol{\theta_{2}})}[\nabla\ell(Z;\boldsymbol{\theta_{1}})]$.
Also, there exist constants $\sigma_{0},\sigma_{1}\geq 0$, such that

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

Note that smoothness (Assumption [3.1](https://arxiv.org/html/2509.17304v2#S3.Thmtheorem1)) and $\epsilon$-sensitivity defined by Def. [3.2](https://arxiv.org/html/2509.17304v2#S3.Thmtheorem2) are standard in the machine learning optimization literature on PP (e.g., (Perdomo et al., 2021; Mendler-Dünner et al., 2020)), while the Lipschitzness of loss function (Assumption [3.3](https://arxiv.org/html/2509.17304v2#S3.Thmtheorem3)) and the bounded variance assumption (Assumption [3.4](https://arxiv.org/html/2509.17304v2#S3.Thmtheorem4)) are specific to non-convex PP settings but are still common in optimization literature.
Under all of the above assumptions, Li and Wai (2024) proved the following convergence result:

###### Lemma 3.5 (Convergence of SGD-GD (Li and Wai, 2024 ) ) .

Consider iterative updates of SGD-GD given in Eqn. [3](https://arxiv.org/html/2509.17304v2#S1.E3), the following holds for any $T>1$:

$$ $\displaystyle\sum_{t=0}^{T-1}\frac{\gamma_{t+1}}{4}\mathbb{E}\left[\|\nabla J(\boldsymbol{\theta}_{t};\boldsymbol{\theta}_{t})\|^{2}\right]\leq\Delta_{0}+L_{0}\epsilon\left(\sigma_{0}+(1+\sigma_{1}^{2})L_{0}\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},$ (4) $$

where $\gamma_{t}$ is the learning rate at iteration $t$, and $\Delta_{0}=\mathcal{J}(\boldsymbol{\theta_{0}};\boldsymbol{\theta_{0}})-\ell_{*}$ is an upper bound of the initial optimality gap between the initial performative risk and the optimal performative risk $\ell_{*}$.
If one sets $\gamma_{t}=1/\sqrt{T}$ for all $t$, the bound in [4](https://arxiv.org/html/2509.17304v2#S3.E4) can be reduced to:

$$ $\displaystyle\mathbb{E}\left[\left\|\nabla J(\boldsymbol{\theta}_{\tau};\boldsymbol{\theta}_{\tau})\right\|^{2}\right]\leq 4\underbrace{\left(\Delta_{0}+\frac{L}{2}\sigma_{0}^{2}\right)}_{\mathrm{Variance-Dependent}}\cdot\frac{1}{\sqrt{T}}+4L_{0}\epsilon\underbrace{\left(\sigma_{0}+(1+\sigma_{1}^{2})L_{0}\epsilon\right)\}}_{\mathrm{Variance-Dependent}}.$ (5) $$

Lemma [3.5](https://arxiv.org/html/2509.17304v2#S3.Thmtheorem5) shows that SGD-GD converges at a rate of $\mathcal{O}(1/\sqrt{T})$ with an error neighborhood of $\mathcal{O}(4L_{0}\epsilon(\sigma_{0}+(1+\sigma_{1}^{2})L_{0}\epsilon)$. However, this term is directly influenced by the variance of the stochastic gradient estimate. Additionally, the convergence rate is affected by the third term on the right-hand side of Eqn. [4](https://arxiv.org/html/2509.17304v2#S3.E4), which is also variance-dependent. These limitations of the state-of-the-art in PP motivate us to propose a new algorithm for non-convex PP that achieves faster convergence to SPS solutions and develop new variance reduction techniques to enable the SGD-GD-Based PP approach to be variance-independent.

Figure: Algorithm 1 Proposed Algorithm: Stochastic Performative Prediction with Variance Reduction ($\mathsf{SPRINT}$)

###### Assumption 3.1 (Smoothness and Lower Bound of the Gradients) .

###### Definition 3.2 ( ϵ \epsilon -Sensitivity) .

###### Assumption 3.3 (Lipschitzness) .

###### Assumption 3.4 (Variance of Stochastic Gradient Estimates) .

###### Lemma 3.5 (Convergence of SGD-GD (Li and Wai, 2024 ) ) .

## 4 The proposed 𝖲𝖯𝖱𝖨𝖭𝖳 \mathsf{SPRINT} algorithm

To address the aforementioned limitations, we propose a new algorithm called $\mathsf{SPRINT}$ for PP, which is inspired by the stochastic variance reduced gradient (SVRG) technique in the VR literature.
As shown in Alg. [1](https://arxiv.org/html/2509.17304v2#alg1), the total number of $T$ update iterations is divided into $S$ epochs, each consisting of $m$ iterations.
At the end of each epoch $s+1$, the algorithm stores a snapshot of the expected full gradient $\mathcal{J}(\boldsymbol{\widetilde{\theta}}^{s};\boldsymbol{\widetilde{\theta}}^{s})$ (Line 4).
Then, during each update iteration $k$ of epoch $s$, the variance reduction step computes $\boldsymbol{v}_{k}^{s+1}$ (Line 8), where the variance of the current stochastic gradient estimate $\nabla\ell(z_{ik};\boldsymbol{\theta}_{k}^{s+1})$ is offset by the variance of the previous stochastic gradient estimate $\nabla\ell(z_{ik};\widetilde{\boldsymbol{\theta}}^{s})$.
Meanwhile, this step adds the expected gradient $\nabla\mathcal{J}(\widetilde{\boldsymbol{\theta}}^{s};\widetilde{\boldsymbol{\theta}}^{s})$ back to ensure that the overall expectation remains stable.

It is important to note that, unlike variance reduction in the non-PP settings, $\boldsymbol{v}^{s}_{k}$ is not an unbiased estimator of $\nabla\mathcal{J}(\boldsymbol{\theta}^{s}_{k};\boldsymbol{\theta}^{s}_{k})=\mathbb{E}_{z\sim\mathcal{D}(\boldsymbol{\theta}^{s}_{k})}[\nabla\ell(z;\boldsymbol{\theta}^{s}_{k})]$ in the PP setting.
This is because $\mathcal{D}(\widetilde{\boldsymbol{\theta}}^{s-1})\neq\mathcal{D}(\boldsymbol{\theta}^{s}_{k})$.
This introduces new bias compared to the non-PP setting and significantly complicates the convergence analysis of the algorithm, as we detail in Sec. [5](https://arxiv.org/html/2509.17304v2#S5).
Table [1](https://arxiv.org/html/2509.17304v2#S4.T1) summarizes a comparison of our algorithm with existing SGD-GD-Based PP methods.
Notably, only Li and Wai (2024) and our method address the non-convex PP setting.
However, our $\mathsf{SPRINT}$ approach achieves both a faster convergence rate and a (non-vanishing) error neighborhood that is independent of the variance of the stochastic gradient estimate, as illustrated in the next section.

**Table 1: $\mathsf{SPRINT}$ compared to previous literature using SGD-GD. Since the main objective of the paper is to improve the existing SGD-GD-Based PP approaches, we do not compare our results with algorithms that focused on PO solutions (Izzo et al., 2021; Miller et al., 2021) or those using repeated risk minimization (RRM) (Perdomo et al., 2021; Mofakhami et al., 2023). For these methods, we refer readers to Table 1 of Li and Wai (2024) for detailed comparisons.**
| Literature | Loss function $\ell$ | Convergence rate | Error neighborhood | Comments |
| --- | --- | --- | --- | --- |
| Mendler-Dünner et al. (2020) | Strongly convex | $\mathcal{O}(\frac{1}{T})$ | None | No error due to strong convexity |
| Li and Wai (2024) | Nonconvex | $\mathcal{O}(\frac{1}{\sqrt{T}})$ | $\mathcal{O}(\sigma_{0}\epsilon+\sigma_{1}^{2}\epsilon^{2})$ | Variance-dependent error neighborhood |
| This work: $\mathsf{SPRINT}$ | Nonconvex | $\mathcal{O}(\frac{1}{T})$ | $\mathcal{O}(\epsilon^{2}+\epsilon^{4})$ | Variance-independent error neighborhood |

## 5 Convergence analysis of our 𝖲𝖯𝖱𝖨𝖭𝖳 \mathsf{SPRINT} algorithm

In this section, we analyze the convergence of the proposed $\mathsf{SPRINT}$ algorithm.

1) Assumptions: Our theoretical analysis only relies on the standard smoothness of the loss function (Assumption [3.1](https://arxiv.org/html/2509.17304v2#S3.Thmtheorem1)), the $\epsilon$-sensitivity of the distribution map (Def. [3.2](https://arxiv.org/html/2509.17304v2#S3.Thmtheorem2)), and the Lipschitzness of loss function (Assumption [3.3](https://arxiv.org/html/2509.17304v2#S3.Thmtheorem3)). Most notably, unlike Li and Wai (2024), we do not assume bounded variance of the stochastic gradient estimate (Assumption [3.4](https://arxiv.org/html/2509.17304v2#S3.Thmtheorem4)).

2) Overview of the Proofs of the Main Theoretical Results: Motivated by the analysis in (Reddi et al., 2016), which solely focused on SVRG in the non-PP settings, we aim to establish the convergence of $\mathsf{SPRINT}$ by constructing a suitable Lyapunov function and showing that the corresponding function sequence decreases over iterations.
However, it turns out that one cannot trivially extend (Reddi et al., 2016) because the loss $\ell(\boldsymbol{\theta};z)$ in the non-PP setting solely depends on the current model parameter $\boldsymbol{\theta}$.
In stark contrast, the decoupled performative risk $\mathcal{J}(\boldsymbol{\theta};\boldsymbol{\theta}^{\prime})=\mathbb{E}_{z\sim\mathcal{D}(\boldsymbol{\theta}^{\prime})}[\ell(\boldsymbol{\theta};z)]$ in the PP setting depends on two variables: 1) $\boldsymbol{\theta}$, which determines the model prediction, and 2) $\boldsymbol{\theta}^{\prime}$, which controls the data distribution. Because updates to the model parameter $\boldsymbol{\theta}$ simultaneously affect the data distribution, bounding the difference between consecutive terms in the Lyapunov sequence becomes significantly more challenging, which necessitates new proof and analysis techniques.

In what follows, we overcome these above challenges and establish the SPS convergence of $\mathsf{SPRINT}$ in three steps: (i) We construct an intermediate function sequence (i.e., $R_{k}^{s+1}(\boldsymbol{\theta})$ in Lemma [5.1](https://arxiv.org/html/2509.17304v2#S5.Thmtheorem1)), where consecutive terms differ only in the second argument of decoupled performative risk $\mathcal{J}$ and show that $R_{k}^{s+1}(\boldsymbol{\theta}_{k}^{s+1})$ is smaller than $R_{k+1}^{s+1}(\boldsymbol{\theta}_{k}^{s+1})$; (ii) We define the final Lyapunov function ($\widetilde{R}_{k}^{s+1}$ in Lemma [5.3](https://arxiv.org/html/2509.17304v2#S5.Thmtheorem3)) where consecutive terms differ in both arguments of $\mathcal{J}$; (iii) by proving the decreasing nature of $\widetilde{R}_{k}^{s+1}$, we establish the convergence of Algorithm [1](https://arxiv.org/html/2509.17304v2#alg1) in Theorem  [5.4](https://arxiv.org/html/2509.17304v2#S5.Thmtheorem4) and derive the IFO complexity bound in Corollary [5.5](https://arxiv.org/html/2509.17304v2#S5.Thmtheorem5).

Step 1) Construct an intermediate function sequence: In the $s+1$st epoch, we identify an intermediate function evaluated at the iterations $k$ and $k+1$ as:

$$ $\displaystyle R_{k}^{s+1}(\boldsymbol{\theta})$ $\displaystyle\overset{\triangle}{=}\mathbb{E}\left[\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta})+c_{k}\|\boldsymbol{\theta}_{k}^{s+1}-\widetilde{\boldsymbol{\theta}}^{s}\|^{2}\right],$ $\displaystyle R_{k+1}^{s+1}(\boldsymbol{\theta})$ $\displaystyle\overset{\triangle}{=}\mathbb{E}\left[\mathcal{J}(\boldsymbol{\theta}_{k+1}^{s+1};\boldsymbol{\theta})+c_{k+1}\|\boldsymbol{\theta}_{k+1}^{s+1}-\widetilde{\boldsymbol{\theta}}^{s}\|^{2}\right].$ $$

Then, the following Lemma [5.1](https://arxiv.org/html/2509.17304v2#S5.Thmtheorem1) characterizes the relation between $R_{k}^{s+1}(\boldsymbol{\theta}^{s+1}_{k})$ and $R_{k+1}^{s+1}(\boldsymbol{\theta}^{s+1}_{k})$:

###### Lemma 5.1 .

Let $\beta_{k}$ be some positive constant, and $c_{k},c_{k+1}$ are some constants defined in intermediate functions $R_{k}^{s+1},R_{k+1}^{s+1}$ that satisfy the following:

$$ $\displaystyle c_{k}-\frac{1}{2}\gamma_{k}=c_{k+1}(1+\gamma_{k}\beta_{k}+2L_{0}\epsilon\gamma_{k})+(2L^{2}+L_{0}^{2}\epsilon^{2})(L\gamma_{k}^{2}+2c_{k+1}\gamma_{k}^{2})+\frac{\beta_{k}}{2}L_{o}\epsilon\gamma_{k}$ $$

Then, we have

$$ $\displaystyle R_{k+1}^{s+1}(\boldsymbol{\theta}_{k}^{s+1})\leq R_{k}^{s+1}(\boldsymbol{\theta}_{k}^{s+1})-\Gamma_{k}\cdot\mathbb{E}[\|\nabla\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})\|^{2}]-\frac{1}{2}\gamma_{k}\|\boldsymbol{\theta}_{k}^{s+1}-\widetilde{\boldsymbol{\theta}}^{s}\|^{2},$ $$

where $\Gamma_{k}=\left(\gamma_{k}-\frac{c_{k+1}\gamma_{k}+\frac{1}{2}L_{0}\epsilon\gamma_{k}}{\beta_{k}}-2L\gamma_{k}^{2}-4c_{k+1}\gamma_{k}^{2}\right)$.

Note that $R_{k+1}^{s+1}(\boldsymbol{\theta}_{k}^{s+1})$ includes a term $\mathbb{E}[\mathcal{J}(\boldsymbol{\theta}_{k+1}^{s+1};\boldsymbol{\theta}_{k}^{s+1})]$, which represents the expected performative risk when the model parameter is $\boldsymbol{\theta}_{k+1}^{s+1}$ but the data distribution is induced by $\boldsymbol{\theta}_{k}^{s+1}$. This enables us to focus on bounding the influence of $\boldsymbol{\theta}$ through the first argument of $\mathcal{J}$, similar to analyses in the non-PP settings.
However, since the snapshot $\mathcal{J}(\widetilde{\boldsymbol{\theta}}^{s};\widetilde{\boldsymbol{\theta}}^{s})$ corresponds to the full gradient computed under a previous distribution, we still need to bound the influence of distribution shifts.
This implies that we need to bound the influence of distribution shifts from $\mathcal{D}(\widetilde{\boldsymbol{\theta}}^{s})$ to $\mathcal{D}(\boldsymbol{\theta}^{s+1}_{k})$.
To this end, we provide the following auxiliary lemma.

###### Lemma 5.2 ( Li and Wai ( 2024 ) ) .

Under Assumption [3.3](https://arxiv.org/html/2509.17304v2#S3.Thmtheorem3) and $\epsilon$-sensitivity in Def. [3.2](https://arxiv.org/html/2509.17304v2#S3.Thmtheorem2), we have:

$$ $\displaystyle|\mathcal{J}(\boldsymbol{\theta};\boldsymbol{\theta}_{1})-\mathcal{J}(\boldsymbol{\theta};\boldsymbol{\theta}_{2})|\leq L_{0}\epsilon\|\boldsymbol{\theta}_{1}-\boldsymbol{\theta}_{2}\|,~~\forall\boldsymbol{\theta}.$ $$

We refer readers to (Li and Wai, 2024, App. B) for the proof of Lemma [5.2](https://arxiv.org/html/2509.17304v2#S5.Thmtheorem2).

To prove Lemma [5.1](https://arxiv.org/html/2509.17304v2#S5.Thmtheorem1), we first apply the $\beta$-smoothness of $\ell$ on $\mathbb{E}[\mathcal{J}(\boldsymbol{\theta}_{k+1}^{s+1};\boldsymbol{\theta}_{k}^{s+1})]$ and leverage Lemma [5.2](https://arxiv.org/html/2509.17304v2#S5.Thmtheorem2) to bound $\|\boldsymbol{v}_{k}^{s+1}\|$ and $\|\boldsymbol{v}_{k}^{s+1}\|^{2}$. The complete proof of Lemma [5.1](https://arxiv.org/html/2509.17304v2#S5.Thmtheorem1) is provided in App. [B](https://arxiv.org/html/2509.17304v2#A2) due to space limitation.

Step 2) Construct the final Lyapunov function: We then construct the final Lyapunov function: at epoch $s+1$, we define Lyapunov function evaluated at iterations $k$ and $k+1$ as:

$$ $\displaystyle\widetilde{R}_{k}^{s+1}$ $\displaystyle\overset{\triangle}{=}\mathbb{E}[\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})]+c_{k}\|\boldsymbol{\theta}_{k}^{s+1}-\widetilde{\boldsymbol{\theta}}^{s}\|^{2},$ $\displaystyle\widetilde{R}_{k+1}^{s+1}$ $\displaystyle\overset{\triangle}{=}\mathbb{E}[\mathcal{J}(\boldsymbol{\theta}_{k+1}^{s+1};\boldsymbol{\theta}_{k+1}^{s+1})]+c_{k+1}\|\boldsymbol{\theta}_{k+1}^{s+1}-\widetilde{\boldsymbol{\theta}}^{s}\|^{2}.$ $$

The following Lemma [5.3](https://arxiv.org/html/2509.17304v2#S5.Thmtheorem3) characterizes the relation between $\widetilde{R}_{k+1}^{s}$ and $\widetilde{R}_{k+1}^{s+1}$:

###### Lemma 5.3 .

Let $c_{k},c_{k+1},\beta_{k},\Gamma_{k}$ be defined the same as in Lemma [5.1](https://arxiv.org/html/2509.17304v2#S5.Thmtheorem1), we have:

$$ $\displaystyle\widetilde{R}_{k+1}^{s+1}\leq\widetilde{R}_{k}^{s+1}+(L_{0}^{4}\epsilon^{4}+2L_{0}^{2}\epsilon^{2}L^{2}+2L_{0}^{2}\epsilon^{2})\gamma_{k}-\left(\Gamma_{k}-\frac{\gamma_{k}}{4}\right)\cdot\mathbb{E}\left[\|\nabla\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})\|^{2}\right].$ (6) $$

Lemma [5.3](https://arxiv.org/html/2509.17304v2#S5.Thmtheorem3) can be proved by leveraging Lemmas [5.1](https://arxiv.org/html/2509.17304v2#S5.Thmtheorem1) and [5.2](https://arxiv.org/html/2509.17304v2#S5.Thmtheorem2).
The proof also establishes a bound on the difference between $\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})$ and $\mathcal{J}(\boldsymbol{\theta}_{k+1}^{s+1};\boldsymbol{\theta}_{k+1}^{s+1})$. The full proof is provided in App. [C](https://arxiv.org/html/2509.17304v2#A3).

Step 3) Establish the convergence of the $\mathsf{SPRINT}$ algorithm: Building on the results of Lemma [5.3](https://arxiv.org/html/2509.17304v2#S5.Thmtheorem3), we are now ready to prove the main convergence theorem:

###### Theorem 5.4 (Convergence to an SPS point) .

Consider $\mathsf{SPRINT}$ in Algorithm [1](https://arxiv.org/html/2509.17304v2#alg1) with $T$ rounds and $S=\lceil T/m\rceil$ epochs, where each epoch consists of $m$ iterations. For the final iteration of each epoch, let the constant $c_{m}$ in the Lyapunov function be set to 0. Then we have:

$$ $\displaystyle\frac{1}{T}\sum_{s=0}^{S-1}\sum_{k=0}^{m-1}\mathbb{E}\|\nabla\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})\|^{2}\leq\frac{\Delta_{0}}{T\Gamma}+\Delta_{1},$ (7) $$

where terms $\Delta_{0}$ and $\Delta_{1}$ are defined as

$$ $\displaystyle\Delta_{0}\overset{\triangle}{=}\mathcal{J}(\boldsymbol{\theta}_{0}^{0};\boldsymbol{\theta}_{0}^{0})-\ell_{*};~~~~~\Delta_{1}\overset{\triangle}{=}\frac{(L_{0}^{4}\epsilon^{4}+2L_{0}^{2}\epsilon^{2}L^{2}+2L_{0}^{2}\epsilon^{2})\sum_{s=0}^{S}\sum_{k=0}^{m}\gamma_{k}}{T\Gamma}$ $$

and $\Gamma>0$ is a lower bound of $\Gamma_{k}-\frac{\gamma_{k}}{4}$ and is guaranteed to exist as a constant.
Thus, Algorithm [1](https://arxiv.org/html/2509.17304v2#alg1) achieves a convergence rate of $\mathcal{O}(1/T)$ and an variance-independent error neighborhood $\Delta_{1}$ of $\mathcal{O}(\epsilon^{2}+\epsilon^{4})$.

Eqn. [7](https://arxiv.org/html/2509.17304v2#S5.E7) is established by summing the terms in Lemma [5.3](https://arxiv.org/html/2509.17304v2#S5.Thmtheorem3), and then lower bounding $\Gamma_{k}-\frac{1}{4}\gamma_{k}$ to guarantee the convergence rate $\mathcal{O}(1/T)$. The complete proof can be found in App. [D](https://arxiv.org/html/2509.17304v2#A4). Thm. [5.4](https://arxiv.org/html/2509.17304v2#S5.Thmtheorem4) improves upon previous results in three key aspects: (i) the $\mathcal{O}(1/T)$ convergence rate is faster than that of Li and Wai (2024);
(ii) we eliminate Assumption [3.4](https://arxiv.org/html/2509.17304v2#S3.Thmtheorem4), and the error neighborhood $\Delta_{1}$ is no longer dependent on the variance of the stochastic gradient; and (iii) the error neighborhood is $\mathcal{O}(\epsilon^{2}+\epsilon^{4})$, which is guaranteed to be smaller than that of the SGD-GD-Based PP algorithm in (Li and Wai, 2024) when the distribution map is not sensitive ($\epsilon<1$). Next, we derive the IFO complexity as defined in Reddi et al. (2016) (i.e., the number of operations required to sample and compute the gradient to achieve a $\mathcal{O}(\delta)$-SPS) for Algorithm [1](https://arxiv.org/html/2509.17304v2#alg1). Compared to the IFO Complexity of SGD-GD ($\mathcal{O}(\frac{(\Delta_{0}+\frac{L}{2}\sigma_{0}^{2})^{2}}{\delta^{2}})$), $\mathsf{SPRINT}$ has a clear advantage when the population size $n$ is not too large or when the variance parameter $\sigma_{0}^{2}$ is large or even does not exist. We provide more details in App. [H](https://arxiv.org/html/2509.17304v2#A8).

###### Corollary 5.5 (IFO Complexity of $\mathsf{SPRINT}$ ) .

Assume a finite distribution setting where the population consists of $n$ samples, and $\alpha\in(0,1)$. For the algorithm with fixed parameters at each round, i.e., $\gamma_{k}=\gamma=\mathcal{O}(\frac{1}{n^{\alpha}}),\beta_{k}=\beta=\mathcal{O}(n^{\frac{\alpha}{2}}),m=\mathcal{O}(n^{\frac{\alpha}{2}})$, we have $\Gamma=\mathcal{O}(\frac{1}{n^{\alpha}})$, and the IFO complexity to achieve $\mathcal{O}(\delta)$-SPS, in addition to the error neighborhood, is $\mathcal{O}(\frac{(n^{\alpha}+n^{1+\frac{\alpha}{2}})\Delta_{0}}{\delta})$.

###### Lemma 5.1 .

###### Lemma 5.2 ( Li and Wai ( 2024 ) ) .

###### Lemma 5.3 .

###### Theorem 5.4 (Convergence to an SPS point) .

###### Corollary 5.5 (IFO Complexity of 𝖲𝖯𝖱𝖨𝖭𝖳 \mathsf{SPRINT} ) .

## 6 Experiments

In this section, we conduct comprehensive experimental studies to verify the effectiveness of $\mathsf{SPRINT}$ in practice.
We compare the convergence behavior of the training loss and training accuracy of $\mathsf{SPRINT}$ with SGD-GD (Li and Wai, 2024) on $3$ real-world datasets with non-convex models and study the performative effects.
The results of Credit dataset (Kaggle, 2012) and CIFAR-10 dataset (Krizhevsky et al., 2009) are included in the main paper, while results of MNIST dataset (Deng, 2012) are relegated in App. [F](https://arxiv.org/html/2509.17304v2#A6) due to space limitation.
We also show the magnitude of the squared gradient norm $\|\nabla\mathcal{J}(\widetilde{\boldsymbol{\theta}}^{s};\widetilde{\boldsymbol{\theta}}^{s})\|^{2}$ in App. [F](https://arxiv.org/html/2509.17304v2#A6) to validate Thm. [5.4](https://arxiv.org/html/2509.17304v2#S5.Thmtheorem4).

Notably, besides the credit and MNIST datasets, we are the first to provide results using CNN models on the CIFAR-10 (Krizhevsky et al., 2009) dataset in the PP settings.
As neural network architectures are widely used, our experiments go one step further for the practical applications of PP algorithms.

### 6.1 Experiments on the Credit dataset

We use the Give me some credit data (Kaggle, 2012) with $10$ features $X\in\mathbb{R}^{10}$ to predict creditworthiness $Y\in\{0,1\}$ (Perdomo et al., 2021; Jin et al., 2024).
We preprocess the dataset by selecting a balanced subset of $5,000$ examples and then normalizing features similar to Perdomo et al. (2021). We then assume a credit scoring agency uses a two-layer MLP model to predict $Y$ given the individuals’ features $X$ with cross entropy loss.
To simulate performative effects, we assume individuals are strategic (Perdomo et al., 2021; Hardt et al., 2015), i.e., there is a subset of features $X_{s}\in X$ which individuals can change to $X_{s}^{\prime}$ based on the current model.
Given the current model $f_{\boldsymbol{\theta}}$, individuals with feature $x_{s}$ will best respond to the model in the direction of $\nabla_{x}f_{\boldsymbol{\theta}}(x_{s})$ and change $x_{s}$ to $x_{s}^{\prime}=x_{s}+\alpha\nabla_{x}f_{\boldsymbol{\theta}}(x_{s})$ to improve their scoring (Rosenfeld et al., 2020; Xie and Zhang, 2024).

In this setting, we train the MLP model using $\mathsf{SPRINT}$ and SGD-GD for 40 epochs and illustrate the accuracy and training loss in Fig. [1](https://arxiv.org/html/2509.17304v2#S6.F1).
We simulate different performative effects using of $\alpha\in\{0.01,0.2,0.4\}$.
In Fig. [1](https://arxiv.org/html/2509.17304v2#S6.F1), we show that the $\mathsf{SPRINT}$ algorithm achieves lower training loss and converges faster compared with the SGD-GD in all $3$ cases and the advantages are quite prominent when performative effects are larger (i.e., when $\alpha=0.2$ and $0.4$).

### 6.2 Experiments on the CIFAR-10 dataset

Furthermore, we use the CIFAR-10 dataset (Krizhevsky et al., 2009) to predict 10 possible class labels $Y\in\{0,1,...,9\}$.
Each label corresponds to a common and well-separated category such as dog and truck.
The dataset contains $50,000$ images in total and each class includes $5,000$ samples.
We design a CNN model with two convolutional layers to train the algorithms with a $0.05$ learning rate and cross-entropy loss.
We simulate performative effects motivated by retention dynamics (Hashimoto et al., 2018; Jin et al., 2024; Zhang et al., 2019) where the class distribution will change based on the performance of the current model.
The fraction in class $c$ at iteration $k$ of $p^{k+1}_{c}=\frac{e^{-\alpha\ell^{k}_{c}}}{\sum_{c^{\prime}}e^{-\alpha\ell^{k}_{c^{\prime}}}}$, where $\ell^{k}_{c}$ is the loss expectation of class $c$ in iteration $k$.
This means that the image class with larger loss now will be less likely to appear in the next iteration (Zhang et al., 2019).
We set $\alpha$ at $20,50,80$ and train 400 epochs in total.
The experimental results are shown in Fig. [2](https://arxiv.org/html/2509.17304v2#S6.F2).
Consistent with the results observed on the Credit dataset, $\mathsf{SPRINT}$ also exhibits a faster convergence on the CIFAR-10 dataset.

Figure: (a) $\alpha=0.01$
Refer to caption: x1.png

Figure: (a) $\alpha=20$
Refer to caption: x7.png

## 7 Conclusion and future work

This paper focused on improving the convergence of stochastic optimization methods in non-convex PP settings by adapting a variance reduction method.
We proposed $\mathsf{SPRINT}$ as a new algorithm to achieve better convergence results without the need for the bounded variance assumption on stochastic gradient descent estimates.
We also presented extensive experiments with different non-convex neural network models to compare $\mathsf{SPRINT}$ with SGD-GD, which demonstrated the superior performance of our algorithm $\mathsf{SPRINT}$. However, $\mathsf{SPRINT}$ only converges to a non-vanishing error neighborhood. Future work may include developing algorithms to converge to performative optimal solutions in the non-convex PP settings, or adapting other more advanced variance reduction methods (SARAH or SPIDER) to the PP setting and trying to completely eliminate the non-vanishing error neighborhood.

## References

- Brown et al. [2022]
Gavin Brown, Shlomi Hod, and Iden Kalemaj.
Performative prediction in a stateful world.
In Gustau Camps-Valls, Francisco J. R. Ruiz, and Isabel Valera, editors, *Proceedings of The 25th International Conference on Artificial Intelligence and Statistics*, volume 151 of *Proceedings of Machine Learning Research*, pages 6045–6061. PMLR, 28–30 Mar 2022.
URL [https://proceedings.mlr.press/v151/brown22a.html](https://proceedings.mlr.press/v151/brown22a.html).
- Chi et al. [2022]
Jianfeng Chi, Jian Shen, Xinyi Dai, Weinan Zhang, Yuan Tian, and Han Zhao.
Towards return parity in markov decision processes.
In *International Conference on Artificial Intelligence and Statistics*, pages 1161–1178. PMLR, 2022.
- Defazio et al. [2014]
Aaron Defazio, Francis Bach, and Simon Lacoste-Julien.
Saga: A fast incremental gradient method with support for non-strongly convex composite objectives.
*Advances in neural information processing systems*, 27, 2014.
- Deng [2012]
Li Deng.
The mnist database of handwritten digit images for machine learning research.
*IEEE Signal Processing Magazine*, 29(6):141–142, 2012.
- Drusvyatskiy and Xiao [2023]
Dmitriy Drusvyatskiy and Lin Xiao.
Stochastic optimization with decision-dependent distributions.
*Mathematics of Operations Research*, 48(2):954–998, 2023.
- Fang et al. [2018]
Cong Fang, Chris Junchi Li, Zhouchen Lin, and Tong Zhang.
Spider: Near-optimal non-convex optimization via stochastic path-integrated differential estimator.
*Advances in neural information processing systems*, 31, 2018.
- Gurbuzbalaban et al. [2021]
Mert Gurbuzbalaban, Umut Simsekli, and Lingjiong Zhu.
The heavy-tail phenomenon in sgd.
In *International Conference on Machine Learning*, pages 3964–3975. PMLR, 2021.
- Hardt et al. [2015]
Moritz Hardt, Nimrod Megiddo, Christos Papadimitriou, and Mary Wootters.
Strategic Classification, November 2015.
URL [http://arxiv.org/abs/1506.06980](http://arxiv.org/abs/1506.06980).
Number: arXiv:1506.06980 arXiv:1506.06980 [cs].
- Hashimoto et al. [2018]
Tatsunori Hashimoto, Megha Srivastava, Hongseok Namkoong, and Percy Liang.
Fairness without demographics in repeated loss minimization.
In *International Conference on Machine Learning*, pages 1929–1938. PMLR, 2018.
- Izzo et al. [2021]
Zachary Izzo, Lexing Ying, and James Zou.
How to learn when data reacts to your model: Performative gradient descent.
*CoRR*, 2021.
arXiv: 2102.07698.
- Jin et al. [2024]
Kun Jin, Tian Xie, Yang Liu, and Xueru Zhang.
Addressing polarization and unfairness in performative prediction, 2024.
URL [https://arxiv.org/abs/2406.16756](https://arxiv.org/abs/2406.16756).
- Jothimurugesan et al. [2018]
Ellango Jothimurugesan, Ashraf Tahmasbi, Phillip Gibbons, and Srikanta Tirthapura.
Variance-reduced stochastic gradient descent on streaming data.
*Advances in neural information processing systems*, 31, 2018.
- Kaggle [2012]
Kaggle.
Give me some credit.
[https://www.kaggle.com/c/GiveMeSomeCredit/data](https://www.kaggle.com/c/GiveMeSomeCredit/data), 2012.
- Kim and Perdomo [2022]
Michael P Kim and Juan C Perdomo.
Making decisions under outcome performativity.
*arXiv preprint arXiv:2210.01745*, 2022.
- Krizhevsky et al. [2009]
Alex Krizhevsky, Geoffrey Hinton, et al.
Learning multiple layers of features from tiny images.
2009.
- 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 and Wai [2024]
Qiang Li and Hoi-To Wai.
Stochastic optimization schemes for performative prediction with nonconvex loss.
*arXiv preprint arXiv:2405.17922*, 2024.
- Li et al. [2024]
Qiang Li, Michal Yemini, and Hoi-To Wai.
Clipped sgd algorithms for performative prediction: Tight bounds for clipping bias and remedies.
*arXiv preprint arXiv:2404.10995*, 2024.
- 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. [2022]
Celestine Mendler-Dünner, Frances Ding, and Yixin Wang.
Anticipating performativity by predicting from predictions.
*Advances in Neural Information Processing Systems*, 35:31171–31185, 2022.
- Mendler-Dünner et al. [2020]
Celestine Mendler-Dünner, Juan Perdomo, Tijana Zrnic, and Moritz Hardt.
Stochastic optimization for performative prediction.
In *Advances in neural information processing systems*, pages 4929–4939. Curran Associates, Inc., 2020.
- Miller et al. [2021]
John Miller, Juan C. Perdomo, and Tijana Zrnic.
Outside the Echo Chamber: Optimizing the Performative Risk.
*arXiv:2102.08570 [cs, stat]*, June 2021.
URL [http://arxiv.org/abs/2102.08570](http://arxiv.org/abs/2102.08570).
arXiv: 2102.08570.
- 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.
- Nguyen et al. [2017]
Lam M Nguyen, Jie Liu, Katya Scheinberg, and Martin Takáč.
Sarah: A novel method for machine learning problems using stochastic recursive gradient.
In *International conference on machine learning*, pages 2613–2621. PMLR, 2017.
- Perdomo et al. [2021]
Juan C. Perdomo, Tijana Zrnic, Celestine Mendler-Dünner, and Moritz Hardt.
Performative Prediction.
*arXiv:2002.06673 [cs, stat]*, February 2021.
URL [http://arxiv.org/abs/2002.06673](http://arxiv.org/abs/2002.06673).
arXiv: 2002.06673.
- 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*, volume 36, pages 8081–8088, 2022.
- Reddi et al. [2016]
Sashank J Reddi, Ahmed Hefny, Suvrit Sra, Barnabas Poczos, and Alex Smola.
Stochastic variance reduction for nonconvex optimization.
In *International conference on machine learning*, pages 314–323. PMLR, 2016.
- Rosenfeld et al. [2020]
Nir Rosenfeld, Anna Hilgard, Sai Srivatsa Ravindranath, and David C Parkes.
From predictions to decisions: Using lookahead regularization.
*Advances in Neural Information Processing Systems*, 33:4115–4126, 2020.
- Roux et al. [2012]
Nicolas Roux, Mark Schmidt, and Francis Bach.
A stochastic gradient method with an exponential convergence _rate for finite training sets.
*Advances in neural information processing systems*, 25, 2012.
- Xie and Zhang [2024]
Tian Xie and Xueru Zhang.
Non-linear welfare-aware strategic learning, 2024.
- Xie et al. [2024]
Tian Xie, Xuwei Tan, and Xueru Zhang.
Algorithmic decision-making under agents with persistent improvement, 2024.
- Zhang et al. [2019]
Xueru Zhang, Mohammad Mahdi Khalili, Cem Tekin, and Mingyan Liu.
*Group Retention When Using Machine Learning in Sequential Decision Making: The Interplay between User Dynamics and Fairness*.
Curran Associates Inc., Red Hook, NY, USA, 2019.
- Zheng et al. [2024]
Xue Zheng, Tian Xie, Xuwei Tan, Aylin Yener, Xueru Zhang, Ali Payani, and Myungjin Lee.
Profl: Performative robust optimal federated learning.
*arXiv preprint arXiv:2410.18075*, 2024.
- 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:47902–47913, 2023.

## Appendix A The Use of Large Language Models (LLMs)

We only use LLMs for polishing the writing without using them to generate any single sentence solely on their own.

## Appendix B Proof of Lemma 5.1

###### Lemma .

Let $\beta_{k}$ be some positive constant, and $c_{k},c_{k+1}$ are some constants defined in intermediate functions $R_{k}^{s+1},R_{k+1}^{s+1}$ that satisfy the following

$$ $\displaystyle c_{k}-\frac{1}{2}\gamma_{k}=c_{k+1}(1+\gamma_{k}\beta_{k}+2L_{0}\epsilon\gamma_{k})+(2L^{2}+L_{0}^{2}\epsilon^{2})(L\gamma_{k}^{2}+2c_{k+1}\gamma_{k}^{2})+\frac{\beta_{k}}{2}L_{0}\epsilon\gamma_{k}$ $$

Then, we have

$$ $\displaystyle R_{k+1}^{s+1}(\boldsymbol{\theta}_{k}^{s+1})\leq R_{k}^{s+1}(\boldsymbol{\theta}_{k}^{s+1})-\Gamma_{k}\cdot\mathbb{E}[\|\nabla\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})\|^{2}]-\frac{1}{2}\gamma_{k}\|\boldsymbol{\theta}_{k}^{s+1}-\widetilde{\boldsymbol{\theta}}^{s}\|^{2}$ $$

where $\Gamma_{k}=\left(\gamma_{k}-\frac{c_{k+1}\gamma_{k}+\frac{1}{2}L_{0}\epsilon\gamma_{k}}{\beta_{k}}-2L\gamma_{k}^{2}-4c_{k+1}\gamma_{k}^{2}\right)$

###### Proof.

Since $\ell(\boldsymbol{\theta};z)$ is $L$-smooth in $\boldsymbol{\theta}$ and $v_{k}^{s+1}=\frac{1}{\gamma_{k}}\bigl(\boldsymbol{\theta}_{k}^{s+1}-\boldsymbol{\theta}_{k+1}^{s+1}\bigr)$, we can apply descent lemma to get the following result:

$$ $\displaystyle\mathbb{E}[\mathcal{J}(\boldsymbol{\theta}_{k+1}^{s+1};\boldsymbol{\theta}_{k}^{s+1})]$ $\displaystyle\leq\mathbb{E}[\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})+\frac{L}{2}\|\boldsymbol{\theta}_{k+1}^{s+1}-\boldsymbol{\theta}_{k}^{s+1}\|^{2}+\nabla\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})^{T}(\boldsymbol{\theta}_{k+1}^{s+1}-\boldsymbol{\theta}_{k}^{s+1})]$ (8) $\displaystyle=\mathbb{E}[\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})]+\frac{L\gamma_{k}^{2}}{2}\|v_{k}^{s+1}\|^{2}-\gamma_{k}\nabla\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})^{T}v_{k}^{s+1}$ $$

By definition of $v_{k}^{s+1}$, we have:

$$ $\displaystyle\mathbb{E}[v_{k}^{s+1}]$ $\displaystyle=\mathbb{E}[\nabla\ell(z_{i_{k}};\boldsymbol{\theta}_{k}^{s+1})-\nabla\ell(z_{i_{k}};\widetilde{\boldsymbol{\theta}}^{s})+\nabla\mathcal{J}(\widetilde{\boldsymbol{\theta}}^{s};\widetilde{\boldsymbol{\theta}}^{s})]$ (9) $\displaystyle=\mathbb{E}[\nabla\ell(z_{i_{k}};\boldsymbol{\theta}_{k}^{s+1})-\nabla\ell(z_{i_{k}};\widetilde{\boldsymbol{\theta}}^{s})+\nabla\mathcal{J}(\widetilde{\boldsymbol{\theta}}^{s};\boldsymbol{\theta}_{k}^{s+1})]-\mathbb{E}[\nabla\mathcal{J}(\widetilde{\boldsymbol{\theta}}^{s};\boldsymbol{\theta}_{k}^{s+1})-\nabla\mathcal{J}(\widetilde{\boldsymbol{\theta}}^{s};\widetilde{\boldsymbol{\theta}}^{s})]$ $\displaystyle=\mathbb{E}[\nabla\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})]-\mathbb{E}[\nabla\mathcal{J}(\widetilde{\boldsymbol{\theta}}^{s};\boldsymbol{\theta}_{k}^{s+1})-\nabla\mathcal{J}(\widetilde{\boldsymbol{\theta}}^{s};\widetilde{\boldsymbol{\theta}}^{s})]$ $$

The last equality holds since $\mathbb{E}[\nabla\ell(z_{i_{k}};\widetilde{\boldsymbol{\theta}}^{s})]=\nabla\mathcal{J}(\widetilde{\boldsymbol{\theta}}^{s},\boldsymbol{\theta}_{k}^{s+1})$. Then applying Lemma [5.2](https://arxiv.org/html/2509.17304v2#S5.Thmtheorem2) on the second term above to derive the following:

$$ $\displaystyle\gamma_{k}\nabla\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})^{T}v_{k}^{s+1}$ $\displaystyle=\gamma_{k}\nabla\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})^{T}\bigl(\nabla\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})-\mathbb{E}[\nabla\mathcal{J}(\widetilde{\boldsymbol{\theta}}^{s};\boldsymbol{\theta}_{k}^{s+1})-\nabla\mathcal{J}(\widetilde{\boldsymbol{\theta}}^{s};\widetilde{\boldsymbol{\theta}}^{s})]\bigr)$ $\displaystyle\geq\gamma_{k}\|\nabla\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})\|^{2}-L_{0}\epsilon\gamma_{k}\|\nabla\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})\|\|\widetilde{\boldsymbol{\theta}}^{s}-\boldsymbol{\theta}_{k}^{s+1}\|$ $$

Next, applying Young’s Inequality, we can get:

$$ 8 $\displaystyle\leq\mathbb{E}[\mathcal{J}(\boldsymbol{\theta}_{k+1}^{s+1};\boldsymbol{\theta}_{k}^{s+1})-\gamma_{k}\|\nabla\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})\|^{2}+L_{0}\epsilon\gamma_{k}(\frac{1}{2\beta_{k}}\|\nabla\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})\|^{2}$ (10) $\displaystyle~~+\frac{\beta_{k}}{2}\|\boldsymbol{\theta}_{k}^{s+1}-\widetilde{\boldsymbol{\theta}}^{s}\|^{2})+\frac{L\gamma_{k}^{2}}{2}\|v_{k}^{s+1}\|^{2}]$ $$

Next, bound $\mathbb{E}[\|\boldsymbol{\theta}^{s+1}_{k+1}-\widetilde{\boldsymbol{\theta}}^{s}\|^{2}]$ by adding and subtracting as follows:

$$ $\displaystyle\mathbb{E}[\|\boldsymbol{\theta}^{s+1}_{k+1}-\widetilde{\boldsymbol{\theta}}^{s}\|^{2}]$ $\displaystyle=\mathbb{E}[\|\boldsymbol{\theta}^{s+1}_{k+1}-\boldsymbol{\theta}_{k}^{s+1}+\boldsymbol{\theta}_{k}^{s+1}-\widetilde{\boldsymbol{\theta}}^{s}\|^{2}]$ $\displaystyle=\mathbb{E}[\gamma_{k}^{2}\|v_{k}^{s+1}\|^{2}+\|\boldsymbol{\theta}_{k}^{s+1}-\widetilde{\boldsymbol{\theta}}^{s}\|^{2}]+2\gamma_{k}\mathbb{E}[\|\boldsymbol{\theta}_{k}^{s+1}-\widetilde{\boldsymbol{\theta}}^{s}\|\|v_{k}^{s+1}\|]$ $\displaystyle\leq\mathbb{E}[\gamma_{k}^{2}\|v_{k}^{s+1}\|^{2}+\|\boldsymbol{\theta}_{k}^{s+1}-\widetilde{\boldsymbol{\theta}}^{s}\|^{2}]+2\gamma_{k}\mathbb{E}[\frac{1}{2\beta_{k}}\|\nabla\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})\|^{2}$ (11) $\displaystyle~~+\frac{\beta_{k}}{2}\|\boldsymbol{\theta}_{k}^{s+1}-\widetilde{\boldsymbol{\theta}}^{s}\|^{2}]+2L_{0}\epsilon\gamma_{k}\mathbb{E}\|\widetilde{\boldsymbol{\theta}}^{s}-\boldsymbol{\theta}_{k}^{s+1}\|^{2}$ $$

Similary, $\mathbb{E}[\|v_{k}^{s+1}\|^{2}]$ can also be bounded:

$$ $\displaystyle\mathbb{E}[\|v_{k}^{s+1}\|^{2}]$ $\displaystyle=\mathbb{E}\left[\nabla\ell(z_{i_{k}};\boldsymbol{\theta}^{s+1}_{k})-\nabla\ell(z_{i_{k}};\widetilde{\boldsymbol{\theta}}^{s})+\nabla\mathcal{J}(\widetilde{\boldsymbol{\theta}}^{s};\widetilde{\boldsymbol{\theta}}^{s})\right]^{2}$ (12) $\displaystyle\leq 2\mathbb{E}\left[\|\nabla\ell(z_{i_{k}};\boldsymbol{\theta}^{s+1}_{k})-\nabla\ell(z_{i_{k}};\widetilde{\boldsymbol{\theta}}^{s})+\nabla\mathcal{J}(\widetilde{\boldsymbol{\theta}}^{s};\boldsymbol{\theta}_{k}^{s+1})\|^{2}\right]$ $\displaystyle~~+2\mathbb{E}\left[\|\nabla\mathcal{J}(\widetilde{\boldsymbol{\theta}}^{s};\boldsymbol{\theta}_{k}^{s+1})-\nabla\mathcal{J}(\widetilde{\boldsymbol{\theta}}^{s};\widetilde{\boldsymbol{\theta}}^{s})\|^{2}\right]$ $$

The above inequality is obtained by adding and subtracting $\nabla\mathcal{J}(\widetilde{\boldsymbol{\theta}}^{s};\boldsymbol{\theta}_{k}^{s+1})$ and then applying Young’s Inequality. For the first term, it is worth noting that the term already gets rid of the performative effects, i.e., the distribution is fixed as $\mathcal{D}(\boldsymbol{\theta}_{k}^{s+1})$. Then bounding this term is equivalent to bound the norm of update in Reddi et al. [2016]. Specifically, let $\delta_{k}^{s+1}=\nabla\ell(z_{i_{k}};\boldsymbol{\theta}^{s+1}_{k})-\nabla\ell(z_{i_{k}};\widetilde{\boldsymbol{\theta}}^{s})$, we know $\mathbb{E}[\delta_{k}^{s+1}]=\mathbb{E}[\ell(z_{i_{k}};\boldsymbol{\theta}^{s+1}_{k})]-\mathbb{E}[\nabla\ell(z_{i_{k}};\widetilde{\boldsymbol{\theta}}^{s})]=\nabla\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})-\nabla\mathcal{J}(\widetilde{\boldsymbol{\theta}}^{s};\boldsymbol{\theta}_{k}^{s+1})$. Thus, we have:

$$ $\displaystyle\mathbb{E}\left[\|\nabla\ell(z_{i_{k}};\boldsymbol{\theta}^{s+1}_{k})-\nabla\ell(z_{i_{k}};\widetilde{\boldsymbol{\theta}}^{s})+\nabla\mathcal{J}(\widetilde{\boldsymbol{\theta}}^{s};\boldsymbol{\theta}_{k}^{s+1})\|^{2}\right]$ (13) $\displaystyle=\mathbb{E}[\|\delta_{k}^{s+1}+\nabla\mathcal{J}(\widetilde{\boldsymbol{\theta}}^{s};\boldsymbol{\theta}_{k}^{s+1})\|^{2}]$ $\displaystyle=\mathbb{E}[\|\delta_{k}^{s+1}+\nabla\mathcal{J}(\widetilde{\boldsymbol{\theta}}^{s};\boldsymbol{\theta}_{k}^{s+1})-\nabla\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})+\nabla\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})\|^{2}]$ $\displaystyle=\mathbb{E}[\|\delta_{k}^{s+1}-\mathbb{E}[\delta_{k}^{s+1}]+\nabla\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})\|^{2}]$ $\displaystyle\leq 2\mathbb{E}[\|\delta_{k}^{s+1}-\mathbb{E}[\delta_{k}^{s+1}]\|^{2}]+2\mathbb{E}[\|\nabla\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})\|^{2}]$ $\displaystyle\leq 2\mathbb{E}[\|\delta_{k}^{s+1}\|^{2}]+2\mathbb{E}[\|\nabla\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})\|^{2}]$ $\displaystyle\leq 2L^{2}\mathbb{E}[\|\widetilde{\boldsymbol{\theta}}_{k}^{s+1}-\widetilde{\boldsymbol{\theta}}^{s}\|^{2}]+2\mathbb{E}[\|\nabla\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})\|^{2}]$ $$

From Lemma [5.2](https://arxiv.org/html/2509.17304v2#S5.Thmtheorem2), we know $\|\nabla\mathcal{J}(\widetilde{\boldsymbol{\theta}}^{s};\boldsymbol{\theta}_{k}^{s+1})-\nabla\mathcal{J}(\widetilde{\boldsymbol{\theta}}^{s};\widetilde{\boldsymbol{\theta}}^{s})\|\leq L_{0}\epsilon\|\widetilde{\boldsymbol{\theta}}_{k}^{s+1}-\widetilde{\boldsymbol{\theta}}^{s}\|$ always holds. So we have:

$$ $\displaystyle\mathbb{E}\left[\|\nabla\mathcal{J}(\widetilde{\boldsymbol{\theta}}^{s};\boldsymbol{\theta}_{k}^{s+1})-\nabla\mathcal{J}(\widetilde{\boldsymbol{\theta}}^{s};\widetilde{\boldsymbol{\theta}}^{s})\|^{2}\right]\leq L_{0}^{2}\epsilon^{2}\mathbb{E}\|\boldsymbol{\theta}_{k}^{s+1}-\widetilde{\boldsymbol{\theta}}^{s}\|^{2}$ (14) $$

Sum them together, we have:

$$ $\displaystyle\mathbb{E}[\|v_{k}^{s+1}\|^{2}]\leq 4\|\nabla\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})\|^{2}+(4L^{2}+2L_{0}^{2}\epsilon^{2})\mathbb{E}\|\boldsymbol{\theta}_{k}^{s+1}-\widetilde{\boldsymbol{\theta}}^{s}\|^{2}$ (15) $$

Now consider $R_{k+1}^{s+1}$ and plug [8](https://arxiv.org/html/2509.17304v2#A2.E8) and [B](https://arxiv.org/html/2509.17304v2#A2.Ex22) into its formula, we have:

$$ $\displaystyle R_{k+1}^{s+1}(\boldsymbol{\theta}_{k}^{s+1})$ $\displaystyle=\mathbb{E}\left[\mathcal{J}(\boldsymbol{\theta}_{k+1}^{s+1};\boldsymbol{\theta}_{k}^{s+1})+C_{k+1}\|\widetilde{\boldsymbol{\theta}}_{k+1}^{s+1}-\widetilde{\boldsymbol{\theta}}^{s}\|^{2}\right]$ (16) $\displaystyle\leq\mathbb{E}[\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})-\gamma_{k}\|\nabla\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1},\boldsymbol{\theta}_{k}^{s+1})\|^{2}$ $\displaystyle~~+L_{0}\epsilon\gamma_{k}(\frac{1}{2\beta_{k}}\|\nabla\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})\|^{2}+\frac{\beta_{k}}{2}\|\boldsymbol{\theta}_{k}^{s+1}-\widetilde{\boldsymbol{\theta}}^{s}\|^{2})$ $\displaystyle~~+\frac{L\gamma_{k}^{2}}{2}\|v_{k}^{s+1}\|^{2}]+\mathbb{E}[c_{k+1}\gamma_{k}^{2}\|v_{k}^{s+1}\|^{2}$ $\displaystyle~~+c_{k+1}\|\boldsymbol{\theta}_{k+1}^{s+1}-\widetilde{\boldsymbol{\theta}}^{s}\|^{2}+2c_{k+1}\gamma_{k}\mathbb{E}(\frac{1}{2\beta_{k}}\|\nabla\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})\|^{2}$ (17) $\displaystyle~~+\frac{\beta_{k}}{2}\|\boldsymbol{\theta}_{k}^{s+1}-\widetilde{\boldsymbol{\theta}}^{s}\|^{2})]+2c_{k+1}L_{0}\epsilon\gamma_{k}\mathbb{E}\left\|\widetilde{\boldsymbol{\theta}}^{s}-\boldsymbol{\theta}_{k}^{s+1}\right\|^{2}$ $\displaystyle=\mathbb{E}\left[\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})\right]+\left(\frac{L\gamma_{k}^{2}}{2}+c_{k+1}\gamma_{k}^{2}\right)\mathbb{E}\left[\|v_{k}^{s+1}\|^{2}\right]$ $\displaystyle~~-\left(\gamma_{k}-\frac{c_{k+1}\gamma_{k}+\frac{1}{2}L_{0}\epsilon\gamma_{k}}{\beta_{k}}\right)\mathbb{E}\left[\|\nabla\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})\|^{2}\right]$ $\displaystyle~~+\left(c_{k+1}(1+\gamma_{k}\beta_{k}+2L_{0}\epsilon\gamma_{k})+\frac{\beta_{k}L_{0}\epsilon\gamma_{k}}{2}\right)\mathbb{E}\left[\|\widetilde{\boldsymbol{\theta}}^{s}-\boldsymbol{\theta}_{k}^{s+1}\|^{2}\right]$ $$

Regarding the $\|v_{k}^{s+1}\|^{2}$, we plug [15](https://arxiv.org/html/2509.17304v2#A2.E15) in [16](https://arxiv.org/html/2509.17304v2#A2.E16) to get the expression we need:

$$ $\displaystyle R_{k+1}^{s+1}(\boldsymbol{\theta}_{k}^{s+1})$ $\displaystyle\leq\mathbb{E}\left[\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})\right]-\Gamma_{k}\mathbb{E}\left[\|\nabla\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})\|^{2}\right]+(c_{k}-\frac{1}{2}\gamma_{k})\mathbb{E}\|\boldsymbol{\theta}_{k}^{s+1}-\widetilde{\boldsymbol{\theta}}^{s}\|^{2}$ $\displaystyle=R_{k}^{s+1}-\Gamma_{k}\mathbb{E}\left[\|\nabla\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})\|^{2}\right]-\frac{1}{2}\gamma_{k}\mathbb{E}\|\boldsymbol{\theta}_{k}^{s+1}-\widetilde{\boldsymbol{\theta}}^{s}\|^{2}$ $$

where

$$ $\displaystyle c_{k}-\frac{1}{2}\gamma_{k}=c_{k+1}(1+\gamma_{k}\beta_{k}+2L_{0}\epsilon\gamma_{k})+(2L^{2}+L_{0}^{2}\epsilon^{2})(L\gamma_{k}^{2}+2c_{k+1}\gamma_{k}^{2})+\frac{\beta_{k}}{2}L_{0}\epsilon\gamma_{k}$ $\displaystyle\Gamma_{k}=\left(\gamma_{k}-\frac{c_{k+1}\gamma_{k}+\frac{1}{2}L_{0}\epsilon\gamma_{k}}{\beta_{k}}-2L\gamma_{k}^{2}-4c_{k+1}\gamma_{k}^{2}\right)$ $$

∎

###### Lemma .

###### Proof.

## Appendix C Proof of Lemma 5.3

###### Lemma .

With $c_{k},c_{k+1},\beta_{k},\Gamma_{k}$ same as Lemma [5.1](https://arxiv.org/html/2509.17304v2#S5.Thmtheorem1), denote

$$ $\displaystyle\widetilde{R}_{k}^{s+1}\overset{\triangle}{=}\mathbb{E}[\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})+c_{k}\|\boldsymbol{\theta}_{k}^{s+1}-\widetilde{\boldsymbol{\theta}}^{s}\|^{2}]$ $\displaystyle\widetilde{R}_{k+1}^{s+1}\overset{\triangle}{=}\mathbb{E}[\mathcal{J}(\boldsymbol{\theta}_{k+1}^{s+1};\boldsymbol{\theta}_{k+1}^{s+1})+c_{k+1}\|\boldsymbol{\theta}_{k+1}^{s+1}-\widetilde{\boldsymbol{\theta}}^{s}\|^{2}]$ $$

Then we have:

$$ $\displaystyle\widetilde{R}_{k+1}^{s+1}$ $\displaystyle\leq\widetilde{R}_{k}^{s+1}+(L_{0}^{4}\epsilon^{4}+2L_{0}^{2}\epsilon^{2}L^{2}+2L_{0}^{2}\epsilon^{2})\gamma_{k}-(\Gamma_{k}-\frac{\gamma_{k}}{4})\cdot\mathbb{E}[\|\nabla\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})\|^{2}]$ (18) $$

###### Proof.

According to Lemma [5.1](https://arxiv.org/html/2509.17304v2#S5.Thmtheorem1) (Eqn. [6](https://arxiv.org/html/2509.17304v2#S5.E6)), we have:

$$ $\displaystyle\mathbb{E}[\mathcal{J}(\boldsymbol{\theta}_{k+1}^{s+1};\boldsymbol{\theta}_{k}^{s+1})+c_{k+1}\|\boldsymbol{\theta}_{k+1}^{s+1}-\widetilde{\boldsymbol{\theta}}^{s}\|^{2}]$ $\displaystyle\leq\mathbb{E}[\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})]+c_{k}\|\boldsymbol{\theta}_{k}^{s+1}-\widetilde{\boldsymbol{\theta}}^{s}\|^{2}$ (19) $\displaystyle-\Gamma_{k}\|\nabla\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})\|^{2}-\frac{\gamma_{k}}{2}\|\boldsymbol{\theta}_{k}^{s+1}-\widetilde{\boldsymbol{\theta}}^{s}\|^{2}]$ $$

By aggregating terms, we will have:

$$ $\displaystyle~~~~c_{k+1}\mathbb{E}\|\boldsymbol{\theta}_{k+1}^{s+1}-\widetilde{\boldsymbol{\theta}}^{s}\|^{2}-\left(c_{k}-\frac{\gamma_{k}}{2}\right)\mathbb{E}\|\boldsymbol{\theta}_{k}^{s+1}-\widetilde{\boldsymbol{\theta}}^{s}\|^{2}+\Gamma_{k}\mathbb{E}\|\nabla\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})\|^{2}$ (20) $\displaystyle\leq\mathbb{E}[\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})-\mathcal{J}(\boldsymbol{\theta}_{k+1}^{s+1};\boldsymbol{\theta}_{k}^{s+1})]$ $$

We can further bound the above expression by:

$$ $\displaystyle\mathbb{E}[\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})-\mathcal{J}(\boldsymbol{\theta}_{k+1}^{s+1};\boldsymbol{\theta}_{k}^{s+1})]$ $\displaystyle=\mathbb{E}\left[\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})-\mathcal{J}(\boldsymbol{\theta}_{k+1}^{s+1};\boldsymbol{\theta}_{k+1}^{s+1})\right]$ (21) $\displaystyle~~+\mathbb{E}\left[\mathcal{J}(\boldsymbol{\theta}_{k+1}^{s+1};\boldsymbol{\theta}_{k+1}^{s+1})-\mathcal{J}(\boldsymbol{\theta}_{k+1}^{s+1};\boldsymbol{\theta}_{k}^{s+1})\right]$ $\displaystyle\leq\mathbb{E}\left[\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})-\mathcal{J}(\boldsymbol{\theta}_{k+1}^{s+1};\boldsymbol{\theta}_{k+1}^{s+1})\right]$ $\displaystyle~~+L_{0}\epsilon\mathbb{E}\|\boldsymbol{\theta}_{k+1}^{s+1}-\boldsymbol{\theta}_{k}^{s+1}\|$ $$

Regarding $\mathbb{E}\|\boldsymbol{\theta}_{k+1}^{s+1}-\boldsymbol{\theta}_{k}^{s+1}\|=\gamma_{k}\mathbb{E}\|v_{k}^{s+1}\|$, we can bound it according to Eqn. [9](https://arxiv.org/html/2509.17304v2#A2.E9) and Eqn. [15](https://arxiv.org/html/2509.17304v2#A2.E15):

$$ $\displaystyle\mathbb{E}[\|v_{k}^{s+1}\|]$ $\displaystyle=\mathbb{E}[\|\nabla\ell(z_{i_{k}};\boldsymbol{\theta}_{k}^{s+1})-\nabla\ell(z_{i_{k}};\widetilde{\boldsymbol{\theta}}^{s})+\nabla\mathcal{J}(\widetilde{\boldsymbol{\theta}}^{s};\widetilde{\boldsymbol{\theta}}^{s})\|]$ (22) $\displaystyle\leq\mathbb{E}[\|\nabla\ell(z_{i_{k}};\boldsymbol{\theta}_{k}^{s+1})-\nabla\ell(z_{i_{k}};\widetilde{\boldsymbol{\theta}}^{s})+\nabla\mathcal{J}(\widetilde{\boldsymbol{\theta}}^{s};\boldsymbol{\theta}_{k}^{s+1})\|]$ $\displaystyle~~~~+\mathbb{E}[\|\nabla\mathcal{J}(\widetilde{\boldsymbol{\theta}}^{s};\boldsymbol{\theta}_{k}^{s+1})-\nabla\mathcal{J}(\widetilde{\boldsymbol{\theta}}^{s};\widetilde{\boldsymbol{\theta}}^{s})\|]$ $$

Let $\nabla\ell(z_{i_{k}};\boldsymbol{\theta}_{k}^{s+1})-\nabla\ell(z_{i_{k}};\widetilde{\boldsymbol{\theta}}^{s})+\nabla\mathcal{J}(\widetilde{\boldsymbol{\theta}}^{s};\boldsymbol{\theta}_{k}^{s+1})$ be $A$ and $\nabla\mathcal{J}(\widetilde{\boldsymbol{\theta}}^{s};\boldsymbol{\theta}_{k}^{s+1})-\nabla\mathcal{J}(\widetilde{\boldsymbol{\theta}}^{s};\widetilde{\boldsymbol{\theta}}^{s})$ be $B$.
According to Eqn. [13](https://arxiv.org/html/2509.17304v2#A2.E13), we have:

$$ $\displaystyle\mathbb{E}[\|A\|^{2}]\leq 2\mathbb{E}\left[\|\nabla\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})\|^{2}\right]+2L^{2}\mathbb{E}\|\boldsymbol{\theta}_{k}^{s+1}-\widetilde{\boldsymbol{\theta}}^{s}\|^{2}$ $$

which means

$$ $\displaystyle\mathbb{E}[\|A\|]\leq\sqrt{\mathbb{E}[\|A\|^{2}]}\leq\sqrt{2\mathbb{E}\left[\|\nabla\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})\|^{2}\right]+2L^{2}\mathbb{E}\|\boldsymbol{\theta}_{k}^{s+1}-\widetilde{\boldsymbol{\theta}}^{s}\|^{2}}$ $$

Applying Young’s Inequality on rhs, for any positive $\lambda$ we have:

$$ $\displaystyle\mathbb{E}\|A\|\leq\frac{1}{\lambda}\left(\mathbb{E}\left[\|\nabla\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})\|^{2}\right]+L^{2}\mathbb{E}\left\|\boldsymbol{\theta}_{k}^{s+1}-\widetilde{\boldsymbol{\theta}}^{s}\right\|^{2}\right)+\frac{\lambda}{2}$ (23) $$

Let $\lambda_{0}=4L_{0}\epsilon(L^{2}+1)$, we have:

$$ $\displaystyle\mathbb{E}\|A\|$ $\displaystyle\leq\frac{1}{\lambda_{0}}\left(\mathbb{E}\left[\|\nabla\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})\|^{2}\right]+L^{2}\mathbb{E}\left\|\boldsymbol{\theta}_{k}^{s+1}-\widetilde{\boldsymbol{\theta}}^{s}\right\|^{2}\right)+\frac{\lambda_{0}}{2}$ (24) $\displaystyle\leq\frac{1}{4L_{0}\epsilon}\mathbb{E}\left[\|\nabla\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})\|^{2}\right]+\frac{1}{4L_{0}\epsilon}\mathbb{E}\left\|\boldsymbol{\theta}_{k}^{s+1}-\widetilde{\boldsymbol{\theta}}^{s}\right\|^{2}+\frac{\lambda_{0}}{2}$ $$

According to the point-wise Lipschitzness in Lemma [5.2](https://arxiv.org/html/2509.17304v2#S5.Thmtheorem2), we have:

$$ $\displaystyle\mathbb{E}[\|B\|]\leq\sqrt{\mathbb{E}[\|B\|^{2}]}\leq\sqrt{L_{0}^{2}\epsilon^{2}\mathbb{E}\left\|\boldsymbol{\theta}_{k}^{s+1}-\boldsymbol{\theta}^{s}\right\|^{2}}\leq\frac{\mathbb{E}\left\|\boldsymbol{\theta}_{k}^{s+1}-\widetilde{\boldsymbol{\theta}}^{s}\right\|^{2}}{4L_{0}\epsilon}+L_{0}^{3}\epsilon^{3}$ (25) $$

Thus, we have

$$ $\displaystyle\mathbb{E}[\|v_{k}^{s+1}\|]$ $\displaystyle\leq\mathbb{E}\|A\|+\mathbb{E}\|B\|$ (26) $\displaystyle\leq\frac{\mathbb{E}\|\nabla\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})\|^{2}}{4L_{0}\epsilon}+(\frac{1}{4L_{0}\epsilon}+\frac{1}{4L_{0}\epsilon})\,\mathbb{E}\left\|\boldsymbol{\theta}_{k}^{s+1}-\widetilde{\boldsymbol{\theta}}^{s}\right\|^{2}+\frac{\lambda_{0}}{2}+L_{0}^{3}\epsilon^{3}$ $\displaystyle=\frac{\mathbb{E}\|\nabla\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})\|^{2}}{4L_{0}\epsilon}+\frac{1}{2L_{0}\epsilon}\mathbb{E}\left\|\boldsymbol{\theta}_{k}^{s+1}-\widetilde{\boldsymbol{\theta}}^{s}\right\|^{2}+(\frac{\lambda_{0}}{2}+L_{0}^{3}\epsilon^{3})$ $$

Finally, we plug Eqn. [26](https://arxiv.org/html/2509.17304v2#A3.E26) into Eqn. [21](https://arxiv.org/html/2509.17304v2#A3.E21), and then plug Eqn. [21](https://arxiv.org/html/2509.17304v2#A3.E21) into Eqn. [20](https://arxiv.org/html/2509.17304v2#A3.E20). This leads to the inequality relationship between $\mathcal{J}(\boldsymbol{\theta}_{k+1}^{s+1};\boldsymbol{\theta}_{k+1}^{s+1})$ and $\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})$.

$$ $\displaystyle\mathbb{E}\left[\mathcal{J}(\boldsymbol{\theta}_{k+1}^{s+1};\boldsymbol{\theta}_{k+1}^{s+1})+c_{k+1}\|\boldsymbol{\theta}_{k+1}^{s+1}-\widetilde{\boldsymbol{\theta}}^{s}\|^{2}\right]$ $\displaystyle\leq\mathbb{E}\left[\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})+c_{k}\|\boldsymbol{\theta}_{k}^{s+1}-\widetilde{\boldsymbol{\theta}}^{s}\|^{2}\right]$ $\displaystyle~~-(\Gamma_{k}-\frac{\gamma_{k}}{4})\cdot\mathbb{E}[\|\nabla\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})\|^{2}]$ $\displaystyle~~+(L_{0}^{4}\epsilon^{4}+2L_{0}^{2}\epsilon^{2}L^{2}+2L_{0}^{2}\epsilon^{2})\gamma_{k}$ $$

which is exactly what we want.
∎

###### Lemma .

###### Proof.

## Appendix D Proof of Thm. 5.4

###### Theorem .

Let the total rounds be $T$ and the epochs be $S=\lceil T/m\rceil$. For the last epoch, $c_{m}=0$. Define $\Delta_{0}$ and $\Delta_{1}$ as follows:

$$ $\displaystyle\Delta_{0}\overset{\triangle}{=}\mathcal{J}(\boldsymbol{\theta}_{0}^{0};\boldsymbol{\theta}_{0}^{0})-\ell_{*};~~~~~\Delta_{1}\overset{\triangle}{=}\frac{(L_{0}^{4}\epsilon^{4}+2L_{0}^{2}\epsilon^{2}L^{2}+2L_{0}^{2}\epsilon^{2})\sum_{s=0}^{S}\sum_{k=0}^{m}\gamma_{k}}{T\Gamma}$ $$

where $\Gamma>0$ is guaranteed to exist as a constant and is the lower bound of $\Gamma_{k}-\frac{\gamma_{k}}{4}$. Then we have:

$$ $\displaystyle\frac{1}{T}\sum_{s=0}^{S-1}\sum_{k=0}^{m-1}\mathbb{E}\|\nabla\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})\|^{2}\leq\frac{\Delta_{0}}{T\Gamma}+\Delta_{1}$ (27) $$

Thus, Alg. [1](https://arxiv.org/html/2509.17304v2#alg1) can achieve $\mathcal{O}(\frac{1}{T})$ convergence rate and the error neighborhood $\Delta_{1}$ is $\mathcal{O}(\epsilon^{2}+\epsilon^{4})$.

###### Proof.

Note the $\gamma_{t}$ is just an abbreviation of the learning rate at global round $s$ and local round $k$. Using Lemma [5.3](https://arxiv.org/html/2509.17304v2#S5.Thmtheorem3) (Eqn. [6](https://arxiv.org/html/2509.17304v2#S5.E6)) and telescoping:

$$ $\displaystyle\sum_{s=0}^{S-1}\sum_{k=0}^{m-1}(\Gamma_{k}-\frac{\gamma_{k}}{4})\cdot\mathbb{E}[\|\nabla\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})\|^{2}]$ $\displaystyle\leq(\widetilde{R}_{0}^{1}-\widetilde{R}_{m}^{S})+(L_{0}^{4}\epsilon^{4}+2L_{0}^{2}\epsilon^{2}L^{2}+2L_{0}^{2}\epsilon^{2})\sum_{t=0}^{T-1}\gamma_{t}$ (28) $$

Thus,

$$ $\displaystyle\frac{1}{T}\cdot\sum_{s=0}^{S-1}\sum_{k=0}^{m-1}\mathbb{E}[\|\nabla\mathcal{J}(\boldsymbol{\theta}_{k}^{s+1};\boldsymbol{\theta}_{k}^{s+1})\|^{2}]\leq\frac{\Delta_{0}}{T\Gamma}+\Delta_{1}$ $$

There is one remaining step: prove there exists values of $\gamma_{k},\beta_{k}$ to let $\Gamma>0$. We first begin with a lemma:

###### Lemma D.1 .

For any $\eta>0$ and $\beta_{k}$ is some constant $\beta$, there exists a sequence of learning rate $\{\gamma_{k}\}_{k=1}^{m}$ to let $c_{k}\leq\eta L_{0}\epsilon$ hold for any $k$.

###### Proof.

We can start from $c_{m}=0$ and are given the relationship between $c_{k}$ and $c_{k+1}$:

$$ $\displaystyle c_{k}-\frac{1}{2}\gamma_{k}=c_{k+1}(1+\gamma_{k}\beta_{k}+2L_{0}\epsilon\gamma_{k})+(2L^{2}+L_{0}^{2}\epsilon^{2})(L\gamma_{k}^{2}+2c_{k+1}\gamma_{k}^{2})+\frac{\beta_{k}}{2}L_{0}\epsilon\gamma_{k}$ $$

For convenience, define $A_{k}=1+\gamma_{k}\beta+2L_{0}\epsilon\,\gamma_{k}+(4L^{2}+2L_{0}^{2}\epsilon^{2})\gamma_{k}^{2}$, $B_{k}=(2L^{2}+L_{0}^{2}\epsilon^{2})L\gamma_{k}^{2}+\frac{\beta}{2}L_{0}\epsilon\,\gamma_{k}$. Then, the recurrence can be rewritten as $c_{k}=\frac{1}{2}\gamma_{k}+A_{k}\,c_{k+1}+B_{k}$. Next, we use induction to prove that for any $\eta>0$, there exists a sequence of $\{\gamma_{k}\}_{k=1}^{m}$ (depending only on $L$, $L_{0}$, $\epsilon$, and $\beta$) such that for all indices $k$ we have $c_{k}\leq\eta L_{0}\epsilon$.

Since we start from $c_{m}=0$, we have $c_{m}=0\leq\eta\,L_{0}\epsilon$ holding for any positive $\eta$. Next, assume for some $k+1$, $c_{k+1}\leq\eta\,L_{0}\epsilon$. Then we can work out $c_{k}$ as:

$$ $\displaystyle c_{k}=\frac{1}{2}\gamma_{k}+A_{k}\,c_{k+1}+B_{k}\leq\frac{1}{2}\gamma_{k}+A_{k}\Bigl(\eta\,L_{0}\epsilon+C\,\gamma_{k}\Bigr)+B_{k}.$ $$

Since we can adjust $\gamma_{k}$ to be small, we can expand $A_{k}=1+\gamma_{k}\beta+2L_{0}\epsilon\,\gamma_{k}+\mathcal{O}(\gamma_{k}^{2})$ and $B_{k}=\frac{\beta}{2}L_{0}\epsilon\,\gamma_{k}+\mathcal{O}(\gamma_{k}^{2})$. Thus,

$$ $\displaystyle c_{k}$ $\displaystyle\leq\frac{1}{2}\gamma_{k}+\Bigl(1+\gamma_{k}\beta+2L_{0}\epsilon\,\gamma_{k}\Bigr)\eta\,L_{0}\epsilon+\frac{\beta}{2}L_{0}\epsilon\,\gamma_{k}+C_{1}\gamma_{k}^{2}$ $\displaystyle=\eta\,L_{0}\epsilon+\Biggl\{\frac{1}{2}+\eta\,L_{0}\epsilon(\beta+2L_{0}\epsilon)+\frac{\beta}{2}L_{0}\epsilon\Biggr\}\gamma_{k}+C_{1}\gamma_{k}^{2}$ $\displaystyle=\eta\,L_{0}\epsilon+C_{0}\gamma_{k}+C_{1}\gamma_{k}^{2}.$ $$

where $C_{0},C_{1}$ are constants represented by $L,L_{0},\epsilon,\beta$. This means for arbitrarily small number $\omega_{k}$, we can simply set $\gamma_{k}$ to be $\min\{1,\frac{\omega_{k}}{C_{0}+C_{1}}\}$ to let $c_{k}\leq(\eta+\omega_{k})L_{0}\epsilon$ . Then from $c_{m}\leq\eta L_{0}\epsilon$ we can derive that $c_{k}\leq(\eta+\Sigma_{i=0}^{m}\omega_{i})L_{0}\epsilon$ holds. Finally, for any $\eta$, we can surely find a $\widetilde{\eta}>0$ to let $\widetilde{\eta}+\Sigma_{i=0}^{m}\omega_{i}<\eta$, so we prove that $c_{k}\leq\eta L_{0}\epsilon$.
∎

Next, we prove that with the learning rates in Lemma [D.1](https://arxiv.org/html/2509.17304v2#A4.Thmtheorem1), then $\Gamma>0$. We first assume $\gamma_{k}<\frac{1}{8L+\frac{8}{3}L_{0}\epsilon}$.

$$ $\displaystyle\Gamma_{k}$ $\displaystyle=\gamma_{k}-\frac{c_{k+1}\gamma_{k}+\frac{1}{2}L_{0}\epsilon\gamma_{k}}{\beta}-2L\gamma^{2}-4c_{k+1}\gamma_{k}^{2}.$ $$

Using the bound $c_{k+1}\leq\eta L_{0}\epsilon$ and setting $\beta=4L+\frac{4}{3}L_{0}\epsilon$, we obtain

$$ $\displaystyle\Gamma_{k}$ $\displaystyle\geq\gamma_{k}-\frac{(\eta L_{0}\epsilon+\frac{1}{2}L_{0}\epsilon)\gamma_{k}}{4L+\frac{4}{3}L_{0}\epsilon}-2L\gamma_{k}^{2}-4\eta L_{0}\epsilon\gamma_{k}^{2}$ $\displaystyle=\gamma_{k}-\frac{(\eta+\frac{1}{2})L_{0}\epsilon\gamma_{k}}{4L+\frac{4}{3}L_{0}\epsilon}-2L\gamma_{k}^{2}-4\eta L_{0}\epsilon\gamma_{k}^{2}.$ $$

Then we can bound the quadratic terms as follows:

- •
$2L\gamma_{k}^{2}\leq 2L\left(\frac{1}{8L+\frac{8}{3}L_{0}\epsilon}\right)\gamma_{k}<2L\cdot\frac{1}{8L}\gamma_{k}=\frac{\gamma_{k}}{4}$
- •
$4\eta L_{0}\epsilon\gamma_{k}^{2}\leq\frac{2\eta L_{0}\epsilon\gamma_{k}}{4L+\frac{4}{3}L_{0}\epsilon}$ (since $\gamma_{k}\leq\frac{1}{8L+\frac{8}{3}L_{0}\epsilon}$).

Thus,

$$ $\displaystyle\Gamma_{k}$ $\displaystyle\geq\gamma_{k}-\frac{(\eta+\frac{1}{2})L_{0}\epsilon\gamma_{k}}{4L+\frac{4}{3}L_{0}\epsilon}-\frac{\gamma_{k}}{4}-\frac{2\eta L_{0}\epsilon\gamma_{k}}{4L+\frac{4}{3}L_{0}\epsilon}.$ $$

Combine the two fractions involving $L_{0}\epsilon$:

$$ $\displaystyle\Gamma_{k}$ $\displaystyle\geq\gamma_{k}\left[1-\frac{1}{4}-\frac{(3\eta+\frac{1}{2})L_{0}\epsilon}{4L+\frac{4}{3}L_{0}\epsilon}\right].$ $$

To bound $\Gamma_{k}\geq\frac{\gamma_{k}}{4}$, we need:

$$ $\displaystyle 1-\frac{1}{4}-\frac{(3\eta+\frac{1}{2})L_{0}\epsilon}{4L+\frac{4}{3}L_{0}\epsilon}\geq\frac{1}{4}.$ $$

This means if $\eta\leq\frac{4L+\frac{1}{3}L_{0}\epsilon}{6L_{0}\epsilon}$, we have $\Gamma\geq\frac{\gamma_{k}}{4}$. Finally, denote $\frac{4L+\frac{1}{3}L_{0}\epsilon}{6L_{0}\epsilon}$ as $\eta_{0}$, then consider $\eta=\eta_{0}$ in Lemma [D.1](https://arxiv.org/html/2509.17304v2#A4.Thmtheorem1), we can let $\omega_{k}=\omega=\frac{\eta_{0}}{2m}$ to ensure $\Sigma_{i=1}^{m}\omega_{k}<\eta_{0}$. Then we can set $\gamma=\min\{1,\frac{1}{8L+\frac{8}{3}L_{0}\epsilon}\frac{\omega}{C_{0}+C_{1}}\}$ and set $\gamma_{k}=\gamma$ to let so that Lemma [D.1](https://arxiv.org/html/2509.17304v2#A4.Thmtheorem1) also holds. Aggregating all these together, we finally prove that $\Gamma\geq\min\frac{\gamma_{k}}{4}=\frac{\gamma}{4}$ which is a constant, resulting in the $\mathcal{O}(\frac{1}{T})$ convergence given that $\gamma$ is a constant.

Finally, consider $\Delta_{1}$ when $\Gamma\geq\frac{\gamma}{4}$, it is obvious it reduces to $L_{0}^{2}\epsilon^{2}+L_{0}^{4}\epsilon^{4}$ multiplied by some constant. When $\epsilon$ is smaller than $1$, the error neighborhood is $\mathcal{O}(\epsilon^{2}+\epsilon^{4})$.
∎

###### Theorem .

###### Proof.

###### Lemma D.1 .

###### Proof.

## Appendix E Proof of Corollary 5.5

###### Corollary .

Assume we have a finite distribution setting where the population consists of $n$ samples in total and $\alpha\in(0,1)$. Then when we have fixed parameters at each round, i.e., $\gamma_{k}=\gamma=\mathcal{O}(\frac{1}{n^{\alpha}}),\beta_{k}=\beta=\mathcal{O}(n^{\frac{\alpha}{2}}),m=\mathcal{O}(n^{\frac{\alpha}{2}})$, then we have $\Gamma=\mathcal{O}(\frac{1}{n^{\alpha}})$, and the IFO complexity is $\mathcal{O}(\frac{(n^{\alpha}+n^{1+\frac{\alpha}{2}})\Delta_{0}}{\delta})$ to achieve $\mathcal{O}(\delta)$-SPS in addition to the error neighborhood.

###### Proof.

From the Lemma [D.1](https://arxiv.org/html/2509.17304v2#A4.Thmtheorem1) and the proof in App. [D](https://arxiv.org/html/2509.17304v2#A4), we already know that when $m=\mathcal{O}(n^{\frac{\alpha}{2}})$, any $\gamma_{k}\leq\gamma=\mathcal{O}(\frac{\eta_{0}}{2m(C_{0}+C_{1})})$ will make $\Gamma>\frac{\gamma_{k}}{2}$. Since $C_{0},C_{1}$ are all $\mathcal{O}(\beta)=\mathcal{O}(n^{\frac{\alpha}{2}})$, then we know $\gamma=\mathcal{O}(\frac{1}{n^{\alpha}})$ is a plausible choice to satisfy Lemma [D.1](https://arxiv.org/html/2509.17304v2#A4.Thmtheorem1) to let $\Gamma>\frac{\gamma}{2}$. Moreover, $\Gamma>\frac{\gamma}{2}=\mathcal{O}(\frac{1}{n^{\alpha}})$
With this fact, we let $\frac{\Delta_{0}}{T\Gamma}\leq\delta$ to get $T$ should be $\mathcal{O}(\frac{n^{\alpha}\Delta_{0}}{\delta})$. Since for each $m$ rounds ($1$ epoch), we would have $2m+n$ evaluations of gradients. So the average sample complexity at each round is $2+\frac{n}{m}$. Then, the final IFO complexity is:

$$ $\displaystyle T(2+\frac{n}{m})=\mathcal{O}\left(\frac{\Delta_{0}n^{\alpha}(2+n^{1-\frac{\alpha}{2}})}{\delta}\right)=\mathcal{O}(\frac{(n^{\alpha}+n^{1+\frac{\alpha}{2}})\Delta_{0}}{\delta})$ $$

∎

###### Corollary .

###### Proof.

## Appendix F Additional experimental results

Additional experimental settings. We conduct all the experiments on a server which has two Intel Xeon 6326 CPUs and a Nvidia A6000 GPU. We implement our code using the pytorch of version 2.4.1 and seed $2024$. In the PP settings, each training requires extensive steps to model distributional shifts, which significantly increases the computational cost, especially on large-scale datasets like CIFAR-10. We are only able to run one experiment on CIFAR-10, but we include experiments using $3$ random seeds for the MNIST dataset as follows.

Figure: (a) $\alpha=20$
Refer to caption: x13.png

Figure: (a) $\alpha=0.01$
Refer to caption: x19.png

Experiments on MNIST dataset. We use the MNIST dataset [Deng, 2012] and train a two-layer MLP to predict the $10$ possible digits. We randomly sample approximately $12000$ images and train the model using a learning rate set at $0.003$. Similar to retention dynamics [Hashimoto et al., 2018, Jin et al., 2024, Zhang et al., 2019], we assume the class distribution will change based on the performance of the current model. The fraction in class $c$ at iteration $k$ of $p^{k+1}_{c}=\frac{e^{-\alpha\ell^{k}_{c}}}{\sum_{c^{\prime}}e^{-\alpha\ell^{k}_{c^{\prime}}}}$, where $\ell^{k}_{c}$ is the loss expectation of class $c$ in iteration $k$. This means the image class with larger loss now will less likely to appear at the next iteration [Zhang et al., 2019]. We train the model for $80$ epochs with cross entropy loss and visualize the accuracy and training loss curve in Fig. [3](https://arxiv.org/html/2509.17304v2#A6.F3). We use three random seeds, which are 2024, 2025, 2026, to run the experiments. The $\alpha$ is set at $20,50,80$ to simulate different magnitudes of performative effects. In Fig. [3](https://arxiv.org/html/2509.17304v2#A6.F3), the training loss of $\mathsf{SPRINT}$ consistently converges faster, which is also reflected in the training accuracy.

The average of squared gradient norm. We visualize the cumulative average of the squared gradient norm at the end of each epoch $\|\nabla\mathcal{J}(\widetilde{\boldsymbol{\theta}}^{s};\widetilde{\boldsymbol{\theta}}^{s})\|^{2}$ of SGD-GD and $\mathsf{SPRINT}$ in all $3$ settings as specified in Sec. [6](https://arxiv.org/html/2509.17304v2#S6) shown in Fig. [4](https://arxiv.org/html/2509.17304v2#A6.F4). Notably, in MNIST setting (Fig. [4(d)](https://arxiv.org/html/2509.17304v2#A6.F4.sf4) to Fig. [4(f)](https://arxiv.org/html/2509.17304v2#A6.F4.sf6)), $\mathsf{SPRINT}$ constantly has much smaller squared gradient norm, demonstrating the effectiveness of our algorithm. In Credit setting (Fig. [4(a)](https://arxiv.org/html/2509.17304v2#A6.F4.sf1) to Fig. [4(c)](https://arxiv.org/html/2509.17304v2#A6.F4.sf3)), $\mathsf{SPRINT}$ constantly has a smaller squared gradient norm when $\alpha=0.01$ and $0.4$. When $\alpha=0.2$, $\mathsf{SPRINT}$ has a much worse initial point than SGD-GD, resulting in the large squared gradient norm at first but similar norm at the end. Similarly, in CIFAR-10 setting (Fig. [4(g)](https://arxiv.org/html/2509.17304v2#A6.F4.sf7) to Fig. [4(i)](https://arxiv.org/html/2509.17304v2#A6.F4.sf9)), $\mathsf{SPRINT}$ has a smaller squared gradient norm when $\alpha=50$ and $80$. When $\alpha=20$, $\mathsf{SPRINT}$ has a much worse initial point than SGD-GD, resulting in the large squared gradient norm at first but similar norm at the end.

Note that Thm. [5.4](https://arxiv.org/html/2509.17304v2#S5.Thmtheorem4) is on the expectation of the squared gradient norm, but we are only able to conduct one experiment for each setting due to the limitation of computational resource. Thus, we can have some bad initialization point of $\mathsf{SPRINT}$ to influence the results.

## Appendix G Extend 𝖲𝖯𝖱𝖨𝖭𝖳 \mathsf{SPRINT} to Infinite Sum Settings

Our current finite-sum analysis avoids bounded variance by using a full gradient. For the population expectation minimization setting, we may have two plausible ways based on our current proof framework:

- 1.
We may use a minibatch of samples as the gradient snapshot, but we do need the bounded variance assumption to obtain SVRG-type rates.
- 2.
We can consider methods similar to Jothimurugesan et al. [2018], where the minibatch is growing and a control variate method is used to keep the gradient variance bounded. In this way, we do not need an explicit bounded variance assumption.

With the above methods and replacing the full gradients with mini-batch gradients, all other derivations should be similar to our current work. However, some nuanced constant parameters tuning may be needed to take care of the additional error produced by the mini-batch gradients.

## Appendix H Comparison of IFO Complexity between 𝖲𝖯𝖱𝖨𝖭𝖳 \mathsf{SPRINT} and SGD-GD

According to Theorem 1 of Li and Wai [2024], we can directly get the IFO Complexity of SGD-GD. When $n$ is not too large and when we require $\delta$ to be small, $\mathsf{SPRINT}$ will have much smaller complexity. Meanwhile, when $\sigma_{0}$ is large or even intractable (e.g., heavy-tailed/fat-tailed noise settings [Gurbuzbalaban et al., 2021]), SGD-GD can also have extremely large or even infinite IFO Complexity, while the complexity of $\mathsf{SPRINT}$ is independent of the variance parameter $\sigma_{0}$.

## Peer Reviews

confidence: 4
contribution: 2
initial_rating: 4
presentation: 3
questions: No further questions.
soundness: 3
strengths: The paper is clearly written and well-structured. The theoretical results appear sound. The proofs are well organized and generally convincing.
summary: This work considers the setting of Performative Prediction (PP), a novel risk minimization framework proposed by [Perdomo et al, 2021]. Unlike the tradietional ERM setting, PP assumes that the underlying data distribution also depends on the model parameters. This line of work gained much attention in the recent years, and much effort is devoted to finding a performative stable solution using variants of SGD, under convex/nonconvex loss function with bounded variance.

The authors observe that the nonconvex analysis in [Li and Wai., 2024] suffers from a critical reliance on a growth condition for the stochastic noise and an $O(1/\sqrt{T}) $ rate. To overcome this limitation, they propose a new algorithm and incorporate the variance-reduction technique. The authors further show that their method achieves an $O(1/T)$ convergence rate without assuming the previous noise condition.
weaknesses: My concern majorly comes from a contribution perspective. The improvements appear to be restricted to a rather narrow setting of performative prediction. Moreover, the incorporation of variance reduction techniques follows a fairly standard approach. As a result, the paper does not seem to offer many new insights for the broader optimization community. Unless the authors can present more decisive evidence on the significance of their work, the paper in its current form falls short of the expectations for an ICLR publication.

---

confidence: 3
contribution: 3
initial_rating: 6
presentation: 3
questions: NA
soundness: 3
strengths: see Summary
summary: Strengths
1. Innovative integration of variance reduction into performative prediction: This paper is the first to systematically introduce variance reduction techniques into the performative prediction framework, extending the convergence theory to the non-convex setting.

2. Strong theoretical guarantees under nonconvex settings: Theoretical analysis establishes an O(1/T) convergence rate without assuming bounded gradient variance, representing a major improvement over prior stochastic PP methods which only achieved O(1/\sqrt{T}) rates with variance-dependent error bounds.

3. Empirical results align with theory: Experiments on both tabular (credit scoring) and vision (CIFAR-10, MNIST) datasets show consistent improvements in convergence speed and stability across varying performativity strengths, supporting the theoretical claims.

Weaknesses

1. Experimental scope is narrow: Evaluations are limited to small-scale networks and datasets. The method’s scalability and robustness under more realistic, high-dimensional or large-data performative environments remain untested.

2. Complexity results lack empirical confirmation: The complexity analysis is solid theoretically, but the paper does not include runtime or gradient-call comparisons to validate practical efficiency gains.
weaknesses: see Summary

---

confidence: 4
contribution: 2
initial_rating: 4
presentation: 2
questions: Your experiments report training accuracy. Is there a corresponding notion of test accuracy here for  performative prediction? I'm guessing there is no such notion. Which make me think that online learning is a better paradigm for this problem.



**Minor issues**

1. Assumption 1: You should state $\forall z \in \mathcal{Z}.$

2. Definition 1. Because you evaluate the full batch gradient on the same $\theta$ as the one from which $z$ is sampled, this condition is exactly the standard $\delta$ stationarity condition, why has it been renamed to "$\delta$-Stationary Performative Stable Solution"? I find this confusing.


3. Formal statement of assumptions. You need to be crystal clear about your assumptions. In a formal Lemma, to be clear, you should always include the assumptions in the statement. For instance, does each Lemma 5.1 and 5.3 use all of your assumptions or only a subset? These things have to be clear for a theory paper.

4.  On a subjective note here, the name "Performative Prediction" is a very unfortunate name. It does not give any hint that this is an optimization problem where the data distribution is parameter dependent. This naming will make it harder for researchers in RL/control/VI working on this exact problem to find these works.
soundness: 2
strengths: One clear strength I see is that, as compared to Li & Tai 2024, these new convergence results do not rely on an explicit bounded variance assumption. Though implicitly the variance is bounded by the remaining assumptions. In any case, they use fewer assumptions.
summary: The paper proposes using the SVRG method for performative prediction, that is, minimizing a loss where the data distribution depends on the parameters being optimized. By adapting the SVRG proof to this setting, the author can establish a 1/T rate to a stationary point. Their proof also drops a previously used "bounded variance assumption".
weaknesses: The main weakness I see, and the one that my current stance hinges on, is the contradiction between the finite sum structure, and having a parameter dependent distribution.

**Finite sum breaks the data dependent assumption**.  How can you assume finite sum and that $z_i$ depends on $\theta$? By assuming the full loss has a finite sum structure, you loose the functional dependency on $\theta$, and you can no longer differentiate through a discrete sample. That is, having
$$ \nabla \mathcal{J}(\theta;\theta)  = \frac{1}{n} \sum_{I=1}^n \nabla \ell (\theta; z_i)$$ you essentially loose the dependence between the samples $z_i$'s and $\theta$. If not, how exactly are you now parametrizing the sampling of $z_t \in \{z_1, \ldots, z_n\}$ from this finite set so that it depends on $\theta_t$? This is why the related fields I mention below all have to contend with true expectations (infinite population). I read your Appendix G that briefly mentions how to avoid this issue, but it is not enough. You must contend with infinite samples, otherwise there is a contradiction between finite sum and samples depending on $\theta$.



**Limited related work.** This problem as stated 

$$ \min\_{\theta \sim \mathcal{D}\_{\theta}} \mathbb{E}\_{z \sim \mathcal{D}\_{\theta}}[\ell(\theta, z)] $$

Is verbatim the same problem stated in numerous fields, including RL, Optimal control, Variational Inference, and more. It has been an ongoing challenge in all of these fields to establish a non-asymptotic converges of SGD and variance reduced versions. These different fields have to deal with the exact same issues you address, such as having biased estimates of the gradient, assuming smoothness over a support of $z$, and analogous variance estimates. Because of this, you should extend your background section beyond predictive performance, and explain to the reader why this hasn't been done in RL/Control/VI before. A few standard references in convergence of SGD in RL below, but you should cast a wider net:

Lin Xiao, On the Convergence Rates of Policy Gradient Methods,  Journal of Machine Learning Research 23 (2022) 1-36

Alekh Agarwal, Sham M. Kakade, Jason D. Lee, and Gaurav Mahajan. On the theory of policy gradient methods: optimality, approximation, and distribution shift. J. Mach. Learn. Res. 22, 1, Article 98 (January 2021), 76 pages.

Rui Yuan, Robert M. Gower Alessandro Lazaric, A general sample complexity analysis of vanilla policy gradient, AISTATS 2022




**Variance independent bounds** This emphasis on not having "variance" appear in your bound is perhaps overplayed. This is because you assume the loss is $L_0$-Lipschitz, which bounds the gradient norm, and thus plays now the role of "variance" in your bounds. To be specific, using Lemma 5.2 from Li & Wai, assuming $L_0$-Lipschitz loss and $\epsilon$--sensitivity, you have immediately that Assumption 3.4 holds with $\sigma_0^2 = 2L_0^2 + 2\epsilon^2 L_0^2$ and $\sigma_1^2 =0.$ As such, I think this emphasis on not having dependency on variance is a bit confusing. A stronger message is that you simply rely on less assumptions (you do not explicitly assume  Assumption 3.4.
*Questions*

---

confidence: 4
contribution: 2
initial_rating: 2
presentation: 2
questions: see  Weaknesses
soundness: 2
strengths: 1. The method converges faster and more stably than SGD-GD, reaching a final neighborhood that does not depend on gradient variance.

2. It operates under weaker assumptions by avoiding the bounded-variance requirement and explicitly addressing deployment-induced distribution shift.

3. It delivers stronger empirical performance, achieving lower loss and steadier training on Credit, MNIST, and CIFAR-10, with especially large gains when performative effects are strong.
summary: This paper proposes SPRINT for performative prediction: an SVRG-style, snapshot-based variance-reduction method that tames noise from deployment-induced distribution shifts, proves fast and stable convergence without bounded-variance assumptions, and empirically outperforms SGD-GD on Credit, MNIST, and CIFAR-10—especially when feedback effects are strong.
weaknesses: Equation numbering is inconsistent (two equations labeled as Eq. 14).

The roles of $c_k$, $\beta_k$, and $\gamma_k$ are not clearly explained. Their definitions and interactions in the update rule should be clarified.

It is unclear whether Assumption 3.3 could be relaxed; this assumption seems restrictive for broader stochastic settings.

The paper should also discuss scenarios where SPRINT might underperform relative to SGD-GD (e.g., small population, low-variance cases) and analyze why.

Although hyperparameters are reported, the lack of details on the hyperparameter search process raises fairness concerns for baseline comparison.

The claimed “variance reduction” effect is not empirically validated—there is no report of gradient variance, convergence variance, or ablation removing the VR component.

 The experiments only report training loss; test accuracy should also be presented, especially given SPRINT’s claimed robustness under data distribution shifts.

Quantitative data on computational and memory costs are missing.

Experiments rely on small MLP and CNN architectures; results on more realistic models such as ResNet or VGG would strengthen the paper.

 The sensitivity of SPRINT to the learning rate remains unclear—this should be evaluated experimentally.

 It is also important to clarify whether SPRINT introduces additional computational or memory overhead compared to SGD-GC.


Decision: Reject