Strategic Distribution Shift of Interacting Agents via Coupled Gradient Flows
Authors: Lauren E Conger
Franca Hoffman
Eric Mazumdar
Lillian J Ratliff
Venue: NeurIPS 2023 poster (2023)
Keywords: distribution shift
partial differential equations

Abstract: We propose a novel framework for analyzing the dynamics of distribution shift in real-world systems that captures the feedback loop between learning algorithms and the distributions on which they are deployed. Prior work largely models feedback-induced distribution shift as adversarial or via an overly simplistic distribution-shift structure. In contrast, we propose a coupled partial differential equation model that captures fine-grained changes in the distribution over time by accounting for complex dynamics that arise due to strategic responses to algorithmic decision-making, non-local endogenous population interactions, and other exogenous sources of distribution shift. We consider two common settings in machine learning: cooperative settings with information asymmetries, and competitive settings where a learner faces strategic users. For both of these settings, when the algorithm retrains via gradient descent, we prove asymptotic convergence of the retraining procedure to a steady-state, both in finite and in infinite dimensions, obtaining explicit rates in terms of the model parameters. To do so we derive new results on the convergence of coupled PDEs that extends what is known on multi-species systems. Empirically, we show that our approach captures well-documented forms of distribution shifts like polarization and disparate impacts that simpler models cannot capture.

Strategic Distribution Shift of Interacting Agents via Coupled Gradient Flows
Lauren Conger California Institute of Technology 1conger@caltech.. edu
Franca Hoffmann California Institute of Technology franca. hoffmann@caltech edu
Eric Mazumdar
Lillian Ratliff
California Institute of Technology
University of Washington ratliffl@ux edu
mazundar@caltech. edu
Abstract
We propose a novel framework for analyzing the dynamics of distribution shift real-world systems that captures the feedback loop between learning algorithms and the distributions on which they are deployed. Prior work largely models feedback- induced distribution shift adversarial via overly simplistic distribution-shift structure. In contrast, we propose a coupled partial differential equation model that captures fine-grained changes in the distribution over time by accounting for com- plex dynamics that arise due to strategic responses to algorithmic decision-making, non-local endogenous population interactions, and other exogenous sources of dis- tribution shift. We consider two common settings in machine learning: cooperative settings with information asymmetries, and competitive settings where learner faces strategic users. For both of these settings, when the algorithm retrains via gradient descent, we prove asymptotic convergence the retraining procedure steady-state, both in finite and in infinite dimensions, obtaining explicit rates in terms of the model parameters. To do so we derive new results on the convergence of coupled PDEs that extends what is known multi-species systems. Empirically, we show that our approach captures well-documented forms of distribution shifts like polarization and disparate impacts that simpler models cannot capture.
Introduction
In many machine learning tasks, there are commonly sources of exogenous and endogenous dis- tribution shift, necessitating that the algorithm be retrained repeatedly over time. Some of these shifts occur without the influence of an algorithm; for example, individuals influence each other to become more or less similar in their attributes, or benign forms of distributional shift occur [Qui+]. Other shifts, however, are in response to algorithmic decision-making. Indeed, the very use of gpeision-making algorithm can incentivize individuals change or mis-report their data achieve desired outcomes- phenomenon known in economics as Goodhart's law. Such phenomena have been empirically observed, a well-known example being in [CCIT], where researchers observed population in Columbia strategically mis-reporting data to game poverty index score used for distributing government assistance. Works such as MMi++20; Wil+2T], which investigate the effects of distribution shift over time machine learning algorithm, point toward the need for evaluating the robustness of algorithms to distribution shifts. Many existing approaches for modeling distribution shift focus simple metrics like optimizing over moments covariates [DYIO, LHL21 BBS09] Other methods consider worst-case scenarios, as in distributionally robust optimization AZ22
LFG22 DN21 Kuh+19] However, when humans respond algorithms, these techniques may not
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

