This paper studies the performative prediction problem which optimizes a stochastic loss function with data distribution that depends on the decision variable. We consider a setting where the agent(s) provides samples adapted to the learner's and agent's previous states. The said samples are used by the learner to optimize a loss function. This closed loop algorithm is studied as a state-dependent stochastic approximation (SA) algorithm, where we show that it finds a fixed point known as the performative stable solution. Our setting models the unforgetful nature and the reliance on past experiences of agent(s). Our contributions are three-fold. First, we demonstrate that the SA algorithm can be modeled with biased stochastic gradients driven by a controlled Markov chain (MC) whose transition probability is adapted to the learner's state. Second, we present a novel finite-time performance analysis of the state-dependent SA algorithm. We show that the expected squared distance to the performative stable solution decreases as ${\cal O}(1/k)$, where $k$ is the iteration number. Third, numerical experiments are conducted to verify our findings.

§ INTRODUCTION

Many supervised learning algorithms are built around the assumption that learners can obtain samples from a static distribution independent of the state of the learner and/or the agent who provides the sample. This assumption is reasonable for static tasks such as image classification. Oftentimes, it simplifies the design and analysis of algorithms such as stochastic gradient methods.

On the other hand, in certain applications agents can be performative where the samples are drawn from a decision-dependent distribution. This is relevant to the framework of strategic classification [Hardt et al., 2016, Cai et al., 2015, Kleinberg and Raghavan, 2020]. For instance, when training a classifier for loan applications, given the classifier published by the learner (bank), as the best response strategy, the agent(s) (loan applicants) may manipulate their profile prior to the submission, e.g., by spending more with credit cards, making unnatural transactions, etc., in order to increase their chance of successful application. The latter manipulation effectively shifts the data distribution and may affect the convergence properties, or even the stability of learning algorithms.

Earlier works  [Bartlett, 1992, Quiñonero-Candela et al., 2009] studied the effects of exogenous changes with shifts in data distribution.
Recently, [Perdomo et al., 2020] considered the convergence properties of learning algorithms when they are agnostic to the shifted distributions. Specifically, the learner is interested in the following performative prediction problem:

min_ θ∈^d    V(θ) = _ z ∼D(θ)  [ ℓ( θ; z ) ],
$\ell( \theta; z )$ is the loss function given the sample $z \in {\sf Z}$. The loss function is strongly-convex with respect to (w.r.t.) the parameter $\theta \in \RR^d$, and the gradient map $\grd_\theta \ell( \theta; z )$ is Lipschitz continuous w.r.t. $z$, $\theta$. In addition, the distribution ${\cal D}(\theta)$ on ${\sf Z}$ is parameterized by the decision vector $\theta$, which captures the distribution shift due to the learner's state. In other words, problem (<ref>) finds a parameter $\theta$ which minimizes the expected loss that takes care of the decision-dependent distribution.

Despite the strong convexity of $\ell(\theta ;z)$, problem (<ref>) is non-convex in general due to the coupling with $\theta$ in the data distribution ${\cal D}(\theta)$. As a remedy, [Perdomo et al., 2020] studied population-based algorithms that converge to a performative stable point, $\theta_{PS}$, which is a fixed point to the system $\theta = \argmin_{ \theta' \in \RR^d } \EE_{ z \sim {\cal D}(\theta) } \big[ \ell( \theta'; z ) \big]$.
Along the same line, [Mendler-Dünner et al., 2020] analyzed stochastic algorithms which deploy minibatches of i.i.d. samples from the shifted distribution at each iteration, [Izzo et al., 2021, Miller et al., 2021] studied gradient estimation techniques and developed algorithms that converge to an optimal solution of (<ref>) through introducing a gradient correction term (also see [Munro, 2020] which has considered a related setting), [Drusvyatskiy and Xiao, 2020] studied the stability of proximal stochastic gradient algorithms (and their variants) in the performative prediction setting, [Brown et al., 2020] studied population-based algorithms where the state-dependent distribution is updated iteratively.

This paper studies the convergence of stochastic algorithms where only one sample (or a minibatch of samples) is required at each iteration. Specifically, we consider a greedy deployment scheme similar to [Mendler-Dünner et al., 2020] where the learner deploys the most recent model after each update round. Moreover, the agent(s) is modeled with a memory property such that the update of his or her state depends on the past state. The closed-loop algorithm can be studied as a state-dependent stochastic approximation (SA) algorithm. In contrast to the setting analyzed in [Mendler-Dünner et al., 2020, Drusvyatskiy and Xiao, 2020] where only the learner's state is incrementally updated and the agent draws i.i.d. samples from the distribution shifted by the learner's state, the agent's state evolves according to a controlled Markov chain (MC) whose stationary distribution is the shifted distribution.

Our study is motivated by the stateful (or unforgetful) nature of the agents who depend on past experiences when adapting to a shifted target data distribution. For example, a loan applicant may take months to build up his/her credit history to adapt to changes in the published classifier. Several questions naturally arise from such dynamical performative prediction problems: will the stochastic algorithm converge to a performative stable point similar to [Mendler-Dünner et al., 2020, Drusvyatskiy and Xiao, 2020]? what is the sample complexity? This paper addresses these questions as we make the following contributions:

* We develop a fully state-dependent performative prediction framework which extends the analysis in [Mendler-Dünner et al., 2020, Drusvyatskiy and Xiao, 2020]. The proposed extension relies on a state-dependent stochastic approximation (SA) algorithm with noise originating from a controlled Markov chain [cf. Algorithm 1].

* Our main result consists of a finite-time convergence analysis of the state-dependent SA algorithm under a setting which does not assume the iterates to be bounded a-priori. Previous works either assumed the latter condition a-priori (e.g., [Benveniste et al., 2012]), or they require a compact constraint set (e.g., [Atchadé et al., 2017]). Using a novel analysis, we show that the mean squared error between the SA iterates and the unique performative stable solution [cf. (<ref>)] converges at a rate of ${\cal O}(1/k)$ in expectation. Additionally, we discuss about the convergence to an approximate stationary point of (<ref>) when the loss function $\ell(\theta;z)$ is not strongly convex (possibly non-convex).

* We demonstrate the efficacy of the SA algorithm with several experiments. We show that it has a comparable performance as in [Mendler-Dünner et al., 2020] which assumes an ideal setting with i.i.d. samples taken from the shifted distribution.

The rest of this paper is organized as follows. <ref> formally describes the performative prediction problem and a state-dependent SA algorithm for tackling the problem, <ref> presents the main theoretical results for the convergence of the state dependent SA algorithm, <ref> gives an overview of the proof strategy, and <ref> presents the numerical experiments.

Related Works Analysis for state-dependent stochastic approximation (SA) algorithms with controlled MC, which extend over the classical SA [Robbins and Monro, 1951], has been considered in a number of works.  Benveniste et al., 2012, Kushner and Yin, 2003 studied the asymptotic convergence of such algorithms, also see [Tadić et al., 2017] which analyzed the case of biased SA.

Recent works have analyzed the finite-time performance of state-dependent SA algorithms that are related to ours. Atchadé et al., 2017 considered a proximal SA algorithm where the proximal function has a compact domain; Karimi et al., 2019 analyzed the plain SA algorithm without projection but assumed that the updates are bounded; Sun et al., 2018, Doan et al., 2020 studied SA algorithms with a static MC not suitable for performative prediction.
Our analysis relaxes the restrictions of these prior works and focuses on convergence to a performative stable solution unique to (<ref>) [cf. see (<ref>)].

Lastly, our analysis technique is related to the recent endeavors on obtaining finite time bounds for reinforcement learning (RL) algorithms. Notice that  (<ref>) can in fact be regarded as a special case of policy optimization [Sutton and Barto, 2018]. To this end, recent works [Wu et al., 2020, Xu et al., 2020, Zhang et al., 2020]  studied the sample complexity of actor critic algorithms in finding a stationary point of an average reward function, where  controlled MCs are considered. In comparison to our analysis, the latter works have only studied the convergence to a stationary point of a simple cost function.

Notations We denote ${\sf Z}$ as the state space of samples provided by the agent(s) and ${\cal Z}$ is a $\sigma$-algebra on ${\sf Z}$. A Markov transition kernel is a map given by $\MK : {\sf Z} \times {\cal Z} \rightarrow \RR_+$. At the state $z \in {\sf Z}$, the next state is drawn as $z' \sim \MK( z, \cdot )$. It holds for any measurable function $f : {\sf Z} \rightarrow \RR$ that $\EE[ f(z') | z ] = \int_{\sf Z} f(z') \MK( z, {\rm d}\!~z') =: \MK f(z)$. Unless otherwise specified, the operator $\grd$ takes the gradient of a function w.r.t. the first argument for a multivariate function, e.g., $\grd \ell( \theta; z )$ denotes the gradient of $\ell(\theta; z)$ taken w.r.t. $\theta$. For $x, y \in \RR^d$, we denote the inner product as $\pscal{x}{y} = x^\top y$.

§ STATE-DEPENDENT PERFORMATIVE PREDICTION

Due to the intractability of Problem (<ref>), we are interested in evaluating the performative stable (PS) solution:

θ_PS = _ θ∈^d   _ z ∼D(θ_PS) [ ℓ( θ; z ) ].
Note that the expectation is taken with respect to $z \sim {\cal D}( \theta_{PS})$. It is known that $\theta_{PS}$ is in general different from an optimal or stationary solution to (<ref>); see <cit.>. Instead, $\theta_{PS}$ is a fixed point solution to the procedure when the learner repeatedly update the parameter $\theta$ with the drifted data distribution provided by the agent.

We consider a state dependent stochastic approximation (SA) algorithm motivated by the stateful nature of agents.
The latter is modeled using a controlled Markov chain.
For any $\theta$, we define a Markov transition kernel $\MK_{\theta}$ which induces a Markov chain with a unique stationary distribution $\pi_{\theta}(\cdot)$ such that $\EE_{z \sim \pi_{\theta}(\cdot)} [ \grd \ell(\theta; z) ] = \EE_{z' \sim {\cal D}(\theta)} [ \grd \ell( \theta ; z' )]$.

The learner and agent interact through the following rules:

Algorithm 1: State-dependent SA

Input: initialization $\theta_0$, step sizes $\{ \gamma_k \}_{k \geq 0}$.

For $k=0,1,2, \ldots$

Agent draws $z_{k+1} = {\sf SAMPLE}( \theta_k, z_k )$.

Learner updates

θ_k+1=θ_k-γ_k+1 ℓ(θ_k; z_k+1),

and deploys $\theta_{k+1}$.

${\sf SAMPLE}(\theta , z )$:

Draw the next sample as

z' ∼_ θ ( z, ·),

where $\MK_{\theta} : {\sf Z} \times {\cal Z} \rightarrow \RR_+$ is a Markov transition kernel.

Output: $z' \in {\sf Z}$.

We observe that (<ref>) is a standard SA recursion based on $\grd \ell(\theta_{k}; z_{k+1})$, where the learner deploys the most recent model $\theta_k$ and takes the sample $z_{k+1}$ directly from the agent.

Specifically, we consider a state-dependent setting where the sampling of $z_{k+1}$ can be affected by both the learner's and agent's states. Formally, the state-dependency is captured by modeling the samples sequence $\{z_k\}_{k \geq 1}$ as a controlled MC in (<ref>).
Notice that the stationary distribution $\pi_{\theta}(\cdot)$ does not need to be the same as ${\cal D}(\theta)$ as long as the former yields an unbiased gradient. However, for simplicity, we assume $\pi_{\theta}(\cdot) \equiv {\cal D}(\theta)$ for any $\theta \in \RR^d$.

The greedy deployment scheme studied by [Mendler-Dünner et al., 2020] assumed that the agents draw $z_{k+1}  {\sim} {\cal D}( \theta_k )$ as independent samples. The latter implies that $\grd \ell( \theta_k; z_{k+1})$ is an unbiased estimator of $\EE_{z' \sim {\cal D}(\theta_k)} [ \grd \ell( \theta_k; z' )]$. As a significant departure, in Algorithm 1, the stochastic gradient $\grd \ell( \theta_k; z_{k+1})$ is a biased estimator for $\EE_{z' \sim {\cal D}(\theta_k)} [ \grd \ell( \theta_k; z' )]$.
For any $\theta \in \RR^d$, we observe

[ ℓ( θ; z_k+1 )  |  z_k ] = _θ_k ℓ( θ; z_k ) = ∫_Z ℓ( θ; z ) _θ_k ( z_k, d   z ) .

Since $\MK_{\theta_k} ( z_k, \cdot ) \neq \pi_{\theta_k}(\cdot)$, we have $\EE[ \grd \ell ( \theta; z_{k+1} ) \!~|\!~ z_k ] \neq \EE_{z' \sim {\cal D}(\theta_k)} [ \grd \ell( \theta_k; z' )]$.
Under the setting (<ref>) with restricted access to the shifted data distribution ${\cal D}(\theta_k)$, one possibility to obtain an unbiased gradient estimate