be sufficient to holistically capture the impact an algorithm has on population. For example, an algorithm that takes into account shifts in a distribution' 's mean might inadvertently drive polarization, rendering portion of the population disadvantaged.
Motivated by the need for more descriptive model, we present alternative perspective which allows us to fully capture complex dynamics that might drive distribution shifts in real-world systems. Our approach general enough to capture various sources of exogenous and endogenous distribution shift including the feedback loop between algorithms and data distributions studied in the literature on performative prediction (0;
R2y+22 Nar+22, MPZ2T], the strategic interactions studied in strategic classification [Har+ 16; Don+18), and also endogenous factors like intra populloooot dynamics and distributional shifts. Indeed, while previous works have studied these phenomena in isolation, our method allows us to capture all of them as well as their interactions. For example, in [mm+TT], the authors investigate the effects of dynamies in strategic classification problems- but the model they analyze does not capture individual interactions in the population. In [[YZ21], the authors model the interaction between a population that repeatedly responds to algorithmic decision-make makiing by shifting its mean. Additionally, [2ay+22] study settings in which the population has both exogenous and endogenous distribution shifts due to feedback, but much like the other cited work, the focus remains average performance. Each of these works fails to account for diffusion or intra- antra-population interactions that can result in important qualitative changes to the distribution.
Contributions, Our approach to this problem relies on a detailed non-local PDE model of the data distribution which captures each of these factors. One term driving the evolution of the distribution over time captures the response of the population to the deployed algorithm, another draws on models used in the PDE literature for describing non-local effects and consensus in biological systems to model intra- poraappoplationn dynamics, and the last captures a background source of distribution shift. This is coupled with an ODE, lifted to a PDE, which describes the training of a machine learning algorithm results in a coupled PDE system which we analyze to better understand the behaviors that can arise among these interactions.
In one subcase, our model exhibits a joint gradient flow structure, where both PDEs can be written as gradients flows with respect the same joint energy, but considering infinite dimensional gradients with respect to the different arguments. This mathematical structure provides powerful tools for analysis and has been an emerging area of study with relatively small body of prior work, none of which related to distribution shifts in societal systems, and general theory for multi-species gradient flows is still lacking. We give a brief overview of the models that are known to exhibit this joint gradient flow structure: in [DS20) the authors consider a two-species tumor model with coupling through Brinkman's is Law, A number of works consider coupling via convolution kernels [F113 Giu+22, JPZ22| CHSI8; DT20, Dou+23] and cross-diffusion (LY22, AB21 MKB14]] with applications in chemotaxis among other areas. In the models we consider here, the way the interaction between the two populations manifests is neither via cross-diffusion, nor via the non-local self-interaction term. A related type of coupling has recently appeared in (HPS22) HPS23), however in the setting of graphs. Recent work [Dom+21] provides particle-basee methods to approximately compute the solution to a minimax problem where the optimization space is over measures; following that work, (WC23) provides another particle- based method using mirror descent-ascent to solve similar problem. Other recent work [Lu23) proves that a mean-field gradient ascent-descent scheme with an entropy annealing schedule converges to the solution of a minimax optimization problem with timescale separation parameter that is also time-varying; in contrast, our work considers fixed timescale separation setting. [GG23] show that the mean-field description of a particle method for solving minimax problems has proveable convergence guarantees in the Wasserstein-Fiisher metric. Each of these references considers an energy functional that is linear in the distribution of each species respectively; our energy includes nonlincarities in the distributions via a self-interaction term as well as diffusion for the population. Moreover, the above works introduce gradient flow dynamic as a tool for obtaining and characterizing the corresponding steady states, whereas in our setting we seek to capture the time-varying behavior that models distributions shifts. In the other subcase, we prove exponential convergence in two competitive, timescale separated settings where the algorithm and strategic population have conflicting objectives. We show numerically that retraining in a competitive setting leads to polarization in the population, illustrating the importance of fine-grained modeling.
2

---

2 Problem Formulation
Machine learning algorithms that deployed into the real world for decision-makes makinng often become part of complex feedback loops with the data distributions and data sources with which they interact. In effort model these interactions, consider machine learning algorithm that has loss given by Lz,, .) where æ E I² are the algorithm parameters and 2 E R² are the population attributes, and the goal is to solve
argmin E L(2,z), PEA zmup
where X' the class of model parameters and p(z) is the population distribution. Individuals have an objective given by [(z, () in response to model parameterized by 2, and they seek to solve
argmin J(, #ER4
When individuals in the population and the algorithm have access to gradients, we model the optimization process as a gradient-descent-lype process. Realistically, individuals in the population will have nonlocal information and influences, as well as external perturbations, the effects of which we seek to capture addition to just minimization. To address this, we propose a partial differential equation (PDE) model for the population, that is able to capture nonlocal interactions between individuals on the level collective population. To analyse how the population evolves over time, notion of derivative in infinite dimensions is needed. A natural, and in this context physically meaningful, way of measuring the dissipation mechanism for probability distributions is the Wasserstein-2 metric (see Definition 4 The following expression appears when computing the gradient of an energy functional with respect to the Wasserstein-2 topology.
Definition 1. [First Variation] For a map G P(R") + and fixed probability distribution E ((RR), the first variation of G at the point p is denoted by 5pG(p) Rdd + R. and is defined via the relation
I 8,Gl)l())())dd lim 1 (G(p -Gp) e
for all vi E C22(TR)) such that dip O, assuming that G is regular enough for all quantities to exist.
Here, P( denotes the space of probability measures on the Borel sigma algebra. Using the first variation, express the gradient in Wasserstein-2 space, for example (Vi103) Exercise 8.8].
Lemma 1. The gradient of an energy G (2R()) + R in the Wasserstein-2 space is given by
Vw, G(p) -iv (pV8,Gloo))
Here, P2 IR denotes the set of probability measures with bounded second moments, also see
Appendix A.2 As consequence, the infinite dimensional steepest descent in Wasserstein-2 space
can be expressed as the PDE
0.p = - Ve. = div (pVá,G((p)
(I)
All the coupled gradient flows considered in this work have this Wasserstein-2 structure. In particular, when considering that individuals minimize their own loss, we can capture these dynamics via gradient flow in the Wasserstein-2 metric the level of the distribution of the population. Then for given algorithm parameters x e Rdd, the evolution for this strategic population given by
0,p = div (DV6, E (2, + (((=a)+E(()) zmup
Blo))
where E(()) is a functional including terms for internal influences and external perturbations. In real-world deployment of algorithmss, decision makers update their algorithm over time. leading to an interaction between the two processes. We also consider the algorithm dynamics over time, which we model
- E ((2,x))
2rip
3

---

In this work, we analyze the behavior of the dynamics under the following model. The algorithm suffers cost (2, x for data point under model parameters x in the strategic population, and cost J2 (2, .)) for data point in fixed, non-strategic population. 'he strategic population is denoted by p E P, and the non-strategic population by p E P. The algorithm aims to minimize
2)) [1(,,))p((=) (2(2,z)d0(a)) + || - =oll² 2mp
E
where the norm is the vector inner product ||=||² = (0,x)) and weights the cost of updating the model parameters from its initial condition.
We consider two settings: (t) aligned objectives, and (ii) competing objectives. Case (i) captures the setting in which the strategic population minimization improves the performance of the algorithm, subject to a cost for deviating from reference distribution required to manipulate features, such as loan applicant adding or closing credit cards. On the other hand, Case (ii) captures the setting in which the strategic population minimization worsens the performance of the algorithm, again incurring cost from distributional changes.
P. This cost stems from effort
2.1 Case (i): Aligned Objectives
In this setting, we consider the case where the strategic population and the algorithm have aligned objectives. This occurs in examples such as recommendation systems, where users and algorithm designers both seek develop accurate recommendations for the users. This corresponds to the population cost
/ʃ
E smppzrup
[2, -
(2, fi(e, + aKL(p)p)
where (-) denotes the Kullback-Leibler divergence. Note that the KL divergence introduces diffusion to the dynamics for p. The weight population. To account for nonlocal information and influence among members of the population, we include a kernel term E(p) }/2 pW * pdz, where (W * p)(() ([ (0==)==pp(pp 2)dp(z) ()) convolution integral and W is suitable interaction potential.
parameterizes the cost of distribution shift the
2.2 Case (ii): Competing Objectives
In settings such as online internet forums, where algorithms and users have used manipulative strategies for marketing [De106), the strategic population may be incentivized to modify or mis-report their attributes. The algorithm has competitive objective, in that it aims to maintain performance against a population whose dynamics cause the algorithm performance suffer. When the strategic population seeks an outcome contrary to the algorithm, we model strategic population cost as
12,x) = // f1 f1(2, + aK p).
E
significant factor in the dynamics for the strategic population is the timescale separation between the two "species"-i.e. the population and the algorithm. In analysis, we will consider two cases: one, where the population responds much faster than the algorithm, and two, where the algorithm responds much faster than the population. We illustrate the intermediate in simulation example.
Results
We are interested characterizing the long-tim asymptotic behavior of the population distribution, as it depends on the deaisioonmakenn action over time. The structure of the population distribution gives us insights about how the decision-man makers actions influences the entire population of users. For instance, as noted in the preceding sections, different behaviors such as bimodal distributions or large tails or variance might emerge, and such effects are not captured in simply looking at average performance. To understand this intricate interplay, one would like to characterize the behavior of both the population and the algorithm over large times. Our main contribution towards this goal is novel analytical framework well analysis of the long-time time asymptotics.

---

A key observation is that the dynamics in 2 and 3 can be re-formulated as a gradient flow; we lift a probability distribution μ by representing the evolution of /a will be governed by PDE, and combined with the PDE for the population, we obtain a system of coupled PDEs,
Dirac delta j sitting at the point z. As result,
3.p div
/V.0pl
pV
E smppatoi
(2, + B(o)))
8,1 div 1 I pmpparne T(2,,))
where 5p and 6, are first variations with respect to p and μ according to Definition The natural candidates for the asymptotic profiles of this coupled system are its steady states, which thanks to the gradient flow structure can characterized as ground states of the corresponding energy functionals. In this work, we show existence and uniqueness of minimizers (maximizers) for the functionals under suitable conditions on the dynamics. We also provide criteria for convergence and explicit convergence rates. We begin with the case where the interests of the population and algorithm are aligned, and follow with analogous results in the competitive setting. We show convergence in energy, which in turn ensures convergence in product Wasserstein metric. For convergence in energy, we use the notion of relative energy and prove that the relative energy converges to zero as time increases.
Definition 2 (Relative Energy). The relative energy of a functional G is given by G(1l7s>) = G(+) G(%s)). where G(% is the energy the steady state.
Since we consider the joint evolution of two probability distributions, define distance metric W on the product space of probability measures with bounded second moment.
Definition 3 (Joint Wasserstein Metric). The metric over P2(R2) x P2(R2) is called W and is given by
((p, p). (p. A)) = W2(p, x)3 + W2(H, p2)²
for all pairs (p. p1, (1) E (R() (RR,, and where W2 denotes the Wasserstein-2 metric (see
Definition 4 We denote by W(R") = (P2(RR) P2( W) the corresponding metric space.
3.1 Gradient Flow Structure
In the case where the objectives of the algorithm and population are aligned, we can write the dynamies as gradient flow by using the same energy functional for both species. Let Ga P(R4) P(Rd) (0, oo] be the energy functional given by
Ga( = /ʃ fi (((xx)dd(s))d(()) + ʃf fa(z, + aKL(p))) + 2 ow
+ 8 / || xall dd(x).
This expression well-defined as the relative entropy KL can be extended to the full set P(Rd) by setting (p. +oo in case p is not absolutely continuous with respect p.
In the competitive case we define Ge(p.z) ((R)) x Rd + [-o0, 20] by fe(z,z')dp(z'))
Ge(p,z) = (222))))((2) - /
-aapnz/)>aK((PP)) - // pW pW*p= * pp
8 zoll
2 |1
settings like recommender systems, the population and algorithm have aligned objectives; they seek to minimize the same cost but are subject to different dynamic constraints and influences. modeled by the regularizer and convolution terms. In the case where the objectives are aligned, the dynamics are given by
0,p = div (pV (o.p..)))l.
(4)
0,p div (nV G.
Note that 44) is a joint gradient flow, because the dynamics can be written in the form
a, = div (7 (7
5

---

where y= (p. jo) and where the gradient and divergence are taken in both variables (2,x). We discuss the structure of the dynamics (4) as well as the meaning of the different terms appearing in the energy functional G. in Appendix A.1
In other settings, such credit score reporting, the objectives of the population competitive with respect to the algorithm. Here we consider two scenarios; one, where the algorithm responds quickly relative to the population, and two, where the population responds quickly relative to the algorithm. In the case where the algorithm can immediately adjust optimally (best-respond) to the distribution, the dynamics are given by
O.p = -iiv (PV (a=b(p))
(5)
b(p) G.(p.
Next can consider the population immediately responding the algorithm, which has dynamics
d
dt x V.G.(p.z))oorrg
(6)
r(x) = argmin x)
PEP
In this time-scale separated setting, model represents dyamic maximization of G. with respect to p W space, and an instantancous minimization of G. with respect to the algorithm parameters x. Model (6) represents an instantaneous maximization of G. with respect to and a dynamic minimization of Ge with respect to the algorithm parameters The key results existence and uniqueness of ground state as well as the convergence behavior of solutions depend on convexity (concavity) of G. and Ge. The notion of convexity that will employ for energy functionals in the Wasserstein-2 geometry is (uniform) displacement convexity, which is analogous to (strong) convexity in Euclidean spaces. One think of displacement convexity for an energy functional defined on P3 as convexity along the shortest path in the Wasserstein-2 metric (linear interpolation in the Wasserstein-2 space) between any two given probability distributions. For detailed definition of (uniform) displacement convexity and concavity, see Section A.2 In fact, suitable convexity properties of the input functions f1, f2, W and ò will ensure (uniform) displacement convexity of the resulting energy functionals appearing in the gradient flow structure, see for instance [vi103 Chapter 5.2].
We make the following assumptions in both the competitive case and aligned interest cases. Here, la denotes the d x d identity matrix, Hess (f) denotes the Hessian of in all variables, while V denotes the Hessian of ff in the variable x only.
Assumption 1 (Convexity of fi and f2). The functions fI, E CR x Rd; [0, oo)) satisfy for all (2,x) Rd Rdd the following:
"Thre exists constants X1, X2 2 0 such that Hess (fi) 2 A1 Ind and Vh.f2 ł A2 Ia:
* TThree exist constants O; O such that vnJi(z,x) 2 for 1,2; Assumption 2 (Reference Distribution Shape). The reference distribution ò E P(Rd) n L (RR) satisfies log P E (22(R2)) and V log p(z) 3 - for some > 0.
Assumption (Convex Interaction Kernel). The interaction kernel W E C [0, oo)) is conver. symmetric W (-=))= (()), and for some D satisfies
(2W(=) 2 DD, [v)((=) < D(1+ Vz VzER" Rd.
We make the following observations regarding the assumptions above:
The convexity in Assumption can be relaxed and without affecting the results outlined below by following more detailed analysis analogous to the approach in [CMV03).
If fi and f2 are strongly convex, the proveable convergence rate increases, but without strict or strong convexity of fi and fz, the regularizers K LL)p)p) and f ||z zo zall dr provide the convexity guarantees necessary for convergence.
For concreteness, one can consider the following classical choices of input functions to the evolution:
Using the log-loss function for fi and f2 satisfies Assumption

---

Taking the reference measure ` to be the normal distribution satisfies Assumption
2
which
ensures the distribution is not too flat.
Taking quadratic interactions (()) = 2|2l² satisfies Assumption 3
Remark 1 (Cauchy-Probleem) To complete the arguments convergence equilibrium, we require sufficient regularity of solutions to the PDEs under consideration. In fact, it is sufficient if we can show that equations 4 ( and (6) can be approximated by equations with smooth solutions. Albeit tedious, these are standard techniques in the regularity theory for partial differential equations, for example [CMV03] Proposition 2.1 and Appendix AJ, [OV000, [Vil03] Chapter 9], and the references therein. Similar arguments as in [DV00) are expected to apply to the coupled gradiens flows considered here, guaranteeing existence of smooth solutions with fast enough decay at infinity, and we leave a detailed proof for future work.
3.2 Analysis of Case (i): Aligned Objectives
The primary technical contribution of this setting consists of lifting the algorithm dynamics from an ODE to a PDE, which allows us to model the system as a joint gradient flow on the product space of probability measures. The coupling occurs in the potential function, rather than as cross- coss-diffusion or non-local interaction as more commonly seen in the literature for multi-species systems.
Theorem 2. Suppose that Assumptions
3 are satisfied and let Xa == À + min(Az + +3, aअ) Consider solutions "/t = (pt: pe) the dynamics (4) with initial conditions satisfying 7o P2 P2((R P2 (2( and G1(r)) < oo. Then the following hold:
(a) There exists unique minimizer "Yoo (Pocr poo) of G., which is also a steady state for equation Moreover, Pro E L' (RR"), has the same support as p. and its density is continuous.
(b) The solution " converges exponentially fast in Gal- (ol-1 and W.
Ga(Te (oo) e -2au'Ga(To (oo) and W(1n: ce- for all t O, where >> 0 is constant only depending on yoo You and the parameter Aa-
Proof. (Sketch) For existence and uniqueness, we leverage classical techniques in the calculus of variations. To obtain convergence to equilibrium in energy, our key result is a new HWI-type inequality, providing as a consequence generalizations the log-Sobolev inequality and the Talagrand inequality. Together, these inequalities relate the energy (classically denoted by H in the case of the Boltzmann entropy), the metric (classically denoted by W in the case of the Wasserstein-2 metric) and the energy dissipation (classically denoted by / in the case of the Fisher information]" Combining these inequalities with Gronwall's inequality allows us to deduce convergence both in energy and in the metric W.
3.3 Analysis Case (ii): Competing Objectives
In this setting, we consider the case where the algorithm and the strategic population have goals in opposition to each other; specifically, the population benefits from being classified incorrectly. First, we will show that when the algorithm instantly best-responds to the population, then the distribution of the population converges exponentially in energy and in W2. Then we will show a similar result for the case where the population instantly best-responds to the algorithm.
In both cases, we begin by proving two Danskin-type results (see [Dan67 Ber7l)) which will be used for the main convergence theorem, including convexity (concavity) results. To this end, we make the following assumption ensuring that the regularizing component in the evolution of p is able control the concavity introduced by fi and fz.
Assumption (Upper bounds for fi and f2). There exists a constant A1 such that
(2,)) 3 Aa for all (z, x) TR.d x Rdd
and for any R 0 there exists constant C2 = c2(R such that
/ f2(2,))dp(2) < c2.
sup =x
'Hence the name HWI inequalities.
7

---

Equipped with Assumption 4 we state the result for a best-responding algorithm.
Theorem 3. Suppose Assumptions 4 are satisfied with at > A1. Let X = ad A1. Define
Go(p) = Ge(p. b()). Consider solution Pt to the dynamics 5
with initial condition Pa E PR(()
such that Go(Pa) co. Then the following hold:
(a) There exists unique maximizer Poo of Go(p), which is also steady state for equation S Moreover, Peo E L (R*). has the same support as p. and its density is continuous.
(b) The solution Pt. converges exponentially fast to Pra. with rate An in Ga(-I Poo) and W2.
Golpe Poo) S e- -22ha GPool Poo) and We(pe. S ce -sst for all 2 0,
where e > O is constant only depending on po, Poo. and the parameter Ab-
Proof. (Sketch) The key addition in this setting as compared with Theorem 2 is proving that Gp is bounded below, uniformly displacement concave and guaranteeing its smoothness via Berge's 's Maximum Theorem. This is non-trivial as uses the properties of the best response b(p). A central observation for our arguments to work is that G. = (6, Gepp. æI) (z=b(p)- We can then conclude using the direct method in the calculus of variations and the HWI method.
Here, the condition that at must be large enough corresponds to the statement that the system must be subjected a strong enough regularizing effect.
In the opposite case, where p instantly best-responds to the algorithm, we show Danskin-like results for derivatives through the best response function and convexity of the resulting energy in which allows to deduce convergence.
Theorem Suppose Assumptions are satisfied with ax > D,, and that r(x) is differentiable (as shown by example conditions in Lemmas 27 and 28 Define Ga(z) = (r)())), Then it holds:
(a) There exists unique minimizer 2 of Gx(z) which is also steady state for 66)
(b) The vector x(t) solving the dynamics
(6) with initial condition x(0) E Rd converges + X2
exponentially fast to Too with rate Xa-
0 in Ga and in the Euclidean norm:
||=(t) - Iooll S e [=x(0) xsoll: Ga(z(t)) Go( Se- 2dat (Ga(z(0)) Ga(x
---
for all O.
These theorems illustrate that, under sufficient convexity conditions the cost functions, expect the distribution and the algorithm x to converge to a steady state. In practice, when the distributions are close enough to the steady state there is no need to retrain the algorithm.
While we have proven results for the extreme timescale cases, we anticipate convergence to the same equilibrium in the intermediate cases. Indeed, it is well known [Bor09] (especially for systems in Euclidean space) that for two-timescale stochastic approximations of dynamical systems, with appropriate stepsize choices, converge asymptotically, and finite-time high probability concentration bounds can also be obtained. These results have been leveraged in strategic classification [Zm++T]] We leave this intricate analysis future work.
and Stackelberg games [FCR20 FR21 Fie+21..
In the following section we show numerical results in the case of best-repponding æ, best-responding p, and in between where and p evolve
similar timescale. Note that in these settings. the dynamics do not have a gradient flow structure due to a sign difference in the energies, requiring conditions to ensure that one species does not dominate the other.
Numerical Examples
We illustrate numerical results for the case of a classifier, which are used scenarios such as loan or government aid applications [CCIT], school admissions (PST3), residency match [Ree18], and recommendation algorithms [LSWIO], of which have some population which incentivized to submit data that will result in desirable classification. For all examples, select classifiers

---

of the form x E R. so that a data point 2 E IR is assigned label of 1 with probability ((2,z)
(1 + exp (-b" x))- where b > 0. Let fI and f2 be given by
==
(2,,))),
fz(2, log (2,
Note that Hess (f1) 0 and V2f22 0, so \ of the functional coming from the regularizers, not the cost functions, with p a scaled normal distribution. We show numerical results for scenarios with additional settings in the appendix. First we illustrate competitive interests under three different timescale settings. Then we simulate the classifier taking were implemented based on the finite volume method from (CCHI5).
O. Here, the strictness of the convexity
more naive strategy than gradient descent and discuss the results. The PDEs
4.1 Competitive Objectives
In the setting with competitive objectives, we utilize Ge(p, with W = 0, and S2 defined above with 6 = 3 fixed as it only changes the steepness of the classifier for 1, and = 0.1 and = 0.05. In Figure D we simulate two extremes of the timescale setting; first when is nearly best-responding and then when æ is best-responding. The simulations have the same initial conditions and end with the same distribution shape; however, the behavior of the strategic population differs in the intermediate stages. When p is nearly best-responding, we see that the distribution quickly shifts
best-responding
grobebi ity
0.75
t-0.0
t-1.0
P
t-177
t-5.0
0.50
0.25
Pa
Po
0.00
7
2
0
-2
4
-2
4
best-respond
0.75
t-0.0
t-15.0
arobabiiit
0.50
0.25
0
2
Figure 1: When æ versus p best-responds, we observe the same final state but different intermediate states. Modes appear the strategic population which simpler models cannot capture.
mass over the classifier threshold. Then the classifier shifts right, correcting for the shift in which then incentivizes p to shift more mass back to the original mode. In contrast, when best-responds, the right-hand mode slowly increases in size until the system converges.
Figure 2 shows simulation results from the setting where p and æ evolve the same timescale. We observe that the distribution shift in p appears to fall between the two extreme timescale cases, which we expect. We highlight two important observations for the competitive case. One, single-mode
and at similar rates
0.75
tw10.0
t=50.0
orobabilit
0.50
0.25
Pu
2
Figure 2: In this experiment the population and classifier have similar rates of change, and the distribution change for p exhibits behaviors from both the fast o and fast æ simulations; the right-hand mode does not peak as high as the fast case but does exceed its final height and return to the cquilibrium.
distribution becomes bimodel, which would not be captured using simplistic metrics such as the mean and variance. This split can seen as polarization in the population, phenomenon that mean-based strategic classification model would not capture. Two, the timescale which the

---

classifier updates significantly impacts the intermediate behavior of the distribution. In our example, when z updated slowly relative to the strategic population, the shifts in the population were greater than in the other two cases. This suggests that understanding the effects of timescale separation are important for minimizing volatility of the coupled dynamies.
4.2
Naïve Behavior
In this example, we explore the results of the classifier adopting a non- con-gradient-foon strategy, where the classifier chooses initially-subopting value for and does not move, allowing the strategic population to respond. All functions and parameters are the same as in the previous example. When
Classifier Loss
Po
Loss
-1.8
-1.9
probability
Population Loss
1.0
Lo33
0.2
0.8
20
40
0.0
time
gradient descent
naive
-2
0
2
4
(a) Both species minimize their respective losses; when the classifier uses naive strategy, the final performance is better for the classifier and uni- formly worse for the population.
(b) The classifier selects suboptimal initial condi- tion x = 2.2, instead of 1.5 which minimizes the initial loss, and then does not move in response to the population.
Figure 3: Although the classifier starts with a larger cost by taking the naive strategy, the final loss is better. This illustrates how model can be used compare robustness of different strategies against a strategic population.
comparing with the gradient descent strategy, we observe that while the initial loss for the classifier worse for the naive strategy, the final cost is better. While this results is not surprising, because one can view this as a general-sum game where the best response to fixed decision may be better than the equilibrium, illustrates how our method provides a framework for evaluating how different training strategies perform in the long run against strategic population.
Future Directions, Limitations, and Broader Impact
Our work presents a method for evaluating the robustness of an algorithm to a strategic population, and investigating a variety of robustness using techniques opens range of future research directions. Our application suggests many questions relevant to the PDE literature, such as: (1) Does convergence still hold with the gradient replaced by estimated gradient? (2) Can prove convergence in between the two timescale extremes? (3) How do multiple dynamic populations respond algorithm, multiple algorithms? In the realm of learning algorithms,, framework can be extended to other learning update strategies and presents a way to model how we can design these update strategies induce desired behaviors in the population.
challenge in our method is that numerically solving nigh-dimensional PDEs is computationally expensive and possibly unfeasible. Here note that in many applications, agents in the population do not alter more than few features due to the cost of manipulation. We are encouraged by the recent progress using deep learning to solve PDEs, which could be used in our application.
Broader Impacts Modeling the full population distribution rather than simple metrics of the distri- bution is important because not all individuals are affected by the algorithm in the same way. For example, if there are tails of the distribution that have poor performance even if on average the model is good, we need to know how that group is advantaged or disadvantaged relative to the rest of the population. Additionally, understanding how people respond to algorithms offers an opportunity to incentivise people to move in direction that increases social welfare.
10

---

Acknowledgments and Disclosure of Funding
LC is supported by an NDSEG fellowhip from the Air Force Office of Scientific Research. FH is supported by start-up funds at the California Institute of Technology. LR is supported by ONR
YIP N00014-22 1--000--------77 1-251- PO0003 and NSF Awards CAREER 1844729, and CPS 1931718. EM acknowledges support from NSF Award 2240110. 10. We are grateful for helpful discussions with José A. Carrillo.
References
[Ala40]
Leon Alaoglu. "Weak Topologies of Normed Linear Spaces". In: The Annals of Mathe- matics 41.1 (Jan. 1940), p. 252. ISSN: 0003486X, 001|110.2 2300.2307111 20072209//998829
[Dan67]
John M. Danskin. The Theory of Max-Min and its Application to Weapons Al- location Problems. Red. by M. Beckmann et al. Vol. 5. Â-konometrie und Un- lernehmensforschunng Econometries and Operations Research. Berlin, Heidelberg: Springer. 1967. ISBN: 9783642460944 9783642460920. DOL: 10.1007/978--3-0000 46092-0
[Ber71]
Dimitri P. Bertsekas. "Control of Uncertain Systems with a Set- Set-Membership Description of Uncertainty". PhD thesis. Cambridge, MA, USA: MIT, 1971.
[Pos75]
Posner. "Random coding strategies for minimum entropy". In: IEEE Transactions on Information Theory 21.4 (July 1975), pp. 388-391, ISSN: 1557-9654. DOI: 10. 1109/
TIT.1975, NII.1975.005541
[Lic83]
Elliott H. Lieb. "Sharp Constants in the Hardy-Littlewon Littlewood-Soboler and Related Inequal- ities". In: Annals of Mathematics 118.2 (1983), pp. 349-374. ISSN: 0003-486X. DOI:
10. 2307 200772200732
[Rud91]
Walter Rudin. Functional analysis. 2nd ed. International series in pure and applied mathematics. New York: McGraw-Hill, 1991. 424 pp. ISBN: 9780070542365.
[McC97]
Robert J. McCann. "A Convexity Principle for Interacting Gases". In: Advances in
Mathematics 128.1 (June 1, 1997), pp. 153-179, ISSN: 0001-8708. DOI: 10. 1006
aima.1997 19997.6634
[BB00]
Jean-David Benamou and Yann Brenier. "A computational fluid mechanics solution the Monge-Kantorove
mass transfer problem". In: Numerische Mathematik 84.3 (Jan. 1, 2000), pp. 375-393, ISSN: 0945-3245. DO1: 10. 0007//002110050000 10050002)
[DV00]
Laurent Desvillettes and Cedric Villani. "On the spatially homogeneous landau equation for hard potentials part existence, uniqueness and smoothness".. In: Communications in Partial Differential Equations 25.1 (Jan. 1, 2000), pp. 179-259. ISSN: 0360-5302.
DOL: 10. 1080/ 0080000605300008821555
[OV00]
F. Otto and C. Villani. "Generalization of an Inequality by Talagrand and Links with the Logarithmic Sobolev Inequality", In: Journal of Functional Analysis 173.2 (June 1, 2000), pp. 361-400. ISSN: 0022- DOT: 10. 1006100//00 10066 1999. 1999.3557
[Gho02]
B. K. Ghosh. "Probability Inequalities Related to Markov's ' Theorem". In: The American Statistician 56.3 (2002), pp. 186-190. ISSN: 00031305.
[CMV03]
José A. Carrillo, Robert J. McCann, and Cédric Villani. "Kinetic cquilibration rates for granular media and related equations: entropy dissipation and mass transportation estimates". In: Revista Matemática Iberoamericana 19.3 (Dec. 2003), pp. 97 971-1018. ISSN: 0213-2230.
[Vi103]
Cédric Villani. Topics optimal transportation. Vol. 58. Graduate Studies in Math- ematics. American Mathematical Society, Providence, RI, 2003, pp. xvi+370. ISBN: 0-8218-3312-X. DOI: 10. 1090/ 058
[Ste04]
J. Michael Steele. The Cauchy- Schwarz Muster Class: An Introduction to the Art of Mathematical Inequalities. Cambridge University Press, Apr. 26, 2004. 320 pp. ISBN: 9780521546775.
[AB06]
"Correspondences.. In: Infinite Dimensional Analysis: A Hitchhiker's Guide. Ed. by Charalambos D. Aliprantis and Kim C. Border. Berlin, Heidelberg: Springer, 2006,
pp. 555-590. ISBN: 9783540295877.
DOI: 10. 10070000///3 / 540-29587-99 29587-9_17

---

[CMV06] José A. Carrillo, Robert J. McCann, and Cedric Villani. "Contractions in the 2- Wasserstein Length Space and Thermalization of Granular Media". In: Archive for Rational Mechanics and Analysis 179.2 (Feb. 1, 2006), pp. 217-263. ISSN: 1432-0673.
DOL: 10. 10070000//// 00000/00020000000006666
[Del06]
Chrysanthos Dellarocas. "Strategic Manipulation Internet Opinion Forums: Implica- tions for Consumers and Firms". In: Management Science 52.10 (Oct. 2006), pp. 1577- 1593. ISSN: 1925-1909, 1526-5501. DOT: |0. 1287 22088/mmss. mnsc. 1060.006)
[De 07]
Camillo De Lellis. "Ordinary Differential Equations with Rough Coefficients and the Renormalization Theorem of Ambrosio". In: Séminaire Bourbaki. Vol. 59. 972. Mar. 2007, pp. 175-204.
[BCL09]
Andrea L. Bertozzi, Jose A. Carrillo, and Thomas Laurent. "Blow-up in multidimen- sional aggregation equations with mildly singular interaction kernels*", In: Nonlinearity 22.3 (Feb. 2009), p. 683. ISSN: 19715, DOL: 10. 1088/095111001 7715/22/ 22/322//02
[BBS09]
Steffen Bickel, Michael Brückner, and Tobias Scheffer. "Discriminative Learning Under Covariate Shift". . In: The Journal of Machine Learning Research 10 (Dec. 1, 2009). pp. 2137-2155. ISSN: 1532-4435.
[Bor09]
Vivek S Borkar. Stochastic approximation: a dynamical systems viewpoint. Vol. 48.
Springer, 2009.
[DY10]
Erick Delage and Yinyu Ye. "Distributionally Robust Optimization Under Moment Uncertainty with Application to Data-Driven Problems". In: Operations Research 58.3
(June 2010), pp. 595-612. ISSN: 0030-364X, 1526-5463. DOI: 10. 0741
1090
[LSWI0]
Juan Lang, Matt Spear, and S. Felix Wu. "Social Manipulation of Online Recommender Systems". In: Social Informatics. Ed. by Leonard Bolc, Marek Makowski, and Adam Wierzbicki. Lecture Notes in Computer Science. Berlin, Heidelberg: Springer, 2010, pp. 125-139. ISBN: 9783642165672. DOT: 10. 1007 0000000000-00000 900-9-78-239 124-2-1656722.00
[CC11]
Adriana Camacho and Emily Conover. "Manipulation of Social Program Eligibility"". In: American Economic Journal: Economic Policy 3.2 (2011), pp. 41-65. ISSN: 1945-7731.
[Car+11]
J. Carrillo et al. "Global-in-time weak measure solutions and finite-time aggregation
for nonlocal interaction equations". In: Duke Mathematical Journal 156.2 (Feb. 1, 2011).
ISSN: 0012-7094. DOI: 10. 12155000227744-2000 221
[BCY12]
D. Balagu's e, J. Carrillo, and Yao. "Confinement for repulsive-attractie kernels". In: Discrete and Continuous Dynamical Systems- B (2012).
[BLL12]
Andrea L. Bertozzi, Thomas Laurent, and Flavien Lager. "Aggregation and spreading via the newtonian potential: the dynamics of patch solutions". In: Mathematical Models and Methods in Applied Sciences 22 (supp01 Apr. 2012), p. 1140005. ISSN: 0218-2025.
DOL: 10. 1 1111//1221100/ 14000057
[FF13]
Marco Di Francesco and Simone Fagioli. "Measure solutions for non-local interaction PDEs with two species". In: Nonlinearity 26. 10 (Oet. 1, 2013), pp. 2777-2808. ISSN: 0000009771///10//2000000
0951-7715, 1361-6544. DOI: 10.
[PS13]
Parag A. Pathak and Tayfun Sönmez. "School Admissions Reform in Chicago and England: Comparing Mechanisms by Their Vulnerability to Manipulation", In: American Economic Review 103.1 (Feb. 2013), pp. 80-106. ISSN: 0002-8282. DOI: 10. 1257
aer 103. 1.800
[CCH14)
José Antonio Carrillo, Michel Chipot, and Yanghong Huang. "On global minimizers of repulsive-attractiv power-law interaction energies". In: Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 372.2028 (Nov. 13, 2014), p. 20130399, ISSN: 1364-503X, 1471-2962. DO1: 10. 1098 reta 2013. 2019.0399
[MKB14]
Alan Mackey, Theodore Kolokolnikov, and Andrea L.. Bertozzi. "Two-species particle aggregation and stability of co-dimension one solutions". In: Discrete and Continuous Dynamical Systems 19.5 (2014). pp. 1411-1436.
[CCH15]
José A. Carrillo, Alina Chertock, and Yanghong Huang. "A Volume- Method for Nonlinear Nonlocal Equations with a Gradient Flow Structure". In: Communications in Computational Physics 17.1 (Jan. 2015), pp. 233-258. ISSN: 1815-2406, 1991-7120.
DOT: 10. 4208/ 4208/cicp. 160214. 010814a
12

---

[San15]
Filippo Santambrogio. "Optimal Transport for Applied Mathematicians. Calculus of Variations, PDEs and Modeling". In: (2015).
[Har+16] Moritz Hardt et al. "Strategic Classification". In: Proceedings of the 2016 ACM Conference on Innovations in Theoretical Computer Science. ITCS '16. New York, NY, USA: Association for Computing Machinery, Jan. 14, 2016, pp. 111-122. ISBN: 97814 9781450340571. DOI: 10. 1145/ 10111111444222 2840728 2840730
[CHS18]
José A. Carrillo, Yanghong Huang, and Markus Schmidtchen. "Zoology of a Nonlocal Cross-Diffusion Model for Two Species". In: SIAM Journal on Applied Mathemat- ics 78.2 (Jan. 2018), pp. 1078-1104 1104. 17M1128782
ISSN: 0036-13300000 1399, 1095-712X. DOE: 10. 1137/
[Don+18
Jinshuo Dong et al. "Strategic classification from revealed preferences". In: Proceedings of the 2018 ACM Conference on Economics and Computation. 2018. pp. 55-70.
[Ree18]
Alex Rees-Jones. "Suboptimal behavior in strategy-proot mechanisms: Evidence from the residency match". In: Games and Economic Behavior. Special Issue in Honor of Lloyd Shapley: Seven Topics in Game Theory 108 (Mar. 1, 2018), pp. 317-330. ISSN:
0899-8256. DOT: 10. 1016/ geb. 2017 .0.011
[Kuh+19] Daniel Kuhn et al. "Wasserstein Distributionally Robust Optimization: Theory and Applications in Machine Learning". In: Operations Research Management Science in the Age of Analytics. Ed. by Serguei Netessine, Douglas Shier, and Harvey J. Greenberg. INFORMS, Oct. 2019, pp. 130-166. ISBN; 9780990615330. DOI: 10. 1287 /duuc 2019 0198
[DS20]
Tomasz Debiec and Markus Schmidtchen. "Incompressible Limit for a Two-Species Tumour Model with Coupling Through Brinkman's 's Law in One Dimension". In: Acta Applicandae Mathematicae 169.1 (Oct. 1, 2020), pp. 593-6111
ISSN: 1572-9036. DOI:
N0.1007/31000 1001107//1140000000000000
[DT20]
Manh Hong Duong and Julian Tugaut. "Coupled McKean- Vlasov diffusions: wellposed- ness, propagation of chaos and invariant measures", In: Stochastics 92.6 (Aug. 17, 2020), pp. 900-943. ISSN: 1744-2508. DOI: 10. 1080/17442506 2019. 1677663
[FCR20]
Tanner Fiez, Benjamin Chasnov, and Lillian Ratliff. "Implicit learning dynamics in stackelberg games: Equilibria characterization, convergence analysis, and empirical study". In: International Conference on Machine Learning. PMLR. 2020, pp. 3133- 3144.
[Mil+20]
John Miller et al. "The Effect of Natural Distribution Shift on Question Answering Models". In: Proceedings of the 37th International Conference Machine Learning. International Conference on Machine Learning. PMLR, Nov. 21, 2020, pp. 6905-691 6916. Juan Perdomo et al. "Performative Prediction". In: Proceedings of the 37th International Conference on Machine Learning. International Conference Machine Learning. PMLR, Nov. 21, 2020, pp. 7599-7609.
[Per+20]
[AB21]
Abdulaziz Alsenafi and Alethea B. T. Barbaro. "A Multispecies Cross-Diffusion Model for Territorial Development". In: Mathematics 9.12 (Jan. 2021), p. 1428. ISSN: 2227- 7390. DOI: /0t::1033990 30990/9
[Dom+21]
Carles Domingo- Enrich et al. A mean-field analysis of two-player zero-sum games. May 6, 2021. DOI: 10. 48550/ arwiv. 2002 2002.06277 stat]
arXiv: 2002, 06277 [es, math
[DN21]
John C. Duchi and Hongseok Namkoong. "Learning models with uniform performance via distributionally robust optimization" " In: The Annals of Statistics 49.3 (June 2021), pp. 3378-1406. ISSN: 0090-5364, 2168-8966. DO1: 10. 1212/22-A022004
[FR21]
Tanner Fiez and Lillian J Ratifff. "Local convergence analysis of gradient descent ascent with finite timescale separation". In: Proceedings of the International Conference on Learning Representation. 2021.
[Fic+21] Tanner Fiez al. "Global convergence to local minmax equilibrium in classes of nonconvex zero-sum games". In: Advances in Neural Information Processing Systems 34 (2021), pp. 29049 29063..
[TYZ21] Zachary Izzo, Lexing Ying, and James Zou. "How to Learn when Data Reacts to Your Model: Performative Gradient Descent". In: Proceedings of the 38th International Conference on Machine Learning. International Conference on Machine Learning. PMLR, July 1, 2021. pp. 4641-4650.
13

---

[LHL21]
Qi Lei, Wei Hu, and Jason D. Lee. "Near-Optimal Linear Regression under Distribution Shift". In: ArXiv (June 23, 2021).
[MPZ21] John P Miller, Juan C Perdomo, and Tijana Zmnic. "Outside the echo chamber: Optimiz- ing the performative risk". In: International Conference Machine Learning. PMLR. 2021. pp. 7710-7720.
[Wil+21]
Olivia Wiles et al. "A Fine-Grained Analysis on Distribution Shift". In: ArXiv (Oct. 21, 2021).
[Zm+21]]
Tijana Zmnc et al. "Who Leads and Who Follows in Strategic Classification?" In: Advances in Neural Information Processing Systems. Vol. 34. Curran Associates, Inc., 2021. pp. 15257-15269.
[AZ22]
Alekh Agarwal and Tong Zhang. "Minimax Regret Optimization for Robust Machine Learning under Distribution Shift". In: Proceedings of Thirty Fifth Conference on Learning Theory. Conference on Learning Theory. PMLR, June 28, 2022, pp. 2704- 2729.
[Giu+22] Valeria Giunta et al. "Local and Global Existence for Nonlocal Multispecies Advection- Diffusion Models". In: SIAM Journal on Applied Dynamical Systems 21.3 (Sept. 2022),
pp. 1686-1708. ISSN: 1536-0040. DOI: 10. 1137 1137/21M1 /215//2111459999
[HPS22] Georg Heinze, Jan-Frederik Pietschmann, and Markus Schmidtchen. Nonlocal cross- interaction systems on graphs: Nonquadratic Finslerian structure and nonlinear mobili- ties. Jan. 11, 2022. DOI: 10. 40550/ arxiv 2107 1 11289 arXiv: 2107 11289 [math]
[JPZ22] Ansgar Jungel, Stefan Portisch, and Antoine Zurek. "Nonlocal cross-diffusion systems for multi-species populations and networks". : In: Nonlinear Analysis 219 (June 1, 2022), 12800. ISSN: 0362-546X. DOT: 10. 1016/j na. 2022.. 1 128800
[LY22]
Guanlin Li and Yao Yao. "Two-species competition model with chemotaxis: well- posedness, stability and dynamics". In: Nonlinearity 35.3 (Feb. 3, 2022), p. 1329. ISSN:
0951-7715. DO1: 10 10. 1000/03611111 1066116544/ac4a00
[LFG22]
Fengming Lin, Xiaolei Fang, and Zheming Gao. "Distributionally Robust Optimization: A review on theory and applications" In: Numerical Algebra, Control and Optimization 12.1 (2022), pp. 159-212.
[Nar+22]
Adhyyan Narang et al. "Learning in Stochastic Monotone Games with Decision- Dependent Data". In: International Conference Artificial Intelligence and Statistics. PMLR.. 2022, pp. 5891-5912.
[Ray+22]
Mitas Ray et al. "Decision-dependent risk minimization in geometrically decaying dy- namic environments" In: Proceedings of the AAAI Conference on Artificial Intelligence. Vol. 36. 7. 2022, pp. 8081-8088.
[CFG23]
J. Carrillo, Alejandro Fernandez-Jimenez, and D. G'omez-Castre "Partial mass concen- tration for fast-diffusions with non-local aggregation terms". In: Apr. 10, 2023.
[Dou+23] Marie Doumic et al. Multispecies cross- difssdiifuss:ns:: from a nonlocal mean-field -feld to a porous medium system without self-diffusion. June 12, 2023. DOT: 10. 40550850/// aarXiv 2306.01777 arXiv: 23066 01777 [math]
[GG23]
Camilo Garcia Trillos and Nicolas Garcia Trillos. On adversarial robustness and the use of Wasserstein ascent-descent dynamics enforce it. Jan. 9, 2023. DO1: 001:|10.a 11.11..8858 48550/ arXiv. 2301 03662 Xivi [3011 03662 [cs, math]]
[HPS23]
Georg Heinze, Jan-Frederik Pietschmann, and Markus Schmidtchen. "Nonlocal cross- interaction systems on graphs: Energy landscape and dynamics". In: Kinetic and Related Models (KRM) (2023).
[Lu23]
Yulong Lu. Two-Scale Gradient Descent Ascent Dynamics Finds Mixed Nash Equilibria of Continuous Games: A Mean-Fiell Perspective. Jan. 8, 2023. DO1: 10. 48550/ arxiv 2212.08791 arXiv: 2212. 08791 [cs, math, stat]
[WC23]
Guillaume Wang and Lenaic Chizat. An Exponentially Converging Particle Method for the Mixed Nash Equilibrium of Continuous Games. Apr. 13, 2023. DOI: 10. 48550/ arXiv. 2211 012800 arxiv: 2211. 01280 [ca, math]]
[Qui+]
Joaquin Siinonero-Candell et al. Dataset Shift in Machine Learning. MIT Press.
14

---

A General structure and preliminaries
In this section, we give more details on the models discussed in the main article, and introduce definitions and notation that needed for the subsequent proofs.
A.1 Structure of the dynamics
For the case of aligned objectives, the full coupled system of PDEs 4 can be written as
0.p- = adpt div (pv- (ʃ fidy log alog W *
(7a)
dep div (uv. (f / sodp+ _ 2 2llee 2 200l" zoll²
(7b)
In other words, the population p in 7a is subject to an isotropic diftusive force with diffusion coefficient 0. a drift force due the time-varye time-varying confining potential f fıdu((t) - log p. and a self-interaction force the interaction potential W. If we consider the measure μ to be given and fixed in time, this corresponds exactly to the type of paraholic equation studied in [CMV03]. Here however the dynamics are more complex due to the coupling of the confining potential with the dynamies (7h) for s(t) via the coupling potential fi. Before presenting the analysis of this model, let give bit more intuition on the meaning and the structure of these dynamics.
the setting where j represents a binary classifier. we can think the distribution P as modelling all those individuals carrying the true label I, say, and the distribution p(t) as modelling all those individuals carrying a true label 0, say. where 0 and 1 denote the labels of two classes of interest. The term I [1(2.x) x)t, dr) represents a penalty for incorrectly classifying an individual at z with true label 0 when using the classifier 2). other words, J f2(2 dx) E 0, oo) increasingly large the more digresses from the correct classification O. Similarly, f f1 ptt, dz) E [0, oo) is increasingly large if the population p shifts mass to locations in where the classification is incorrect. The terminology aligned objectives refers to the fact that 7 both the population and the classifier are trying to evolve in way as to maximize correct classification. Analogously, the term [ f2(2, (2
large if a would incorrectly classify the population f that carries the label 1. A natural extension of the model (7) would be a setting where also the population carrying labels I evolves over time, which is simulated in Section E.2 Most elements of the framework presented here would likely carry over the setting of three coupled PDES: one for the evolution of p(t), one for the evolution of p(t) and one for the classifier (t)..
The term
aAp- - adiv (pV log - adiv (pV
forces the evolution of p(t) to approach p. In other words, it penalizes (in energy) deviations from given reference measure 5. In the context of the application hand, we take ý be the initial distribution p(t 0). The solution p(t) then evolves away from e over time due to the other forces that are present. Therefore, the term LL p) in the energy both provides smoothing the flow and a penalization for deviations away from the reference measure p.
The self-interaction term W +p introduces non-locality into the dynamics, as the decision for any given individual to move in a certain direction is influenced by the behavior of all other individuals in the population. The choice of W is application dependent. Very often, the interaction between two individuals only depends on the distance between them. This suggests a choice of W as radial function, i.e. (z) = w([=]), A choice of w - such that (r) > O corresponds to an attractive force between individuals, whereas w (r) < O conresponds to a repulsive force. The statement [2|=> ([=l) = z V. ((2) 2 - D in Assumption therefore corresponds to a requirement that the self-in self-interaction force is not too repulsive. Neglecting all other forces in 7a). we obtain the
non-local interaction equation
0,p = div (PV * p)
which appears in many instances in mathematical biology, mathematical physics, and material science, and it is an equation that has been extensively studied over the past few decades, see for example [Car+1l] BCY12 CCH14 BCL09 BLL12 CMV06 CFG23 and references therein. Using the results from these works, assumptions the interaction potential W can he relaxed in many ways, for example hy allowing discontinuous derivatives at zero for W, or by allowing W to be negative.
The dynamics /bb for the algorithm y is a non-autonomous transport equation,
dip = div (1v),
where the time -ene-dependence in the velocity field
vtt x)) V= (f [f((a.)ddo()))) + falz, / - 01== xoll?)
15

---

comes through the evolving population p(t). This structure allows to obtain explicit solution for p(() in terms of the initial condition μi and the solution p(() to 7a using the method of characteristics.
Proposition Assume that there exists constant such that
|/ (=fi(e,z)dp(2))+ (afele, / s <l |=]l))
4p P2(R*) and Var E Rd..
(8)
Then the unique distributional solution p(f) to 7b is given by
p(t) @tt, +)
(9)
where (f,, 8, x) solves the characteristic equation
0,0(s,t, + v(s, (8,, t,)) O,
(tt, æ.
(10)
Proof. Thanks to Assumption 1 we have that v E C (R x RE RR).. and by 8 we have
[]v(t, æ)1l s c(1 + [|=]]) for all 2 0, R
By classical Cauchy-Lipschitg theory for ODEs, this guarantces the existence of unique global solution
(t s, solving (T0). Then can be checked directly that p(t) as defined in ( is distributional solution to (7b).
In the characteristic equation 10
f, is parametrization of all trajectories: if particle was at location at time t, then it is at location f, at time Our assumptions fi, and also ensure that (ss t) Rdd + ied a C -diffeomorphism for all s, E IR. For more details on transport equations, see for example [De (7).
Remark 2. Consider the special case where μ4 6xg for some initial position so Rd. Then by Proposition the solution to (bb given by jutt) 52(1). where x(t) = (t,, 2o) solves the ODE
*(f) (t), x00) = 2o:
which is precisely of type
For the case competing objectives, the models consider written a,p- = div (p (V(f((z)b)) a log(p/ - W pl)
bp) = argmin / fi (z.2)))((2) + f2(2, /a(o,==)dp(=') +
B || rol
for S and
dt x = v (( + f2(U, ('(x,=')dp())) + 3 || -|ج--xal 2oll?
r(a) = argmax / - aKL(ppp)
1/2 pw *
for 6
A.2 Definitions and notation
Here, and in what follows, Ia denotes the x d identity matrix, and id denotes the identity map. The energy functionals considering are usually defined on the set of probability measures Rd denoted by ) of probability measures with bounded second moment,
we consider the subset (RR
Pa(R = (o E P(R)) f ||æl|=dp(x) < bo),
then we can endow this space with the Wasserstein-2 metric. Definition 4 (Wasserstein-22 Metrie), The Wasserstein-2 metric between two probability measures p, (R given by
Wa(p, inf rET(ie))
-1==
where I' the set of all joint probability distributions with marginals μ and v, i.e. j(dz) = f (dzz, z')dz' and v(dzz) 7/2, dz' ddz.
16

---

The restriction to P2(R") ensures that W2 is always finite. Then the space (Pr(R"), Wa) is indeed metric space. We will make of the fact that W2 metrizes statement precise, us introduce two common notions convergence probability measures, which are subset of the finite signed Radon measures ((R)).
convergence probability measures. To make this
Definition Consider a sequence (pn) E (M(R)) and limit μa E EM(R)..
(Narrow topology) The sequence (1a)) converges narrowly to denoted by pnn - if for all continuous bounded functions f: R R.
J fx))du(())) f())dp()).
(Weak-* topology) The sequence (Ha) converges weakly- to μ. denoted by pmn I p.. if for all continuous functions vanishing at infinity (ie. f. RR" + such that for all >> 0 there exists compact set K, C. g such that (f(z) < con I8" K.). we have
J f(2)dpm(2)
J f())dp(z).
Let us denote the set of continuous functions on IRd vanishing at infinity by Co(R",, and the set of continuous bounded functions CL (.B.). Note that narrow convergence immediately implies that pa (IR")
=(ll²) as the constant function is in Cs (TR)). This is not necessarily true for weak-s convergence. We will later make use of the Banach-Alaoge llaollu theorem (Al44). which gives weak-+ compactness of the unit hall in dual space. Note that M(R") indeed the dual of (I endowed with the sup-norm, and ((xd)) is the unit ball in M(R") using the dual norm. Moreover, we can ensure that mass does not escape to infinity, the two notions convergence in Definition are in fact equivalent.
Lemma 6. Consider sequence (mn) E MM(") and measure M(R"). Then phn - if and only if M4 - and pt, (R.) ((R").
This follows directly from Definition Here, the condition pan (I2d) - ((x)) and follows from Markov's 's inequality Ghoo22 if we can establish uniform bounds on the second moments, i.e. we want to show that there exists a constant C >> 0 independent of n such that
) is equivalent to tightness of (po)
f [=]))up(=)><CC
<
Vn E N.
(II)
Definition (Tightness probability measures). A collection of measures (pn) E M(RR)) is tight if for all 0 there exists compact set K. RRd such that [R,](R" K.) < for all N, where |μ] denotes the total variation of μ.
Another classical result is that the Wasserstein-2 metric metrizes narrow convergence and weak-x convergence of probability measures, see for example [San15] Theorem 5.11] or [Vil03] Theorem 7.12].
Lemma 7. Let pnny J Pa(RR"). Then W2(pn, p)) and only
Lf || |1=1)dd((z).
Pn.
and
gd
|l=l)d
Remark 3. Note that + can be replaced by uh in the above statement thanks to the fact that the limit is probability measure with mass I, see Lemma
Next, we consider two measures (4, P P( that are atomless, u((+}) for all z Rd. By Brenier's theorem [BO0) (also see [Vi003 Theorem 2.32]) there exists a unique measurable map T' : Rd Bd such that I Tw and T = V for some convex function Rd + IR. Here, the push-forward operator V
defined as
[ f(=)d = f ((V)(=))dpp(s)
for all Beasl-mmasurable functions f:R R+. If Velw we denote by P. [(1 -s) +VV
the discplacement interpolant between po and P1- We convexity, which the same as geodesic convexity in the geodesic space (Pz()))) W2). We will state
ready introduce the notion of displacement
definition here for atomless measures, but can be relaxed to any pair of measures in P2 using optimal transport plans instead of transport maps. In what follows, we will use to denote the interpolation parameter for geodesics, and to denote time related to solutions of and (6)
17

---

Definition 7 (Displacement Convexity). A fuuntiinall G P -> IR is displacement convex for all po. that atomlesy have
G(p.) s (1- (1-8))((p) + sG(pr),
where = [11 - s) id + is the displacement interpolant berween po and p1- Further, G:
P P
uniformly displacement convex with constant 1) 0 if
G(a.) S s)C(()) sG(pi) s(1 - (a
where = [11 - +d Po is the displacement interpolant between Po and pr.
Remark 4. In other words, is displacement convex (concave) if the function Gp.)) is convex (concave) with
P. [[1 - s sid +VVV) # Po being the displacement interpolant between Pa and p1. Contrast this with the
classical notion of convexity (concavity) for G, where we require that the function G -) + spt) is convex (concave).
In fact, if the energy G is twice differentiable along geodesics, then the condition d.a any geodesic (p=)ee(o,1) between pu and pr is sufficient to obtain displacement convexity. Similarly, when
d2 G(?+)) 0 along
dez G(..) 2 n/22(Po then G is uniformly displacement convex with constant 1) >> 0. For more details, McC97 and [Vi003 Chapter 5.2].
A.3
Steady states
The main goal in our theoretical analysis characterize the asymptotic behavior for the models 4 and (6) as time goes infinity. The steady states of these equations are the natural candidates to be asymptotic profiles for the corresponding equations. Thanks to the gradient flow structure, we expect to be able to make a connection between ground states of the energy functionals, and the steady states of the corresponding gradient flow dynamics. More precisely, any minimizer or maximizer is in particular critical point of the energy, and therefore satisfies that the first variation constant on disconnected components of its support. this ground state also has enough regularity (weak differentiability) to a solution to the equation, it immediately follows that it is in fact a steady state.
To make this connection precise, we first introduce what exactly we mean by a steady state.
Definition (Steady states for (.) Given Poo E (R*) Liee (R") with Il 1 and uoo E P2(R")
then (Poc, Poo) steady state for the system 4 if Poo E Wio. (R"), VW Poo (RR",, Poo is absolutely
continuous with respect to P. and (peas [oo) sarisfy
V. (f (1 + alog a
(Po(2) p(3)
+ W pz)
=
Vz E supp(p>>);
(12a)
V. (ʃ ((Ll + f2(2,z)dp(2)+ (ale + +
8
111x 1
201l3
=0
Ve supp(use)
(12b)
in the sense of distributions.
Here, L), (R()) := tp L' (IR2) p= 0}.
Definition (Steady states for (5). Let Poo E I.), (RR)) Lio (RR) with || 1. Then pee is steady
state for the system (5) if Poo E 1h (R.), VW poo E (B"),, Poo is absolutely continuous with respect and pso satisfies
V. ((((E.b(po))) a log
W pro(3)) =0 =0 Vz R"
(13)
p(z)
in the sense of distributions, where b(poe)== == argmin, Ge(.oo,
Definition 10 (Steady states for (6). The vector Tos E Rd is a steady state for the system ④ if satisfies
V=Ga(p=)=00) 0.
In fact, with the above notions of steady state, we obtain improved regularity for Poo- Lemma 8. Let Assumptions) 1 hold. Then the steady states Pao for and are contimuous.
Proof. We present here the argument for equation 5 only. The result for follows in exactly the same way by replacing fi(z, b(px)) with - fi(z, dd) (x).
18

---

Thanks to our assumptions. we have (i ((,))) + a log p(:) E C², which implies that V(i(( (- b b(pt)) +
log p()) Line. By the definition of steady state, Pro L Line and thanks Assumption have W C which implies that VW poo E L Let
h(2) P (2) (2) + + p(z) (W pre)
Then by the aforementioned regularity, obtain h E Loc n Liee By interpolation, it follows that h Lio for all 1 < < co. This implics that div (h) E Wine." Since Poo is weak Ioc -oluutio of (13) we have
Ap
div (h))
and so by classic elliptic regularity theory we conclude Pre. E Whe Finally, applying Morre''s 's inequality. we have Peo C where and for any < p < op. Therefore Poo C(R)) (after possibly being redefined set of measure zero). p
B Proof of Theorem 2
For case of notation, we write G. : x + 10, co]
Ga(le: = akL(() + V(a.)) W(p).
where we define
(( = // Si [((2,))d + / ((o)du(a),
W(p) = 2 // W 22)dp(21)dp(zz).
with potential given by ((x) f fz(z,=)dp(2) + 1 =ool"
In order to prove the existence of unique ground state for G.. natural approach consider the corresponding Lagrane equations
p(2)
p(z) f [((,,a)dp(())) + (W ()() = ei jp. 14]
a log
for all = E supp(p).
(14a)
[ + V(z) = calp. 14] for all x E supp(u),
(14b)
where Cı, C2 are constants that may differ different connected components of supp(p) and Supp()() These equations are easy to solve explicitly, and calculus of variations. We first show continuity and convexity properties the functional G. (Lemma and Proposition TO), essential properties that will allow us to deduce existence and uniqueness of ground states using the direct method in the calculus of variations (Proposition MT]. Using the Euler-Lagrange equation [4] we then prove properties on the support of the ground state (Corollary 12). To obtain convergence results, apply the HWI method: we first show general 'interpolation" inequality between the energy, the energy dissipation the metric (Proposition[1) [33]; this fundamental incquality will then imply generalized logarithmic Sobolev inequality (Corollary (4) relating the energy to the energy dissipation, and a generalized Talagrand inequality (Corollary 15) that allows to translate convergence in energy into convergence in metric. Putting all these conclude for the statements in Theorem
therefore using general non-constructive techniques from
ingrediends together will then allow
Lemma 9 (Lower semi-continuity). Let Assumptions hold. Then the functional G.. P P R is lower semi-comtinuous with respect to the weak- topology.
Proof. We split the energy G. into three parts: (i) K KL p). (ii) the interaction energy W, and (iii) the potential energy V. For (i), lower semi-continuity has been shown in [Pos75]. For (ii), we can directly apply [San Proposition 7.2] using Assumption 3 For (iii), note that V and fi are lower semi-continuous and bounded below thanks to Assumption [ and so the result follows from San Proposition 7.1].
Proposition 10 (Uniform displacement convexity). Let a, 9 > 0. Fix 70,70 71 E P2xx Pa and let Assumptions hold. Along any geodesic (7,)=εlu.j E P2 x Pa connecting To to 7h, we have for all 8 [0, 1]
d²
a G.(T%) 2 A, A. W(no.))²
A= min( + B,a
(15)
As result, the functional Ga : P R is uniformly displacement convex with constant Ao 0.
19

---

Proef. Let To and yi be two probability measures with bounded second moments. Denote by :. RR the optimal Kantorovich potentials pushing po pl, and Ho μr, respectively:
W2(po, pz)² J ||2 Vo(()))
P1 Vo
such that
Vd such that
W2(Ho: 11)² = 12d || V
The now classical results in BB00 guarantee that there always exists convex functions o, & that satisfy the conditions above. Then the path (7x (T)uu] = (ps, ) defined by
= [1 8) id id+
H, [(1 iddd # #/o
is a W -geodesscc from "Yo to ?
The first derivative of V along geodesics in the Wasserstein metric is given by
ds V(7a) ds [/ [1((1 8)2 + vVo(z), - s)x VVb(z)) dpoo))dde)))
+ +/v(1 v((1-
5)x s)x +) duo(z)
// V. vV
(Ve(<z) - x) (po())duo(T)
// V.fi((1 - 5))= -)) (1-))))((1))))))0000 (Vo(z) 2) (po(s)dpooz))
+ +/ V.V((1- =)x+ (Vb(2)- dp)(z),
and taking another derivative we have
(V((x)) - =77"" dsa V(Ta) = //
d²
(Vo(2) - ] -,(2,
[V)())) =3]] dpo(=)dpo(z))
[Vo((=)
+ // (Vb(x) =)" viV((1 - 8)L veb(()) (V+(() 2) dpo(s)dppo(a
n² (1) where we denoted D, x) = Hess(f))(( -)=--==--- s)2 + vVo(z), (1 s)x + +Vv(z)) and the last incquality follows from Assumption and the optimality of the potentials & and v.
2 Au W (10,
(Au + W)
Following CMV03 Vi103 and using Assumption the second derivatives of the diffusion term and the
interaction term along geadesics given by
d²
d²
dsa KLppppp) 2 at p1)
dsa W(p)) 20.
(16)
Putting the above estimates together, we obtain
Remark 5. Alternatively, one could assune strong convexity of W, which would improve the lower-bound on the second derivative along geodesics.
Proposition 11. (Ground state) Let Assumptions 3 hold for a, 3 > Then the functional G. : P(RR) x
P(RR) + [0, oo] admits unique minimizer 7 P2(R"), ), and is absolutely continuous with respect to p.
(pe, J.), and it satisfies p. P2((R)) L (Rd). p
Proof. We show existence of minimizer of G. using the direct method in the calculus of variations. Denote by y= (p, E PxP C MMMM a pair probability measures a point in the product space of Radon measures. Since G. 2 on P x P (see Assumption []] and not identically 4oo0 everywhere there exists minimizing sequence (7m) E x P. Note that (xx) is in the closed unit ball of the dual space of continuous functions
vanishing at infinity (Rd(Rd) x Co(R)) endowed with the dual norm l/ = sup
[f -
gdgeel over [f
f.g E (R() with [I 00 == |l + ||g]loc
By the Banach-Alasn Alaoglu thcorem Rud91 Thm 3.15]
there cxists a limit 7. = (p-, nm) E M x M = (Co x Co) and convergent subsequence (not relabelled)
such that 7oe I Jv. In fact, since K LLpp p) < oo it follows that p. absolutely continuous with respect to p.
implying P. L (RR) thanks Assumption 2
Further, μaa has bounded second moment, else would have
20

---

xpx Ga(7)= = which yields a contradiction. It remains to show that f f djr to conclude that 7, E P x P. To this aim, it sufficient show tightness of (px)) and (po), preventing the escape of mass infinity as we have [ ['pp. = f dyn
all n 2 Tightness follows from Markov's 's inequality [Gho02] we can establish uniform bounds on the second moments, i.c. we want to show that there exists a constant C > independent of n such that
/ [I=IPapo(s) + ||=l) <C
Vn N.
(17)
To establish (77. observe that thanks Assumption
2 there exists constant co E (possibly negative) such
that log p()) 2 co |ll=|/ for = E R2. Then
at
|= d= -cc a / log p(z)dpm
Therefore, using S dpp djn, and writing Ç min aĀ bound on the second moments of the minimizing sequence,
2) > O, we obtain the desired uniform upper #
sʃʃ (Il=ll [1x)l() dpn.um. I -acs
log logg(2)do +
|| zoll + 31
S -aco + 01|=o| 2| + Ge(Tn) S -aco + =| G.(m) <
This concludes the proof that the limit satisfies indeed 7, E x.P, and indeed E (PR" as well. Finally, 7 is indeed minimizer of C. thanks to weak-* lower-semicontinuity of C. following Lemma 9
Next we show uniqueness using contradiction argument. Suppose 7 (p., p.) and 7/4 (pl. pd) minimizers of Ga. For s E [0, 1J. define % : ((1 8) id +sT, (1 - s) id + where T, S R + R are the optimal transport maps such that p. p- and μ.. Syp.. By Proposition 10 the energy G. is uniformly displacement convex, and
have
Ga(T.) S (1 - s)Ga(7)) + sGa(T.) = G.(1-).
If # 7 and E (0, 1), then strict inequality holds by applying similar arguments in McC97 Proposition 1.2]. However, the strict inequality (T(7)) < G..))) for 7 * 7 is contradiction to the minimality of 7 Hence, the minimizer is unique.
Remark 6. A O, then the strict convexity of fi can be used to deduce uniqueness, and the axsumptions on log can be weakened from strict convexity convexity.
Corollary 12. Let Assumptions 7 hold. Any minimizer 7, = (p-,)e) of Ga is a steady state for equation according to Definition 8 and satishes supp(p.)) = supp(p).
Proof. By Proposition we have pa, Hr, P2, as as p. E 4., ||p. = continuous with respect to p. Since W E (R²),, it follows that VW P. E Lloc- In order to show that r,. is steady state for equation ( it remains to prove that p. E w n Liee- As As is a minimizer, is in particular a critical point, and therefore satisfies equations 14). Rearranging, we obtain (for possible different constant (p., [.] o 0) from (1a)) that
that p. absolutely
[-d (f ((lt,))a(o)+WWpp))
(=) ca calp,, n. (l))
supp(p.).
(18)
Then for any compact set K C Rd.
(-2 ( ((le))))o)))
(- W p.
sup (x(x) calp-, u-] sup p(2) sup exp =ek
sup exp =EK
zek
sek
As fi 2 by Assumption I and W 2 by Assumption
3 the last terms finite. The first supremum is finite thanks to contimuity of p. Therefore p. Lice- To show that E wh12 note that for any compact
set K C Rd. we have Jx (== < oo as consequence of p E Lofe- Moreover, defining Th)(=)) = (ʃ f1(2 u(x) + W p(2)) < O, we have
/ K =cle ' / |V [V (xp(2
S 2201lp. / /Dp exp(2T + 2es 2e1 m.. / ee
ex.(22
21

---

which is bounded noting that exp(2T(n-]) I and that T.))(.). <+]((() and Vp are in Lioe where we used that f1 (..z), r), (()) (() E (R*) by Assumptions 3 We conclude that pa E Wie and indeed (p-, (co) solves (22 in the sense of distributions as a consequence of (14).
supp(p) using again the relation (8) Firstly, note that supp(p.) _ supp(p) since P. is absolutely continuous with respect to P. Secondly, we claim that exp I-a (ʃ (L + W p=()))] > for all z E IR.. In other words, we claim that fi 2) (x) < oo and W p2) < oo for all & E Rt. Indeed, for the first term, fix any & R
Next, we show that supp(p-)
and choose R > O large enough such that E (0). Then, thanks to continuity of fi according to Assump-
tion T we have
/ f1(2,)) < oc. =EB
sup
the second term, note that by Assumption
we have for any 2 R and e
W(2) W(0) VW(x)
W(0 + 2f [VW(2)]]
+ 2 11=1
S W + D
(1+ + < W(0) + D²
(
12
|=|
2:
D
( + + =|1
v2
where the last equality follows by choosing the optimal = V2D. We conclude that
W * px(2) < (() + + V2D / ||= 2 (x())
D
V2
D
S (0) +
+ 22/2D11=11 + 2V2D |=|l p.(3),
(19)
V2
which is finite for any fixed E Rd thanks to the fact that p. E P2(Rd). Hence, supp(p.) = supp(p).
Remark If we have in addition that " E L (RC),, then the minimizer p. of Go is I (R*) as well. This follows directly by bounding the right hand side of (18).
The following inequality referred to as HWI inequality and represents the key result obtain convergence cquilibrium.
Proposition 13 (HWI inequality). Define the dissipation functional
Da(r) = // Ga.
Assume ax, B and let da defined in 15 Let 70, 72 E Pa Pa such that Ga(zu), Ga(m)) co. Then
Au W (1)
Ga(ro) Ga(m) & W(70, m) v Da(70)
(20)
2
Prouf. For simplicity, consider quo 72 that have smooth Lebesgue densities of compact support. The general case can be recovered using approximation arguments. Let (ys))ejo,1) denote W W-geodesic between foo n Following similar arguments as in [CMV03) and (OVO0) Section 5) and making use of the calculations in the proof of Proposition 10 we have
=Iʃ [z) [2(x)
(Vo(z) - =3] d70(2, (Ve((x) -
Ga(Ts)
ds
===
where
po(2)
(=fi(z,,))d((0o) + / V.W((= V.W((==2)dpp(1 z')dpo(='), p(z)
(i(To[(=) V. log
(2/10(() / v. (2,)))po) v.V(s).
Note that the dissipation functional can then be written
Da(70)= = // (I& + (ea(z))?) d70(2,
22

---

Using the double integral Schwarz inequality (Ste04). we obtain
d
2 --
Ei Ea
Vo(z)= = (v//
Ga(T))
1
dho
Ve((x)) - d70
ds
==0
/ [[b((2) 2| + [IV((z))) x ddo
v (70)
VD. W ()).
Next, compute Taylor expansion of Ga(T)) when considered as a function in and use the bound on da
de? G. from [5p
= = + d
Ga(xs)
+ ==0 +%o f.-u-
(
dp ds
di
ds
2 Ga(ro) - /D (((o, d W(10
Remark 8. The HWI inequality in Proposition L3 immediately implies uniqueness of minimizers for Ga in the set {ry E P P Da(r)) < t. Indeed, if to such that D. = 0, then for any other "n in the above set we have Ga(70) ) S Ga(r)) with equality if and only if W ("oo m1) = 0.
Corollary 14 (Generalized Log- Log-Sobolev inequality). Denote by 7/. the unique minimizer of Ga. With Xa as defined in 1S any product measure E P2 x P2 such that G(r)), Du < satisfies
D.(7) 2h, Gu(7l7))-
(21)
Proof. This statement follows immediately from Proposition Indeed, "/1 Y and 7o = Then
Aa W 71)
Ga(rl W(> Da(t)
2
S max 120
Da((t)t
A
Aa
{²
= Da(1)
2
2Aa
Corollary 15 (Talagrand inequality). Denote by 7t. the unique minimizer of Ga. With defined in holds
W(r
No Ga(xln))
for any E P2 x P2 such that G1(T) <
Proof. This is also direct consequence of Proposition 13 by setting To and ?1 = Then Go(7)) and Da(7+) and the result follows.
Proof of Theorem 2 The entropy term f plog P produces diffusion in p for the corresponding PDE consequence, solutions P: to 4 and minimizers for G. have to be L' functions. As there is no diffusion for the evolution of He, solutions may have singular part. In fact, for initial condition μo 8xg. the corresponding solution will be from the fact that the evolution for fle is tramsport equation (also see Section A. for more details). Results (a)
form f 62(1). where s=() solves
ODE ( with initial condition æo. This follows
and (b) are the statements in Proposition [T] Corollary 12] and Corollary [5 To obtain (c), we differentiate the energy G. along solutions ?/2 the equation 4
d Ca(a) f 8, (,GG(ll())0pppdd dt
= - |V (o. o.
1ial--IIvvea (V du
= -D.(m)) < -
2X, Ga(Tu (7)))
where the last bound follows from Corollary 14 Applying Gronwall's inequality, immediately obtain decay energy.
Ga(% (x) 2dat Ga 1x).
Finally, applying Talagrand's inequality (Corollary 15). the decay in energy implies decay in the product Wasserstein metric,
W W(r,, y-) ce Aat
where 0 is constant only depending 70. 74 and the parameter Aa-
23

---

C Proof of Theorem 3
In the case of competing objectives, we rewrite the energy Ge(p,z) P(RR)) x R - [-00, follows:
+ 3
G.(p,x)= f1(2,))d0(=)+ fz(2,z)dd((2++
2 1| - xall - P(p))
where
P(p) = aKL(plT) +
pW
Note that any fixed p P. the energy Gelp, -) is strictly convex æ, therefore has unique minimizer. Define the best response by
b(p) = argmin Ge(p,
denote Gp(p) = b(p)). begin with auxiliary results computing variations of the best. response b and then the different terms in (b(p) using Definitio [
Lemma 16 (First variation of the best response). The first variation of the best response of the classifier at fif it exists) is
8, b,b((())) = V..o)-VV.. for almost every E B
where (p)) (8 A2) symmetric matrix, constant in and x. defined
@(p) BIa V=.li(eb V..a((, (..
In particular, we then have for any vi E C (Rd) with J wd = that
b(p evll blpl | / 6,b10)())0())d
lim
Proof. Let di E C2 (Ied) with f de and fix > 0. Any minimizer of Gelp + esb, x)) for tixed must satisfy
V. bb, + eg)) =0
Differentiating in e, we obtain
f 6, V.Gepp oò, b(p ct)) VGGopp cbb b(p evi) / 6oblp + (c)()))( (22)
Next, we explicitly compute all terms involved in
Computing the derivatives yields
V=Gelpxz= Vdd(() V=.f(e, + B(x - xo)
6, V. V.G.)),)(>) -.f((z,a)
V) Vifi(zz
(2 BLa.
Note that VaGa is invertible by Assumption which states that fi and fa have positive-defiite Hessians. Inverting this term and substituting these expressions into (22) for e = gives
f (,)b))()))(()) BTa (2)((e
= =- Vv.
Finally, the lower bound on ((p) follows thanks to Assumption
Remark 9. we include the additional assumption that E Ca Rd: ox)) for 1,2, then the Hessian of blp) dell-defined. More precisely, the Hessian is given by
d²
dez blo blp + euj)c== Q(p)
de Qlp ev)l V(f((sbb)))))))d Q(p) "u(p)
where
= = V
and bbt / 0r, V.fa(z, (2D)) D - ulp. v).
Qs(p+ - / 0r, ,z, (is. / Volo: -(o uupp.
de
Therefore, can Taylor expand b(p) up to second order and control the remainder term of order ?
24

---

Lemma 17 (First variation of Gb). The first variation of Go. is given by (,Golp](z) (1)= hz(2)) + Blu(2) 6,L(ol(2),
where
8 (/ fn (EE
(E) ( Vz b()))dp(2). 5b
30 [p](2)
(1(=) 5p
+ +(((,,(())).
8p
bz(z) 8p (ʃ [E(E.((o))dp(E)))
(2) ( Vala(E, (()))d(()) [p((s)
ōb
hu(=) = 2 öp [b(p) =all? ((()) mo. 6p [+](z)
and
(,1()((2) = a llog(p())/(())) + (W * p)))
Proof. We begin with general expressions for Taylor expansions P(RR) - R and f+(z,b()) P(RR)
R for i 1,2 around p. Let 1h E T with T = { ((2)d= 0}. Then
b(p ew) b(p) /
86
(= +
(23)
8p
and
f((2, 6(p+ = f((,,(())))) + (V.Silz, k())),
8b ( + +0()).
(24)
We compute explicitly each the first variations:
i) Using
24 we have
/ (()))1(2)dz
/ + lim
++)=))d) +V) -
-// Va bo)))d(())
/ (p(ce'm(s))a (((;;)))))))))
5b(0)
8p
= =/ (/ Vali(z. (o)))
5b(0)
(l(z') "(=')dE'+ + (1(2,
8p
5b = (=) V.,iEE, (o))) op [1(z)) + (=,bb(()))
(ii Similarly, using again 24
/ (c))ba())d
lim 1 [ fe(2, b(p b0 e0))dp(2) / falt,
--d
= - Vufa(z, (()))
6p (D((e) u(z)dz
ha(z) V.fa(E, (o)))p()),
30 [)(()))
8b
(ii) Finally, from 23 it follows that
1
(1))ho(z)dz
lim -+0 2
(blp evd) zo, b(p eth) mo) (B(p) zu, b(p) xz))
= ( (b(0) - xo, 5p (0((z) d(())dz
5b
ha(2)) = b(p) vo, 5p lp(())
Finally, the expression for P[p] follows by direct computation
Lemma 18. Denote Gp(p) Go b(p)) with ((p) given by (5) Then 8,Golpl 6,Golp) (
25

---

Proef. We start by computing 8, x,G.(t,)))p(( for 2, E R"
(.Ga(t.w))ppe = (1,, Plp((()
(25)
Next, we compute 5,Gr. Using Lemma 17 first variation G. given by
()G()p((z) = h1 (z) hz(2) (z) (z)()((2)
= (1/ V=file, (o)))p())
V=Ta(E + +Bb()))))B)))))))))))))))))))))) 200] (5b())(2)
+ (((2, - 5ppp((()).
Note that
v.G.lp.)) V=f((E2))d(()) + V=
f(x - -20).
(26)
and by the definition of the best response b(p). we have (p = O. Substituting into the expression for Go and using (25). obtain
6, Gololl())) = ((p) -
-,00)))-_,,,0000 P(o((( = 5, G.G((s)))l((a))
b=b(())
This concludes the proof.
Lemma 19 (Uniform boundedness of the best response). Assumption hold. Then for any p E P(RR)) we have
2(a1 a2)
||b()]]1 < ||=o|]2
B
Proef. By definition of the best response b(p), we have
f (e. + V=.f(z,b(p))ddp + B(KP)) - -x0) = =0.
To show that that b(p) is uniformly bounded, we take the inner product of the above expression with b(p) itself
(1(((p)]1 = 32n b(p) V=fi(e,b(p)) b())pp(z) V..f(2,, b(p)) b(p)dp(z).
Using Assumption on the right-hand side, we obtain
I to bound the two integrals, together with using Young's 's inequality to bound the first term
8 B|b(p))1 S |=oo| 1 [[b(e]] aj + a2.
which concludes the proof after rearranging terms.
Lemma (Upper semi- coniicontinuity,. Let Assumptiios I holld. The functional G. P(R)) (-00 +oo is upper semi-continen continuous when P(R") x B is endowed with the product topology of the weak- topology and the Euclidean topology. Moreover, the fanctional Gs P(RR) > [-o0, t is upper semi-continuous with respect to the weak-# topology.
Proof. The functional G. P(RR) Rd 1-00, is continuous in the second variable thanks to Assumption I Similarly, f f2 x)dp(2) + zd7() continuous in p thanks [San15 Proposition 7.11 using the continuity fi and fz. Further, -P is upper semi-continuous using [Pos75] and [San15] Proposition 7.2] thanks to Assumptions(]andll
This concludes the continuity properties for G..
The upper semi-continuity of Go then follows from a direct application of a version of Berge's maximum theorem AB06 Lemma 16.30]. Let = ||a0||| 2(a:+az)
0. We define y (P(RR)) Wz) - IR the correspondence that maps any P(RR") to the closed ball Bx( C Rd. Then the graph of is = P(R") x (B.(O)) With this definition of the range of is compact and is continuous with respect weak-+ convergence, and so it in particular upper hemicontinuous. Thanks to Lemma [9 the best response function b(p) is always contained in (0) for any choice of E P(RR). As a result, maximizing -Ge(p, x) x over Rd for fixed p P(Id)) reduces to maximizing it over Bn Using the notation introduced above, we can restrict G. to Ge Gr p + and write
Go(p) = max #EP(P)
(G(
Because (p. .r) is upper www.--oomimmonnooo continuous when ((R)) x I is endowed with the product topology of the weak- * topology and the Euclidean topology, [AB06 Lemma 16.30] guarantees that G.(-) is upper continuusu in the weak- topology.
L
26

---

Proposition 21. Let 3 > 0 and assume Assumptions 1
4 hold with the parumeters satisfying ax >> A1. Fix
pa, ((R)). Along any geodesic (p=)sE (2
connecting po pt. have all E [0.
d² ds
Golps) S -2o
Aa ==
(27)
As result, the functional G1 P2(R 1-00, tool is uniformly displacement concave with constant A 0.
Pruof. Consider any Po p E (R Then any Wz-geodesic (pa)sEju.aj connecting with solves the following system of geodesic equations:
Osps + div (pav.) = 0, v.(,.v) div (p.v. @ vs) =0, 0,
(28)
where R² + 18 and vu
Rd + 1e The first derivative of Gb along geodesics can be computed explicidly
Gold.)= - / (x)) d d ds
P(o.)
ds
([/ x)dp.(;) + Vela(E, + 3(x 20)
d
b(px)
ds
x=b(o,))
The left-hand side of the inner product is zero by definition of the best response b(p=) to Pa. see 26 Therefore
Go(p.) V=,i(z,b(p.)) (=(:), ds P(p.).
ds
Differentiating a second time, using 28 and integration by parts, we obtain dt
=o(p.)= Li(p)) L2(ps)
do² P(p.),
where
Li L.(p)) = / V.,i(zzb(p.)) (v. 8 v.. p.dz= / ((mm v2hi(z,((p.) v.) p,dz,
Lz(p.) = f
d b(..) VeVasi(s,b(o.) vo(2) pz(2)z..
ds
From 16 have that
d² P(p.) 2 aā aaa (z
ds²
and thanks to Assumption Dit follows that
Li(s) S A.W
This leaves La bound; we first consider the bp.))
ds b(p.) / 8,blos((3)0,pp = = (,bl (p
/ V=o,blo-)(3) (3))))))
Defining u(ps) E Ie by
u(p.) = V=V=fi(2, - v. vz(2)dps(=),
using the results from Lemma have
16 Vz5,b(pa), Assumptioo D and the fact that Q(p) is constant in and x, we
L2(p=) = /f [1 (3)] V. V.V.f( dp
=- (u(ps). ((p.)) 5 0
Combining all terms together, obtain
d²
dsa Golp.) (ai-A1) 1 Aa) Wz(oo.pa)
27

---

Remark 10. Under some additional assumptions on the functions f1 and fz. we can obtain an improved convergence rale. In particular, assume that for all 2, R..
there exists a constant A2 2 X2 2 such that V,f2(2, 3 A2
there exists constant 2 0 such that VV 20.
Then we have - -(ps))++ -- -1/(3 + A1 + A2) la. Using Lemma 16 we then obiain a stronger bound on L as follows:
|a(e)))1 < 8+ /
IVVaafle dp.(s) / [us(=)]l
L2(Ps)
B + N1 + A2
0²
1 3 A + A2 Ws(po.
This means we can improve the convergence rate in to As == ai +
B+di+A2 An.
o²
Proposition 22 (Ground state), Let Assumptions hold for aZ > A 2 and > 0. Then there exists a unique maximizer pa for the functional Gb over P(R.), and it satisfies pa E P2(RR)) o (RR)) and p. is absolutely continuous with respect p.
Proof. Uniqueness of the maximizer (if it exists) is guaranteed by the uniform concavity provided by Lemma 21 To show existence of maximizer, following key properties for Gs: (1) boundedness from above, (2) upper contininn,,,,,, and (3) tightness
the direct method in the calculus of variations, requiring the
of any minimizing sequence. To show (I), note that (2,,( + +a (()) I -(aX A1) Ia for all 2, x E Rd x Rd by Assumptions and and so
f1(2,0) + t co(x) (a) A1) 12|
((2 E x Rd
(29)
with co(x) == =0 + aolog - V. (0. a log ((0] 11.. Therefore,
Go(p) / bz))) b(p)) + a + log d(()) + fa(2,b()))d7(=)) + _ (lb(p) =zoll
-af plugp closs-fove DW>P
2 (1b))) xoll?
S co(((())) + (2(2,(()))7((s) +
To estimate each of the remaining terms the right-hand side, denote It. || 26r) and recall that
[[b(p) || < R for any p E PPRRd) thanks to Lemma 19 By continuity of fa and log p. there exists a constant
Ci E R such that
sup co(x) sup (0)
[1(0,z)) + a log p(0) + ai - A
|/ (fi 2) log (())
ci (30)
(0)
The second term is controlled by ez thanks Assumption obtain
And the third term be bounded directly
Go(p) S + + B(R² + |
This concludes the proof of (1), Statement (2) was shown in Lemma 20 Then we obtain minimizing sequence (px) P(R)) which in the closed unit ball of (IR")" and the Banach-Anaoglu theorem [Rud91 Theorem
3.15] there exists limit p. in the Radon measures and subsequence (not relabeled) such that Pr: : p.. In fact, P. is absolutely continuous with respect to ò as otherwise Go(p.) -oo, which contradicts that G.(-) > -00 somewhere. We conclude that p. E L' (R") since E L' (R") by Assumption 2 To ensure p. E P(R)). require (3) tightness of the minimizing sequence (p.)). By Markov's inequality Gho02] it is sufficient establish a uniform bound on the second moments:
f [l]]pdp(<))< <C
E.
(31)
To see this we proceed a similar way as in the proof of Proposition I0 Defining
K(p) = - (fi b()) aog p())) a plogpde+ + oW * pdz.
28

---

we have (p) = -G((p) (2(2,b())ddp + 1b(p)) woll Then using again the bound on b(p) from Lemma
19
K(p) - GG(()) + sup / 22(2,2))dp(2++ 8 (1² + xe
-GG(p) + + B (1² + |
where the inequality thanks to Assumption Hence, using estimates noting that the sequence (pr) minimizing
29 and 30 from above, and
Gb), we have
(ax- An) / |= [==1/2ppn33) s co(b(pn))
/ [[((z, log p(=)] (p)
I c + k(p)) _ c1 - Go(pa) + + B ( + || c Gb(pr) + + 3 (R² + ||ab|| #) < x0.
which uniformly bounds the second moments of (pa). This concludes the proof for the estimate 31 and also ensures that p. E (2(
Corollary 23. Any maximizer p. of G. a steady state for equation according Definition and satisfies supp(p.) = supp(p).
Prouf. To show that p. is steady state follow exactly the argument as in the proof of Corollary [2] just replacing - fi (2, (x) with f1 b(p.). It remains show that supp(p.) supp(p).
X² is maximizer, it is in particular critical point, and therefore satisfics that 6,Go [e-](2) constant on all
connected components of supp(p.). Thanks to Lemma [8 this means there exists constant c(p.] (which may be different on different components of supp(p.)) such that
p=(=)
f1(2,b(p.) log
W p.(2) clpa]
on supp(p.).
((2)
Rearranging, we obtain (for a possible different constant c/p.] +
p. (z) = clp-1 exp
(fi(2,bb - W pp)) p.))))
on supp(e.).
(32)
Firstly, note that sp.) supp(p) since p. is absolutely continuous with respect p. Secondly, note that exp _ fi b(p.)) 2 1 for all E Rd since fi 2 0. Finally, we claim that exp (- W pa (z)) > 0 for all = E Rd In other words, we claim that W * p=(2) < do for all E Rd. This follows by exactly the same argument as in Corollary 12 see equation (1).. We conclude that supp(p.) = supp(p).
Remark 11. If have in addition that (KE) and fi(:,)) L*((() for all E Bd. then the maximizer of Gy is in I (R") as well. This follows directly by bounding the right-hand side of (32)
With the above preliminary results, inequality and generalized logarithmic Sobolev inequality.
show the HWI inequality, which implies again Talagrand-typ
Proposition 24 (HWI inequalities). Define the dissipation functional
Delm) = /ʃ (V,,GGolllaa)aattt))
Assume a, B such that at As + o², and X as defined (2 Denote by pa the unique maximizer of Gb.
(HWI) Let Pa, P2(RR) such that Golpo), Do(po) xc. Then
Golpo) Go(pi) ((o. Do(oo) An
pi))
(33)
2
(logSobolev) Any E P2(R') such that G(p), Ds(p) < satisfies
D.(p) 2 2A, G.(pp))).
(34)
(Talagrand) For any p E P2(RR) such that Gv(p) < 0o, we have
Wa(p.p.) Golpl p-).
(35)
29

---

Proef. The proof for this result follows analogously to the arguments presented in the proofs of Proposition 13
Corollary 14d Corollary 13
using the preliminary results established Proposition 21
and Proposition 22
Proof of Theorem 3
Following the same approach as in the proof of Theorem 2
the results in Theorem
immediately follow by combining Proposition 22 Corollary 23 and Proposition 24
applied to solutions of the
PDE
D Proof of Theorem
The proof for this theorem uses similar strategics that of Theorem but considers the cvolution of an ODE
rather than a PDE. Recall that for any x E IR" the best response (-))(-) E P(R") in (6 is defined as
r(x) = argmax Ge(,, x),
PEP
where the energy Ge(p,)): P(RR) x Rd + 1-00, oo) is given by
- f1(2, Gc(p.z) = (
1
f2(2, + || - voll" aKK((KK))))) - pW
hold and assume ad > A1. Then for each a E IR4 there exists a unique maximizer == r(x) solving argmaxper; Gelā,z), x). Further, r(x) E L (B"), supp(r(z)) supp(p), and there exists function c : Rd + I such that the best response pa(2) = rz))(2) solves the Fuler-Lagrange equation
Lemma 25. Let Assumptions
Gel, =](2 log (al (()))aaglooo())))))))
c(x) for all (2, E supp(p) x xR=
(36)
Proof. Equivalently, consider the minimization problem for F(p) = - (1(2, + aKL(p) I pW p with some fixed Note that we can rewrite F(p) as
F(P) = p plogpdz + V(2, + 1 / pW x
where V(2, x) == -(i x) + a log f(z)) is strictly convex in for fixed x by Assumptions and Together with Assumption 3 we can directly apply the uniqueness and existence result from [CMV03 Theorem 2.1 (i)).
The result on the support of r(z) and the expression for the Euler-Lagrange equation follows by exactly the same arguments as in Corollary and Corollary 23
Lemma 26. The density of the best response r(x) is continuous on Ie for any fixed E Rdd
Proof. Instead of solving the Euler-Lagrange equation 36), we can also obtain best response r(x) as long-time asymptotics for the following gradient flow:
0.p = div (pVo,
(37)
Following Definitions 8 and we characterize the steady states pso of the PDE 37) by requiring that
E L (R") n Lice (R") with ll ll = such that Poo E w.Le VW Poo LRo((") Poo is absolutely continuous with respect and pro satisfies
V. ( (f)(2z aog
(oc p(z)
+ W +
Vz E Edd
(38)
the sense of distributions. Noting that because energy functional F(p) differs from Ga μ1) only in
sign of fi f1( x)) if viewing G.(p. as function of only. Note that F()) is still uniformly displacement
convex p due to Assumption 4 Then the argument to obtain that Poo C(R)) follows exactly that of
Lemma
Lemma
Let E (L...., d}. If the energy H. P(R)) Rd given by
p(2) Hdlp,)) = 2 r(x)(2) dz 1 2 pW p- f Oz,
(39)
admits a critical point x E RR,, then the best response r(x) P(R))) is differentiable the ith coordinate direction E 1dd Further, the critical point of H; is in the subdifferential rerr(x).
30

---

Proef. First, note that (=(r())((=)(u) = 0 for all directions E C2 (IR*) and for all E R thanks to optimality of rax) Here, DF denotes the Fréchet derivative F, associating to every p E P(RR) the bounded linear operator DF(p) C + R
DF(p)(u)== := 6,
and we note that F(p) depends through the potential V. Fixing an index E (1, d], and diffcrentiating the optimality condition with respect to x, we obtain
O2 DF(r())((<)))))) + D2 ((r())](<))v, 02, =0 Vu C. (R*).
(40)
Both terms can he made more explicit using the expressions for the Fréchet derivative of F:
Oz, (F))(=)((=)())) - (12,)))(()))
and for the second term note that the second Fréchet derivative of F at E P(RR)) along directions u, C (R*) such that (supp(u) U supp(u)) C supp(p) given by
() = [u())(())
/f
W(2 - ())(
p(z)
In other words, assuming supp(r(z:) spp = R relation 40 can be written as
Oz a (ootte u(z) r())(z)
/ (W * o., rz)) (z)u(z) da - / 22,00 (z.
For case of notation, given r(x) P(R"), ) we define the function 9: P(R")) Lioc (R*) by
p(2)
( = (=) + W p 8r, (=,
The question whether the partial derivative 0z, r(x) exists then reduces to the question whether there exists some P. E P(RR)) such that / p. solves the equation
d (2)
for almost every E R
and for some constant O. This is preciscly the Euler-Lagrange condition for the functional H, defined in which has solution thanks the assumption of Lemma
39
We observe that first term in H; is precisely (up to constant) the ddiveergee with respect to r(z),
/( 0 ((x)) (-be)1) ((x)((
r(z)dz= dz (())))=
0² d r(x)
Depending on the shape the best response r(x), the - -divergencce may displacement convex. Similarly, the last term - J 0x f2, (2,22))d) x))p(z) ) in the energy H: is in fact displacement concave due to the convexity properties of fi in 2. The interaction term is displacement convex thanks to Assumption 3 As a result, the overall convexity properties of H: are not known in general. Proving the existence of critical for H1, under functional
assumptions f1, f2, L and W would be interesting result in own right, providing inequality that expands on the literature of related functional inequalities such the related Hardy-Littlewood Sobolev inequality [Lie83].
It remains to show that H, indeed admits a critical point. Next, we provide examples of additional assumptions that would guarantce for Lemma 27 to apply.
Lemma 28. If either C sup.Erd [W(x) 00, or
C= sup la (ogg)( + (2, =eed
<
then for each Rd and for large enough coordinate 0.,, r(x) given by the unique coordinate-wir solutions of the Euler- Lagrann condition for H..
the best response r(z) is differentiable with the gradient
Proof. We will show this result using the Banach Fixed Point Theorem for the mapping TІ T. (RR.) - (R.))
for each fixed E (L,.., d) given by
r(=)(z) IW p)(=) (2, +c.
(p)p)
31

---

noting that T(pp)) is the Euler- BagranLagrange condition for a critical point of Hs. It remains to show that T; is contractive mapping. For the first assumption, note that
- (())lW (0 (0-p')ldz
[T((p) -T(d')l
_ /f r(x)(2)W( - 2)lo()) (')))dādd à
| (fro)ca))) (ʃ (o(3) (8)aaE s a 1|p- 1'11
C
Similarly, for the second assumption we estimate
[T((p) T((D')l1 = _ r(o)lW (0- p')
// r(x)()) 2)1p(2) - (
< a
(W)())))))) d'(=)ld=
x | 1p -P'l1
which requires bound |W r(=)) Using
|W r( sup la log(r())) + (z,2z)
sERd
we conclude that T; is contraction map for large enough In both cases, we can then apply the Banach Fixed-Point Theorem to conclude that (xr(x) exists and unique.
Lemma 29. Let r(x) as defined in 6 (V.G.lp,)) p=r(e)
If r(x) is differentiable in then we have VzGa(x)
Proof. We start by computing V.Ga(x). have
V.Ga(r) = v (G.(r(x).x))
(V=Gdpp
=c(xVV, / ((s)(addz+ (V..G.(p.))) (V. (=r(())
where we used that r(x) solves the Euler Lagrange equation 36 and that r(z) E P(R") for any x Rd so that f (z)) is independent of x.
Lemma 30. Let Assumption hold. Then Ga 12d U (+00) is strongly convex with constant Xa A1 + + 3
Proef. The energy Golp, .)) is strongly convex in æ due to our assumptions on fi. fz. and the regularizing term || - 20 This means that for any E P.
Gelp, Go(p, VeG.lpp))) - x0) + Ad Au 1| 22112
Selecting p = r(x)), ) we have
x') ) da 1 1
Ge(r(z'). Gelr(el)r)))))
Since
by definition of r(x). we obtain the required convexity condition:
Ga(z)= = +
x')
A/
Ad
2|
x 11
2
Proof of Theorem 4 For any reference measure Pe P, we have
Ga(x) 2 Gelpo, x) 2 LaKL(pu (0) /
3 -11 xall
poW * pu +
32

---

and therefore, Ga is coercive. Together with the strong convexity provided by Lemma 30 we obtain the existence of unique minimizer aho R Convergence in norm immediately follows also using Lemma 30 for solutions zx() (6. we have
( x I (Ga(z(tt) Ga(xx))) (a(t) x=0) S ->al[a(t) x=oll"
2 dt
similar result holds for convergence in entropy using the Dolyak-Lojssiewiizzcoom convexity incquality
[VCa( 2 >e(Ga(s)) Ca(aoo))
which is itself a direct consequence strong convexity provided Lemma 30 Then
dt (Ga(x(t)) Ga(zs))) V.GG(a(t)) =(t) =-IV=Ga(a())) 22aa (Ga(z(t)) Ga(x>>))
so result in Theorem 4 follows.
E Additional Simulation Results
We simulate number additional scenarios to illustrate extensions beyond the setting with provable guarantees and in the settings for which we have results but no numerical implementations in the main paper. First, we simulate the aligned objectives setting in one dimension, corresponding to 4 Then we consider two settings which not covered in theory: (1) the previously-fixed distribution ş also time varying, and (2) the algorithm not have access the full distributions p and o and instead samples from them to update. Lastly, we illustrate classifier with the population attributes in two dimensions, which requires a different volume implementation CCHI5 Section 2.2] than the one dimension version of the PDE due to flux in two dimensions.
E.1 Aligned Objectives
Here we show numerical simulation results for the aligned objectives case, where the population and distribution have the same cost function. In this setting, the dynamics are of the form
p.p- div (pV.dh L.V.&,Gall S,pGalo nl)
= div pV: (f ((ll,a))d(()) + log(e//) + W+ xx)))
dt -z (f ((le.m)dpp))) / (2(2,)))7(s)) 1| moll?
where fi and fa are as defined section 4.1 and W= 1 1) consensus kereel. Note that does not satisfy Assumption 3 but we still observe convergence in the simulation. This is expected: in other works such as (CMV03). the assumptions on W are relaxed and convergence results proven given sufficient convexity of other terms. The regularizer is set pu. which models penalty for the effort required of individuals alter their attributes. The coefficient weights are a 0.1 and with discretization parameters dz 0.1. dt = 0.01.
0.8
0.6
P
probability
0.4
4
-2
2
6
2
Figure 4: The dynamics include a consensus kernel, which draws neighbors in z-space closer together. We see that the population moves to make the classifier performer better, as the two distributions become more easily separable by the linear classifier.
33

---

In Figure we observe the strategic distribution separating itself from the stationary distribution, improving the performance of the classifier and also improving the performance of the population itself. The strategic distribution classifier appear to stationary by time = 40.
E.2 Multiple Dynamical Populations
We also want to understand the dynamics when both populations are strategic and respond to the classifier. In this example, we numerically simulate this and in future work we hope to prove additional results regarding convergence. This corresponds to modeling the previously-fixet distribution distribution be P2. We consider the case where is competitive with and is aligned with 2, with dynamics given by
time-dependent; let this
0,p = - div (pV, (f)(2, llog(ppp) W * p))
al T. = div (rV. (f2(2,)) +
Φ V (/ [)(p())))(s)+ / /2(2,x)d7(() +
x
di
We use W = 0 and f1, f2 as in section
4.1 and the same discretization parameters as in Scction
E.1
In Figure
1.0
probabiliyy
0.8
0.4
0.0
-2
Figure 5: The population p aims to classified with the r population, while the classifier moves to delineate between the two. We observe that adjusts to improve the performance of the classifier while P competes against it. The distributions are plotted time O, corresponding to ý and 7, and time t 20, corresponding to and
we observe that the population moves to the right, assisting the classifier in maintaining accurate scoring. In contrast, also moves the right, rendering the right tail to classified incorrectly, which is desirable for individuals the P population but not desirable for the classifier. While we leave analyzing long-term bchavior mathematically for future work, the distributions and classifier appear to converge by time t 20.
E.3 Sampled Gradients
In real-world applications of classifiers, the algorithm may not know the exact distribution of the population, relying sampling to estimate it. In this section explore the effects of the classifier updating based an approximated gradient, which computed hy sampling the underlying distributions p and We use the same parameters for the population dynamics as section 4.1 and for classifier we use the approximate gradient
V. L( ~
72 Σ (Vef( (1) + Vafa(R, + 3(x, 2o), 21 ~ Pts
2, ~ Pr.
First, we simulate the dynamics with the classifier and the strategic population updating at the same rate, using 0.05, B = 1, and the same consensus kemel as used previously, with the same discretization parameters as In Figure
observe visual difference between the results with 4 versus 40 samples, which suggests that not many samples are needed to estimate the gradient.
Next, we consider the setting where classifier is best-responding to strategic population.
Unlike the first setting, we observe in Figure 7 noticeable difference between the evolution of pt with n versus n 40 samples. This is not surprising because optimizing with a very poor estimate of the cost function at each time step would cause to vary wildly, and this method fails to take advantage of correct "average" behavior that gradient descent provides.
34

---

sampled gradient under equal rate, n
0.75
1-0.0
T-6.0
t- 10.0
t- 30.0
P
probab biitt
0.50
P
0.25
Po
Po
Pi
0.00
-2
2
2
sampled gradient under equal rate, 40
0.75
t-0.0
t-6.0
t-10.0
t- C-30.0
probabilty
0.50
0.25
Po
Po
0.00
2
Figure When the classifier updating at the same rate as the population, we do not see a significant change in the evolution of both species, suggesting that as long as the gradient estimate for the classifier is correct on average, the estimate itself does not need be particularly accurate.
sampled gradient equal rate,
0.75
t=0.0
t=6.0
t=10.0
t= C=30.0
probabiity
0.50
0.25
Po
0.00
sampled gradient under rate, 40
0.75
t=6.0
t-30.0
probablity
0.50
0.25
X
Po
Pu
0.00
-2
2 -2
2
2
Figure 7: When the classifier is best-responding to the population, we observe that using n 4 samples leads to different behavior for both the classifier and population, compared with a more accurate estimate using 40 samples.
E.4 Two-dimensional Distributions
practice, individuals may alter more that one their attributes in response to an algorithm, for example, both cancelling credit card and also reporting different income in an effort to change a credit score. We model this case with R and RR?, and simulate the results for the setting where the classifier and the population evolving at the same While this setting is covered in theory, it interpolates between the two timescale extremes.
We consider the following classifier:
(2,22) (
exp Txp
(41)
1 1 fe(2,2) = eeexxxx
with W Again, the reference distribution á corresponds the initial shape of the distribution, instituting a penalty for deviating from initial distribution. We use and / for t = with 0.005 and dr = dy = 0.2 for the discretization. In this case, the strategic population is competing with the classifier, with dynamics given by
penalty weights, run
O.p - div (pV: (fi(z, log(p/p)))
x V (/ f1(2,))dp(=)+ / fa(=,))d7(=))
xoll?)
dt
In Figure 8 we observe the strategic population increasing mass toward the region of higher probability of being labeled "I" while the true underlying label is zero, with the probability plotted at time 4. This illustrates
35

---

p
-5
-5
0.04
0.04
0
0.02
0.02
55
5
-5
21
Pt
probability
-5
5
0.04
0.4
0
0.02
-5
0.0
-5
0
5
-5
0
Z1
Z1
Figure 8: We (41) for the classifier functions, using Gaussian initial condition and regularizer for p. We see the distribution moving toward the region with higher probability of misclassification.
similar behavior to the mie-dimensional case, including the distribution splitting into two modes, which another example of polarization induced by the classifier. Note that while in this example, æ E R and we use linear classifier; we could have x E IRd' with d > 2 and different functions for fi and fa which yield nonlinear classifier; our theory in the timescale-separatedd case holds long the convexity and smoothness assumptions on fi and f2 are satisfied.
36


Reviews:
[{"review_id": "gzFXTrdoYr", "reviewer": "Reviewer_P1kC", "summary": "The paper proposes a new framework for analyzing the interaction and feedback loop between learning algorithms and the distributions on which they are deployed. The framework in the paper is more general than prior models, allowing to capture more complex interactions and shifts in the algorithm and distribution. This is done using a partial differential equation that modifies the algorithm\\distribution in a way that optimizes its cost with respect to the choice of distribution\\algorithm. Under several assumptions which include strong convexity, the authors were able to prove that gradient flow converges and provide fast convergence rates. The authors are provide experimental results for a simple 1-dimensional setting which demonstrate that their framework can explain general phenomenon which are not explained by prior models, such as polarization.", "questions": "See above.", "limitations": "See above.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper is well-written and easy to follow. Moreover, as mentioned above, the framework in the paper offers several advantages compared to prior frameworks, allowing to capture more complex interactions and shifts of the distribution. This model is therefore able (via simple examples) to theoretically explain forms of distribution shifts which we see in practice such as polarization. \n", "weaknesses": "\n1. While the framework allows for complex forms of shifts of the model and distribution, the analysis assumes that the model\\distribution changes to be the exact minimizer of a certain cost function. This may be a little strict in general. Can the same analysis work when the shifts are approximate minimizers for the cost function?\n\n2. Extreme timescale setting: the analysis in the paper only considers two extreme timescale settings: (i) the algorithm responds much faster than the population, and (ii) the population responds much faster than the algorithm. This is somewhat strict and doesn’t allow for intermediate changes in the algorithm\\distribution. \n\n3. Assumptions: the paper has many assumptions, some of which are less common. The authors should give some motivation and justification for their assumptions, and provide examples where all of these assumptions hold.\n\n4. The experiments are a little simple (1-dimensional classifiers) and it would be more interesting to see an example which is more relevant in practice. This is perhaps one of the limitations of this framework: even though it is able to capture more complex structures of shifts, it is harder to simulate for more sophisticated examples and models (e.g. high dimensional classifiers, other timescale settings, and so on). The authors should therefore also discuss the limitation of their framework compared to existing ones. \n\n5. Experimental discussion: the experimental section jumps too fast to the experiments without explaining the objectives of these experiments or discussing its connections to the previous theoretical results. The authors should explain their main objectives in this setting, discuss connection to the theory in the previous section, and also compare to existing prior work (e.g. show how the example in section 4.1 cannot be captured by prior frameworks).\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your detailed response and helpful comments. \n\n**Model Optimization:** The setting where the dynamics are not exactly minimizing a cost functional is an interesting idea for future work; our work can be seen as a necessary first step for building more general frameworks in this direction. While we considered the case where the distribution evolves towards the exact minimizer of a given cost function, in practice, this cost function would be approximated from data resulting in a shift that is not exactly solving the true minimization problem. We also note that for the gradient flow formulation, the minimizer is reached as time goes to infinity. In practice, the evolution would be solved up to finite time only, therefore providing an approximation of the true minimizer. By choosing this final time large enough, arbitrary accuracy can be achieved.  The convergence rates we obtain provide quantitative error bounds both in energy and in Wasserstein distance that measure how close the shift obtained at any chosen final time is to the true minimizer.\n\n**Extreme timescale setting:** While we provide numerical examples of the intermediate timescale separation case, we leave rigorous analysis for future work. Different analytical tools are likely needed for a rigorous treatment of this case.\n\n**Assumptions:** We will add more details on motivation and justification of these assumptions, extending the top of page 6 in our paper where we have examples of functions which satisfy the assumptions, including the log-loss function which is commonly used in machine learning problems.\n\n**Regarding Limitations:** We will add a limitations section highlighting scaling with dimension as that is a key challenge for our framework. Besides using SDEs and deep learning to solve these systems of PDEs in higher dimensions, we would like to investigate settings where agents are resource bound and can only perturb one or two features in a larger set of features, since in reality humans have limited resources for feature manipulation.\n\nIn comparison with prior work that addresses learning in the face of distribution shift, consider the methods used in [MPZ21], wherein the learning algorithm fits a linear mean-shift model to the algorithm parameters. Our framework provides a method for comparing the outcome of using a mean-shift model versus continually retraining the algorithm via gradient descent, which is shown in Figure 3. This simulation shows that state-of-the-art techniques for responding to distribution shift can perform poorly compared to naive techniques, emphasizing the need for better modeling. The framework itself is novel as PDEs have not been used to model distribution shift in the context of people interacting with algorithms, and takes a similar form as PDEs which model populations of animals and bacteria.\n\nTo illustrate the relevance of this new framework, we show how our model can be used to fit real-world data. In Figure 1 [see the attached PDF], poverty index scores collected in Colombia [CC11] are shown before and after the threshold for qualifying for government aid (shown in red) was released. We observe that many people shifted their data to be just under the threshold. We can model this with our proposed PDE framework, as shown in Figure 2. We will add both of these examples to our paper.\n\n**Experimental Discussion:** We will provide more context and explain our experiments more clearly; thank you for the suggestion. Other than the 2d classifier we simulated, we leave implementations in higher dimensions for future work. Using stochastic differential equations (SDEs) would be a reasonable strategy for simulations in higher dimensions. Additionally, a number of recent works such as [Li+22,SS18] show promise in using deep learning to solve PDEs in high dimensions.\n\n[CC11] Adriana Camacho and Emily Conover. “Manipulation of Social Program Eligibility”. In: American Economic\nJournal: Economic Policy 3.2 (2011), pp. 41–65. issn: 1945-7731.\n\n[SS18] Justin Sirignano and Konstantinos Spiliopoulos. “DGM: A deep learning algorithm for solving partial differential\nequations”. In: Journal of Computational Physics 375 (Dec. 2018), pp. 1339–1364.\n\n[MPZ21] John P. Miller, Juan C. Perdomo, and Tijana Zrnic. “Outside the Echo Chamber: Optimizing the Performa-\ntive Risk”. In: Proceedings of the 38th International Conference on Machine Learning. International Conference\non Machine Learning. \n\n[Li+22] Zongyi Li et al. Fourier Neural Operator with Learned Deformations for PDEs on General Geometries.\n"}, {"review_id": "JpT7OnT8pJ", "reviewer": "Reviewer_HH7Z", "summary": "This paper presents a framework for analyzing dynamics of distribution shift that captures feedback loop between learning algorithms and the distributions that they are deployed on using a coupled differential equations model. The paper considers both co-operative and competitive settings and suggests that a learner updated with gradient descent converges to steady state both in finite and infinite dimensions (in terms of model parameters). ", "questions": "One way to make the paper appeal to a general audience is to have a discussion surrounding claims and assumptions comparing this work against prior results in the literature.", "limitations": "n/a", "rating": 5, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "The paper presents a rigorous treatment of the problem and presents general results, although I am not an expert in this area and cannot present a rating with any meaningful confidence.", "weaknesses": "The paper works with several assumptions that I do not see a clear cut comparison of how these assumption compares against prior work, so it is unclear how general these results are.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your review and feedback.\n\n**Questions:** In comparison with prior work that addresses learning in the face of distribution shift, we can analyze with the methods used in [MPZ21], wherein the learning algorithm fits a linear mean-shift model to the algorithm parameters. Our framework provides a method for comparing the outcome of using a mean-shift model versus continually retraining the algorithm via gradient descent, which is shown in Figure 3. This simulation shows that state-of-the-art techniques for responding to distribution shift can perform poorly compared to naive techniques, emphasizing the need for better modeling. The framework itself is novel as PDEs have not been used to model distribution shift in the context of people interacting with algorithms, and takes a similar form as PDEs which model populations of animals and bacteria. We will add this to the paper to provide more context for our framework.\n\nWe will add more details on motivation and justification of our assumptions, extending the top of page 6 in our paper where we have examples of functions which satisfy the assumptions, including the log-loss function which is commonly used in machine learning problems. More broadly, the assumptions on the input functions are commonly used and well-studied for related classes of PDE models. There are several ways how they could be relaxed, see for example the remarks in [CMV03]; we did not opt for the most general presentation of assumptions in this work as it would increase the technicality of the proofs without much conceptual gain for the presented frameworks, given that the current assumptions already cover large classes functions relevant in practice, including the prototypical example considered in Section 4.\n\nTo illustrate the relevance of this new framework, we show how our model can be used to fit real-world data. In Figure 1 [see the attached PDF], poverty index scores collected in Colombia [CC11] are shown before and after the threshold for qualifying for government aid (shown in red) was released. We observe that many people shifted their data to be just under the threshold. We can model this with our PDE, as shown in Figure 2. We will also add this to the paper.\n\nIs there any other information we could provide to help improve our paper?\n\n[CMV03] Jose A. Carrillo, Robert J. McCann, and Cedric Villani. “Kinetic equilibration rates for granular media and related equations: entropy dissipation and mass transportation estimates”. In: Revista Matematica Iberoamericana 19.3 (Dec. 2003), pp. 971–1018. issn: 0213-2230.\n\n[MPZ21] John P. Miller, Juan C. Perdomo, and Tijana Zrnic. “Outside the Echo Chamber: Optimizing the Performative Risk”. In: Proceedings of the 38th International Conference on Machine Learning. International Conference\non Machine Learning. \n\n[CC11] Adriana Camacho and Emily Conover. “Manipulation of Social Program Eligibility”. In: American Economic\nJournal: Economic Policy 3.2 (2011), pp. 41–65. issn: 1945-7731."}, {"review_id": "UkKryAMt9E", "reviewer": "Reviewer_vBit", "summary": "This paper studies the long-term behavior of a dynamical system where there are distribution shifts in response to algorithmic decision-making. The authors derive PDEs to capture distributional changes and describe the result of the algorithm for both the coorperative setting and the competitive setting. Under certain regularization conditions, it is shown that the PDEs converge to unique steady-state distribution. Finally, some numerical simulations are conducted to iilustrate the problem.", "questions": "1. Many real-world applications involve discrete rather than continuous time. Is it possible to generalize the analysis in this paper to discrete time setting?\n\n2. Is it reasonable to assume that $f_1,f_2$ are convex?", "limitations": "The authors have adequately addressed the limitations and potential negative societal impact of their work.", "rating": 7, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "1. The paper is well-written and not hard to read. All definitions, theoretical results and experiments are accompanied by explanations when necessary. Moreover, the authors motivate the problem of interest in a very convincing way.\n\n2. The paper is technically solid. The process of interaction between algorithm and population is rigorously modelled and its asymptotic behavior is characterized.", "weaknesses": "The authors do not include many discussions on the insights provided by the theoretical findings. This is probably because the limit distribution of the PDEs are difficult to characterize. Moreover, the setting considered is over-simplified since it only contains one agent. Given that this paper is kind of a first step towards mathematical analysis of algorithm interaction, I think these flaws are totally acceptable and can be left for future work.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your helpful feedback and questions.\n\n**Weaknesses:** We agree that more of a discussion about insights would improve our paper and we can add this in the next version. One key insight is the identification of a large classes of input functions that lead to exponential convergence with explicit rates that can be computed directly from model parameters, telling us how fast we can expect the dynamics to converge. Another key insight is the potential of this PDE formulation to allow for methodology transfer from other application fields such as math-biology and physics in terms of the analysis tools used for predicting the behavior of solutions. \nWe are excited about extending this framework beyond one type of agent (note that currently we have an infinite number agents, same type) in our next iteration, where we would have different \"species\". For instance, in our current setting $\\bar \\rho$ is a fixed reference measure; we expect our framework to be able to capture multi-species settings as well, such as the case when $\\bar \\rho$ also changes dynamically. We see no methodological limitations to extend our results to such multi-species systems in the future.\n\n**Question 1:** We could also consider these dynamics in discrete time and that is a good future direction. The goal of this paper is to analyze the precise distribution shift from a continuous optimization perspective. Any results achieved for the continuum dynamics directly translate into approximate statements on the level of the time discretized analogues.\n\n**Question 2:** We assumed that $f_1$ and $f_2$ are convex for the purposes of the analysis, which covers many reasonable scenarios such as the log-loss, a commonly-used ML loss function. There are several ways how the assumptions could be relaxed, see for example the remarks in [CMV03]; we did not opt for the most general presentation of assumptions in this work as it would increase the technicality of the proofs without much conceptual gain for the presented frameworks, given that the current assumptions already cover large classes functions relevant in practice, including the prototypical example considered in Section 4. For instance, in [CMV03] the convexity of the potential term (corresponding to $f_1$ and $f_2$ in our case) can be relaxed and $\\tilde\\rho$ and/or $W$ can compensate to ensure convergence; similar generalizations could be applied in our setting.\n\n**Note on additional simulations:** In comparison with prior work that addresses learning in the face of distribution shift, we consider the methods used in [MPZ21], wherein the learning algorithm fits a linear mean-shift model to the algorithm parameters. Our framework provides a method for comparing the outcome of using a mean-shift model versus continually retraining the algorithm via gradient descent, which is shown in Figure 3 [see attached PDF]. This simulation shows that state-of-the-art techniques for responding to distribution shift can perform poorly compared to naive techniques, emphasizing the need for better modeling. The framework itself is novel as PDEs have not been used to model distribution shift in the context of people interacting with algorithms, and takes a similar form as PDEs which model populations of animals and bacteria.\n\nTo illustrate the relevance of this new framework, we show how our model can be used to fit real-world data. In Figure 1, poverty index scores collected in Colombia [CC11] are shown before and after the threshold for qualifying for government aid (shown in red) was released. We observe that many people shifted their data to be just under the threshold. We can model this with our proposed PDE framework, as shown in Figure 2.\n\n[CMV03] Jose A. Carrillo, Robert J. McCann, and Cedric Villani. “Kinetic equilibration rates for granular media and related\nequations: entropy dissipation and mass transportation estimates”. In: Revista Matematica Iberoamericana 19.3\n(Dec. 2003), pp. 971–1018. issn: 0213-2230.\n\n[CC11] Adriana Camacho and Emily Conover. “Manipulation of Social Program Eligibility”. In: American Economic\nJournal: Economic Policy 3.2 (2011), pp. 41–65. issn: 1945-7731.\n\n[MPZ21] John P. Miller, Juan C. Perdomo, and Tijana Zrnic. “Outside the Echo Chamber: Optimizing the Performative Risk”. In: Proceedings of the 38th International Conference on Machine Learning. International Conference\non Machine Learning. "}, {"review_id": "c1TlfTqSjx", "reviewer": "Reviewer_nK6f", "summary": "The authors consider a the problem where the change in a distribution for an objective can be modeled as a set of coupled nonlinear parabolic PDEs. These methods have nonlocal interactions and describe a model of the influence of the model on the population and vice-versa. Additionally, these correspond to the evolution of the measures associated with the Wasserstein gradient flow that minimizes a set of defined energy functionals, which relate to the original optimization problem. Such equations have been well studied in terms of the granular media equation and a wide theory has been devoted to existence and uniqueness of their solution. Two scenarios are considered: one with cooperation and another with an adversarial behavior of the population. The main results of the paper are convergence rates to the steady state for the coupled PDEs. The behavior of both the optimization algorithm and the population are studied in some numerical experiments at the end of the paper. The numerical experiments illustrate the importance of considering such a model rather than relying on more simple summary statistics for notions of distributional shift. ", "questions": "The second equation should be the argmin of x \\in R^d?\n\nIs there intuition on how the different functionals could be estimated e.g. for a particular application? Or is it assumed that these are already known? \n\nRelated to the previous question, how difficult would it be to apply to a real scenario, beyond the simulations the authors provided? This may not be a focus of the paper, but seems like something worth discussing. \n\nCould particle methods be used to solve the problem in high dimensions as is sometimes done in high dimensional cases? For example, this was explored in [1] and I was wondering if it could be applied in this scenario by approximating the associated McKean-Vlasov process (there might be an issue with some of the terms in the PDEs)? \n\n[1] Crucinio et al, Solving Fredholm Integral Equations of the First Kind via Wasserstein Gradient Flows, 2022", "limitations": "The authors sufficiently discuss limitations of their work.\n\n\n", "rating": 7, "confidence": 3, "soundness": 4, "presentation": 3, "contribution": 4, "strengths": "The proposed model is a very elegant model for describing distribution shift and the feedback mechanisms between the shift in population and objective functions. It provides a much richer class of perturbations than what is generally seen in the literature in cases such as distributionally robust optimization or in adversarial optimization. The model also has implications in providing additional distribution information (beyond low order moments) at equilibrium. The authors show the importance of this where different geometries of the distribution can imply different effects on the population distribution. Overall, the paper is well written, provides nice ideas, and describes a more unique take on the problem of distribution shifts. ", "weaknesses": "The numerical evaluation is a bit limited, but I think that’s not a problem since most of the results are theoretical. Still, it would be nice to see some more concrete applications of the theory by simulating the PDEs described, where, for example, all components are not simulated. It also appears that this may be difficult to apply to real world scenarios as the authors alluded in their limitations. A real application where all components need to be estimated was not discussed, but that's beyond the scope of the paper. \n\nSome of the assumptions are somewhat strong, but these are usually made to provide convergence and existence of solutions of these PDEs. In that regard, processes that satisfy these assumptions may not be entirely general, but that’s not a big deal since the authors clearly constrain their analysis based on the assumptions they make. \n\n\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your detailed, specific, and very helpful feedback on our work.\n\n**Question 1:** Thank you for catching this; we had a typo. The equation should be $$\\underset{z\\in\\mathbb{R}^d}{\\operatorname{argmin}}J(z,x)$$ because this is modeling an individual agent, whose data is given by $z$, optimizing their own utility with respect to given algorithm parameters $x\\in\\mathbb{R}^d$, whereas the first equation represents the algorithm optimizing a loss over the entire population of agents. \n\n**Question 2:** This is a great idea; we are really interested and excited about investigating this in our future work. We have been considering  kernel methods and Gaussian processes to accomplish this task in addition to machine learning approaches. \n\n**Question 3:** To illustrate the relevance of this new framework, we show how our model can be used to fit real-world data. In Figure 1 [see attached PDF], poverty index scores collected in Colombia [CC11] are shown before and after the threshold for qualifying for government aid (shown in red) was released. We observe that many people shifted their data to be just under the threshold. We can model this with our PDE, as shown in Figure 2, using basic intuition to fit the PDE. We will add this to our paper.\n\n**Question 4:** Using stochastic differential equations (SDEs) as in the paper that was referenced would be a good way to simulate our PDEs in higher dimensions. Additionally, a number of recent works such as [Li+22,SS18] show promise in using deep learning to solve PDEs in high dimensions. As a subcase of the high-dimensional setting, we would like to investigate settings where agents are resource bound and can only perturb one or two features in a larger set of features, since in reality humans have limited resources for feature manipulation.\n\n**Note on additional figure:**\nIn comparison with prior work that addresses learning in the face of distribution shift, we consider the methods used in [MPZ21], wherein the learning algorithm fits a linear mean-shift model to the algorithm parameters. Our framework provides a method for comparing the outcome of using a mean-shift model versus continually retraining the algorithm via gradient descent, which is shown in Figure 3. This simulation shows that state-of-the-art techniques for responding to distribution shift can perform poorly compared to naive techniques, emphasizing the need for better modeling. The framework itself is novel as PDEs have not been used to model distribution shift in the context of people interacting with algorithms, and takes a similar form as PDEs which model populations of animals and bacteria.\n\n[CC11] Adriana Camacho and Emily Conover. “Manipulation of Social Program Eligibility”. In: American Economic\nJournal: Economic Policy 3.2 (2011), pp. 41–65. issn: 1945-7731.\n\n[SS18] Justin Sirignano and Konstantinos Spiliopoulos. “DGM: A deep learning algorithm for solving partial differential\nequations”. In: Journal of Computational Physics 375 (Dec. 2018), pp. 1339–1364.\n\n[MPZ21] John P. Miller, Juan C. Perdomo, and Tijana Zrnic. “Outside the Echo Chamber: Optimizing the Performa-\ntive Risk”. In: Proceedings of the 38th International Conference on Machine Learning. International Conference\non Machine Learning. \n\n[Li+22] Zongyi Li et al. Fourier Neural Operator with Learned Deformations for PDEs on General Geometries."}]

Author rebuttal:
See figures in attachment.

Meta-review: The paper proposes a novel framework for analyzing distribution shift and feedback mechanisms between learning algorithms and the distributions they operate on using coupled nonlinear parabolic partial differential equations (PDEs). The reviews generally found the paper technically sound and appreciated its contribution. The reviewers noted the elegance of the model and its potential to capture complex interactions in distribution shifts as well as a richer class of perturbations compared to existing literature.

However, there are some concerns and suggestions raised by the reviewers:

1. Reviewers suggest that the paper could benefit from a more extensive discussion of insights derived from the theoretical findings. They recommend adding more context about the practical implications and applications of the model.
2. While the paper is primarily theoretical, some reviewers mention that additional numerical experiments could enhance its relevance. This could include simulating the PDEs described in various scenarios or providing more concrete examples of how the framework can be applied.
3. Some reviewers note that the assumptions made in the paper are strong, but they acknowledge the need for such assumptions to ensure convergence and existence of solutions for the PDEs. The authors should further clarify and justify these assumptions.
4. Reviewers express interest in understanding how challenging it would be to apply the framework in real-world scenarios. While the limitations are discussed, providing more insights into the practical challenges and potential applications could be valuable.

Overall, the paper was well-received, and the reviewers commended the authors for their detailed responses to the feedback. Addressing the suggestions and adding more context and insights to the paper could further enhance its impact and practical relevance. Therefore, I recommend accepting the paper with revisions based on the feedback provided.
Decision: Accept (poster)