is to hold $\theta_k$ as fixed and repeat the sampling process $z' = {\sf SAMPLE}(\theta_k , z)$ indefinitely. In this case, we have the unbiased estimate $\lim_{ n \rightarrow \infty } \MK_{\theta_k}^n {\grd \ell}( \theta ; z ) = \EE_{z' \sim {\cal D}(\theta_k)} [ \grd \ell(\theta; z') ]$ for any initial agent state $z \in {\sf Z}$.
Instead of searching for an unbiased gradient estimator through repeated MC updates, the state-dependent SA algorithm uses the instantaneous samples from (<ref>) that are co-evolving with the learner's iterate in (<ref>).

§.§ Example of Controlled Markov Chain (<ref>)

The model (<ref>) captures the stateful and stochastic nature of the agent as the sample $z_{k+1}$ depends on the previous one $z_k$. Concretely, our study is motivated by the following application example which satisfies $\pi_{\theta}(\cdot) \equiv {\cal D}(\theta)$ (a case where $\pi_{\theta}(\cdot) \neq {\cal D}(\theta)$ will be discussed in Appendix <ref>):

[Strategic Classification with Adapted Best Response]  We consider the problem of strategic classification [Hardt et al., 2016] involving some agents and a learner. In a typical scenario, the agent provides the best-response (i.e., optimized) samples upon knowing the current learner's state $\theta_k$. Ideally, the sample $z_{k+1} \sim {\cal D}(\theta_k)$ shall be drawn as

z_k+1 ∈_ z' ∈Z  U( z'; z̃_k+1 , θ_k ),  z̃_k+1 ∼D_0,
where ${\cal D}_0$ is the base distribution and $U( z'; z, \theta )$ is strongly concave in $z'$ for any $(z,\theta)$. The best response $\max_{ z' \in {\sf Z} }U( z'; z , \theta )$ perturbs the base sample $z$ in favor of the agent.

In practice, the exact maximization in (<ref>) is not obtained unless the agent(s) are given sufficient time to respond to the learner's state. Instead, we consider a setting where the agent(s) improve their responses via a gradient ascent dynamics evolving simultaneously with the learner.

Concretely, consider a setting where the learning problem (<ref>) utilizes data provided by $m$ agents.
Let ${\cal D}_0$ be the empirical distribution of $m$ data points $\{ \bar{d}_1,..., \bar{d}_m \}$, where $\bar{d}_i \in {\sf Z}$ is the initial data held by agent $i$.

At iteration $k$, a subset of agents ${\cal I}_k \subset \{1,...,m\}$ (selected uniformly with $|{\cal I}_k| = pm$, $p \in (0,1]$) becomes aware of the learner's state and they search for the best response through a gradient descent update. Then, the learner selects uniformly an agent $i_k \in \{1,...,m\}$ and requests the data sample $z_{k+1}$ from him/her.
In summary, the inexact best response dynamics follows:

Step 1:  d_i^k+1 =

d_i^k + αU( d_i^k; d̅_i, θ_k ),     i ∈I_k,

d_i^k,     i ∉I_k,

Step 2:  z_k+1 = d_i_k^k+1 ,

where the gradient is taken w.r.t. the first argument of $U( \cdot )$, $\alpha > 0$ is the agents' response rate (stepsize).
For the initialization, we set $d_i^0 = \bar{d}_i$ for all $i=1,...,m$.
The above dynamics highlights the stateful nature of the agents as they improve the responses based on their past experiences.

The best response dynamics executed by the agent(s) in (<ref>) leads naturally to a controlled MC (<ref>). Specifically, the MC's state is given by the tuple $\hat{z}_k = \{ d_1^k, ..., d_m^k , z_k \}$ and the application of the Markov kernel $\MK_{\theta_k}$ to $\hat{z}_k$ yields the inexact best response dynamics (<ref>). Furthermore, the latter admits a stationary distribution where
$\lim_{n \rightarrow \infty} \MK_{\theta_k}^n \grd \ell( \theta; \hat{z} ) = \EE_{ z \sim {\cal D}(\theta_k)} [ \grd \ell( \theta; z ) ]$.

We provide detailed properties about the controlled MC in Appendix <ref>. $\square$

The analysis of our state-dependent SA algorithm  (<ref>), (<ref>) entails unique challenges not found in the literature. First, the agent's states $\{ z_k \}_{k \geq 1}$ form a controlled MC whose transition probabilities are changing according to the learner's states $\{ \theta_k \}_{k \geq 0}$. Second, a plain SA update is used in (<ref>) which does not require a compact constraint set as in the prior works such as [Atchadé et al., 2017]. In fact, $\theta_k$ can be unbounded when the step size is not carefully selected.

§ MAIN RESULTS

Let us define the following shorthand notations
f( θ_1 ; θ_2 ) = _z ∼D(θ_2)  [ ℓ( θ_1; z) ],  f(θ_1;θ_2) = _ z ∼D(θ_2)  [ ℓ( θ_1 ; z) ],
where the first argument $\theta_1$ controls the loss function value and the second argument $\theta_2$ controls the distribution shift.
Notice that $\grd f(\theta_{PS}; \theta_{PS}) = 0$. We consider the following assumptions. First, the learner's loss is strongly convex in $\theta$, and its gradient map is Lipschitz continuous in $(\theta,z)$, i.e.,

For each $z \in {\sf Z}$, there exists $\mu>0$ such that
ℓ( θ; z ) ≥ℓ( θ'; z ) +  ℓ( θ'; z )  θ- θ'  + (μ/2)  θ' - θ^2,   ∀ θ, θ' ∈^d.

There exists $L \geq 0$ such that

ℓ( θ; z ) -  ℓ( θ' ; z' )  ≤L {  θ- θ'  +  z - z'  },   ∀ θ, θ' ∈^d, z,z' ∈Z.

Notice that as a consequence, the expected objective function $f(\theta_1, \theta_2)$ and gradient $\grd f(\theta_1; \theta_2)$ are $\mu$-strongly convex in $\theta_1$, and $L$-Lipschitz in $\theta_1$, respectively. These are standard assumptions in the optimization literature. As indicated by [Drusvyatskiy and Xiao, 2020], these conditions are necessary for finding a performative stable solution in (<ref>).

Second, we have the following assumption on the oscillation of the stochastic gradient $\grd \ell(\theta;z)$:
There exists $\sigma \geq 0$ such that
sup_ z ∈Z   ℓ( θ; z ) - f( θ; θ_PS )  ≤σ( 1 +  θ- θ_PS  ),   ∀ θ∈^d.

The above is slightly stronger than the assumptions on second order moments typically found in the stochastic gradient literature, e.g., [Bottou et al., 2018], as we require a uniform bound on the gradient noise. This condition is common for the algorithms using Markovian samples [Sun et al., 2018, Srikant and Ying, 2019, Karimi et al., 2019], which requires that the oscillation of stochastic gradient is controlled. For strategic classicaition problems, it is satisfied for the finite dataset setting in <Ref>.
Moreover, similar to [Doan et al., 2020], this bound is adapted to the growth of $\| \theta - \theta_{PS} \|$ which is compatible with the strong convexity of the loss function $\ell(\theta; z)$.

Our next set of assumptions pertain to the Markov kernels $\MK_{\theta} $ that generate $\{ z_k \}_{k \geq 1}$:

There exists a solution $\widehat{\grd \ell}: \RR^d \times {\sf Z} \rightarrow \RR^d$ to the Poisson equation:

ℓ( θ'; z) - f( θ'; θ) = ℓ (θ'; z) - _θ ℓ ( θ'; z),   ∀  θ, θ' ∈^d, z ∈Z.

Consider the Poisson equation's solution $\widehat{\grd \ell} (\cdot; \cdot )$ defined in <Ref>. There exists $\Lph \geq 0$ such that

sup_ z ∈Z   _θ ℓ( θ; z) - _θ' ℓ (θ'; z)  ≤ θ- θ' ,   ∀ θ, θ' ∈^d.

For Assumption <ref>, the existence of $\widehat{\grd \ell}$ in (<ref>) holds under mild assumptions on the Markov chains (MCs). For instance, it holds when the MC is irreducible and aperiodic, and satisfying a Lyapunov drift condition, or in the simpler case, when the MC is uniform geometrically ergodic, see <cit.>.
In the case for performative prediction with stateful agents, the above condition holds when repeated applications of the iterative map (<ref>) where agents adopt their data to the current learner's model $\theta$ converges linearly to the best response.
Assumption <ref> is a smoothness condition on the kernel $\MK_\theta$ with respect to $\theta$. The condition can be satisfied when the Markov kernel is only slightly modified when $\theta$ is perturbed.

<Ref> is also linked to our next assumption which is central to the study of performative prediction. Particularly, we require the distribution map $\cal D(\theta)$ to be $\epsilon$-sensitive w.r.t. $\theta$:
There exists $\epsilon \geq 0$ such that
W_1( D(θ), D(θ') ) =
inf_ J ∈J( D(θ), D(θ') )
_ (z,z') ∼J  [  z - z' _1 ]

≤ϵ θ- θ' ,
for any $\theta, \theta' \in \RR^d$. Notice that $W_1(\cdot)$ denotes the Wasserstein-1 distance and ${\cal J}( {\cal D}(\theta), {\cal D}(\theta') )$ is the set of all joint distributions on ${\sf Z} \times {\sf Z}$ with ${\cal D}(\theta), {\cal D}(\theta')$ as its marginal distribution.

The above is a common condition for performative prediction, e.g., [Perdomo et al., 2020]. Intuitively, it allows for performative prediction algorithm to behave stably as the perturbation to ${\cal D}(\theta)$ is under control. In the subsequent analysis, we demonstrate that carefully controlling the step size in relation to $\mu, \epsilon, L$ is crucial to the convergence of Algorithm 1.

Before presenting our main result, we notice that Assumptions <ref>, <ref>, <ref> imply that there exists constants $\overline{L}$, $\widehat{L} > 0$ such that for any $z \in {\sf Z}$, $\theta \in \RR^d$,

ℓ( θ; z)  ≤L (1 +  θ- θ_PS  ),   max{  ℓ( θ; z)  , _θ ℓ( θ; z )  } ≤(1 +  θ- θ_PS  ),
In other words, $\grd \ell(\theta; z)$, $\widehat{\grd \ell}(\theta; z)$, $\MK_\theta \widehat{\grd\ell}( \theta; z )$ are all locally bounded functions.
Notice that $\overline{L}$ is proportional to $\sigma$ in Assumption <ref>, while $\widehat{L}$ is proportional to the maximum mixing time of the Markov chain induced by the kernel $\MK_\theta$ over all $\theta \in \RR^d$.

Our main result for the state-dependent SA algorithm is summarized as follows:

Under Assumptions <ref>–<ref>. Suppose that the problem parameters satisfy $\epsilon < \frac{\mu}{L}$, the step sizes  $\{\gamma_{k}\}_{k\geq 1}$ are non-increasing and satisfy for any $k \geq 1$,

γ_k-1/γ_k≤1+  γ_k ( μ - L ϵ)/ 4 ,    γ_k ≤min{ μ- L ϵ/2L^2,  μ- L ϵ/2 , min{ (μ- L ϵ)/3, 3 } / + 3 ( μ- L ϵ) , 1/6  }.
Then for any $k \geq 1$, the expected distance between $\theta_{k}$ and the performative stable solution $\theta_{PS}$ satisfies

[  θ_k - θ_PS ^2 ] ≤∏_i=1^k ( 1 - γ_i μ- L ϵ/2 )  θ_0 - θ_PS ^2 + C   γ_k ,
where $\EE[\cdot]$ is the expectation taken over all the randomness in (<ref>), (<ref>), and we have defined:

C :=
3 Δ  + 4 ς/ μ- L ϵ  ( 2 (2 σ^2 + ) + (μ- L ϵ) + ( + 3 (μ- L ϵ) ) Δ ) ,

with $\varsigma := 1 + \gamma_1 ( \mu - L \epsilon ) / 4$, and $\Cone, \Ctwo, \Cthree$, $\overline{\Delta}$ are constants defined in (<ref>), (<ref>), respectively.

The above establishes the finite-time convergence of the studied state-dependent SA algorithm (<ref>), (<ref>). To understand this result, we observe that the step size conditions in (<ref>) can be satisfied by a variety of step size schedules. For instance, it can be satisfied by the constant step size $\gamma_k \equiv \gamma$; and the diminishing step size $\gamma_{k}=a_{0}/(a_{1}+k) = {\cal O}(1/k)$ with appropriate $a_0, a_1 > 0$.
Moreover, we require the SA algorithm to work in the regime when $\epsilon < \mu / L$. This is similar to [Perdomo et al., 2020] which ensures that the solution $\theta_{PS}$ is a stable fixed point to (<ref>).

The main result is stated on the expected squared distance of $\| \theta_k - \theta_{PS} \|^2$ in (<ref>). Here, the bound consists of a transient term and a fluctuation term. The transient term decays sub-exponentially as ${\cal O}( \exp( - \frac{\mu-L\epsilon}{2} \sum_{i=1}^k \gamma_i ) )$ and is scaled by the initial gap $\| \theta_0 - \theta_{PS} \|^2$. The fluctuation term is in the order of ${\cal O}(\gamma_k)$ and is scaled by $\mathds{C}$ which depends on the oscillation of stochastic gradient (via $\sigma$) and the mixing time of the controlled Markov chain (via $\LZ$).
With a diminishing step size schedule such as $\gamma_{k}= c_{0}/( c_{1}+k)$, <Ref> shows that the state dependent SA algorithm finds the performative stable solution $\theta_{PS}$ at the rate of ${\cal O}(1/k)$ in expectation.

On Non-strongly-convex Loss Function An obvious shortcoming with <Ref> is the requirement of strongly convex loss functions [cf. <Ref>]. Below, we comment on the convergence of the state-dependent SA algorithm (<ref>), (<ref>) when the loss function is not strongly convex (possibly non-convex). Notice that in the absence of <Ref>, the performative stable solution $\theta_{PS}$ may not be well defined. We resort to finding a stationary point to the performative prediction problem (<ref>).

Our idea is to view (<ref>), (<ref>) as a biased SA algorithm with mean field $h(\theta) = \EE_{z \sim {\cal D}(\theta)} [ \grd \ell( \theta; z )] = \grd f(\theta; \theta)$. This mean field is correlated with the gradient for the performative loss in (<ref>). Under additional assumptions on $\ell(\theta;z)$, ${\cal D}(\theta)$, in Appendix <ref> we show

h(θ)  _ z ∼D(θ)  [ ℓ( θ; z ) ]  ≥ h(θ) ^2 / 2 - c_0 ,  ∀ θ∈^d,
where ${\rm c}_0 > 0$ is a bias dependent on the sensitivity of distribution shift [cf. <Ref>].
We define the constant ${\rm c}_0$ in (<ref>) of Appendix <ref>, which is shown to be dependent on $\ell(\theta; z)$,  $\grd_{\theta} \log(p_{\mathcal{D}(\theta)}(z))$.

It should be pointed that the algorithm (<ref>) may not provide a `good' solution to the non-convex performative learning problem (<ref>). However, if the  state-dependent distribution is not too sensitive to the change of state, the following corollary shows that (<ref>) would still converge to a ${\rm c}_{0}$-neighborhood of a stationary solution. Before we discuss the main statement, we need two additional assumptions:
The function $V(\theta) = \EE_{ z \sim {\cal D}(\theta) } [ \ell( \theta; z ) ]$ is continuously differentiable and there exists $L_V \geq 0$ such that $\| \grd V(\theta) - \grd V(\theta') \| \leq L_V \| \theta - \theta' \|$ for any $\theta, \theta' \in \RR^d$.

There exists $\sigma \geq 0$ such that $\| \grd \ell(\theta; z) - \grd f(\theta; \theta) \| \leq \sigma$ for any $\theta \in \RR^d, z \in {\sf Z}$.

The above are stronger conditions than Assumptions <ref>, <ref>, yet are reasonable settings for certain non-convex loss functions.
For instance, as we will show in Appendix <ref>, Assumption <ref> holds if $\ell(\theta ; z), \nabla_{\theta} \log \left(p_{D(\theta)}(z)\right)$ are bounded, and $\nabla_{\theta} \log \left(p_{D(\theta)}(z)\right)$ to be Lipschitz w.r.t. $\theta$, e.g., when $\mathcal{D}(\theta)$ is `smooth', e.g., it is Gaussian or softmax distribution.
Assumption <ref> holds under similar condition as Assumption <ref>, e.g., ${\sf Z}$ is compact.

We obtain the following as a corollary of <cit.>:

Under Assumptions <ref>, <ref>, <ref>, <ref>, and let  (<ref>) holds. With a step size sequence that decays in the order of $\gamma_k = {\cal O}(1/\sqrt{k})$, it holds for any $K \geq 1$ that
[  V( θ_K ) ^2 ] = O ( logK / √(K) + c_0 ),  ( K = k ) = γ_k / ∑_j=1^K γ_j, ∀ k ∈{1,...,K}.
Note that ${\sf K} \in \{1,...,K\}$ is a discrete r.v. independent of the randomness in the SA algorithm and $\EE[\cdot]$ denotes the total expectation.

See the details in Appendix <ref>.
The above corollary shows that even without the strong convexity assumption, the state dependent SA algorithm finds an ${\cal O}(\log K / \sqrt{K} + {\rm c}_0)$-stationary solution to the performative prediction problem (<ref>) in at most $K$ iterations.

§ PROOF OUTLINE OF <REF>

We outline the main steps in proving <Ref>. Our proof strategy consists in tracking the progress of the mean squared error $\Delta_k \eqdef \EE[ \| \theta_k - \theta_{PS} \|^2 ]$.
To simplify notations, we define $\tilde{\mu} \eqdef \mu - L \epsilon$, the scalar product $G_{m:n} = \prod_{i=m}^n (1 - \gamma_i \tilde{\mu})$, for $n > m \geq 1$, and $G_{m:n} = 1$ if $n \leq m$.

The following lemma describes the one-step progress of the SA algorithm.

Under Assumptions <ref>, <ref>, <ref>, <ref>. For any $k \geq 0$, it holds

θ_k+1 - θ_PS ^2     ≤( 1 - 2 γ_k+1 μ̃ + 2L^2 γ_k+1^2 )  θ_k - θ_PS ^2

+ 2 σ^2 γ_k+1^2 - 2 γ_k+1  θ_k - θ_PS  ℓ( θ_k; z_k+1 ) - f( θ_k; θ_k )  .

The proof can be found in Appendix <ref>, which involves a simple expansion of the squared error.
The above lemma suggests that the sensitivity parameter shall satisfy $\epsilon < {\mu} / {L}$ to ensure $\tilde{\mu} > 0$. Furthermore, the step size condition
$\sup_{k \geq 1} \gamma_{k} \leq {\tilde{\mu}} / (2 L^2)$ in (<ref>) leads to $1 - 2 \gamma_{k+1} \tilde{\mu} + 2L^2 \gamma_{k+1}^2 \leq 1 - \gamma_{k+1} \tilde{\mu}$ such that the first term in the r.h.s. of (<ref>) is a contraction.
Under the above premises and suppose $z_{k+1} \overset{\text{i.i.d.}}{\sim} {\cal D}( \theta_k )$ as in the case of [Mendler-Dünner et al., 2020, Drusvyatskiy and Xiao, 2020], the stochastic gradient in (<ref>) is conditionally unbiased. As such, <Ref> leads to the recursion $\Delta_{k+1} \leq ( 1 - \gamma_{k+1} \tilde{\mu} ) \Delta_k + 2 \sigma^2 \gamma_{k+1}^2$, implying $\Delta_k = {\cal O}(\gamma_k)$.

However, for the state-dependent SA algorithm  (<ref>), (<ref>),  the stochastic gradient $\grd \ell( \theta_k; z_{k+1} )$ is conditionally biased and is driven by a controlled MC. Under the stepsize condition $\sup_{k \geq 1} \gamma_{k} \leq {\tilde{\mu}} / (2 L^2)$, taking the total expectation and solving the recursion in (<ref>) yield

Δ_k     ≤G_1:k Δ_0 + 2 σ^2 ∑_s=0^k-1 G_s+2:k γ_s+1^2

+ 2 ∑_s=0^k-1 G_s+2:k γ_s+1 [  θ_PS - θ_s  ℓ( θ_s; z_s+1) - f(θ_s; θ_s)   ] .

It can be shown that the first two terms are bounded by ${\cal O}(\gamma_k)$. We are interested in analyzing the last term when the samples $\{ z_{k} \}_{k \geq 1}$ are drawn according to (<ref>). Observe that:

Under Assumptions <ref>–<ref> and the stepsize conditions in (<ref>). For any $k \geq 1$, it holds

2 ∑_s=0^k-1 G_s+2:k γ_s+1 [  θ_PS - θ_s  ℓ( θ_s; z_s+1) - f(θ_s; θ_s)   ] ≤

∑_s=2^k γ_s^2 G_s+1:k ( + Δ_s-1 + Δ_s-2 ) + γ_1 G_2:k { (1 + 3 Δ_0 ) + γ_1 } + γ_k { 1 + 3 Δ_k-1 },

where we have defined the constants:

:= ςL + 4 ςL + ( 1 + μ̃ ) ς,   := 2 ςL,
:= + 2 ( 1 + μ̃ )ς.

The analysis is inspired by [Benveniste et al., 2012] and has been adopted in recent works such as [Atchadé et al., 2017, Karimi et al., 2019]; see the details in Appendix <ref>. To handle the controlled MC, our technique involves applying the Poisson equation in <Ref> and decomposing the gradient error $\grd \ell( \theta_s; z_{s+1}) - \grd f(\theta_s; \theta_s)$ into Martingale and finite difference terms.

A key difference between <Ref> and analysis in the previous works is that the latter assumed that the stochastic gradients $\grd \ell(\theta_k; z_{k+1} )$ are bounded which greatly simplifies the proof. Our assumptions are significantly weaker as the latter actually grows as ${\cal O}( 1 + \| \theta_k - \theta_{PS} \|)$ [cf. <Ref>].

From the analysis standpoint, this demands a new proof technique as we present next.

To proceed, observe that substituting <Ref> into (<ref>) yields:

Δ_k     ≤G_1:k Δ_0 + ∑_s=1^k-1 γ_s+1^2 G_s+2:k ( 2 σ^2 + + Δ_s + Δ_s-1 )

+ γ_1 G_2:k { (1 + 3 Δ_0 ) + γ_1 ( 2 σ^2 + ) } + γ_k { 1 + 3 Δ_k-1 }.

With the first step size condition in (<ref>), we can apply the auxiliary result in <Ref> from the appendix, which simplifies the upper bound as

Δ_k     ≤G_1:k Δ_0 + (  2/ μ̃  (2 σ^2 + ) + ) γ_k + ∑_s=1^k-1 γ_s+1^2 G_s+2:k ( Δ_s + Δ_s-1 )

+ γ_1 G_2:k { (1 + 3 Δ_0 ) + γ_1 ( 2 σ^2 + ) } + 3 γ_k Δ_k-1 .

Observe that the first row in (<ref>) is already in a similar form to the bound presented in the theorem. The key issue lies with the last term $3 \gamma_k \LZ \Delta_{k-1}$ which may be unbounded. We show that our choice of step sizes in (<ref>) ensures the convergence of $\Delta_k$ to ${\cal O}(\gamma_k)$:

Suppose that $\{ \Delta_k \}_{k \geq 0}$ satisfy (<ref>) and the step sizes $\{ \gamma_k \}_{k \geq 1}$ satisfy (<ref>). It holds (i)

sup_k ≥0 Δ_k ≤Δ := 3 Δ_0 + ς/ 9 ^2  ( 2 (2 σ^2 + ) + (μ- L ϵ) ),

and (ii) the following inequality holds for any $k \geq 1$:

Δ_k ≤∏_i=1^k ( 1 - γ_i μ̃/2) Δ_0 + { 3 Δ  + 4 ς/ μ̃  ( 2 (2 σ^2 + ) + μ̃ + ( + 3 μ̃ ) Δ ) } γ_k.

Proving the above lemma requires one to establish the stability of the system (<ref>), which demands a sufficiently small $\gamma_k$ to control the remainder term $3 \LZ \gamma_k \Delta_{k-1}$. Our analysis relies on the special structure of this inequality system; see the proof details in Appendix <ref>.
The convergence bound (<ref>) follows from the boundedness of $\Delta_k$.

Finally, we obtain <Ref> by applying <Ref>.

§ NUMERICAL EXPERIMENTS

This section considers two performative prediction problems to corroborate with our theories. All the experiments are performed with Python on a server using a single thread of an Intel Xeon 6138 CPU. Further details about the experiments below can be found in Appendix <ref>.

Gaussian Mean Estimation
The first problem is concerned with Gaussian mean estimation using synthetic data. Our aim is to validate Theorem <ref> using a simple experiment. Here, (<ref>) is specified as $\min_{ \theta \in \RR } \EE_{ z \sim {\cal D}(\theta)} [ (z-\theta)^2/2 ]$ with ${\cal D}(\theta) \equiv {\cal N}( \bar{z} + \epsilon \theta ; \sigma^2 )$. For $0<\epsilon<1$, the performative stable solution has a closed form $\theta_{PS} = \frac{ \bar{z} }{ 1 - \epsilon }$. For the state-dependent SA, the agent follows an autoregressive (AR) model $z_{k+1} = (1 - \rho) z_k + \rho \tilde{z}_{k+1}$ with independent $\tilde{z}_{k+1} \sim {\cal N}( \bar{z} + \epsilon \theta_k; \sigma^2 )$ and regression parameter $\rho \in (0,1)$. This AR recursion is a controlled MC with a stationary distribution that yields the unbiased gradient of (<ref>), details about the MC are in Appendix <ref>.

We consider a large variance setting with $\bar{z} = 10$, $\sigma = 50$, $\epsilon = 0.1$. The step size is $\gamma_k = \frac{c_0}{c_1+k}$, $c_0 = \frac{500}{\tilde{\mu}}, c_1 = \frac{800}{\tilde{\mu}^2}$. In Fig. <ref> (left), we compare $|\theta_k - \theta_{PS}|^2$ against the iteration number $k$ for the Gaussian estimation problem using our state-dependent SA and greedy deploy [Mendler-Dünner et al., 2020] algorithms. As observed, both algorithms have an asymptotic convergence rate of ${\cal O}(1/k)$ towards $\theta_{PS}$ which is predicted by <Ref>. As $\rho \downarrow 0$, the state-dependent SA algorithm delivers a smaller error as the AR model has a stationary distribution with lower variance.

Gaussian mean estimation – (Left) Under different regression parameter $\rho$; Strategic Classification – (Middle) Under Linear BR $U_{\sf q}(\cdot)$ and different agent response rate $\alpha$ [cf. (<ref>)]; (Right) Under Logistics BR $U_{\sf lg}(\cdot)$. The shaded region shows the 90% confidence interval over 20 trials.

Strategic Classification
The second problem is a strategic classification (SC) problem similar to  [Perdomo et al., 2020]

for a credit scoring classifier with  dataset[Available at <https://www.kaggle.com/c/GiveMeSomeCredit/data>.].
Our aim is to showcase the effects when the agent adapts slowly to the shifted distribution. In particular, our theory implies that while the algorithm will still converge to $\theta_{PS}$, a slower convergence rate will be observed. To specify (<ref>), let $z \equiv (x,y)$ where $x \in \RR^d$ is feature vector, $y \in \{0,1\}$ is label. The learner finds $\theta \in \RR^d$ that minimizes:

_z ∼D(θ)[ ℓ( θ; z) ],   where   ℓ( θ; z ) = β/2  θ^2 + log( 1 + exp(  θx ) ) - y θx.
Observe that when $\beta > 0$, $\ell(\theta;z)$ is a $\beta$-strongly convex function w.r.t. $\theta$ satisfying <Ref>.
For any $\theta \in \RR^d$, the shifted data distribution ${\cal D}(\theta)$ is obtained through evaluating the best response (BR) in (<ref>) of Example <ref>. We consider two types of utility functions adopted by the agents:
U_q(z^';z, θ) = θx'- x'-x^2/2ϵ,   U_lg(z^';z, θ) = y θx' - log(1+exp(θx'))- x'-x^2/2ϵ,
where $z \equiv (x,y)$ is the original unshifted data. The label $y \in \{0,1\}$ is unchanged in the BR. Notice that $U_{\sf q}(\cdot)$, $U_{\sf lg}(\cdot)$ have respectively linear and logistics costs. Both utility functions include a quadratic regularizer where $\epsilon$ controls the sensitivity of the distribution shift [cf. <Ref>].

With a published $\theta_k$, the agent(s) maximize the utility function prior to giving data to the learner for the next round.
For both $U_{\sf q}(\cdot)$ and $U_{\sf lg}(\cdot)$, the BR obtained steers the classifier in favor of the agent(s).
Furthermore,  $U_{\sf lg}(\cdot)$ is motivated by logistic regression which favors towards samples with label `1'.

For details of the numerical experiments, we set $\beta = 1000/m$ in (<ref>), $\epsilon = 0.01$ in the utility functions, and  in (<ref>), we set number of selected agents as $|{\cal I}_k| = 5$, agents' response rate (stepsize) as $\alpha = 0.5 \epsilon$ unless otherwise specified.
The step size for (<ref>) is $\gamma_k = {c_0} / {(c_1+k)}$, $c_0 = {100} / {\tilde{\mu}}, c_1 = { 8 L^2 } / {\tilde{\mu}^2}$.

We first consider when $\theta_{PS}$ is computed with ${\cal D}(\theta)$ defined by the linear BR function $U_{\sf q}(\cdot)$ and compare our state-dependent SA with greedy deploy scheme [Mendler-Dünner et al., 2020] and repeated risk minimization [Perdomo et al., 2020].  As shown in Fig. <ref> (middle), all algorithms converge to $\theta_{PS}$. As $\alpha \downarrow 0$, the state-dependent SA converges at slightly slower rates as the agents adapt to the distribution shift with increased mixing time of the MC.
The result corroborates with <Ref> which established the ${\cal O}(1/k)$ convergence rate with state-dependent SA.

Our last experiment pertains to the same SC problem as before, but we consider a different setting where $\theta_{PS}$ is computed  with ${\cal D}(\theta)$ defined by the logistics BR function $U_{\sf lg}(\cdot)$.
The agents follow a more complicated dynamics since the BR does not admit a closed form solution. Again, we aim to confirm Theorem <ref> that Algorithm 1 converges under the stateful agent setting. We compare the distance $\| \theta_k - \theta_{PS}\|^2$ versus iteration $k$.

In addition to showing the convergence result as predicted in Theorem <ref>, we aim to observe the effects when the learner adopts a lazy deployment scheme [cf. [Mendler-Dünner et al., 2020]], by varying the number of $\theta$-update by the learner per the agents' update.
Fig. <ref> (right) shows the error $\| \theta_k - \theta_{PS} \|^2$ against the number of adaptation steps performed at the agents via (<ref>) as we illustrate the convergence rate from the perspectives of the agents.
We observe that the error decreases at a faster rate when the number of learner's iteration increases. More details about our numerical experiments can be found in Appendix <ref>.

Conclusions & Limitations    We consider a state-dependent SA algorithm for performative prediction. We showed a convergence rate of ${\cal O}(1/k)$ in mean-squared error towards the performative stable solution when the agents provide data drawn from a controlled MC. Our study paved the first step towards understanding and applying performative prediction in a dynamical setting.

There are several limitations. First, despite being a reasonable fixed point, $\theta_{PS}$ can be different from the optimal solution to (<ref>). An open problem is whether the state dependent SA converges when gradient correction, e.g., [Izzo et al., 2021], is used. Second, our convergence analysis of state-dependent SA algorithms relies on smoothness conditions on the controlled MC and distributions [cf. <Ref>, <ref>]. An open problem is to verify if these are necessary for our results to hold.

[Atchadé et al., 2017]
Y. F. Atchadé, G. Fort, and E. Moulines.
On perturbed proximal gradient algorithms.
The Journal of Machine Learning Research, 180
(1):0 310–342, 2017.

[Bartlett, 1992]
P. L. Bartlett.
Learning with a slowly changing distribution.
In Proceedings of the fifth annual workshop on Computational
learning theory, pages 243–252, 1992.

[Benveniste et al., 2012]
A. Benveniste, M. Métivier, and P. Priouret.
Adaptive algorithms and stochastic approximations, volume 22.
Springer Science & Business Media, 2012.

[Borkar, 2021]
V. S. Borkar.
A concentration bound for contractive stochastic approximation.
Systems & Control Letters, 153:0 104947, 2021.

[Bottou et al., 2018]
L. Bottou, F. E. Curtis, and J. Nocedal.
Optimization methods for large-scale machine learning.
Siam Review, 600 (2):0 223–311, 2018.

[Brown et al., 2020]
G. Brown, S. Hod, and I. Kalemaj.
Performative prediction in a stateful world.
arXiv preprint arXiv:2011.03885, 2020.

[Cai et al., 2015]
Y. Cai, C. Daskalakis, and C. Papadimitriou.
Optimum statistical estimation with strategic data sources.
In Conference on Learning Theory, pages 280–296. PMLR, 2015.

[Doan et al., 2020]
T. T. Doan, L. M. Nguyen, N. H. Pham, and J. Romberg.
Finite-time analysis of stochastic gradient descent under markov
arXiv preprint arXiv:2003.10973, 2020.

[Douc et al., 2018]
R. Douc, E. Moulines, P. Priouret, and P. Soulier.
Markov chains.
Springer, 2018.

[Drusvyatskiy and Xiao, 2020]
D. Drusvyatskiy and L. Xiao.
Stochastic optimization with decision-dependent distributions.
arXiv preprint arXiv:2011.11173, 2020.

[Ginart et al., 2021]
T. Ginart, E. Zhang, Y. Kwon, and J. Zou.
Competing ai: How does competition feedback affect machine learning?
In International Conference on Artificial Intelligence and
Statistics, pages 1693–1701. PMLR, 2021.

[Hardt et al., 2016]
M. Hardt, N. Megiddo, C. Papadimitriou, and M. Wootters.
Strategic classification.
In Proceedings of the 2016 ACM conference on innovations in
theoretical computer science, pages 111–122, 2016.

[Izzo et al., 2021]
Z. Izzo, L. Ying, and J. Zou.
How to learn when data reacts to your model: Performative gradient
arXiv preprint arXiv:2102.07698, 2021.

[Kaggle, 2012]
Give me some credit, 2012.
URL <https://www.kaggle.com/c/GiveMeSomeCrdit/data>.

[Kaledin et al., 2020]
M. Kaledin, E. Moulines, A. Naumov, V. Tadic, and H.-T. Wai.
Finite time analysis of linear two-timescale stochastic approximation
with markovian noise.
In Conference on Learning Theory, 2020.

[Karimi et al., 2019]
B. Karimi, B. Miasojedow, E. Moulines, and H.-T. Wai.
Non-asymptotic analysis of biased stochastic approximation scheme.
In Conference on Learning Theory, pages 1944–1974. PMLR,

[Kleinberg and Raghavan, 2020]
J. Kleinberg and M. Raghavan.
How do classifiers induce agents to invest effort strategically?
ACM Transactions on Economics and Computation (TEAC),
80 (4):0 1–23, 2020.

[Kushner and Yin, 2003]
H. Kushner and G. G. Yin.
Stochastic approximation and recursive algorithms and
applications, volume 35.
Springer Science & Business Media, 2003.

[Mendler-Dünner et al., 2020]
C. Mendler-Dünner, J. Perdomo, T. Zrnic, and M. Hardt.
Stochastic optimization for performative prediction.
Advances in Neural Information Processing Systems, 33, 2020.

[Miller et al., 2021]
J. Miller, J. C. Perdomo, and T. Zrnic.
Outside the echo chamber: Optimizing the performative risk.
arXiv preprint arXiv:2102.08570, 2021.

[Munro, 2020]
E. Munro.
Learning to personalize treatments when agents are strategic.
arXiv preprint arXiv:2011.06528, 2020.

[Nesterov, 2012]
Y. Nesterov.
Efficiency of coordinate descent methods on huge-scale optimization
SIAM Journal on Optimization, 220 (2):0
341–362, 2012.

[Patrascu and Necoara, 2015]
A. Patrascu and I. Necoara.
Efficient random coordinate descent algorithms for large-scale
structured nonconvex optimization.
Journal of Global Optimization, 610 (1):0
19–46, 2015.

[Perdomo et al., 2020]
J. C. Perdomo, C. Zrnic, Tijana aCnd Mendler-Dünner, and M. Hardt.
Performative prediction.
In ICML, 2020.

[Qu and Wierman, 2020]
G. Qu and A. Wierman.
Finite-time analysis of asynchronous stochastic approximation and $ q
In Conference on Learning Theory, pages 3185–3205. PMLR,

[Quiñonero-Candela et al., 2009]
J. Quiñonero-Candela, M. Sugiyama, N. D. Lawrence, and A. Schwaighofer.
Dataset shift in machine learning.
Mit Press, 2009.

[Richtárik and Takáč, 2014]
P. Richtárik and M. Takáč.
Iteration complexity of randomized block-coordinate descent methods
for minimizing a composite function.
Mathematical Programming, 1440 (1):0 1–38,

[Robbins and Monro, 1951]
H. Robbins and S. Monro.
A stochastic approximation method.
The annals of mathematical statistics, pages 400–407, 1951.

[Robert and Casella, 2013]
C. Robert and G. Casella.
Monte Carlo statistical methods.
Springer Science & Business Media, 2013.

[Srikant and Ying, 2019]
R. Srikant and L. Ying.
Finite-time error bounds for linear stochastic approximation and TD
In Conference on Learning Theory, pages 2803–2830. PMLR,

[Sun et al., 2018]
T. Sun, Y. Sun, and W. Yin.
On markov chain gradient descent.
In NeurIPS, 2018.

[Sutton and Barto, 2018]
R. S. Sutton and A. G. Barto.
Reinforcement learning: An introduction.
MIT press, 2018.

[Tadić et al., 2017]
V. B. Tadić, A. Doucet, et al.
Asymptotic bias of stochastic gradient search.
Annals of Applied Probability, 270 (6):0
3255–3304, 2017.

[Wu et al., 2020]
Y. Wu, W. Zhang, P. Xu, and Q. Gu.
A finite time analysis of two time-scale actor critic methods.
In NeurIPS, 2020.

[Xu et al., 2020]
T. Xu, Z. Wang, and Y. Liang.
Non-asymptotic convergence analysis of two time-scale (natural)
actor-critic algorithms.
arXiv preprint arXiv:2005.03557, 2020.

[Zhang et al., 2020]
S. Zhang, B. Liu, H. Yao, and S. Whiteson.
Provably convergent two-timescale off-policy actor-critic with
function approximation.
In International Conference on Machine Learning, pages
11204–11213. PMLR, 2020.

§ SUPPLEMENTARY INFORMATION FOR SECTION <REF>

§.§ Example on Gaussian Estimation (A Case where $\pi_{\theta}(\cdot) \neq {\cal D}(\theta)$)

Consider the following instance of (<ref>) with:

min_ θ∈  _ z ∼D(θ)  [ ( z - θ)^2 / 2 ]   where   D(θ) ≡N( z̅ + ϵθ; σ^2 ).
Following (<ref>), the state-dependent SA algorithm reads

θ_k+1 = θ_k - γ_k+1 ℓ(θ_k; z_k+1) = θ_k - γ_k+1 ( θ_k - z_k+1 ).
where the sequence $\{ z_k \}_{k \geq 1}$ is generated by an autoregressive (AR) model, with $\rho \in (0,1]$,

z_k+1 = (1 - ρ) z_k + ρz̃_k+1   where   z̃_k+1 ∼D( θ_k ) = N( z̅ + ϵθ_k ; σ^2 ),
such that the draw of $\tilde{z}_{k+1}$ are independent.
We show that the algorithm (<ref>), (<ref>) can be analyzed as a state-dependent SA (<ref>), (<ref>) considered in our framework. Precisely, we show that the controlled MC in (<ref>) admits a stationary distribution $\pi_{\theta_k}(\cdot)$ such that $\EE_{z \sim \pi_{\theta_k}(\cdot)} [ \grd \ell(\theta; z) ] = \EE_{z' \sim {\cal D}(\theta_k)} [ \grd \ell( \theta ; z' )]$.

Observe that (<ref>) defines a controlled MC with a transition kernel denoted by $\MK_{\theta_k}: \RR \times \RR \rightarrow \RR_+$ at the $k$th iteration, as in (<ref>). For every $\theta \in \RR$, $z \in \RR$, the kernel $\MK_{\theta}$ has a unique stationary distribution given by

lim_ n →∞ _θ^n (z, ·) = π_θ(·) ≡N ( z̅ + ϵθ  ;   ρ/2-ρ σ^2 ).
Notice that the above is different from the distribution ${\cal D}(\theta)$ desired in (<ref>) unless $\rho = 1$. In the latter case, the AR model (<ref>) reduces to drawing i.i.d. samples from ${\cal D}(\theta)$. For general $\rho < 1$, it still satisfies the asymptotically unbiasedness of the stochastic gradient estimate. In particular,
_ z ∼π_θ(·)  [ ℓ( θ; z ) ] = _ z ∼π_θ(·)  [ θ- z ] = _ z ∼D(θ)  [ ℓ( θ; z ) ] .
The key observation is that for this particular performative prediction problem (<ref>), the gradient of the loss function is linear in the sample $z$. As such, with (<ref>) yielding a stationary distribution that has the same mean as ${\cal D}(\theta)$, the asymptotic unbiasedness property is unaffected.
In fact, the stationary distribution in (<ref>) has a reduced variance compared to ${\cal D}(\theta)$. Therefore, we expect the estimation error of $\theta_{PS}$ to be more stable using (<ref>), (<ref>) than [Mendler-Dünner et al., 2020].

§.§ Details on Example <ref> for Adapted Best Response

We continue the discussions in the paper with the procedure (<ref>).
When $\theta \in \RR^d$ is fixed, the procedure in (<ref>) is modelled as a Markov Chain (MC) with unique stationary distribution that corresponds to the best response distribution ${\cal D}(\theta)$ described in (<ref>).

To this end, we model the state of the MC by the tuple $\hat{z} \equiv ( d_1, ..., d_m, z )$. Consider the state space given by ${\sf Z}^{m+1}$ and
denote $\MK_{\theta} : {\sf Z}^{m+1} \times {\cal Z}^{m+1} \rightarrow \RR_+$ as the Markov transition kernel.
We remark that there is a slight abuse of notation here as the stochastic gradient $\grd \ell( \theta; \hat{z} )$ used by the learner depends only on the last term, $z$, in the agents' state variable $\hat{z}$. We have decided to use the current notation in the main paper to avoid introducing complicated notation for the implementation focused readers. Nevertheless, the SC example fits our proposed model.

Turning back on the MC. Observe when the current state is $\hat{z}$, under the action of kernel $\MK_{\theta}$, by following the description in (<ref>), we obtain the next state $\hat{z}' = ( d_1', ..., d_m', z' )$ as

d_j' = d_j + α1_ j ∈I U( d_j; d̅_j, θ), j=1,...,m,   z' = d_i',
with probability

\[
\frac{1}{ { m \choose pm } } \times \frac{1}{m} ,
\]

for any ${\cal I} \subseteq \{1,...,m\}$, $|{\cal I}| = pm$ and $i \in \{1,...,m \}$.

At each transition, the data points $\{ d_1,..., d_m \}$ are updated by the first equation in (<ref>). The latter can be treated as a special case of the random block coordinate gradient descent (RBCD) algorithm to the separable problem:

max_ d_i , i =1,...,m ∑_i=1^m U( d_i ; d̅_i, θ).
Note that the optimal solution to the above, $\{ d_1^\star, ..., d_m^\star \}$, is a set of data points that forms the empirical distribution ${\cal D}(\theta)$.
Furthermore, it is known that the RBCD algorithm converges linearly with high probability and almost surely to the optimal solution for strongly concave maximization; see [Richtárik and Takáč, 2014, Patrascu and Necoara, 2015].

With the above observations, the MC induced by $\MK_\theta$ has a stationary distribution $\pi_\theta (\cdot)$ where for any measurable function $f : {\sf Z}^{m+1} \rightarrow \RR^n$, it holds
lim_k →∞ _θ^k f( d_1,..., d_m, z ) = 1/m ∑_i=1^m f( d_1^⋆, ..., d_m^⋆, d_i^⋆) = _z' ∼D(θ) [ f( d_1^⋆, ..., d_m^⋆, z') ],
for any initial state $\hat{z} = ( d_1, ..., d_m, z )$. The above identity can be derived from the fact that the RBCD algorithm converges almost surely to the optimal solution to (<ref>) and the random variable $z^k$ is uniformly drawn from $\{ d_1^k, ..., d_m^k \}$.

Furthermore, for any $L$-Lipschitz continuous $f$, it holds

_θ^k f( d_1,..., d_m, z ) - _z' ∼D(θ) [ f( d_1^⋆, ..., d_m^⋆, z') ]

(a)≤ 1/m ∑_i=1^m  [ f( d_1^k, ..., d_m^k, d_i^k )] - f( d_1^⋆, ..., d_m^⋆, d_i^⋆)

(b)≤ L (1 + 1/m ) ∑_i=1^m [  d_i^k - d_i^⋆ ] ≤C ρ^k,

where (a) uses $\MK_\theta^k f( d_1,..., d_m, z ) = \sum_{i=1}^m \EE[ f( d_1^k , ..., d_m^k , d_i^k )] / m$ and the expectation is taken with respect to the random subset selection of ${\cal I}_k$ in (<ref>). In the expression that follows (b), the constants $\overline{\rm C}$, $\rho \in [0,1)$ depend on the initial value $\hat{z}$ and the strong concavity property of $U(\cdot)$. The above property is important for establishing the existence of the solution to Poisson equation in <Ref>.

§ CONVERGENCE ANALYSIS WITH NON-CONVEX LOSS FUNCTION

We first verify the inequality (<ref>) by observing the following expression for the gradient of performative loss:

V(θ) = ∫_Z ℓ( θ; z) p_ D(θ) (z) d z = _ z ∼D(θ)  [ ℓ(θ; z) ] + _z ∼D(θ) [ ℓ(θ;z) _θlog( p_ D(θ) ( z) ) ],
where we have denoted $p_{ {\cal D}(\theta) }(z)$ as the probability distribution function for ${\cal D}(\theta)$. The above identity is derived using chain rule and the property $\grd_\theta \log p_{ {\cal D}(\theta) } (z) = \frac{ \grd_\theta p_{ {\cal D}(\theta) } (z) }{ p_{ {\cal D}(\theta) } (z) }$ similar to the policy gradient theorem; see <cit.>.

Observe that
V( θ)  h (θ)  =  h(θ) ^2 +  _z ∼D(θ) [ ℓ(θ;z) _θlog( p_ D(θ) ( z) ) ]  h(θ)
We note

|  _z ∼D(θ) [ ℓ(θ;z) _θlog( p_ D(θ) ( z) ) ]  h(θ) | ≤1/2  h(θ)^2 + 1/2  _z ∼D(θ) [ ℓ(θ;z) _θlog( p_ D(θ) ( z) ) ] ^2

≤1/2  h(θ)^2 + 1/2 _z ∼D(θ) [ |ℓ(θ;z)|^2  _θlog( p_ D(θ) ( z) ) ^2 ] ≤1/2  h(θ) ^2 + c_0

where we have used the Jensen's inequality and set

c_0 := sup_ θ∈^d  1/2 _z ∼D(θ) [ |ℓ(θ;z)|^2  _θlog( p_ D(θ) ( z) ) ^2 ].
The above can be shown to be bounded when the loss function is bounded (e.g., a sigmoid loss), and the state-dependent distribution has bounded gradient w.r.t. $\theta$ (e.g., a soft-max distribution).
Together, we obtain the desired inequality:

V( θ)  h (θ)  ≥1/2  h(θ) ^2 - c_0, ∀ θ∈^d.
Notice that (<ref>) also implies

V( θ)      ≤ h( θ)  +  _z ∼D(θ) [ ℓ(θ;z) _θlog( p_ D(θ) ( z) ) ]

≤ h( θ)  + _z ∼D(θ) [ |ℓ(θ;z)|  _θlog( p_ D(θ) ( z) )  ] ≤ h( θ)  + √(2 c_0).

Proof of Corollary <ref> Notice that (<ref>), (<ref>) imply A1, A2 of [Karimi et al., 2019], respectively. Moreover, the stated assumptions in the corollary imply A3, A5-A7 of [Karimi et al., 2019]. Applying Theorem 2 from [Karimi et al., 2019] shows that
[  V( θ_ K  ) ^2 ] ≲[  h( θ_ K  ) ^2 ] + c_0 ≲1 + ∑_k=1^K γ_k^2 / ∑_k=1^K γ_k  + c_0,
where we have omitted the constants from [Karimi et al., 2019]. Note that ${\sf K} \in \{1,...,K\}$ is a discrete r.v. selected independently with the probability $\PP ( {\sf K} = k ) = \gamma_k / \sum_{j=1}^K \gamma_j$.
Setting the step sizes as $\gamma_k = {\cal O}(1/\sqrt{k})$ shows the desired bound in the corollary.

§ MISSING PROOFS IN SECTION <REF>

Below, we present the detailed proof for the lemmas presented in <ref>.

For <Ref>, if $b=0$, then the smoothness parameter of the function $\MK_{\theta} \widehat{\grd \ell} ( \theta ; z)$ is invariant with the optimality gap $\| \theta - \theta_{PS} \|$.
If $b = 1$, then the smoothness of the composite function $\MK_{\theta} \widehat{\grd \ell} ( \theta ; z)$ grows with the optimality gap $\| \theta - \theta_{PS} \|$, which can be justified provided that the Markov kernel is smooth w.r.t. $\theta$, e.g., $\| \MK_{\theta} - \MK_{\theta'} \|_{TV} \leq {\rm C}_P \| \theta - \theta' \|$ for some constant ${\rm C}_P$. However, as we shall reveal later in the analysis, the convergence of the SA scheme with $b=1$ would require the following assumption on the fourth moment of $\| \theta_k - \theta_{PS}\|$, i.e.,

For any $t \geq 1$, there exists a constant ${\rm C}_m$ such that
√( [  θ_t - θ_PS ^4 ] ) ≤C_m [  θ_t - θ_PS ^2 ].

The above depends on the tail of the distribution of $\| \theta_t - \theta_{PS} \|$ which may be difficult to verify. Our following analysis will be derived on two cases: (a) <Ref> with $b = 0$, or (b) <Ref> with $b = 1$ and <Ref>.

§.§ Proof of Lemma <ref>

We begin our analysis by observing that as $\grd f(\theta_{PS}; \theta_{PS} ) = 0$, we have:

\begin{align*}
& \left\|\theta_{k+1}-\theta_{PS}\right\|^{2} = \left\|\theta_{k}-\gamma_{k+1} \nabla \ell( \theta_{k}; z_{k+1} )-\theta_{PS}\right\|^{2} \\
% & = \left\|\theta_{k}-\theta_{PS}+\nabla f(\theta_{PS}; \theta_{PS})-\gamma_{k+1} \nabla \ell ( \theta_{k}; z_{k+1} ) \right\|^{2} \\
& = \underbrace{ \left\|\theta_{k}-\theta_{PS}\right\|^{2}}_{ =: B_1 } - \, 2 \gamma_{k+1} \underbrace{ \pscal{ \theta_k -\theta_{PS} }{\nabla \ell( \theta_{k}; z_{k+1} ) - \nabla f(\theta_{PS}; \theta_{PS}) } }_{=:B_2} \\
& \qquad + \gamma_{k+1}^2 \underbrace{ \left\|  \nabla f(\theta_{PS}; \theta_{PS})- \nabla \ell ( \theta_{k}; z_{k+1} ) \right\|^{2} }_{ =: B_3 }
\end{align*}

The inner product can be lower bounded as

\begin{eqnarray}
\begin{aligned}
B_{2}&= \pscal{ \theta_k -\theta_{PS}}{ \nabla \ell ( \theta_{k}; z_{k+1} ) -\nabla f(\theta_{PS}; \theta_{PS}) }\\
% &=\EE\left<\theta_k -\theta_{PS}|\nabla \ell(z_{k+1};\theta_{k})-\nabla f(\theta_{PS}; \theta_{PS})-\nabla f(\theta_{k}; \theta_{k})+\nabla f(\theta_{k}; \theta_{k})-\nabla f(\theta_{k}; \theta_{PS})+\nabla f(\theta_{k}; \theta_{PS})\right>\\
&= \pscal{  \theta_k -\theta_{PS} } { \nabla \ell ( \theta_{k}; z_{k+1} ) -\nabla f(\theta_{k}; \theta_{k}) } + \pscal { \theta_k -\theta_{PS} }{\nabla  f(\theta_{k};\theta_{k})-\nabla f(\theta_{k}; \theta_{PS}) } \\
&\quad + \pscal{ \theta_k -\theta_{PS} }{ \nabla f(\theta_{k};\theta_{PS})-\nabla f(\theta_{PS}; \theta_{PS}) }\\
&\overset{(a)}{\geq} \pscal{ \theta_k -\theta_{PS}} {\nabla \ell( \theta_{k}; z_{k+1} )-\nabla f(\theta_{k}; \theta_{k}) } \\
& \quad - \left\|\theta_{k}-\theta_{PS} \right\| \left\|\nabla f(\theta_{k}; \theta_{k})- \nabla f(\theta_{k}; \theta_{PS}) \right\|+ \mu \left\|\theta_{k}-\theta_{PS} \right\|^{2}\\
&\overset{(b)}{\geq} \left<\theta_k -\theta_{PS}|\nabla \ell( \theta_{k}; z_{k+1} ) -\nabla f(\theta_{k}; \theta_{k})\right> +(\mu-L\varepsilon) \left\|\theta_{k}-\theta_{PS} \right\|^{2}
\end{aligned}
\end{eqnarray}

where (a) is due to the Cauchy-schwarz inequality and the $\mu$-strong convexity of $\nabla f(\cdot; \cdot)$; (b) is due to the $L$-smoothness of $f$ and the $\epsilon$-sensitivity of the distribution [c.f Assumption <ref>]; also see [Perdomo et al., 2020]. Furthermore,

\begin{eqnarray}
\begin{aligned}
B_{3}&= \left\|  \nabla \ell ( \theta_{k}; z_{k+1} )-\nabla f(\theta_{PS}; \theta_{PS})+\nabla \ell ( \theta_{PS}; z_{k+1} ) -\nabla \ell ( \theta_{PS}; z_{k+1} ) \right\|^{2} \\
&\leq 2 \left( \left\| \nabla\ell ( \theta_{PS}; z_{k+1} )-\nabla\ell ( \theta_{k}; z_{k+1} ) \right\|^2+\left\|\nabla f(\theta_{PS}; \theta_{PS})-\nabla\ell ( \theta_{PS}; z_{k+1} ) \right\|^2  \right)\\
&\leq 2L^2 \left\|\theta_{k}-\theta_{PS} \right\|^2 + 2 \sigma^2
\end{aligned}
\end{eqnarray}

where the third inequality is due to Assumptions <ref>, <ref>.
Combing the bounds for $B_1$, $B_{2}$ and $B_3$, we can get the desired inequality.

\begin{eqnarray}\label{eq:recursion}
\begin{aligned}
& \left\|\theta_{k+1}-\theta_{PS}\right\|^2 \\
&\leq \left\|\theta_{k}-\theta_{PS}\right\|^2+ 2 \gamma_{k+1}^2\cdot \left(\sigma^2+L^2 \left\|\theta_{k}-\theta_{PS}\right\|^2\right)\\
&\quad -2\gamma_{k+1}\left(\left<\theta_k -\theta_{PS}|\nabla \ell( \theta_{k}; z_{k+1} )-\nabla f(\theta_{k}; \theta_{k})\right> +(\mu-L\varepsilon) \left\|\theta_{k}-\theta_{PS} \right\|^{2}\right) \\
& = \big( 1 - 2 \gamma_{k+1} ( \mu - L \epsilon ) + 2 \gamma_{k+1}^2 L^2 \big) \| \theta_k - \theta_{PS} \|^2 \\
& \quad  + 2 \gamma_{k+1}^2 \sigma^2 - 2\gamma_{k+1} \pscal{ \theta_k - \theta_{PS} }{ \grd \ell( \theta_k; z_{k+1} ) - \grd f( \theta_k; \theta_k ) }.
\end{aligned}
\end{eqnarray}

It is noted that if we consider a case when the SA scheme (<ref>) is non-state-dependent, e.g., $z_{k+1}$ is drawn from $\mathcal{D}(\theta_{k})$  independently, then proving Lemma <ref> suffices to show our desired Theorem <ref> since the last term in equation (<ref>) is zero mean when conditioned on the previous iterates [cf. (<ref>)].

§.§ Proof of Lemma <ref>

Applying Assumption <ref> shows that the sum of inner product can be evaluated as

\begin{eqnarray*}
\begin{aligned}
& \sum_{s=1}^{k}\gamma_{s}G_{s+1:k}\EE\pscal{ \theta_{PS} - \theta_{s-1} } {\nabla \ell( \theta_{s-1}; z_s )-\nabla f(\theta_{s-1}; \theta_{s-1}) }  \\
&= \sum_{s=1}^{k}\gamma_{s}G_{s+1:k}\EE\pscal{ \theta_{PS} - \theta_{s-1} } {\widehat{\nabla\ell}( \theta_{s-1}; z_s )-\MK_{\theta_{s-1}} \widehat{\nabla\ell}( \theta_{s-1}; z_s )} \equiv \EE\left(A_1+A_2+A_3+A_4+A_5\right),
\end{aligned}
\end{eqnarray*}

where we decomposed the sum of inner product into five sub-terms $A_{1}$, $A_{2}$, $A_{3}$, $A_{4}$, $A_{5}$ such that

\begin{align*}
&A_{1}:=-\sum_{s=2}^{k}\gamma_{s}G_{s+1:k}\pscal{ \theta_{s-1}-\theta_{PS}} { \widehat{\nabla\ell}( \theta_{s-1}; z_s )-\MK_{\theta_{s-1}}\widehat{\nabla \ell}( \theta_{s-1}; z_{s-1} ) } \\
&A_{2}:=-\sum_{s=2}^{k}\gamma_{s}G_{s+1:k} \pscal{ \theta_{s-1}-\theta_{PS} }{ \MK_{\theta_{s-1}} \widehat{\nabla\ell}( \theta_{s-1}; z_{s-1} )-\MK_{\theta_{s-2}}\widehat{\nabla\ell}( \theta_{s-2}; z_{s-1} ) }\\
&A_{3}:=-\sum_{s=2}^{k}\gamma_{s}G_{s+1:k} \pscal{ \theta_{s-1}-\theta_{s-2} } { \MK_{\theta_{s-2}}\widehat{\nabla\ell}( \theta_{s-2}; z_{s-1} ) }\\
&A_{4}:=-\sum_{s=2}^{k}(\gamma_{s}G_{s+1:k}-\gamma_{s-1}G_{s:k})\pscal{ \theta_{s-2}-\theta_{PS} }{ \MK_{\theta_{s-2}}\widehat{\nabla\ell}(\theta_{s-2}; z_{s-1}) }\\
&A_{5}:=-\gamma_{1}G_{2:k}\pscal{ \theta_{0}-\theta_{PS} } {\widehat{\nabla\ell}( \theta_{0}; z_1 ) } +
\gamma_{k}  \pscal{ \theta_{k-1}-\theta_{PS} }{ \MK_{\theta_{k-1}}\widehat{\nabla\ell}( \theta_{k-1}; z_k ) } .
\end{align*}

We remark that a similar decomposition can be found in [Benveniste et al., 2012]. However, [Benveniste et al., 2012] proceeded with the analysis by assuming that $\theta_k$ stays in the compact set for all $k \geq 0$. We do not make such assumption in this work.

For $A_{1},$ we note that $\widehat{\nabla\ell}( \theta_{s-1}; z_s )-\MK_{\theta_{s-1}}\widehat{\nabla\ell}( \theta_{s-1}; z_{s-1} )$ is a martingale difference sequence and therefore we have $\EE\left[A_{1}\right]=0$ by taking the total expectation.

For $A_{2}$, as $\theta_{k+1}=\theta_{k}-\gamma_{k+1} \nabla {\ell}\left(  \theta_{k}; z_{k+1} \right) $, we get $\theta_{s-1}-\theta_{s-2}=-\gamma_{s-1} \nabla {\ell} ( \theta_{s-2}; z_{s-1})$. Applying the smoothness condition  <Ref> shows that

A_2   =-∑_s=2^kγ_sG_s+1:kθ_s-1-θ_PS _θ_s-1 ∇ℓ( θ_s-1; z_s-1 )-_θ_s-2∇ℓ( θ_s-2; z_s-1)

≤∑_s=2^k γ_s G_s+1:k θ_s-1-θ_PS θ_s-1-θ_s-2

≤∑_s=2^k  γ_s-1 γ_s G_s+1:k θ_s-1-θ_PS  ℓ( θ_s-2; z_s-1 ).

Combining with the implied bound (<ref>) from the assumptions as well as (<ref>) yield

\begin{align*}
&%\leq \varsigma \Lph \overline{L} \, \sum_{s=2}^k \gamma_{s}^2  G_{s+1:k} \|\theta_{s-1}-\theta_{PS} \|  \left( 1 + b \norm{\theta_{s-1}-\theta_{PS}} \right) \left( 1 + \| \theta_{s-2} - \theta_{PS} \| \right) \\
\leq \varsigma \Lph \overline{L} \, \sum_{s=2}^k \gamma_{s}^2  G_{s+1:k} \|\theta_{s-1}-\theta_{PS} \|  \left( 1 + \| \theta_{s-2} - \theta_{PS} \| \right)
\\
&%\leq \varsigma \Lph \overline{L} \, \sum_{s=2}^k \gamma_{s}^2 G_{s+1:k}  \left\{ \frac{1}{2} + \frac{1}{2} \| \theta_{s-2} - \theta_{PS} \|^2 + (1+b) \| \theta_{s-1} - \theta_{PS} \|^2 + b \| \theta_{s-1} - \theta_{PS} \|^2 \| \theta_{s-2} - \theta_{PS} \| \right\}
\leq \varsigma \Lph \overline{L} \, \sum_{s=2}^k \gamma_{s}^2 G_{s+1:k}  \left\{ \frac{1}{2} + \frac{1}{2} \| \theta_{s-2} - \theta_{PS} \|^2 +  \| \theta_{s-1} - \theta_{PS} \|^2
%+ \| \theta_{s-2} - \theta_{PS} \|
\right\}
\\
& %\leq  \varsigma \Lph \overline{L} \, \Big\{ \frac{1}{2} \sum_{s=2}^k \gamma_s^2 G_{s+1:k} + \frac{1}{2} \sum_{s=2}^k \gamma_s^2 G_{s+1:k} \| \theta_{s-2} - \theta_{PS} \|^2 +  (1+b)\sum_{s=2}^k \gamma_s^2 G_{s+1:k} \| \theta_{s-1} - \theta_{PS} \|^2 \Big\} \\
\leq  \varsigma \Lph \overline{L} \, \Big\{ \frac{1}{2} \sum_{s=2}^k \gamma_s^2 G_{s+1:k} + \frac{1}{2} \sum_{s=2}^k \gamma_s^2 G_{s+1:k} \| \theta_{s-2} - \theta_{PS} \|^2 +  \sum_{s=2}^k \gamma_s^2 G_{s+1:k} \| \theta_{s-1} - \theta_{PS} \|^2 \Big\} ,
%& \qquad + b \, \varsigma \Lph \overline{L} \, \sum_{s=2}^k \gamma_s^2 G_{s+1:k} \| \theta_{s-1} - \theta_{PS} \|^2 \| \theta_{s-2} - \theta_{PS} \|
% \\
% &\leq  \sum_{s=2}^k \Lph (\sigma+1)\gamma_{s}\gamma_{s-1} G_{s+1:k}  \left( 1+2\norm{\theta_{s-1}-\theta_{PS}}^2\right)\cdot \left( 1 + \norm{\grd f(\theta_{s-2};\theta_{PS})-\grd f(\theta_{PS};\theta_{PS})}\right)
% \\
% &\leq  \sum_{s=2}^k \Lph (\sigma+1)\gamma_{s}\gamma_{s-1} G_{s+1:k}  \left( 1+2\norm{\theta_{s-1}-\theta_{PS}}^2\right)\cdot \left( 1 + L\norm{\theta_{s-2}-\theta_{PS}}\right)
% \\
% &\leq  \sum_{s=2}^k \Lph (\sigma+1)\gamma_{s}\gamma_{s-1} G_{s+1:k}  \left( 1+L\norm{\theta_{s-2}-\theta_{PS}} +2\norm{\theta_{s-1}-\theta_{PS}}^{2}+2L \norm{\theta_{s-1}-\theta_{PS}}^{2}\norm{\theta_{s-2}-\theta_{PS}}\right)
% \\
% &\leq  \sum_{s=2}^k \Lph (\sigma+1)\gamma_{s}\gamma_{s-1} G_{s+1:k}  \left( 1+L+L\norm{\theta_{s-2}-\theta_{PS}}^{2} +2\norm{\theta_{s-1}-\theta_{PS}}^{2}+2L \norm{\theta_{s-1}-\theta_{PS}}^{2}\norm{\theta_{s-2}-\theta_{PS}}\right)
% \\
% &\leq  \Lph (\sigma+1)\left( 1+L\right)\sum_{s=2}^k \gamma_{s}\gamma_{s-1} G_{s+1:k}
% + \Lph L(\sigma+1)\sum_{s=2}^k \gamma_{s}\gamma_{s-1} G_{s+1:k}\Delta_{s-2}\\
% &\quad +2\Lph (\sigma+1)\sum_{s=2}^k \gamma_{s}\gamma_{s-1} G_{s+1:k}\Delta_{s-1}
% +2L\Lph (\sigma+1)\sum_{s=2}^k \gamma_{s}\gamma_{s-1} G_{s+1:k}\Delta_{s-1}\sqrt{\Delta_{s-2}}
\end{align*}

where the second inequality applies $a(1 + c) \leq \frac{1}{2} + \frac{1}{2} c^2 +  a^2 $ for any $a,c\in \RR$.

For $A_{3}$, again using (<ref>), we observe that

A_3   =-∑_s=2^kγ_sG_s+1:k θ_s-1-θ_s-2  _θ_s-2∇ℓ( θ_s-2; z_s-1 )

≤∑_s=2^kγ_s G_s+1:kθ_s-1-θ_s-2·_θ_s-2∇ℓ( θ_s-2; z_s-1 )

≤∑_s=2^kγ_sγ_s-1 G_s+1:k ℓ( θ_s-2; z_s-1 )·( 1 + θ_s-2 - θ_PS  )

≤ςL ∑_s=2^kγ_s^2 G_s+1:k (1 +  θ_s-2 - θ_PS  )^2

≤2 ςL ∑_s=2^kγ_s^2 G_s+1:k { 1 +  θ_s-2 - θ_PS ^2 } .

For $A_{4}$, we notice that

A_4   =-∑_s=2^k(γ_sG_s+1:k-γ_s-1G_s:k)  θ_s-2-θ_PS  _θ_s-2∇ℓ( θ_s-2; z_s-1)

≤∑_s=2^k | γ_sG_s+1:k-γ_s-1G_s:k | θ_s-2-θ_PS ·_θ_s-2ℓ(θ_s-2; z_s-1 ) .

It can be shown that $| \gamma_{s}G_{s+1:k}-\gamma_{s-1}G_{s:k} | \leq ( 1 + \tilde{\mu} ) \varsigma \gamma_s^2 G_{s+1:k} $, therefore

A_4    ≤( 1 + μ̃ ) ς∑_s=2^k γ_s^2 G_s+1:k θ_s-2-θ_PS  ( 1 +  θ_s-2 - θ_PS  )

≤( 1 + μ̃ ) ς∑_s=2^k γ_s^2 G_s+1:k { 1/2 + 3/2  θ_s-2 - θ_PS ^2 }

≤( 1 + μ̃ ) ς{ 1/2 ∑_s=2^k γ_s^2 G_s+1:k + 3/2 ∑_s=2^k γ_s^2 G_s+1:k  θ_s-2 - θ_PS ^2 }.

Finally, for $A_{5}$, we have

\begin{align*}
A_{5}&=-\gamma_{1}G_{2:k}\pscal{ \theta_{0}-\theta_{PS}} { \widehat{\nabla\ell}(\theta_{0}; z_{1} ) } +
\gamma_{k} \pscal{ \theta_{k-1}-\theta_{PS}} {\MK_{\theta_{k-1}}\widehat{\nabla\ell} ( \theta_{k-1}; z_{k} ) }\\
&\leq \gamma_{1} G_{2:k} \left\|\theta_{0}-\theta_{PS} \right\| \left\|\widehat{\grd\ell}( \theta_{0}; z_{1} )\right\|+
\gamma_{k} \left\|\theta_{k-1}-\theta_{PS}\right\|\left\|\MK_{\theta_{k-1}}\widehat{\grd\ell}( \theta_{k-1}; z_{k}) \right\|  \\
& \leq \gamma_1 \widehat{L} \, G_{2:k} \| \theta_0 - \theta_{PS} \| \big( 1 + \| \theta_0 - \theta_{PS} \| \big) + \gamma_k \LZ \, \| \theta_{k-1} - \theta_{PS} \| \big( 1 + \| \theta_{k-1} - \theta_{PS} \| \big) \\
& \leq \frac{\gamma_1 \widehat{L} \, G_{2:k}}{2} + \frac{ \gamma_k \LZ }{2} + \frac{3 \gamma_1 \widehat{L}}{2} \, G_{2:k} \| \theta_0 - \theta_{PS} \|^2 + \frac{3 \gamma_k \LZ}{2} \, \| \theta_{k-1} - \theta_{PS} \|^2
% \\
% &\leq \gamma_{1} G_{2:k} \norm{\theta_{0}-\theta_{PS}}\cdot \hat{\sigma}\left(1+\norm{\grd f(\theta_{0}; \theta_{PS}}\right) + 	\gamma_{k}\left\|\theta_{k-1}-\theta_{PS}\right\| \cdot {\LZ}\left( 1+ \norm{\grd f(\theta_{k-1}; \theta_{PS})} \right)\\
% & {\color{red}\leq \gamma_{1} G_{2:k} \norm{\theta_{0}-\theta_{PS}}\cdot \hat{\sigma}\left(1+\norm{\grd f(\theta_{0}; \theta_{PS})-\grd f(\theta_{PS}; \theta_{PS})}\right)} \\
% &\quad +\gamma_{k}\left\|\theta_{k-1}-\theta_{PS}\right\| \cdot {\LZ}\left( 1+ \norm{\grd f(\theta_{k-1}; \theta_{PS})-\grd f(\theta_{PS}; \theta_{PS})} \right)\\
% & \leq \gamma_{1} G_{2:k} \norm{\theta_{0}-\theta_{PS}}\cdot \hat{\sigma}\left(1+L\norm{\theta_{0}-\theta_{PS}}\right)
%  +\gamma_{k}\left\|\theta_{k-1}-\theta_{PS}\right\| \cdot {\LZ}\left( 1+ L\norm{\theta_{k-1}- \theta_{PS}} \right)\\
% & \leq \gamma_{1} G_{2:k} \cdot \hat{\sigma}\left(\norm{\theta_{0}-\theta_{PS}}+L\norm{\theta_{0}-\theta_{PS}}^{2}\right)
% +\gamma_{k} {\LZ}\left( \left\|\theta_{k-1}-\theta_{PS}\right\|+ L\norm{\theta_{k-1}- \theta_{PS}}^{2} \right)\\
% & \leq \gamma_{1} G_{2:k} \cdot \hat{\sigma}\left(1+\norm{\theta_{0}-\theta_{PS}}^{2}+L\norm{\theta-\theta_{PS}}^{2}\right)  +\gamma_{k} {\LZ}\left( 1+\left\|\theta_{k-1}-\theta_{PS}\right\|^{2}+ L\norm{\theta_{k-1}- \theta_{PS}}^{2} \right)\\
% & \leq \gamma_{1}\hat{\sigma}G_{2:k}+\gamma_{k}\LZ+ \gamma_{1} \hat{\sigma}(1+L)G_{2:k}\norm{\theta_{0}-\theta_{PS}}^{2}+\gamma_{k} {\LZ}(1+L)\norm{\theta_{k-1}-\theta_{PS}}^{2}\\
% & \leq \gamma_{1}\hat{\sigma}G_{2:k}+\gamma_{k}\LZ+  \gamma_{1} \hat{\sigma}(1+L)G_{2:k}\Delta_{0}+\gamma_{k}{\LZ}(1+L)\Delta_{k-1} \\
% & \leq \gamma_{1}\hat{\sigma}G_{2:k}+\gamma_{k}\LZ+  \gamma_{1} \hat{\sigma}(1+L)G_{2:k}\Delta_{0}+\gamma_{k} {\LZ}(1+L)\Delta_{k-1}
\end{align*}

Summing up $A_{1}$ to $A_{5}$ and taking the full expectation yield:

\begin{align*}
& 2 \big| \EE\left[A_1+A_2+A_3+A_4+A_5\right] \big| \\
%& \leq \varsigma \Lph \overline{L} \, \Big\{ \sum_{s=2}^k \gamma_s^2 G_{s+1:k} + \sum_{s=2}^k \gamma_s^2 G_{s+1:k} \Delta_{s-2} + 2(1+b) \sum_{s=2}^k \gamma_s^2 G_{s+1:k} \Delta_{s-1} + 2 b {\rm C}_m \, \sum_{s=2}^k \gamma_s^2 G_{s+1:k} \Delta_{s-1} \sqrt{\Delta_{s-2}}  \Big\} \\
%& \quad + 4 \varsigma \overline{L} \LZ \sum_{s=2}^{k}\gamma_{s}^2 G_{s+1:k} \big\{ 1 + \Delta_{s-2} \big\} + \{ 1 + \mu - L \epsilon \} \varsigma \LZ \Big\{ \sum_{s=2}^{k} \gamma_{s}^{2} G_{s+1:k} + 3 \sum_{s=2}^{k} \gamma_{s}^{2} G_{s+1:k} \Delta_{s-2} \Big\} \\
%& \quad + \gamma_1 \widehat{L} \, G_{2:k} + \gamma_k \LZ  + 3 \gamma_1 \widehat{L} \, G_{2:k} \Delta_0 + 3 \gamma_k \LZ \, \Delta_{k-1}
& \leq \varsigma \Lph \overline{L} \, \Big\{ \sum_{s=2}^k \gamma_s^2 G_{s+1:k} + \sum_{s=2}^k \gamma_s^2 G_{s+1:k} \Delta_{s-2} + 2\sum_{s=2}^k \gamma_s^2 G_{s+1:k} \Delta_{s-1}  \Big\} \\
& \quad + 4 \varsigma \overline{L} \LZ \sum_{s=2}^{k}\gamma_{s}^2 G_{s+1:k} \big\{ 1 + \Delta_{s-2} \big\} + ( 1 + \tilde{\mu} ) \varsigma \LZ \Big\{ \sum_{s=2}^{k} \gamma_{s}^{2} G_{s+1:k} + 3 \sum_{s=2}^{k} \gamma_{s}^{2} G_{s+1:k} \Delta_{s-2} \Big\} \\
&\quad + \gamma_1 \widehat{L} \, G_{2:k} + \gamma_k \LZ  + 3 \gamma_1 \widehat{L} \, G_{2:k} \Delta_0 + 3 \gamma_k \LZ \, \Delta_{k-1} .
\end{align*}

Recall the following constants:

:= ςL + 4 ςL + ( 1 + μ̃ ) ς,    := 2 ςL,    := ςL + 4 ςL + 3 ( 1 + μ̃ )ς.

We obtain the desirable bound for the lemma:

\begin{align*}
& 2 \big| \EE\left[A_1+A_2+A_3+A_4+A_5\right] \big| \\
%& \leq \varsigma \Lph \overline{L} \, \Big\{ \sum_{s=2}^k \gamma_s^2 G_{s+1:k} + \sum_{s=2}^k \gamma_s^2 G_{s+1:k} \Delta_{s-2} + 2(1+b) \sum_{s=2}^k \gamma_s^2 G_{s+1:k} \Delta_{s-1} + 2 b {\rm C}_m \, \sum_{s=2}^k \gamma_s^2 G_{s+1:k} \Delta_{s-1} \sqrt{\Delta_{s-2}}  \Big\} \\
%& \quad + 4 \varsigma \overline{L} \LZ \sum_{s=2}^{k}\gamma_{s}^2 G_{s+1:k} \big\{ 1 + \Delta_{s-2} \big\} + \{ 1 + \mu - L \epsilon \} \varsigma \LZ \Big\{ \sum_{s=2}^{k} \gamma_{s}^{2} G_{s+1:k} + 3 \sum_{s=2}^{k} \gamma_{s}^{2} G_{s+1:k} \Delta_{s-2} \Big\} \\
%& \quad + \gamma_1 \widehat{L} \, G_{2:k} + \gamma_k \LZ  + 3 \gamma_1 \widehat{L} \, G_{2:k} \Delta_0 + 3 \gamma_k \LZ \, \Delta_{k-1}
& \leq \sum_{s=2}^k \gamma_s^2 G_{s+1:k} \big( \Cone + \Ctwo \Delta_{s-1} + \Cthree \Delta_{s-2} \big) + \LZ \gamma_k \big\{ 1 + 3 \Delta_{k-1} \big\} + \gamma_1 G_{2:k} \big( \LZ (1 + 3 \Delta_0 ) + \gamma_1 \Cone \big).
\end{align*}

This concludes the proof.

§.§ Proof of Lemma <ref>

Consider the inequality in (<ref>). We consider a non-negative upper bound sequence $\{ {\rm U}_k \}_{k \geq 0}$ defined by the recursion:

U_k     = G_1:k U_0 + (  2/ μ̃  (2 σ^2 + ) + ) γ_k + ∑_s=1^k-1 γ_s+1^2 G_s+2:k ( U_s + U_s-1 )

+ γ_1 G_2:k { (1 + 3 U_0 ) + γ_1 ( 2 σ^2 + ) } + 3 γ_k U_k-1 ,

for any $k \geq 1$, and we have defined ${\rm U}_0 = \Delta_0$.

Notice that by construction, we have $\Delta_k \leq {\rm U}_k$ for any $k \geq 0$.

Using the convention that ${\rm U}_{-1} = 0$, we observe that for any $k \geq 1$,

U_k     = ( 1 - γ_k μ̃ ) U_k-1 + (  2/ μ̃  (2 σ^2 + ) + ) ( γ_k - (1 - γ_k μ̃ ) γ_k-1 )

+ γ_k^2 ( U_k-1 + U_k-2 ) + 3 ( γ_k U_k-1 - (1 - γ_k μ̃) γ_k-1 U_k-2 )

≤( 1 - γ_k μ̃ + γ_k^2 ) U_k-1 + (  2/ μ̃  (2 σ^2 + ) + ) μ̃ ςγ_k^2 + ( + 3 μ̃ )γ_k γ_k-1 U_k-2

+ 3 ( γ_k U_k-1 - γ_k-1 U_k-2 )

≤( 1 - γ_k μ̃ / 2 ) U_k-1 + (  2/ μ̃  (2 σ^2 + ) + ) μ̃ ςγ_k^2 + ( + 3 μ̃ )γ_k γ_k-1 U_k-2

+ 3 ( γ_k U_k-1 - γ_k-1 U_k-2 ),

where the last inequality is due to $\gamma_k \leq \tilde{\mu}/2 \Ctwo$.

We prove part (i) of the lemma.
From (<ref>), we consider an upper bound sequence $\{ \oU_k \}_{k \geq -1}$ defined by the recursion:

_k     = ( 1 - γ_k μ̃ / 2 ) _k-1 + (  2/ μ̃  (2 σ^2 + ) + ) μ̃ ςγ_k^2 + ( + 3 μ̃ )γ_k γ_k-1 _k-2

+ 3 ( γ_k _k-1 - γ_k-1 _k-2 ),   ∀k ≥1.

We have also defined $\oU_0 = {\rm U}_0$, $\oU_{-1} = 0$.

For any $t \geq 1$, summing up the equation (<ref>) from $k=1$ to $k=t$ yields

∑_k=1^t _k     = ∑_k=1^t { ( 1 - γ_k μ̃ / 2 ) _k-1 + (  2/ μ̃  (2 σ^2 + ) + ) μ̃ ςγ_k^2 + ( + 3 μ̃ )γ_k γ_k-1 _k-2 }

+ 3 ∑_k=1^t ( γ_k _k-1 - γ_k-1 _k-2 ),

Rearranging terms leads to
_t = _0 +  ∑_k=1^t {  ( + 3 μ̃ ) γ_k γ_k-1 _k-2 + (  2/ μ̃  (2 σ^2 + ) + ) μ̃ ςγ_k^2 -  μ̃/2 γ_k _k-1 } + 3 γ_t _t-1
Using the step size conditions $\gamma_k \leq \gamma_{k-1}$, $\gamma_k \leq \big( \Cthree + 3 \LZ \tilde{\mu} \big)^{-1} \min\{  \tilde{\mu} / 2 , 3 \LZ \}$, $\gamma_k \leq (6 \LZ)^{-1}$,

_t     ≤_0 + 3 γ_t _t-1

+ ∑_k=1^t { [ ( + 3 μ̃ ) γ_k^2 - μ̃/2 γ_k ] _k-1 + (  2/ μ̃  (2 σ^2 + ) + ) μ̃ ςγ_k^2 }

≤3 γ_t _t-1 + _0 + (  2/ μ̃  (2 σ^2 + ) + ) μ̃  ς∑_k=1^t γ_k^2

≤1/2 _t-1 + { _0 + (  2/ μ̃  (2 σ^2 + ) + ) μ̃ ς∑_k=1^t γ_k^2 },

where we obtain the first inequality after  shifting the summation's index and it is noted that $\oU_{-1}=0.$ Rearranging terms and solving the recursion lead to

_t     ≤( 1/2 )^t _0 + ∑_s=1^t ( 1/2 )^t-s { _0 + (  2/ μ̃  (2 σ^2 + ) + ) μ̃ ς∑_k=1^s γ_k^2 }

≤3 _0 + 2 μ̃ ς(  2/ μ̃  (2 σ^2 + ) + ) ∑_ℓ=1^t γ_ℓ^2 ( 1/2 )^t-ℓ

≤3 _0 +  μ̃ ς/ 9 ^2  (  2/ μ̃  (2 σ^2 + ) + )

Recall that $\overline{\Delta} \eqdef 3 \oU_0 + \frac{ \varsigma}{ 9 \LZ^2 } \Big( 2 (2 \sigma^2 + \Cone) +  \tilde{\mu} \LZ \Big)$, the above shows $\Delta_t \leq {\rm U}_t \leq \oU_t \leq \overline{\Delta}$ for any $t \geq 1$, thus establishing part (i).

We now proceed to proving part (ii) of the lemma. We define $\oG_{m:n} = \prod_{\ell=m}^n (1 - \gamma_\ell \tilde{\mu}/2)$ and observe from (<ref>) that

U_k     ≤_1:k U_0 + ∑_s=1^k _s+1:k { (  2/ μ̃  (2 σ^2 + ) + ) μ̃ ςγ_s^2 + ( + 3 μ̃ )γ_s γ_s-1 U_s-2 }

+ 3 ∑_s=1^k _s+1:k {  ( γ_s U_s-1 - γ_s-1 U_s-2 ) }.

Notice that

∑_s=1^k _s+1:k ( γ_s U_s-1 - γ_s-1 U_s-2 )

= ∑_s=1^k _s+1:k ( γ_s U_s-1 + (1 - γ_s μ̃ / 2) ( γ_s-1 U_s-2 - γ_s-1 U_s-2 ) - γ_s-1 U_s-2 )

= ∑_s=1^k { ( _s+1:k γ_s U_s-1 - _s:k γ_s-1 U_s-2 ) - γ_s γ_s-1 μ̃ U_s-2 / 2 } ≤γ_k U_k-1 ≤γ_k Δ.

By Lemma <ref>, we have $\sum_{s=1}^k \oG_{s+1:k} \gamma_s^2 \leq 4 \gamma_k / \tilde{\mu}$ and the following is obtained
U_k ≤_1:k U_0 + { 4 ς/μ̃ ( 2 (2 σ^2 + ) + μ̃ ) + 3 Δ  + 4 ς/μ̃ ( + 3 μ̃ ) Δ } γ_k .

The proof is completed.

§.§ Auxiliary Lemmas

Let $a>0$ and $\left(\gamma_{k}\right)_{k \geq 1}$ be a non-increasing sequence such that $\gamma_{1}< 2 / a$. If $\gamma_{k-1} / \gamma_k \leq 1 + (a/2)\gamma_k$ for any $k \geq 1$, then for any $k \geq 2$,
∑_j=1^k γ_j^2 ∏_ℓ=j+1^k(1-γ_ℓ a) ≤2/a γ_k.

The proof is elementary. Observe that:

∑_j=1^k γ_j^2 ∏_ℓ=j+1^k(1-γ_ℓ a)      = γ_k ∑_j=1^k γ_j ∏_ℓ=j+1^k  γ_ℓ-1 / γ_ℓ (1-γ_ℓ a)

≤γ_k ∑_j=1^k γ_j ∏_ℓ=j+1^k (1 + (a/2) γ_ℓ) (1-γ_ℓ a)

≤γ_k ∑_j=1^k γ_j ∏_ℓ=j+1^k (1-γ_ℓ (a/2) )

= 2 γ_k/a ∑_j=1^k ( ∏_ℓ=j+1^k (1 - γ_ℓa/2 ) - ∏_ℓ'=j^k (1 - γ_ℓ' a/2 ) )

= 2 γ_k/a ( 1 - ∏_ℓ'=1^k (1 - γ_ℓ' a/2 ) ) ≤2 γ_k/a.

The proof is concluded.

§ DETAILS OF THE NUMERICAL EXPERIMENTS

This section provides details about the numerical experiments on the second problem of strategic classification (SC) in <ref>. Moreover, we provide additional experiment results to better illustrate the performance of the state dependent SA algorithm for this problem.

The experiments conducted in this section are based on the Credit simulator provided at <https://github.com/zykls/performative-prediction>.
Our experiments are conducted on a server with Intel Xeon Gold 6138 CPU. The Python codes are executed in a single-thread environment.

There are two roles in the SC problem – learner and agents. The learner utilizes agents' information to obtain a classifier $f_{\theta}$. Meanwhile, individual agents hope to be assigned to a favorable class. To do so, they modify their features and thereby shifting the data distribution towards the target ${\cal D}(\theta)$.
Specifically, our experiments are done on the  dataset with $m=18357$ samples as we select $d=3$ features to build the classifier. Each (original) data sample is given by $\bar{z}_i = (\bar{x}_i, \bar{y}_i)$ with the label $\bar{y}_i \in \{0,1\}$ and selected feature $\bar{x}_i \in \RR^3$. We associate each data sample to an agent. The task for the learner (bank) is to design a classifier that distinguishes whether the application of an individual (agent) who want to default a loan should be granted or not.

We simulate the adapted best response presented in Example 1 of the main paper. In this setting, the agents rely on their past experience to present data to the learner that is favorable to to agents. The latter is achieved by a gradient descent step that depends on the current learner's state ($\theta_k$), past agent's state ($z_k$) and the original data (${\cal D}_0$).
As the dynamics is coupled between the agents' and learner's update, we present the overall algorithm based on (<ref>), (<ref>) as follows:

Algorithm 2: State-dependent SA with Adapted Best Response.

Input: initial iterate $\theta_0 \in \RR^d$, agents' state $x_i^0 = \bar{x}_i$, $i \in \{1,...,m\}$ such that $\bar{x}_i$ is the $i$th original feature vector, step sizes $\{ \gamma_k \}_{k \geq 0}$, agents' response rate $\alpha > 0$, update parameter ${\sf b}$.

For $k=0,1,2, \ldots$

* A subset of agents, ${\cal I}_k$ with $|{\cal I}_k| = {\sf b}$, is selected uniformly from $\{1,...,m\}$. They adapt their feature vectors based on past experience and $\theta_k$ as:

x_i^k+1 = x_i^k + αU( x_i^k ; z̅_i, θ_k ), ∀ i ∈I_k,   x_i^k+1 = x_i^k, ∀ i ∉I_k.

* An agent $i_k \in \{1,...,m\}$ is drawn uniformly to present data. Set $z_{k+1} = ( x_{i_k}^{k+1} , y_{i_k})$.

* The learner computes the $k+1$th iterate by:

θ_k+1=θ_k-γ_k+1 ℓ(θ_k; z_k+1).
The most recent iterate $\theta_{k+1}$ is deployed and made available to the agent(s).

Steps 1 & 2 in Algorithm 2 resemble the adaptive best response update in (<ref>). We emphasize that these two steps are agnostic to the learner as the latter only sees $z_{k+1}$ at iteration $k$, similarly, the last step is not known to the agents as the latter only sees the classifier given as $\theta_{k+1}$.

Furthermore, we recall that the following two types of utility functions are considered as $U(\cdot)$:

U_q(x^';z, θ)     = θx'- x'-x^2/2ϵ,

U_lg(x^';z, θ)     = y θx' - log(1+exp(θx'))- x'-x^2/2ϵ.

In step 1, the agents' response rate $\alpha$ and parameter ${\sf b}$ control the speed of adaptation among the group of $m$ agents. These parameters will affect the mixing time of the MC which determines the bounds in <Ref>.
Overall, we observe that the agents' states and learner's iterates are evolving simultaneously, highlighting the coupled nature in the analysis of the state-dependent SA algorithm.

In cases such as $U_{\sf lg}(\cdot)$ where the ideal best response $\argmax_{x'} U(x'; z, \theta)$ must be obtained via an iterative algorithm. From an algorithmic standpoint, the stateful nature for the agent is necessary for the performative prediction algorithm to converge to $\theta_{PS}$.

Additional Experiments Next, we provide additional experiments to illustrate the performance of the state-dependent SA algorithm from a few additional perspectives. Unless otherwise specified, we adopt the same parameters set in the experiments presented in the main paper.
In particular, we set $\beta = 1000/m$ in (<ref>), $\epsilon = 0.01$ in the utility functions, and  in (<ref>), we set number of selected agents $|{\cal I}_k| = 5$, agents' response rate $\alpha = 0.5 \epsilon$.
The step size for (<ref>) is $\gamma_k = {c_0} / {(c_1+k)}$, $c_0 = {100} / {\tilde{\mu}}, c_1 = { 8 L^2 } / {\tilde{\mu}^2}$, where $L, \tilde{\mu}$ are estimated as $\sqrt{2 \beta m + \| X \|_F^2 / 2}$, $(1-\epsilon) \beta - \epsilon \| X \|_F^2 / 4m$, respectively.  By default, the SA algorithm is executed as presented in Algorithm 2 with a batch size of ${\sf batch} = 1$ and the agents perform only ${\sf BR}=1$ best response update per SA update in step 3 of Algorithm 2.

Besides, we compare the convergence rates of the algorithms from the perspective of the agents – measured by the number of BR updates performed by the agents. This is the setting used in the plot of Fig. <ref> (right) and is denoted with the $x$-axis label of `no. of agent update'. We also compare the convergence from the perspective of the learner – measured by the number of samples requested from the agents by the learner. This setting is denoted with the $x$-axis label of `no. of samples drawn by learner'.

Convergence of SA algorithm with varying number of learner's updates per iteration.

Effects of Stateful Updates at Agents Notice that the comparison has been made in Fig. <ref> (right). Here, we again plot the convergence of the SA algorithm

to illustrate the convergence rates from the learner's perspective as well. We observe that the SA algorithms with stateful update converges as $k$ increases.

We vary the `learner's iteration' parameter to observe the effects on convergence when the learner is adapting at faster rate than the agents. This is achieved by repeating steps 2 and 3 in Algorithm 2 for multiple times. Notice that this setting is similar to the lazy deploy scheme in [Mendler-Dünner et al., 2020]. From the figure, we observe that doing so improves the convergence from the agents' perspective, while the sample efficiency (from the learner's perspective) is unaffected.

Convergence of SA algorithm with varying number of minibatch size.

Effects of Minibatch Size In this experiment, we consider the setting of $U_{\sf lg}(\cdot)$ and we draw different batch size of samples ($\hat{\sf b} \in \{1,5,10\}$) per iteration. To implement this, at step 2 of Algorithm 2, the learner draws $\hat{\sf b}$ agents uniformly as $\hat{\cal I}_k$, and at step 3, we update the iterate through:

\[
\theta_{k+1} = \theta_k - \gamma_{k+1} \frac{1}{\hat{\sf b}} \sum_{ j \in \hat{\cal I}_k } \grd \ell( \theta_k ; z_{k+1, j}).
\]

In Fig. <ref>, we compare the error $\| \theta_k - \theta_{PS} \|^2$ in terms of the number of agents' best response update[Since the agents only perform one best response update per iteration, the $x$-axis here is equivalent to the iteration number $k$.] and the number of samples drawn by the learner.
We find that increasing the minibatch reduces the variance of the gradient estimate, yet it can be less sample efficient from the learner's perspective.

Convergence of SA algorithm with varying number of best responses.

Effects of Number of Adaptive Best Responses In this experiment, we consider the setting of $U_{\sf lg}(\cdot)$ and at each iteration, the agents execute multiple rounds of adapted best response (${\sf BR} \in \{1,2,4\}$) to simulate the scenario when the agents are allowed with more time to respond to the published classifier $\theta_k$. To implement this, we repeat the update in (<ref>) of step 1 in Algorithm 2 for ${\sf BR}$ times.
Notice that this is reverse of Fig. <ref> (right) where the learner performs multiple iterations per agents' best response update.

In Fig. <ref>, we compare the error $\| \theta_k - \theta_{PS} \|^2$ in terms of the number of agent update and the number of samples drawn by the learner. We observe that increasing the number of best responses improves the performance slightly. However, as a drawback, it requires more computations/updates at the agents to reach the same performance.