# Strategic Distribution Shift of Interacting Agents
via Coupled Gradient Flows

Lauren Conger
California Institute of Technology
<EMAIL_ADDRESS>
&Franca Hoffmann
California Institute of Technology
<EMAIL_ADDRESS>
&Eric Mazumdar
California Institute of Technology
<EMAIL_ADDRESS>
&Lillian Ratliff
University of Washington
<EMAIL_ADDRESS>
###### 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.

##  1 Introduction

In many machine learning tasks, there are commonly sources of exogenous and
endogenous distribution 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 a decision-making algorithm can incentivize
individuals to change or mis-report their data to achieve desired outcomes— a
phenomenon known in economics as Goodhart’s law. Such phenomena have been
empirically observed, a well-known example being in [CC11], where researchers
observed a population in Columbia strategically mis-reporting data to game a
poverty index score used for distributing government assistance. Works such as
[Mil+20, Wil+21], which investigate the effects of distribution shift over
time on a machine learning algorithm, point toward the need for evaluating the
robustness of algorithms to distribution shifts.

Many existing approaches for modeling distribution shift focus on simple
metrics like optimizing over moments or covariates [DY10, LHL21, BBS09]. Other
methods consider worst-case scenarios, as in distributionally robust
optimization [AZ22, LFG22, DN21, Kuh+19]. However, when humans respond to
algorithms, these techniques may not be sufficient to holistically capture the
impact an algorithm has on a population. For example, an algorithm that takes
into account shifts in a distribution’s mean might inadvertently drive
polarization, rendering a portion of the population disadvantaged.

Motivated by the need for a more descriptive model, we present an alternative
perspective which allows us to fully capture complex dynamics that might drive
distribution shifts in real-world systems. Our approach is 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 [Per+20, IYZ21, Ray+22, Nar+22,
MPZ21], the strategic interactions studied in strategic classification
[Har+16, Don+18], and also endogenous factors like intra-population 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 [Zrn+21], the authors investigate the
effects of dynamics in strategic classification problems— but the model they
analyze does not capture individual interactions in the population. In
[IYZ21], the authors model the interaction between a population that
repeatedly responds to algorithmic decision-making by shifting its mean.
Additionally, [Ray+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 on average performance. Each of these
works fails to account for diffusion or intra-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-population 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 to 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 a relatively small body of prior
work, none of which related to distribution shifts in societal systems, and a
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 Law. A number of works consider coupling via
convolution kernels [FF13, Giu+22, JPZ22, CHS18, 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-based 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 a 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 a 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-Fisher-Rao metric. Each of these
references considers an energy functional that is linear in the distribution
of each species respectively; our energy includes nonlinearities in the
distributions via a self-interaction term as well as diffusion for the
population. Moreover, the above works introduce a 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 Problem Formulation

Machine learning algorithms that are deployed into the real world for
decision-making often become part of complex feedback loops with the data
distributions and data sources with which they interact. In an effort to model
these interactions, consider a machine learning algorithm that has loss given
by $L(z,x)$ where $x\in\mathbb{R}^{d}$ are the algorithm parameters and
$z\in\mathbb{R}^{d}$ are the population attributes, and the goal is to solve

$\displaystyle\operatorname*{argmin}_{x\in\mathcal{X}}\operatorname*{\mathbb{E}}_{z\sim\rho}L(z,x),$

where $\mathcal{X}$ is the class of model parameters and $\rho(z)$ is the
population distribution. Individuals have an objective given by $J(z,x)$ in
response to a model parameterized by $x$, and they seek to solve

$\displaystyle\operatorname*{argmin}_{z\in\mathbb{R}^{d}}J(z,x).$

When individuals in the population and the algorithm have access to gradients,
we model the optimization process as a gradient-descent-type 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 in 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 of a
collective population. To analyse how the population evolves over time, a
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:\mathcal{P}(\mathbb{R}^{d})\mapsto\mathbb{R}$
and fixed probability distribution $\rho\in\mathcal{P}(\mathbb{R}^{d})$, the
_first variation of $G$ at the point $\rho$_ is denoted by
$\delta_{\rho}G[\rho]:\mathbb{R}^{d}\to\mathbb{R}$, and is defined via the
relation

$\displaystyle\int\delta_{\rho}G[\rho](z)\psi(z)\mathrm{d}z=\lim_{\epsilon\rightarrow
0}\frac{1}{\epsilon}(G(\rho+\epsilon\psi)-G(\rho))$

for all $\psi\in C_{c}^{\infty}(\mathbb{R}^{d})$ such that
$\int\mathrm{d}\psi=0$, assuming that $G$ is regular enough for all quantities
to exist.

Here, $\mathcal{P}(\mathbb{R}^{d})$ denotes the space of probability measures
on the Borel sigma algebra. Using the first variation, we can express the
gradient in Wasserstein-2 space, see for example [Vil03, Exercise 8.8].

######  Lemma 1.

The gradient of an energy $G:\mathcal{P}_{2}(\mathbb{R}^{d})\to\mathbb{R}$ in
the Wasserstein-2 space is given by

$\nabla_{W_{2}}G(\rho)=-\mathrm{div}\left(\rho\nabla\delta_{\rho}G[\rho]\right)\,.$

Here, $\mathcal{P}_{2}(\mathbb{R}^{d})$ denotes the set of probability
measures with bounded second moments, also see Appendix A.2. As a consequence,
the infinite dimensional steepest descent in Wasserstein-2 space can be
expressed as the PDE

$\displaystyle\partial_{t}\rho=-\nabla_{W_{2}}G(\rho)=\mathrm{div}\left(\rho\nabla\delta_{\rho}G[\rho]\right)\,.$
(1)

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 a gradient flow in the Wasserstein-2
metric on the level of the distribution of the population. Then for given
algorithm parameters $x\in\mathbb{R}^{d}$, the evolution for this strategic
population is given by

$\displaystyle\partial_{t}\rho=\mathrm{div}\left(\rho\nabla\delta_{\rho}\Big{[}\operatorname*{\mathbb{E}}_{z\sim\rho}J(z,x)+E(\rho)\Big{]}\right),$
(2)

where $E(\rho)$ is a functional including terms for internal influences and
external perturbations. In real-world deployment of algorithms, 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 as

$\displaystyle\dot{x}$
$\displaystyle=-\nabla_{x}\big{[}\operatorname*{\mathbb{E}}_{z\sim\rho}L(z,x)\big{]}.$
(3)

In this work, we analyze the behavior of the dynamics under the following
model. The algorithm suffers a cost $f_{1}(z,x)$ for a data point $z$ under
model parameters $x$ in the strategic population, and a cost $f_{2}(z,x)$ for
a data point in a fixed, non-strategic population. The strategic population is
denoted by $\rho\in\mathcal{P}$, and the non-strategic population by
$\bar{\rho}\in\mathcal{P}$. The algorithm aims to minimize

$\displaystyle\operatorname*{\mathbb{E}}_{z\sim\rho}L(z,x)=\int
f_{1}(z,x)\mathrm{d}\rho(z)+\int
f_{2}(z,x)\mathrm{d}\bar{\rho}(z)+\frac{\beta}{2}\left\lVert
x-x_{0}\right\rVert^{2}\,,$

where the norm is the vector inner product $\left\lVert
x\right\rVert^{2}=\langle x,x\rangle$ and $\beta>0$ weights the cost of
updating the model parameters from its initial condition.

We consider two settings: $(i)$ 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 a reference distribution $\tilde{\rho}\in\mathcal{P}$. This
cost stems from effort required to manipulate features, such as a 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.

###  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 to
develop accurate recommendations for the users. This corresponds to the
population cost

$\displaystyle\operatorname*{\mathbb{E}}_{z\sim\rho,x\sim\mu}J(z,x)=\iint
f_{1}(z,x)\mathrm{d}\rho(z)\mathrm{d}\mu(x)+\alpha KL(\rho\,|\,\tilde{\rho}),$

where $KL(\cdot\,|\,\cdot)$ denotes the Kullback-Leibler divergence. Note that
the KL divergence introduces diffusion to the dynamics for $\rho$. The weight
$\alpha>0$ parameterizes the cost of distribution shift to the population. To
account for nonlocal information and influence among members of the
population, we include a kernel term $E(\rho)=\frac{1}{2}\int\rho
W\ast\rho\,\mathrm{d}z$, where $(W\ast\rho)(z)=\int
W(z-\bar{z})\mathrm{d}\rho(\bar{z})$ is a convolution integral and $W$ is a
suitable interaction potential.

###  2.2 Case (ii): Competing Objectives

In settings such as online internet forums, where algorithms and users have
used manipulative strategies for marketing [Del06], the strategic population
may be incentivized to modify or mis-report their attributes. The algorithm
has a competitive objective, in that it aims to maintain performance against a
population whose dynamics cause the algorithm performance to suffer. When the
strategic population seeks an outcome contrary to the algorithm, we model
strategic population cost as

$\displaystyle\operatorname*{\mathbb{E}}_{z\sim\rho,x\sim\mu}J(z,x)=-\iint
f_{1}(z,x)\mathrm{d}\rho(z)\mathrm{d}\mu(x)+\alpha KL(\rho\,|\,\tilde{\rho}).$

A 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 our 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 case in a simulation example.

##  3 Results

We are interested in characterizing the long-time asymptotic behavior of the
population distribution, as it depends on the decision-makers action over
time. The structure of the population distribution gives us insights about how
the decision-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 a novel analytical framework as well as
analysis of the long-time asymptotics.

A key observation is that the dynamics in (2) and (3) can be re-formulated as
a gradient flow; we lift $x$ to a probability distribution $\mu$ by
representing it as a Dirac delta $\mu$ sitting at the point $x$. As a result,
the evolution of $\mu$ will be governed by a PDE, and combined with the PDE
for the population, we obtain a system of coupled PDEs,

$\displaystyle\partial_{t}\rho$
$\displaystyle=\mathrm{div}\left(\rho\nabla_{z}\delta_{\rho}\big{[}\operatorname*{\mathbb{E}}_{z\sim\rho,x\sim\mu}J(z,x)+E(\rho)\big{]}\right)$
$\displaystyle\partial_{t}\mu$
$\displaystyle=\mathrm{div}\left(\mu\nabla_{x}\delta_{\mu}\big{[}\operatorname*{\mathbb{E}}_{z\sim\rho,x\sim\mu}L(z,x)\big{]}\right),$

where $\delta_{\rho}$ and $\delta_{\mu}$ are first variations with respect to
$\rho$ and $\mu$ according to Definition 1. The natural candidates for the
asymptotic profiles of this coupled system are its steady states, which -
thanks to the gradient flow structure - can be 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 a 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(\gamma|\gamma_{\infty})=G(\gamma)-G(\gamma_{\infty})$, where
$G(\gamma_{\infty})$ is the energy at the steady state.

Since we consider the joint evolution of two probability distributions, we
define a distance metric $\overline{W}$ on the product space of probability
measures with bounded second moment.

######  Definition 3 (Joint Wasserstein Metric).

The metric over
$\mathcal{P}_{2}(\mathbb{R}^{d})\times\mathcal{P}_{2}(\mathbb{R}^{d})$ is
called $\overline{W}$ and is given by

$\displaystyle\overline{W}((\rho,\mu),(\tilde{\rho},\tilde{\mu}))^{2}=W_{2}(\rho,\tilde{\rho})^{2}+W_{2}(\mu,\tilde{\mu})^{2}$

for all pairs
$(\rho,\mu),(\tilde{\rho},\tilde{\mu})\in\mathcal{P}_{2}(\mathbb{R}^{d})\times\mathcal{P}_{2}(\mathbb{R}^{d})$,
and where $W_{2}$ denotes the Wasserstein-2 metric (see Definition 4). We
denote by
$\mathcal{\overline{W}}(\mathbb{R}^{d}):=(\mathcal{P}_{2}(\mathbb{R}^{d})\times\mathcal{P}_{2}(\mathbb{R}^{d}),\overline{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 dynamics as a gradient flow by using the same energy
functional for both species. Let
$G_{a}(\rho,\mu):\mathcal{P}(\mathbb{R}^{d})\times\mathcal{P}(\mathbb{R}^{d})\mapsto[0,\infty]$
be the energy functional given by

$\displaystyle G_{a}(\rho,\mu)$ $\displaystyle=\iint
f_{1}(z,x)\mathrm{d}\rho(z)\mathrm{d}\mu(x)+\iint
f_{2}(z,x)\mathrm{d}\bar{\rho}(z)\mathrm{d}\mu(x)+\alpha
KL(\rho|\tilde{\rho})+\frac{1}{2}\int\rho W\ast\rho$
$\displaystyle\quad+\frac{\beta}{2}\int\left\lVert
x-x_{0}\right\rVert^{2}\mathrm{d}\mu(x).$

This expression is well-defined as the relative entropy
$KL(\rho\,|\,\tilde{\rho})$ can be extended to the full set
$\mathcal{P}(\mathbb{R}^{d})$ by setting $G_{a}(\rho,\mu)=+\infty$ in case
$\rho$ is not absolutely continuous with respect to $\tilde{\rho}$.

In the _competitive_ case we define
$G_{c}(\rho,x):\mathcal{P}(\mathbb{R}^{d})\times\mathbb{R}^{d}\mapsto[-\infty,\infty]$
by

$\displaystyle G_{c}(\rho,x)$ $\displaystyle=\int
f_{1}(z,x)\mathrm{d}\rho(z)+\int
f_{2}(x,z^{\prime})\mathrm{d}\bar{\rho}(z^{\prime})-\alpha
KL(\rho|\tilde{\rho})-\frac{1}{2}\int\rho W\ast\rho+\frac{\beta}{2}\left\lVert
x-x_{0}\right\rVert^{2}.$

In 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

$\displaystyle\partial_{t}\rho$
$\displaystyle=\mathrm{div}\left(\rho\nabla_{z}\delta_{\rho}G_{a}[\rho,\mu]\right)$
(4) $\displaystyle\partial_{t}\mu$
$\displaystyle=\mathrm{div}\left(\mu\nabla_{x}\delta_{\mu}G_{a}[\rho,\mu]\right).$

Note that (4) is a joint gradient flow, because the dynamics can be written in
the form

$\partial_{t}\gamma=\mathrm{div}\left(\gamma\nabla\delta_{\gamma}G_{a}(\gamma)\right)\,,$

where $\gamma=(\rho,\mu)$ and where the gradient and divergence are taken in
both variables $(z,x)$. We discuss the structure of the dynamics (4) as well
as the meaning of the different terms appearing in the energy functional
$G_{a}$ in Appendix A.1.

In other settings, such as credit score reporting, the objectives of the
population are 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

$\displaystyle\partial_{t}\rho$
$\displaystyle=-\mathrm{div}\left(\rho\left(\nabla_{z}\delta_{\rho}G_{c}[\rho,x]\right)|_{x=b(\rho)}\right)\,,$
(5) $\displaystyle
b(\rho)\coloneqq\operatorname*{argmin}_{\bar{x}}G_{c}(\rho,\bar{x})\,.$

Next we can consider the population immediately responding to the algorithm,
which has dynamics

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}x$
$\displaystyle=-\nabla_{x}G_{c}(\rho,x)|_{\rho=r(x)}\,,$ (6) $\displaystyle
r(x)\coloneqq\operatorname*{argmin}_{\hat{\rho}\in\mathcal{P}}-G_{c}(\hat{\rho},x)\,.$

In this time-scale separated setting, model (5) represents a dyamic
maximization of $G_{c}$ with respect to $\rho$ in Wasserstein-2 space, and an
instantaneous minimization of $G_{c}$ with respect to the algorithm parameters
$x$. Model (6) represents an instantaneous maximization of $G_{c}$ with
respect to $\rho$ and a dynamic minimization of $G_{c}$ with respect to the
algorithm parameters $x$. The key results on existence and uniqueness of a
ground state as well as the convergence behavior of solutions depend on
convexity (concavity) of $G_{a}$ and $G_{c}$. The notion of convexity that we
will employ for energy functionals in the Wasserstein-2 geometry is (uniform)
displacement convexity, which is analogous to (strong) convexity in Euclidean
spaces. One can think of displacement convexity for an energy functional
defined on $\mathcal{P}_{2}$ 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 a detailed definition of
(uniform) displacement convexity and concavity, see Section A.2. In fact,
suitable convexity properties of the input functions $f_{1},f_{2},W$ and
$\tilde{\rho}$ will ensure (uniform) displacement convexity of the resulting
energy functionals appearing in the gradient flow structure, see for instance
[Vil03, Chapter 5.2].

We make the following assumptions in both the competitive case and aligned
interest cases. Here, $\operatorname{I}_{d}$ denotes the $d\times d$ identity
matrix, $\text{Hess}\left(f\right)$ denotes the Hessian of $f$ in all
variables, while $\nabla^{2}_{x}{f}$ denotes the Hessian of $f$ in the
variable $x$ only.

######  Assumption 1 (Convexity of $f_{1}$ and $f_{2}$).

The functions $f_{1},f_{2}\in
C^{2}(\mathbb{R}^{d}\times\mathbb{R}^{d};[0,\infty))$ satisfy for all
$(z,x)\in\mathbb{R}^{d}\times\mathbb{R}^{d}$ the following:

  * •

There exists constants $\lambda_{1},\lambda_{2}\geq 0$ such that
$\text{Hess}\left(f_{1}\right)\succeq\lambda_{1}\operatorname{I}_{2d}$ and
$\nabla^{2}_{x}{f_{2}}\succeq\lambda_{2}\operatorname{I}_{d}$;

  * •

There exist constants $a_{i}>0$ such that $x\cdot\nabla_{x}f_{i}(z,x)\geq-
a_{i}$ for $i=1,2$;

######  Assumption 2 (Reference Distribution Shape).

The reference distribution $\tilde{\rho}\in\mathcal{P}(\mathbb{R}^{d})\cap
L^{1}(\mathbb{R}^{d})$ satisfies $\log\tilde{\rho}\in C^{2}(\mathbb{R}^{d})$
and
$\nabla^{2}_{z}\log\tilde{\rho}(z)\preceq-\tilde{\lambda}\operatorname{I}_{d}$
for some $\tilde{\lambda}>0$.

######  Assumption 3 (Convex Interaction Kernel).

The interaction kernel $W\in C^{2}(\mathbb{R}^{d};[0,\infty))$ is convex,
symmetric $W(-z)=W(z)$, and for some $D>0$ satisfies

$\displaystyle z\cdot\nabla_{z}W(z)\geq-D,\quad|\nabla_{z}W(z)|\leq
D(1+|z|)\quad\forall\,z\in\mathbb{R}^{d}\,.$

We make the following observations regarding the assumptions above:

  * •

The convexity in Assumption 3 can be relaxed and without affecting the results
outlined below by following a more detailed analysis analogous to the approach
in [CMV03].

  * •

If $f_{1}$ and $f_{2}$ are strongly convex, the proveable convergence rate
increases, but without strict or strong convexity of $f_{1}$ and $f_{2}$, the
regularizers $KL(\rho|\tilde{\rho})$ and $\int\left\lVert
x-x_{0}\right\rVert_{2}^{2}\mathrm{d}x$ 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 $f_{1}$ and $f_{2}$ satisfies Assumption 1.

  * •

Taking the reference measure $\tilde{\rho}$ to be the normal distribution
satisfies Assumption 2, which ensures the distribution is not too flat.

  * •

Taking quadratic interactions $W(z)=\tfrac{1}{2}|z|^{2}$ satisfies Assumption
3.

######  Remark 1 (Cauchy-Problem).

To complete the arguments on convergence to equilibrium, we require sufficient
regularity of solutions to the PDEs under consideration. In fact, it is
sufficient if we can show that equations (4), (5), and (6) can be approximated
by equations with smooth solutions. Albeit tedious, these are standard
techniques in the regularity theory for partial differential equations, see
for example [CMV03, Proposition 2.1 and Appendix A], [OV00], [Vil03, Chapter
9], and the references therein. Similar arguments as in [DV00] are expected to
apply to the coupled gradient 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-diffusion or
non-local interaction as more commonly seen in the literature for multi-
species systems.

######  Theorem 2.

Suppose that Assumptions 1-3 are satisfied and let
$\lambda_{a}:=\lambda_{1}+\min(\lambda_{2}+\beta,\alpha\tilde{\lambda})>0$.
Consider solutions $\gamma_{t}:=(\rho_{t},\mu_{t})$ to the dynamics (4) with
initial conditions satisfying
$\gamma_{0}\in\mathcal{P}_{2}(\mathbb{R}^{d})\times\mathcal{P}_{2}(\mathbb{R}^{d})$
and $G_{a}(\gamma_{0})<\infty.$ Then the following hold:

  1. (a)

There exists a unique minimizer $\gamma_{\infty}=(\rho_{\infty},\mu_{\infty})$
of $G_{a}$, which is also a steady state for equation (4). Moreover,
$\rho_{\infty}\in L^{1}(\mathbb{R}^{d})$, has the same support as
$\tilde{\rho}$, and its density is continuous.

  2. (b)

The solution $\gamma_{t}$ converges exponentially fast in
$G_{a}(\cdot\,|\,\gamma_{\infty})$ and $\overline{W}$,

$G_{a}(\gamma_{t}\,|\,\gamma_{\infty})\leq
e^{-2\lambda_{a}t}G_{a}(\gamma_{0}\,|\,\gamma_{\infty})\,\quad\text{ and
}\quad\overline{W}(\gamma_{t},\gamma_{\infty})\leq
ce^{-\lambda_{a}t}\quad\text{ for all }t\geq 0\,,$

where $c>0$ is a constant only depending on $\gamma_{0}$, $\gamma_{\infty}$
and the parameter $\lambda_{a}$.

###### 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 of 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 $I$ in the case of the Fisher information)111Hence the name HWI
inequalities.. Combining these inequalities with Gronwall’s inequality allows
us to deduce convergence both in energy and in the metric $\overline{W}$. ∎

###  3.3 Analysis of 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
$W_{2}$. 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,
Ber71]) 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 $\rho$ is able to
control the concavity introduced by $f_{1}$ and $f_{2}$.

######  Assumption 4 (Upper bounds for $f_{1}$ and $f_{2}$).

There exists a constant $\Lambda_{1}>0$ such that

$\displaystyle\nabla_{z}^{2}f_{1}(z,x)\preceq\Lambda_{1}I_{d}\qquad\text{ for
all }(z,x)\in\mathbb{R}^{d}\times\mathbb{R}^{d}\,,$

and for any $R>0$ there exists a constant $c_{2}=c_{2}(R)\in\mathbb{R}$ such
that

$\sup_{x\in B_{R}(0)}\int f_{2}(z,x)\mathrm{d}\bar{\rho}(z)<c_{2}\,.$

Equipped with Assumption 4, we state the result for a best-responding
algorithm.

######  Theorem 3.

Suppose Assumptions 1-4 are satisfied with
$\alpha\tilde{\lambda}>\Lambda_{1}$. Let
$\lambda_{b}\coloneqq\alpha\tilde{\lambda}-\Lambda_{1}$. Define
$G_{b}(\rho)\coloneqq G_{c}(\rho,b(\rho))$. Consider a solution $\rho_{t}$ to
the dynamics (5) with initial condition
$\rho_{0}\in\mathcal{P}_{2}(\mathbb{R}^{d})$ such that
$G_{b}(\rho_{0})<\infty$. Then the following hold:

  1. (a)

There exists a unique maximizer $\rho_{\infty}$ of $G_{b}(\rho)$, which is
also a steady state for equation (5). Moreover, $\rho_{\infty}\in
L^{1}(\mathbb{R}^{d})$, has the same support as $\tilde{\rho}$, and its
density is continuous.

  2. (b)

The solution $\rho_{t}$ converges exponentially fast to $\rho_{\infty}$ with
rate $\lambda_{b}$ in $G_{b}(\cdot\,|\,\rho_{\infty})$ and $W_{2}$,

$G_{b}(\rho_{t}\,|\,\rho_{\infty})\leq
e^{-2\lambda_{b}t}G_{a}(\rho_{0}\,|\,\rho_{\infty})\,\quad\text{ and }\quad
W_{2}(\rho_{t},\rho_{\infty})\leq ce^{-\lambda_{b}t}\quad\text{ for all }t\geq
0\,,$

where $c>0$ is a constant only depending on $\rho_{0}$, $\rho_{\infty}$ and
the parameter $\lambda_{b}$.

###### Proof.

(Sketch) The key addition in this setting as compared with Theorem 2 is
proving that $G_{b}(\rho)$ is bounded below, uniformly displacement concave
and guaranteeing its smoothness via Berge’s Maximum Theorem. This is non-
trivial as it uses the properties of the best response $b(\rho)$. A central
observation for our arguments to work is that
$\delta_{\rho}G_{b}[\rho]=\left(\delta_{\rho}G_{c}[\rho,x]\right)|_{x=b(\rho)}$.
We can then conclude using the direct method in the calculus of variations and
the HWI method. ∎

Here, the condition that $\alpha\tilde{\lambda}$ must be large enough
corresponds to the statement that the system must be subjected to a strong
enough regularizing effect.

In the opposite case, where $\rho$ 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 $x$ which allows to deduce
convergence.

######  Theorem 4.

Suppose Assumptions 1-4 are satisfied with
$\alpha\tilde{\lambda}>\Lambda_{1}$, and that $r(x)$ is differentiable (as
shown by example conditions in Lemmas 27 and 28). Define $G_{d}(x)\coloneqq
G_{c}(r(x),x)$. Then it holds:

  1. (a)

There exists a unique minimizer $x_{\infty}$ of $G_{d}(x)$ which is also a
steady state for (6).

  2. (b)

The vector $x(t)$ solving the dynamics (6) with initial condition
$x(0)\in\mathbb{R}^{d}$ converges exponentially fast to $x_{\infty}$ with rate
$\lambda_{d}:=\lambda_{1}+\lambda_{2}+\beta>0$ in $G_{d}$ and in the Euclidean
norm:

$\displaystyle\|x(t)-x_{\infty}\|$ $\displaystyle\leq
e^{-\lambda_{d}t}\|x(0)-x_{\infty}\|\,,$ $\displaystyle
G_{d}(x(t))-G_{d}(x_{\infty})$ $\displaystyle\leq
e^{-2\lambda_{d}t}\left(G_{d}(x(0))-G_{d}(x_{\infty})\right)$

for all $t\geq 0$.

These two theorems illustrate that, under sufficient convexity conditions on
the cost functions, we expect the distribution $\rho$ 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 [Zrn+21] and Stackelberg games [FCR20, FR21,
Fie+21]. We leave this intricate analysis to future work.

In the following section we show numerical results in the case of a best-
responding $x$, best-responding $\rho$, and in between where $x$ and $\rho$
evolve on a 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.

##  4 Numerical Examples

We illustrate numerical results for the case of a classifier, which are used
in scenarios such as loan or government aid applications [CC11], school
admissions [PS13], residency match [Ree18], and recommendation algorithms
[LSW10], all of which have some population which is incentivized to submit
data that will result in a desirable classification. For all examples, we
select classifiers of the form $x\in\mathbb{R}$, so that a data point
$z\in\mathbb{R}$ is assigned a label of $1$ with probability
$q(z,x)=(1+\exp{(-b^{\top}z+x)})^{-1}$ where $b>0$. Let $f_{1}$ and $f_{2}$ be
given by

$\displaystyle f_{1}(z,x)$ $\displaystyle=-\log(1-q(z,x))\,,\qquad
f_{2}(z,x)=-\log q(z,x).$

Note that $\text{Hess}\left(f_{1}\right)\succeq 0$ and
$\nabla_{x}^{2}{f_{2}}\succeq 0$, so $\lambda_{1}=\lambda_{2}=0$. Here, the
strictness of the convexity of the functional is coming from the regularizers,
not the cost functions, with $\tilde{\rho}$ a scaled normal distribution. We
show numerical results for two scenarios with additional settings in the
appendix. First we illustrate competitive interests under three different
timescale settings. Then we simulate the classifier taking an even more naïve
strategy than gradient descent and discuss the results. The PDEs were
implemented based on the finite volume method from [CCH15].

###  4.1 Competitive Objectives

In the setting with competitive objectives, we utilize $G_{c}(\rho,x)$ with
$W=0$, $f_{1}$ and $f_{2}$ as defined above with $b=3$ fixed as it only
changes the steepness of the classifier for $d=1$, and $\alpha=0.1$ and
$\beta=0.05$. In Figure 1, we simulate two extremes of the timescale setting;
first when $\rho$ is nearly best-responding and then when $x$ 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.

Figure 1: When $x$ versus $\rho$ best-responds, we observe the same final
state but different intermediate states. Modes appear in the strategic
population which simpler models cannot capture.

When $\rho$ is nearly best-responding, we see that the distribution quickly
shifts mass over the classifier threshold. Then the classifier shifts right,
correcting for the shift in $\rho$, which then incentivizes $\rho$ to shift
more mass back to the original mode. In contrast, when $x$ best-responds, the
right-hand mode slowly increases in size until the system converges.

Figure 2 shows simulation results from the setting where $\rho$ and $x$ evolve
on the same timescale. We observe that the distribution shift in $\rho$
appears to fall between the two extreme timescale cases, which we expect.

Figure 2: In this experiment the population and classifier have similar rates
of change, and the distribution change for $\rho$ exhibits behaviors from both
the fast $\rho$ and fast $x$ simulations; the right-hand mode does not peak as
high as the fast $\rho$ case but does exceed its final height and return to
the equilibrium.

We highlight two important observations for the competitive case. One, a
single-mode distribution becomes bimodel, which would not be captured using
simplistic metrics such as the mean and variance. This split can be seen as
polarization in the population, a phenomenon that a mean-based strategic
classification model would not capture. Two, the timescale on which the
classifier updates significantly impacts the intermediate behavior of the
distribution. In our example, when $x$ 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 dynamics.

###  4.2 Naïve Behavior

In this example, we explore the results of the classifier adopting a non-
gradient-flow strategy, where the classifier chooses an initially-suboptimal
value for $x$ and does not move, allowing the strategic population to respond.
All functions and parameters are the same as in the previous example.

(a) Both species minimize their respective losses; when the classifier uses a
naïve strategy, the final performance is better for the classifier and
uniformly worse for the population.

(b) The classifier selects a suboptimal initial condition $x=2.2$, instead of
$x=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 our model can
be used to compare robustness of different strategies against a strategic
population.

When comparing with the gradient descent strategy, we observe that while the
initial loss for the classifier is 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 a fixed decision may be
better than the equilibrium, it illustrates how our method provides a
framework for evaluating how different training strategies perform in the long
run against a strategic population.

##  5 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 our
techniques opens a 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 an estimated gradient?
(2) Can we prove convergence in between the two timescale extremes? (3) How do
multiple dynamic populations respond to an algorithm, or multiple algorithms?
In the realm of learning algorithms, our framework can be extended to other
learning update strategies and presents a way to model how we can design these
update strategies to induce desired behaviors in the population.

A challenge in our method is that numerically solving high-dimensional PDEs is
computationally expensive and possibly unfeasible. Here we note that in many
applications, agents in the population do not alter more than a 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 distribution 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 a
direction that increases social welfare.

## 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-20-1-2571 P00003 and NSF Awards
CAREER 1844729, and CPS 1931718. EM acknowledges support from NSF Award
2240110. We are grateful for helpful discussions with José A. Carrillo.

## References

  * [Ala40] Leon Alaoglu “Weak Topologies of Normed Linear Spaces” In _The Annals of Mathematics_ 41.1, 1940, pp. 252 DOI: 10.2307/1968829
  * [Dan67] John M. Danskin “The Theory of Max-Min and its Application to Weapons Allocation Problems” 5, Ã–konometrie und Unternehmensforschung / Econometrics and Operations Research Berlin, Heidelberg: Springer, 1967 DOI: 10.1007/978-3-642-46092-0
  * [Ber71] Dimitri P. Bertsekas “Control of Uncertain Systems with a Set-Membership Description of Uncertainty”, 1971
  * [Pos75] E. Posner “Random coding strategies for minimum entropy” In _IEEE Transactions on Information Theory_ 21.4, 1975, pp. 388–391 DOI: 10.1109/TIT.1975.1055416
  * [Lie83] Elliott H. Lieb “Sharp Constants in the Hardy-Littlewood-Sobolev and Related Inequalities” In _Annals of Mathematics_ 118.2, 1983, pp. 349–374 DOI: 10.2307/2007032
  * [Rud91] Walter Rudin “Functional analysis”, International series in pure and applied mathematics New York: McGraw-Hill, 1991
  * [McC97] Robert J. McCann “A Convexity Principle for Interacting Gases” In _Advances in Mathematics_ 128.1, 1997, pp. 153–179 DOI: 10.1006/aima.1997.1634
  * [BB00] Jean-David Benamou and Yann Brenier “A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem” In _Numerische Mathematik_ 84.3, 2000, pp. 375–393 DOI: 10.1007/s002110050002
  * [DV00] Laurent Desvillettes and Cedric Villani “On the spatially homogeneous landau equation for hard potentials part i : existence, uniqueness and smoothness” In _Communications in Partial Differential Equations_ 25.1, 2000, pp. 179–259 DOI: 10.1080/03605300008821512
  * [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, 2000, pp. 361–400 DOI: 10.1006/jfan.1999.3557
  * [Gho02] B.. Ghosh “Probability Inequalities Related to Markov’s Theorem” In _The American Statistician_ 56.3 [American Statistical Association, Taylor & Francis, Ltd.], 2002, pp. 186–190
  * [CMV03] José A. Carrillo, Robert J. McCann and Cédric Villani “Kinetic equilibration rates for granular media and related equations: entropy dissipation and mass transportation estimates” In _Revista Matemática Iberoamericana_ 19.3, 2003, pp. 971–1018
  * [Vil03] Cédric Villani “Topics in optimal transportation” 58, Graduate Studies in Mathematics American Mathematical Society, Providence, RI, 2003, pp. xvi+370 DOI: 10.1090/gsm/058
  * [Ste04] J. Steele “The Cauchy-Schwarz Master Class: An Introduction to the Art of Mathematical Inequalities” Cambridge University Press, 2004
  * [AB06] “Correspondences” In _Infinite Dimensional Analysis: A Hitchhiker’s Guide_ Berlin, Heidelberg: Springer, 2006, pp. 555–590 DOI: 10.1007/3-540-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, 2006, pp. 217–263 DOI: 10.1007/s00205-005-0386-1
  * [Del06] Chrysanthos Dellarocas “Strategic Manipulation of Internet Opinion Forums: Implications for Consumers and Firms” In _Management Science_ 52.10, 2006, pp. 1577–1593 DOI: 10.1287/mnsc.1060.0567
  * [De ̵07] Camillo De Lellis “Ordinary Differential Equations with Rough Coefficients and the Renormalization Theorem of Ambrosio” In _Séminaire Bourbaki_ 59, 972, 2007, pp. 175–204
  * [BCL09] Andrea L. Bertozzi, Jose A. Carrillo and Thomas Laurent “Blow-up in multidimensional aggregation equations with mildly singular interaction kernels*” In _Nonlinearity_ 22.3, 2009, pp. 683 DOI: 10.1088/0951-7715/22/3/009
  * [BBS09] Steffen Bickel, Michael Brückner and Tobias Scheffer “Discriminative Learning Under Covariate Shift” In _The Journal of Machine Learning Research_ 10, 2009, pp. 2137–2155
  * [Bor09] Vivek S Borkar “Stochastic approximation: a dynamical systems viewpoint” 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, 2010, pp. 595–612 DOI: 10.1287/opre.1090.0741
  * [LSW10] Juan Lang, Matt Spear and S. Wu “Social Manipulation of Online Recommender Systems” In _Social Informatics_ , Lecture Notes in Computer Science Berlin, Heidelberg: Springer, 2010, pp. 125–139 DOI: 10.1007/978-3-642-16567-2_10
  * [CC11] Adriana Camacho and Emily Conover “Manipulation of Social Program Eligibility” In _American Economic Journal: Economic Policy_ 3.2, 2011, pp. 41–65
  * [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, 2011 DOI: 10.1215/00127094-2010-211
  * [BCY12] D. Balagu’e, J. Carrillo and Y. Yao “Confinement for repulsive-attractive kernels” In _Discrete and Continuous Dynamical Systems-series 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, 2012, pp. 1140005 DOI: 10.1142/S0218202511400057
  * [FF13] Marco Di Francesco and Simone Fagioli “Measure solutions for non-local interaction PDEs with two species” In _Nonlinearity_ 26.10, 2013, pp. 2777–2808 DOI: 10.1088/0951-7715/26/10/2777
  * [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, 2013, pp. 80–106 DOI: 10.1257/aer.103.1.80
  * [CCH14] José Antonio Carrillo, Michel Chipot and Yanghong Huang “On global minimizers of repulsive-attractive power-law interaction energies” In _Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences_ 372.2028, 2014, pp. 20130399 DOI: 10.1098/rsta.2013.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 Finite-Volume Method for Nonlinear Nonlocal Equations with a Gradient Flow Structure” In _Communications in Computational Physics_ 17.1, 2015, pp. 233–258 DOI: 10.4208/cicp.160214.010814a
  * [San15] Filippo Santambrogio “Optimal Transport for Applied Mathematicians. Calculus of Variations, PDEs and Modeling”, 2015
  * [Har+16] Moritz Hardt, Nimrod Megiddo, Christos Papadimitriou and Mary Wootters “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, 2016, pp. 111–122 DOI: 10.1145/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 Mathematics_ 78.2, 2018, pp. 1078–1104 DOI: 10.1137/17M1128782
  * [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-proof mechanisms: Evidence from the residency match” In _Games and Economic Behavior_ 108, Special Issue in Honor of Lloyd Shapley: Seven Topics in Game Theory, 2018, pp. 317–330 DOI: 10.1016/j.geb.2017.04.011
  * [Kuh+19] Daniel Kuhn, Peyman Mohajerin Esfahani, Viet Anh Nguyen and Soroosh Shafieezadeh-Abadeh “Wasserstein Distributionally Robust Optimization: Theory and Applications in Machine Learning” In _Operations Research & Management Science in the Age of Analytics_ INFORMS, 2019, pp. 130–166 DOI: 10.1287/educ.2019.0198
  * [DS20] Tomasz Debiec and Markus Schmidtchen “Incompressible Limit for a Two-Species Tumour Model with Coupling Through Brinkman’s Law in One Dimension” In _Acta Applicandae Mathematicae_ 169.1, 2020, pp. 593–611 DOI: 10.1007/s10440-020-00313-1
  * [DT20] Manh Hong Duong and Julian Tugaut “Coupled McKean-Vlasov diffusions: wellposedness, propagation of chaos and invariant measures” In _Stochastics_ 92.6, 2020, pp. 900–943 DOI: 10.1080/17442508.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_ , 2020, pp. 3133–3144 PMLR
  * [Mil+20] John Miller, Karl Krauth, Benjamin Recht and Ludwig Schmidt “The Effect of Natural Distribution Shift on Question Answering Models” In _Proceedings of the 37th International Conference on Machine Learning_ PMLR, 2020, pp. 6905–6916
  * [Per+20] Juan Perdomo, Tijana Zrnic, Celestine Mendler-Dünner and Moritz Hardt “Performative Prediction” In _Proceedings of the 37th International Conference on Machine Learning_ PMLR, 2020, pp. 7599–7609
  * [AB21] Abdulaziz Alsenafi and Alethea B.. Barbaro “A Multispecies Cross-Diffusion Model for Territorial Development” In _Mathematics_ 9.12, 2021, pp. 1428 DOI: 10.3390/math9121428
  * [Dom+21] Carles Domingo-Enrich et al. “A mean-field analysis of two-player zero-sum games” arXiv, 2021 DOI: 10.48550/arXiv.2002.06277
  * [DN21] John C. Duchi and Hongseok Namkoong “Learning models with uniform performance via distributionally robust optimization” In _The Annals of Statistics_ 49.3, 2021, pp. 1378–1406 DOI: 10.1214/20-AOS2004
  * [FR21] Tanner Fiez and Lillian J Ratliff “Local convergence analysis of gradient descent ascent with finite timescale separation” In _Proceedings of the International Conference on Learning Representation_ , 2021
  * [Fie+21] Tanner Fiez et 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
  * [IYZ21] 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_ PMLR, 2021, pp. 4641–4650
  * [LHL21] Qi Lei, Wei Hu and Jason D. Lee “Near-Optimal Linear Regression under Distribution Shift” In _ArXiv_ , 2021
  * [MPZ21] John P Miller, Juan C Perdomo and Tijana Zrnic “Outside the echo chamber: Optimizing the performative risk” In _International Conference on Machine Learning_ , 2021, pp. 7710–7720 PMLR
  * [Wil+21] Olivia Wiles et al. “A Fine-Grained Analysis on Distribution Shift” In _ArXiv_ , 2021
  * [Zrn+21] Tijana Zrnic, Eric Mazumdar, Shankar Sastry and Michael Jordan “Who Leads and Who Follows in Strategic Classification?” In _Advances in Neural Information Processing Systems_ 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_ PMLR, 2022, pp. 2704–2729
  * [Giu+22] Valeria Giunta, Thomas Hillen, Mark Lewis and Jonathan R. Potts “Local and Global Existence for Nonlocal Multispecies Advection-Diffusion Models” In _SIAM Journal on Applied Dynamical Systems_ 21.3, 2022, pp. 1686–1708 DOI: 10.1137/21M1425992
  * [HPS22] Georg Heinze, Jan-Frederik Pietschmann and Markus Schmidtchen “Nonlocal cross-interaction systems on graphs: Nonquadratic Finslerian structure and nonlinear mobilities” arXiv, 2022 DOI: 10.48550/arXiv.2107.11289
  * [JPZ22] Ansgar Jungel, Stefan Portisch and Antoine Zurek “Nonlocal cross-diffusion systems for multi-species populations and networks” In _Nonlinear Analysis_ 219, 2022, pp. 112800 DOI: 10.1016/j.na.2022.112800
  * [LY22] Guanlin Li and Yao Yao “Two-species competition model with chemotaxis: well-posedness, stability and dynamics” In _Nonlinearity_ 35.3, 2022, pp. 1329 DOI: 10.1088/1361-6544/ac4a8d
  * [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 on Artificial Intelligence and Statistics_ , 2022, pp. 5891–5912 PMLR
  * [Ray+22] Mitas Ray, Lillian J Ratliff, Dmitriy Drusvyatskiy and Maryam Fazel “Decision-dependent risk minimization in geometrically decaying dynamic environments” In _Proceedings of the AAAI Conference on Artificial Intelligence_ 36.7, 2022, pp. 8081–8088
  * [CFG23] J. Carrillo, Alejandro Fernandez-Jimenez and D. G’omez-Castro “Partial mass concentration for fast-diffusions with non-local aggregation terms”, 2023
  * [Dou+23] Marie Doumic, Sophie Hecht, Benoit Perthame and Diane Peurichard “Multispecies cross-diffusions: from a nonlocal mean-field to a porous medium system without self-diffusion” arXiv, 2023 DOI: 10.48550/arXiv.2306.01777
  * [GG23] Camilo Garcia Trillos and Nicolas Garcia Trillos “On adversarial robustness and the use of Wasserstein ascent-descent dynamics to enforce it” arXiv, 2023 DOI: 10.48550/arXiv.2301.03662
  * [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-Field Perspective” arXiv, 2023 DOI: 10.48550/arXiv.2212.08791
  * [WC23] Guillaume Wang and Lenaic Chizat “An Exponentially Converging Particle Method for the Mixed Nash Equilibrium of Continuous Games” arXiv, 2023 DOI: 10.48550/arXiv.2211.01280
  * [Qui+] Joaquin Quinonero-Candela, Masashi Sugiyama, Anton Schwaighofer and Neil D. Lawrence “Dataset Shift in Machine Learning”

##  Appendix 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 are 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

$\displaystyle\partial_{t}\rho$
$\displaystyle=\alpha\Delta\rho+\mathrm{div}\left(\rho\nabla_{z}\left(\int
f_{1}\mathrm{d}\mu-\alpha\log\tilde{\rho}+W\ast\rho\right)\right)\,,$ (7a)
$\displaystyle\partial_{t}\mu$
$\displaystyle=\mathrm{div}\left(\mu\nabla_{x}\left(\int
f_{1}\mathrm{d}\rho+\int
f_{2}\mathrm{d}\bar{\rho}+\frac{\beta}{2}\|x-x_{0}\|^{2}\right)\right)\,.$
(7b)

In other words, the population $\rho$ in (7a) is subject to an isotropic
diffusive force with diffusion coefficient $\alpha>0$, a drift force due to
the time-varying confining potential $\int
f_{1}\mathrm{d}\mu(t)-\alpha\log\tilde{\rho}$, and a self-interaction force
via the interaction potential $W$. If we consider the measure $\mu$ to be
given and fixed in time, this corresponds exactly to the type of parabolic
equation studied in [CMV03]. Here however the dynamics are more complex due to
the coupling of the confining potential with the dynamics (7b) for $\mu(t)$
via the coupling potential $f_{1}$. Before presenting the analysis of this
model, let us give a bit more intuition on the meaning and the structure of
these dynamics.

In the setting where $\mu$ represents a binary classifier, we can think of the
distribution $\bar{\rho}$ as modelling all those individuals carrying the true
label 1, say, and the distribution $\rho(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 $\int f_{1}(z,x)\mu(t,\mathrm{d}x)$
represents a penalty for incorrectly classifying an individual at $z$ with
true label 0 when using the classifier $\mu(t,x)$. In other words, $\int
f_{1}(z,x)\mu(t,\mathrm{d}x)\in[0,\infty)$ is increasingly large the more $z$
digresses from the correct classification 0. Similarly, $\int
f_{1}(z,x)\rho(t,\mathrm{d}z)\in[0,\infty)$ is increasingly large if the
population $\rho$ shifts mass to locations in $z$ where the classification is
incorrect. The terminology _aligned objectives_ refers to the fact that in (7)
both the population and the classifier are trying to evolve in a way as to
maximize correct classification. Analogously, the term $\int
f_{2}(z,x)\bar{\rho}(\mathrm{d}z)$ is large if $x$ would incorrectly classify
the population $\bar{\rho}$ that carries the label 1. A natural extension of
the model (7) would be a setting where also the population carrying labels 1
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 $\rho(t)$, one for the evolution of
$\bar{\rho}(t)$ and one for the classifier $\mu(t)$.

The term

$\alpha\Delta\rho-\alpha\mathrm{div}\left(\rho\nabla\log\tilde{\rho}\right)=\alpha\mathrm{div}\left(\rho\nabla\delta_{\rho}KL(\rho\,|\,\tilde{\rho})\right)$

forces the evolution of $\rho(t)$ to approach $\tilde{\rho}$. In other words,
it penalizes (in energy) deviations from a given reference measure
$\tilde{\rho}$. In the context of the application at hand, we take
$\tilde{\rho}$ to be the initial distribution $\rho(t=0)$. The solution
$\rho(t)$ then evolves away from $\tilde{\rho}$ over time due to the other
forces that are present. Therefore, the term $KL(\rho\,|\,\tilde{\rho})$ in
the energy both provides smoothing of the flow and a penalization for
deviations away from the reference measure $\tilde{\rho}$.

The self-interaction term $W\ast\rho$ 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 a radial function, i.e. $W(z)=\omega(|z|)$. A
choice of $\omega:\mathbb{R}\to\mathbb{R}$ such that $\omega^{\prime}(r)>0$
corresponds to an _attractive_ force between individuals, whereas
$\omega^{\prime}(r)<0$ corresponds to a _repulsive_ force. The statement
$|z|\omega^{\prime}(|z|)=z\cdot\nabla_{z}W(z)\geq-D$ in Assumption 3 therefore
corresponds to a requirement that the self-interaction force is not too
repulsive. Neglecting all other forces in (7a), we obtain the non-local
interaction equation

$\partial_{t}\rho=\mathrm{div}\left(\rho\nabla W\ast\rho\right)$

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+11, BCY12, CCH14, BCL09,
BLL12, CMV06, CFG23] and references therein.

Using the results from these works, our assumptions on the interaction
potential $W$ can be relaxed in many ways, for example by allowing
discontinuous derivatives at zero for $W$, or by allowing $W$ to be negative.

The dynamics (7b) for the algorithm $\mu$ is a non-autonomous transport
equation,

$\partial_{t}\mu=\mathrm{div}\left(\mu v\right)\,,$

where the time-dependence in the velocity field

$v(t,x):=\nabla_{x}\left(\int f_{1}(z,x)\mathrm{d}\rho(t,z)+\int
f_{2}(z,x)\mathrm{d}\tilde{\rho}(z)+\frac{\beta}{2}\|x-x_{0}\|^{2}\right)\,,$

comes through the evolving population $\rho(t)$. This structure allows to
obtain an explicit solution for $\mu(t)$ in terms of the initial condition
$\mu_{0}$ and the solution $\rho(t)$ to (7a) using the method of
characteristics.

######  Proposition 5.

Assume that there exists a constant $c>0$ such that

$\displaystyle\left\|\int\nabla_{x}f_{1}(z,x)\mathrm{d}\rho(z)+\int\nabla_{x}f_{2}(z,x)\mathrm{d}\bar{\rho}(z)\right\|$
$\displaystyle\leq
c(1+\|x\|)\quad\forall\rho\in\mathcal{P}_{2}(\mathbb{R}^{d})\text{ and
}\forall x\in\mathbb{R}^{d}\,.$ (8)

Then the unique distributional solution $\mu(t)$ to (7b) is given by

$\mu(t)=\Phi(t,0,\cdot)_{\\#}\mu_{0}\,,$ (9)

where $\Phi(t,s,x)$ solves the characteristic equation

$\displaystyle\partial_{s}\Phi(s,t,x)+v(s,\Phi(s,t,x))=0\,,\qquad\Phi(t,t,x)=x\,.$
(10)

###### Proof.

Thanks to Assumption 1, we have that $v\in
C^{1}(\mathbb{R}\times\mathbb{R}^{d};\mathbb{R}^{d})$, and by (8), we have

$\|v(t,x)\|\leq c(1+\|x\|)\quad\text{ for all }t\geq 0,x\in\mathbb{R}^{d}\,.$

By classical Cauchy-Lipschitz theory for ODEs, this guarantees the existence
of a unique global solution $\Phi(t,s,x)$ solving (10). Then it can be checked
directly that $\mu(t)$ as defined in (9) is a distributional solution to (7b).
∎

In the characteristic equation (10), $\Phi(s,t,x)$ is a parametrization of all
trajectories: if a particle was at location $x$ at time $t$, then it is at
location $\Phi(s,t,x)$ at time $s$. Our assumptions on $f_{1},f_{2}$ and
$\bar{\rho}$ also ensure that
$\Phi(s,t,\cdot):\mathbb{R}^{d}\to\mathbb{R}^{d}$ is a $C^{1}$-diffeomorphism
for all $s,t\in\mathbb{R}$. For more details on transport equations, see for
example [De ̵07].

######  Remark 2.

Consider the special case where $\mu_{0}=\delta_{x_{0}}$ for some initial
position $x_{0}\in\mathbb{R}^{d}$. Then by Proposition 5, the solution to (7b)
is given by $\mu(t)=\delta_{x(t)}$, where $x(t):=\Phi(t,0,x_{0})$ solves the
ODE

$\dot{x}(t)=-v(t,x(t))\,,\qquad x(0)=x_{0}\,,$

which is precisely of type (3).

For the case of competing objectives, the two models we consider can be
written as

$\displaystyle\partial_{t}\rho$
$\displaystyle=-\mathrm{div}\left(\rho\left[\nabla(f_{1}(z,b(\rho))-\alpha\log(\rho/\tilde{\rho})-W\ast\rho\right]\right)\,,$
$\displaystyle b(\rho)$ $\displaystyle:=\operatorname*{argmin}_{\bar{x}}\int
f_{1}(z,\bar{x})\mathrm{d}\rho(z)+\int
f_{2}(\bar{x},z^{\prime})\mathrm{d}\bar{\rho}(z^{\prime})+\frac{\beta}{2}\left\lVert\bar{x}-x_{0}\right\rVert^{2}\,$

for (5), and

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}x$
$\displaystyle=-\nabla_{x}\left(\int f_{1}(z,x)\,r(x)(\mathrm{d}z)+\int
f_{2}(x,z^{\prime})\mathrm{d}\bar{\rho}(z^{\prime})+\frac{\beta}{2}\left\lVert
x-x_{0}\right\rVert^{2}\right)\,,$ $\displaystyle
r(x)\coloneqq\operatorname*{argmax}_{\hat{\rho}\in\mathcal{P}}\int
f_{1}(z,x)\mathrm{d}\hat{\rho}(z)-\alpha
KL(\hat{\rho}|\tilde{\rho})-\frac{1}{2}\int\hat{\rho}W\ast\hat{\rho}\,.$

for (6).

###  A.2 Definitions and notation

Here, and in what follows, $\operatorname{I}_{d}$ denotes the $d\times d$
identity matrix, and $\operatorname{id}$ denotes the identity map. The energy
functionals we are considering are usually defined on the set of probability
measures on $\mathbb{R}^{d}$, denoted by $\mathcal{P}(\mathbb{R}^{d})$. If we
consider the subset $\mathcal{P}_{2}(\mathbb{R}^{d})$ of probability measures
with bounded second moment,

$\mathcal{P}_{2}(\mathbb{R}^{d}):=\\{\rho\in\mathcal{P}(\mathbb{R}^{d})\,:\,\int_{\mathbb{R}^{d}}\|z\|^{2}\mathrm{d}\rho(z)<\infty\\}\,,$

then we can endow this space with the Wasserstein-2 metric.

######  Definition 4 (Wasserstein-2 Metric).

The Wasserstein-2 metric between two probability measures
$\mu,\nu\in\mathcal{P}_{2}(\mathbb{R}^{d})$ is given by

$\displaystyle
W_{2}(\mu,\nu)^{2}=\inf_{\gamma\in\Gamma(\mu,\nu)}\int\left\lVert
z-z^{\prime}\right\rVert_{2}^{2}\mathrm{d}\gamma(z,z^{\prime})$

where $\Gamma$ is the set of all joint probability distributions with
marginals $\mu$ and $\nu$, i.e.
$\mu(\mathrm{d}z)=\int\gamma(\mathrm{d}z,z^{\prime})\mathrm{d}z^{\prime}$ and
$\nu(\mathrm{d}z^{\prime})=\int\gamma(z,\mathrm{d}z^{\prime})\mathrm{d}z$.

The restriction to $\mathcal{P}_{2}(\mathbb{R}^{d})$ ensures that $W_{2}$ is
always finite. Then the space $(\mathcal{P}_{2}(\mathbb{R}^{d}),W_{2})$ is
indeed a metric space. We will make use of the fact that $W_{2}$ metrizes
narrow convergence of probability measures. To make this statement precise,
let us introduce two common notions of convergence for probability measures,
which are a subset of the finite signed Radon measures
$\mathcal{M}(\mathbb{R}^{d})$.

######  Definition 5.

Consider a sequence $(\mu_{n})\in\mathcal{M}(\mathbb{R}^{d})$ and a limit
$\mu\in\mathcal{M}(\mathbb{R}^{d})$.

  * •

(Narrow topology) The sequence $(\mu_{n})$ converges _narrowly_ to $\mu$,
denoted by $\mu_{n}\rightharpoonup\mu$, if for all continuous bounded
functions $f:\mathbb{R}^{d}\to\mathbb{R}$,

$\int_{\mathbb{R}^{d}}f(z)\mathrm{d}\mu_{n}(z)\to\int_{\mathbb{R}^{d}}f(z)\mathrm{d}\mu(z)\,.$

  * •

(Weak-$*$ topology) The sequence $(\mu_{n})$ converges _weakly- $*$_ to $\mu$,
denoted by $\mu_{n}\overset{\ast}{\rightharpoonup}\mu$, if for all continuous
functions vanishing at infinity (i.e. $f:\mathbb{R}^{d}\to\mathbb{R}$ such
that for all $\epsilon>0$ there exists a compact set
$K_{\epsilon}\subset\mathbb{R}^{d}$ such that $|f(z)|<\epsilon$ on
$\mathbb{R}^{d}\setminus K_{\epsilon}$), we have

$\int_{\mathbb{R}^{d}}f(z)\mathrm{d}\mu_{n}(z)\to\int_{\mathbb{R}^{d}}f(z)\mathrm{d}\mu(z)\,.$

Let us denote the set of continuous functions on $\mathbb{R}^{d}$ vanishing at
infinity by $C_{0}(\mathbb{R}^{d})$, and the set of continuous bounded
functions by $C_{b}(\mathbb{R}^{d})$. Note that narrow convergence immediately
implies that $\mu_{n}(\mathbb{R}^{d})\to\mu(\mathbb{R}^{d})$ as the constant
function is in $C_{b}(\mathbb{R}^{d})$. This is not necessarily true for
weak-$*$ convergence. We will later make use of the Banach-Alaoglu theorem
[Ala40], which gives weak-$*$ compactness of the unit ball in a dual space.
Note that $\mathcal{M}(\mathbb{R}^{d})$ is indeed the dual of
$C_{0}(\mathbb{R}^{d})$ endowed with the sup-norm, and
$\mathcal{P}(\mathbb{R}^{d})$ is the unit ball in
$\mathcal{M}(\mathbb{R}^{d})$ using the dual norm. Moreover, if we can ensure
that mass does not escape to infinity, the two notions of convergence in
Definition 5 are in fact equivalent.

######  Lemma 6.

Consider a sequence $(\mu_{n})\in\mathcal{M}(\mathbb{R}^{d})$ and a measure
$\mu\in\mathcal{M}(\mathbb{R}^{d})$. Then $\mu_{n}\rightharpoonup\mu$ if and
only if $\mu_{n}\overset{\ast}{\rightharpoonup}\mu$ and
$\mu_{n}(\mathbb{R}^{d})\to\mu(\mathbb{R}^{d})$.

This follows directly from Definition 5. Here, the condition
$\mu_{n}(\mathbb{R}^{d})\to\mu(\mathbb{R}^{d})$ is equivalent to tightness of
$(\mu_{n})$, and follows from Markov’s inequality [Gho02] 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

$\int\|z\|^{2}\mathrm{d}\mu_{n}(z)<C\qquad\forall n\in\mathrm{N}\,.$ (11)

######  Definition 6 (Tightness of probability measures).

A collection of measures $(\mu_{n})\in\mathcal{M}(\mathbb{R}^{d})$ is _tight_
if for all $\epsilon>0$ there exists a compact set
$K_{\epsilon}\subset\mathbb{R}^{d}$ such that
$|\mu_{n}|(\mathbb{R}^{d}\setminus K_{\epsilon})<\epsilon$ for all
$n\in\mathrm{N}$, where $|\mu|$ denotes the total variation of $\mu$.

Another classical result is that the Wasserstein-2 metric metrizes narrow
convergence and weak-$*$ convergence of probability measures, see for example
[San15, Theorem 5.11] or [Vil03, Theorem 7.12].

######  Lemma 7.

Let $\mu_{n},\mu\in\mathcal{P}_{2}(\mathbb{R}^{d})$. Then
$W_{2}(\mu_{n},\mu)\to 0$ if and only if

$\displaystyle\mu_{n}\rightharpoonup\mu\quad\text{ and
}\quad\int_{\mathbb{R}^{d}}\|z\|^{2}\mathrm{d}\mu_{n}(z)\to\int_{\mathbb{R}^{d}}\|z\|^{2}\mathrm{d}\mu(z)\,.$

######  Remark 3.

Note that $\mu_{n}\rightharpoonup\mu$ can be replaced by
$\mu_{n}\overset{\ast}{\rightharpoonup}\mu$ in the above statement thanks to
the fact that the limit $\mu$ is a probability measure with mass 1, see Lemma
6.

Next, we consider two measures $\mu,\nu\in\mathcal{P}(\mathbb{R}^{d})$ that
are _atomless_ , i.e $\mu(\\{z\\})=0$ for all $z\in\mathbb{R}^{d}$. By
Brenier’s theorem [BB00] (also see [Vil03, Theorem 2.32]) there exists a
unique measurable map $T:\mathbb{R}^{d}\to\mathbb{R}^{d}$ such that
$T_{\\#}\mu=\nu$, and $T=\nabla\psi$ for some convex function
$\psi:\mathbb{R}^{d}\to\mathbb{R}$. Here, the _push-forward_ operator
$\nabla\psi_{\\#}$ is defined as

$\displaystyle\int_{\mathbb{R}^{d}}f(z)\mathrm{d}\nabla\psi_{\\#}\rho_{0}(z)=\int_{\mathbb{R}^{d}}f(\nabla\psi(z))\mathrm{d}\rho_{0}(z)$

for all Borel-measurable functions $f:\mathbb{R}^{d}\mapsto\mathbb{R}_{+}$. If
$\rho_{1}=\nabla\psi_{\\#}\rho_{0}$, we denote by
$\rho_{s}=[(1-s)\operatorname{id}+s\nabla\psi]_{\\#}\rho_{0}$ the
_discplacement interpolant_ between $\rho_{0}$ and $\rho_{1}$. We are now
ready to introduce the notion of displacement convexity, which is the same as
geodesic convexity in the geodesic space
$(\mathcal{P}_{2}(\mathbb{R}^{d}),W_{2})$. We will state the definition here
for atomless measures, but it can be relaxed to any pair of measures in
$\mathcal{P}_{2}$ using optimal transport plans instead of transport maps. In
what follows, we will use $s$ to denote the interpolation parameter for
geodesics, and $t$ to denote time related to solutions of (4), (5) and (6).

######  Definition 7 (Displacement Convexity).

A functional $G:\mathcal{P}\mapsto\mathbb{R}$ is _displacement convex_ if for
all $\rho_{0},\rho_{1}$ that are atomless we have

$\displaystyle G(\rho_{s})\leq(1-s)G(\rho_{0})+sG(\rho_{1})\,,$

where $\rho_{s}=[(1-s)\operatorname{id}+s\nabla\psi]_{\\#}\rho_{0}$ is the
displacement interpolant between $\rho_{0}$ and $\rho_{1}$. Further,
$G:\mathcal{P}\mapsto\mathbb{R}$ is _uniformly displacement convex with
constant $\eta>0$_ if

$\displaystyle
G(\rho_{s})\leq(1-s)G(\rho_{0})+sG(\rho_{1})-s(1-s)\frac{\eta}{2}W_{2}(\rho_{0},\rho_{1})^{2}\,,$

where $\rho_{s}=[(1-s)\operatorname{id}+s\nabla\psi]_{\\#}\rho_{0}$ is the
displacement interpolant between $\rho_{0}$ and $\rho_{1}$.

######  Remark 4.

In other words, $G$ is displacement convex (concave) if the function
$G(\rho_{s})$ is convex (concave) with
$\rho_{s}=[(1-s\operatorname{id}+s\nabla\psi]_{\\#}\rho_{0}$ being the
displacement interpolant between $\rho_{0}$ and $\rho_{1}$. Contrast this with
the classical notion of convexity (concavity) for $G$, where we require that
the function $G((1-s)\rho_{0}+s\rho_{1})$ is convex (concave).

In fact, if the energy $G$ is twice differentiable along geodesics, then the
condition $\frac{\mathrm{d}^{2}}{\mathrm{d}s^{2}}G(\gamma_{s})\geq 0$ along
any geodesic $(\rho_{s})_{s\in[0,1]}$ between $\rho_{0}$ and $\rho_{1}$ is
sufficient to obtain displacement convexity. Similarly, when
$\frac{\mathrm{d}^{2}}{\mathrm{d}s^{2}}G(\rho_{s})\geq\eta
W_{2}(\rho_{0},\rho_{1})^{2}$, then $G$ is uniformly displacement convex with
constant $\eta>0$. For more details, see [McC97] and [Vil03, Chapter 5.2].

###  A.3 Steady states

The main goal in our theoretical analysis is to characterize the asymptotic
behavior for the models (4), (5) and (6) as time goes to 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 a
critical point of the energy, and therefore satisfies that the first variation
is constant on disconnected components of its support. If this ground state
also has enough regularity (weak differentiability) to be 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 8 (Steady states for (4)).

Given $\rho_{\infty}\in L^{1}_{+}(\mathbb{R}^{d})\cap
L^{\infty}_{loc}(\mathbb{R}^{d})$ with $\|\rho_{\infty}\|_{1}=1$ and
$\mu_{\infty}\in\mathcal{P}_{2}(\mathbb{R}^{d})$, then
$(\rho_{\infty},\mu_{\infty})$ is a steady state for the system (4) if
$\rho_{\infty}\in W^{1,2}_{loc}(\mathbb{R}^{d})$, $\nabla
W\ast\rho_{\infty}\in L^{1}_{loc}(\mathbb{R}^{d})$, $\rho_{\infty}$ is
absolutely continuous with respect to $\tilde{\rho}$, and
$(\rho_{\infty},\mu_{\infty})$ satisfy

$\displaystyle\nabla_{z}\left(\int
f_{1}(z,x)\mathrm{d}\mu_{\infty}(x)+\alpha\log\left(\frac{\rho_{\infty}(z)}{\tilde{\rho}(z)}\right)+W\ast\rho_{\infty}(z)\right)=0\qquad\forall
z\in\operatorname{supp}(\rho_{\infty})\,,$ (12a)
$\displaystyle\nabla_{x}\left(\int f_{1}(z,x)\mathrm{d}\rho_{\infty}(z)+\int
f_{2}(z,x)\mathrm{d}\tilde{\rho}(z)+\frac{\beta}{2}\|x-x_{0}\|^{2}\right)=0\qquad\forall
x\in\operatorname{supp}(\mu_{\infty})$ (12b)

in the sense of distributions.

Here, $L^{1}_{+}(\mathbb{R}^{d}):=\\{\rho\in
L^{1}(\mathbb{R}^{d})\,:\,\rho\geq 0\\}$.

######  Definition 9 (Steady states for (5)).

Let $\rho_{\infty}\in L^{1}_{+}(\mathbb{R}^{d})\cap
L^{\infty}_{loc}(\mathbb{R}^{d})$ with $\|\rho_{\infty}\|_{1}=1$. Then
$\rho_{\infty}$ is a steady state for the system (5) if $\rho_{\infty}\in
W^{1,2}_{loc}(\mathbb{R}^{d})$, $\nabla W\ast\rho_{\infty}\in
L^{1}_{loc}(\mathbb{R}^{d})$, $\rho_{\infty}$ is absolutely continuous with
respect to $\tilde{\rho}$, and $\rho_{\infty}$ satisfies

$\displaystyle\nabla_{z}\left(f_{1}(z,b(\rho_{\infty}))-\alpha\log\left(\frac{\rho_{\infty}(z)}{\tilde{\rho}(z)}\right)-W\ast\rho_{\infty}(z)\right)=0\qquad\forall
z\in\mathbb{R}^{d}\,,$ (13)

in the sense of distributions, where
$b(\rho_{\infty})\coloneqq\operatorname*{argmin}_{x}G_{c}(\rho_{\infty},x)$.

######  Definition 10 (Steady states for (6)).

The vector $x_{\infty}\in\mathbb{R}^{d}$ is a steady state for the system (6)
if it satisfies

$\displaystyle\nabla_{x}G_{d}(x_{\infty})=0\,.$

In fact, with the above notions of steady state, we can obtain improved
regularity for $\rho_{\infty}$.

######  Lemma 8.

Let Assumptions 1-3 hold. Then the steady states $\rho_{\infty}$ for (4) and
(5) are continuous.

###### Proof.

We present here the argument for equation (5) only. The result for (4) follows
in exactly the same way by replacing $f_{1}(z,b(\rho_{\infty}))$ with $-\int
f_{1}(z,x)\mathrm{d}\mu_{\infty}(x)$.

Thanks to our assumptions, we have
$f_{1}(\cdot,b(\rho_{\infty}))+\alpha\log\tilde{\rho}(\cdot)\in C^{1}$, which
implies that
$\nabla(f_{1}(\cdot,b(\rho_{\infty}))+\alpha\log\tilde{\rho}(\cdot))\in
L_{loc}^{\infty}$. By the definition of a steady state, $\rho_{\infty}\in
L^{1}\cap L_{loc}^{\infty}$ and thanks to Assumption 3 we have $W\in C^{2}$,
which implies that $\nabla W\ast\rho_{\infty}\in L_{loc}^{\infty}$. Let

$\displaystyle
h(z)\coloneqq\rho_{\infty}(z)\nabla\left[f_{1}(z,b(\rho_{\infty}))+\alpha\log\tilde{\rho}(z)-(W\ast\rho_{\infty})(z)\right]\,.$

Then by the aforementioned regularity, we obtain $h\in L_{loc}^{1}\cap
L_{loc}^{\infty}$. By interpolation, it follows that $h\in L_{loc}^{p}$ for
all $1<p<\infty$. This implies that $\mathrm{div}\left(h\right)\in
W_{loc}^{-1,p}$. Since $\rho_{\infty}$ is a weak $W_{loc}^{1,2}$-solution of
(13), we have

$\displaystyle\Delta\rho_{\infty}=\mathrm{div}\left(h\right)\,,$

and so by classic elliptic regularity theory we conclude $\rho_{\infty}\in
W_{loc}^{1,p}$. Finally, applying Morrey’s inequality, we have
$\rho_{\infty}\in C^{0,k}$ where $k=\frac{p-d}{p}$ for any $d<p<\infty$.
Therefore $\rho_{\infty}\in C(\mathbb{R}^{d})$ (after possibly being redefined
on a set of measure zero). ∎

##  Appendix B Proof of Theorem 2

For ease of notation, we write
$G_{a}:\mathcal{P}(\mathbb{R}^{d})\times\mathcal{P}(\mathbb{R}^{d})\mapsto[0,\infty]$
as

$G_{a}((\rho,\mu))=\alpha
KL(\rho|\tilde{\rho})+\mathcal{V}(\rho,\mu)+\mathcal{W}(\rho)\,,$

where we define

$\displaystyle\mathcal{V}(\rho,\mu)=\iint
f_{1}(z,x)\mathrm{d}\rho(z)\mathrm{d}\mu(x)+\int V(x)\mathrm{d}\mu(x)\,,$
$\displaystyle\mathcal{W}(\rho)=\frac{1}{2}\iint
W(z_{1}-z_{2})\mathrm{d}\rho(z_{1})\mathrm{d}\rho(z_{2})\,,$

with potential given by $V(x):=\int
f_{2}(z,x)\mathrm{d}\bar{\rho}(z)+\frac{\beta}{2}\|x-x_{0}\|^{2}$.

In order to prove the existence of a unique ground state for $G_{a}$, a
natural approach is to consider the corresponding Euler-Lagrange equations

$\displaystyle\alpha\log\frac{\rho(z)}{\tilde{\rho}(z)}+\int
f_{1}(z,x)\mathrm{d}\mu(x)+(W\ast\rho)(z)$
$\displaystyle=c_{1}[\rho,\mu]\quad\text{ for all
}z\in\operatorname{supp}(\rho)\,,$ (14a) $\displaystyle\int
f_{1}(z,x)\mathrm{d}\rho(z)+V(x)$ $\displaystyle=c_{2}[\rho,\mu]\quad\text{
for all }x\in\operatorname{supp}(\mu)\,,$ (14b)

where $c_{1},c_{2}$ are constants that may differ on different connected
components of $\operatorname{supp}(\rho)$ and $\operatorname{supp}(\mu)$.
These equations are not easy to solve explicitly, and we are therefore using
general non-constructive techniques from calculus of variations. We first show
continuity and convexity properties for the functional $G_{a}$ (Lemma 9 and
Proposition 10), essential properties that will allow us to deduce existence
and uniqueness of ground states using the direct method in the calculus of
variations (Proposition 11). Using the Euler-Lagrange equation 14, we then
prove properties on the support of the ground state (Corollary 12). To obtain
convergence results, we apply the HWI method: we first show a general
’interpolation’ inequality between the energy, the energy dissipation and the
metric (Proposition 13); this fundamental inequality will then imply a
generalized logarithmic Sobolev inequality (Corollary 14) 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 ingrediends together will then allow us to conclude for the
statements in Theorem 2.

######  Lemma 9 (Lower semi-continuity).

Let Assumptions 1-3 hold. Then the functional
$G_{a}:\mathcal{P}\times\mathcal{P}\to\mathbb{R}$ is lower semi-continuous
with respect to the weak-$*$ topology.

###### Proof.

We split the energy $G_{a}$ into three parts: (i) $KL(\rho|\tilde{\rho})$,
(ii) the interaction energy $\mathcal{W}$, and (iii) the potential energy
$\mathcal{V}$. For (i), lower semi-continuity has been shown in [Pos75]. For
(ii), we can directly apply [San15, Proposition 7.2] using Assumption 3. For
(iii), note that $V$ and $f_{1}$ are lower semi-continuous and bounded below
thanks to Assumption 1, and so the result follows from [San15, Proposition
7.1]. ∎

######  Proposition 10 (Uniform displacement convexity).

Let $\alpha,\beta>0$. Fix
$\gamma_{0},\gamma_{1}\in\mathcal{P}_{2}\times\mathcal{P}_{2}$ and let
Assumptions 1-3 hold. Along any geodesic
$(\gamma_{s})_{s\in[0,1]}\in\mathcal{P}_{2}\times\mathcal{P}_{2}$ connecting
$\gamma_{0}$ to $\gamma_{1}$, we have for all $s\in[0,1]$

$\displaystyle\frac{\mathrm{d}^{2}}{\mathrm{d}s^{2}}G_{a}(\gamma_{s})\geq\lambda_{a}\overline{W}(\gamma_{0},\gamma_{1})^{2}\,,\qquad\lambda_{a}:=\lambda_{1}+\min(\lambda_{2}+\beta,\alpha\tilde{\lambda})\,.$
(15)

As a result, the functional $G_{a}:\mathcal{P}\times\mathcal{P}\to\mathbb{R}$
is uniformly displacement convex with constant $\lambda_{a}>0$.

###### Proof.

Let $\gamma_{0}$ and $\gamma_{1}$ be two probability measures with bounded
second moments. Denote by $\phi,\psi:\mathbb{R}^{d}\to\mathbb{R}$ the optimal
Kantorovich potentials pushing $\rho_{0}$ onto $\rho_{1}$, and $\mu_{0}$ onto
$\mu_{1}$, respectively:

$\displaystyle\rho_{1}=\nabla\phi_{\\#}\rho_{0}\quad\text{ such that }\quad
W_{2}(\rho_{0},\rho_{1})^{2}=\int_{\mathbb{R}^{d}}\|z-\nabla\phi(z)\|^{2}\mathrm{d}\rho_{0}(z)\,,$
$\displaystyle\mu_{1}=\nabla\psi_{\\#}\mu_{0}\quad\text{ such that }\quad
W_{2}(\mu_{0},\mu_{1})^{2}=\int_{\mathbb{R}^{d}}\|x-\nabla\psi(x)\|^{2}\mathrm{d}\mu_{0}(x)\,.$

The now classical results in [BB00] guarantee that there always exists convex
functions $\phi,\psi$ that satisfy the conditions above. Then the path
$(\gamma_{s})_{s\in[0,1]}=(\rho_{s},\mu_{s})_{s\in[0,1]}$ defined by

$\displaystyle\rho_{s}$
$\displaystyle=[(1-s)\operatorname{id}+s\nabla_{z}\phi]_{\\#}\rho_{0}\,,$
$\displaystyle\mu_{s}$
$\displaystyle=[(1-s)\operatorname{id}+s\nabla_{x}\psi]_{\\#}\mu_{0}$

is a $\overline{W}$-geodesic from $\gamma_{0}$ to $\gamma_{1}$.

The first derivative of $\mathcal{V}$ along geodesics in the Wasserstein
metric is given by

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}s}\mathcal{V}(\gamma_{s})=$
$\displaystyle\frac{\mathrm{d}}{\mathrm{d}s}\left[\iint
f_{1}((1-s)z+s\nabla\phi(z),(1-s)x+s\nabla\psi(x))\,\mathrm{d}\rho_{0}(z)\mathrm{d}\mu_{0}(x)\right.$
$\displaystyle\hskip 19.91684pt\left.+\int
V((1-s)x+s\nabla\psi(x))\,\mathrm{d}\mu_{0}(x)\right]$ $\displaystyle=$
$\displaystyle\iint\nabla_{x}f_{1}((1-s)z+s\nabla\phi(z),(1-s)x+s\nabla\psi(x))\cdot(\nabla\psi(x)-x)\,\mathrm{d}\rho_{0}(z)\mathrm{d}\mu_{0}(x)$
$\displaystyle\iint\nabla_{z}f_{1}((1-s)z+s\nabla\phi(z),(1-s)x+s\nabla\psi(x))\cdot(\nabla\phi(z)-z)\,\mathrm{d}\rho_{0}(z)\mathrm{d}\mu_{0}(x)$
$\displaystyle+\int\nabla_{x}V((1-s)x+s\nabla\psi(x))\cdot(\nabla\psi(x)-x)\,\mathrm{d}\mu_{0}(x)\,,$

and taking another derivative we have

$\displaystyle\frac{\mathrm{d}^{2}}{\mathrm{d}s^{2}}\mathcal{V}(\gamma_{s})=$
$\displaystyle-\iint\begin{bmatrix}(\nabla\psi(x)-x)\\\
(\nabla\phi(z)-z)\end{bmatrix}^{T}\cdot
D_{s}(z,x)\cdot\begin{bmatrix}(\nabla\psi(x)-x)\\\
(\nabla\phi(z)-z)\end{bmatrix}\,\mathrm{d}\rho_{0}(z)\mathrm{d}\mu_{0}(x)$
$\displaystyle+\iint(\nabla\psi(x)-x)^{T}\cdot\nabla^{2}_{x}V((1-s)x+s\nabla\psi(x))\cdot(\nabla\psi(x)-x)\,\mathrm{d}\rho_{0}(z)\mathrm{d}\mu_{0}(x)$
$\displaystyle\geq\lambda_{1}\overline{W}(\gamma_{0},\gamma_{1})^{2}+(\lambda_{2}+\beta)W_{2}(\mu_{0},\mu_{1})^{2}\,,$

where we denoted
$D_{s}(z,x):=\text{Hess}(f_{1})((1-s)z+s\nabla\phi(z),(1-s)x+s\nabla\psi(x))$,
and the last inequality follows from Assumption 1 and the optimality of the
potentials $\phi$ and $\psi$.

Following [CMV03, Vil03] and using Assumption 2, the second derivatives of the
diffusion term and the interaction term along geodesics are given by

$\displaystyle\frac{\mathrm{d}^{2}}{\mathrm{d}s^{2}}KL(\rho_{s}|\tilde{\rho})\geq\alpha\tilde{\lambda}\,W_{2}(\rho_{0},\rho_{1})^{2}\,,\qquad\frac{\mathrm{d}^{2}}{\mathrm{d}s^{2}}\mathcal{W}(\rho_{s})\geq
0.$ (16)

Putting the above estimates together, we obtain (15). ∎

######  Remark 5.

Alternatively, one could assume strong convexity of $W$, which would improve
the lower-bound on the second derivative along geodesics.

######  Proposition 11.

(Ground state) Let Assumptions 1-3 hold for $\alpha,\beta>0$. Then the
functional
$G_{a}:\mathcal{P}(\mathbb{R}^{d})\times\mathcal{P}(\mathbb{R}^{d})\to[0,\infty]$
admits a unique minimizer $\gamma_{*}=(\rho_{*},\mu_{*})$, and it satisfies
$\rho_{*}\in\mathcal{P}_{2}(\mathbb{R}^{d})\cap L^{1}(\mathbb{R}^{d})$,
$\mu_{*}\in\mathcal{P}_{2}(\mathbb{R}^{d})$, and $\rho_{*}$ is absolutely
continuous with respect to $\tilde{\rho}$.

###### Proof.

We show existence of a minimizer of $G_{a}$ using the direct method in the
calculus of variations. Denote by
$\gamma=(\rho,\mu)\in\mathcal{P}\times\mathcal{P}\subset\mathcal{M}\times\mathcal{M}$
a pair of probability measures as a point in the product space of Radon
measures. Since $G_{a}\geq 0$ on $\mathcal{P}\times\mathcal{P}$ (see
Assumption 1) and not identically $+\infty$ everywhere, there exists a
minimizing sequence $(\gamma_{n})\in\mathcal{P}\times\mathcal{P}$. Note that
$(\gamma_{n})$ is in the closed unit ball of the dual space of continuous
functions vanishing at infinity $(C_{0}(\mathbb{R}^{d})\times
C_{0}(\mathbb{R}^{d}))^{*}$ endowed with the dual norm
$\left\lVert\gamma_{n}\right\rVert_{*}=\sup\frac{|\int
f\mathrm{d}\rho_{n}+\int
g\mathrm{d}\mu_{n}|}{\left\lVert(f,g)\right\rVert_{\infty}}$ over $f,g\in
C_{0}(\mathbb{R}^{d})$ with
$\|(f,g)\|_{\infty}:=\|f\|_{\infty}+\|g\|_{\infty}\neq 0$. By the Banach-
Alaoglu theorem [Rud91, Thm 3.15] there exists a limit
$\gamma_{*}=(\rho_{*},\mu_{*})\in\mathcal{M}\times\mathcal{M}=(C_{0}\times
C_{0})^{*}$ and a convergent subsequence (not relabelled) such that
$\gamma_{n}\overset{\ast}{\rightharpoonup}\gamma_{*}$. In fact, since
$KL(\rho_{*}\,|\,\tilde{\rho})<\infty$ it follows that $\rho_{*}$ is
absolutely continuous with respect to $\tilde{\rho}$, implying $\rho_{*}\in
L^{1}(\mathbb{R}^{d})$ thanks to Assumption 2. Further, $\mu_{*}$ has bounded
second moment, else we would have
$\inf_{\gamma\in\mathcal{P}\times\mathcal{P}}G_{a}(\gamma)=\infty$ which
yields a contradiction. It remains to show that
$\int\mathrm{d}\rho_{*}=\int\mathrm{d}\mu_{*}=1$ to conclude that
$\gamma_{*}\in\mathcal{P}\times\mathcal{P}$. To this aim, it is sufficient to
show tightness of $(\rho_{n})$ and $(\mu_{n})$, preventing the escape of mass
to infinity as we have $\int\mathrm{d}\rho_{n}=\int\mathrm{d}\mu_{n}=1$ for
all $n\geq 1$. Tightness follows from Markov’s inequality [Gho02] 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

$\int\|z\|^{2}\mathrm{d}\rho_{n}(z)+\int\|x\|^{2}\mathrm{d}\mu_{n}(x)<C\qquad\forall
n\in\mathrm{N}\,.$ (17)

To establish (17), observe that thanks to Assumption 2, there exists a
constant $c_{0}\in\mathbb{R}$ (possibly negative) such that
$-\log\tilde{\rho}(z)\geq c_{0}+\frac{\tilde{\lambda}}{4}\|z\|^{2}$ for all
$z\in\mathbb{R}^{d}$. Then

$\displaystyle\frac{\alpha\tilde{\lambda}}{4}\int\|z\|^{2}\mathrm{d}\rho_{n}$
$\displaystyle\leq-\alpha
c_{0}-\alpha\int\log\tilde{\rho}(z)\mathrm{d}\rho_{n}$

Therefore, using $\int\mathrm{d}\rho_{n}=\int\mathrm{d}\mu_{n}=1$ and writing
$\zeta:=\min\\{\frac{\alpha\tilde{\lambda}}{4},\frac{\beta}{2}\\}>0$, we
obtain the desired uniform upper bound on the second moments of the minimizing
sequence,

$\displaystyle\zeta\iint\left(\|z\|^{2}+\|x\|^{2}\right)\mathrm{d}\rho_{n}\mathrm{d}\mu_{n}$
$\displaystyle\leq-\alpha
c_{0}-\alpha\int\log\tilde{\rho}(z)\mathrm{d}\rho_{n}+\beta\int\|x-x_{0}\|^{2}\mathrm{d}\mu_{n}+\beta\|x_{0}\|^{2}$
$\displaystyle\leq-\alpha c_{0}+\beta\|x_{0}\|^{2}+G_{a}(\gamma_{n})$
$\displaystyle\leq-\alpha
c_{0}+\beta\|x_{0}\|^{2}+G_{a}(\gamma_{1})<\infty\,.$

This concludes the proof that the limit $\gamma_{*}$ satisfies indeed
$\gamma_{*}\in\mathcal{P}\times\mathcal{P}$, and indeed
$\rho_{*}\in\mathcal{P}_{2}(\mathbb{R}^{d})$ as well. Finally, $\gamma_{*}$ is
indeed a minimizer of $G_{a}$ thanks to weak-* lower-semicontinuity of $G_{a}$
following Lemma 9.

Next we show uniqueness using a contradiction argument. Suppose
$\gamma_{*}=(\rho_{*},\mu_{*})$ and
$\gamma_{*}^{\prime}=(\rho_{*}^{\prime},\mu_{*}^{\prime})$ are minimizers of
$G_{a}$. For $s\in[0,1]$, define
$\gamma_{s}:=((1-s)\operatorname{id}+sT,(1-s)\operatorname{id}+sS)_{\\#}\gamma_{*}$,
where $T,S:\mathbb{R}^{d}\mapsto\mathbb{R}^{d}$ are the optimal transport maps
such that $\rho_{*}^{\prime}=T_{\\#}\rho_{*}$ and
$\mu_{*}^{\prime}=S_{\\#}\mu_{*}$. By Proposition 10 the energy $G_{a}$ is
uniformly displacement convex, and so we have

$\displaystyle
G_{a}(\gamma_{s})\leq(1-s)G_{a}(\gamma_{*})+sG_{a}(\gamma_{*}^{\prime})=G_{a}(\gamma_{*}).$

If $\gamma_{*}\neq\gamma_{*}^{\prime}$ and $s\in(0,1)$, then strict inequality
holds by applying similar arguments as in [McC97, Proposition 1.2]. However,
the strict inequality $G_{a}(\gamma_{s})<G_{a}(\gamma_{*})$ for
$\gamma_{*}\neq\gamma_{*}^{\prime}$ is a contradiction to the minimality of
$\gamma_{*}$. Hence, the minimizer is unique. ∎

######  Remark 6.

If $\lambda_{1}>0$, then the strict convexity of $f_{1}$ can be used to deduce
uniqueness, and the assumptions on $-\log\tilde{\rho}$ can be weakened from
strict convexity to convexity.

######  Corollary 12.

Let Assumptions 1-3 hold. Any minimizer $\gamma_{*}=(\rho_{*},\mu_{*})$ of
$G_{a}$ is a steady state for equation (4) according to Definition 8 and
satisfies $\operatorname{supp}(\rho_{*})=\operatorname{supp}(\tilde{\rho})$.

###### Proof.

By Proposition 11, we have $\rho_{*},\mu_{*}\in\mathcal{P}_{2}$, as well as
$\rho_{*}\in L^{1}_{+}$, $\|\rho_{*}\|_{1}=1$, and that $\rho_{*}$ is
absolutely continuous with respect to $\tilde{\rho}$. Since $W\in
C^{2}(\mathbb{R}^{d})$, it follows that $\nabla W\ast\rho_{*}\in L_{loc}^{1}$.
In order to show that $\gamma_{*}$ is a steady state for equation (4), it
remains to prove that $\rho_{*}\in W^{1,2}_{loc}\cap L^{\infty}_{loc}$. As
$\gamma^{*}$ is a minimizer, it is in particular a critical point, and
therefore satisfies equations (14). Rearranging, we obtain (for a possible
different constant $c_{1}[\rho_{*},\mu_{*}]\neq 0$) from (14a) that

$\displaystyle\rho_{*}(z)=c_{1}[\rho_{*},\mu_{*}]\tilde{\rho}(z)\exp{\left[-\frac{1}{\alpha}\left(\int
f_{1}(z,x)\,\mu_{*}(x)+W\ast\rho_{*}(z)\right)\right]}\qquad\text{ on
}\operatorname{supp}(\rho_{*})\,.$ (18)

Then for any compact set $K\subset\mathbb{R}^{d}$,

$\displaystyle\sup_{z\in K}\rho_{*}(z)\leq c_{1}[\rho_{*},\mu_{*}]\sup_{z\in
K}\tilde{\rho}(z)\sup_{z\in K}\exp\left(-\frac{1}{\alpha}\left(\int
f_{1}(z,x)\,\mu_{*}(x)\right)\right)\sup_{z\in
K}\exp\left(-\frac{1}{\alpha}W\ast\rho_{*}\right).$

As $f_{1}\geq 0$ by Assumption 1 and $W\geq 0$ by Assumption 3, the last two
terms are finite. The first supremum is finite thanks to continuity of
$\tilde{\rho}$. Therefore $\rho_{*}\in L_{loc}^{\infty}$. To show that
$\rho_{*}\in W_{loc}^{1,2}$, note that for any compact set
$K\subset\mathbb{R}^{d}$, we have
$\int_{K}|\rho_{*}(z)|^{2}\mathrm{d}z<\infty$ as a consequence of $\rho_{*}\in
L_{loc}^{\infty}$. Moreover, defining
$T[\gamma](z)\coloneqq-\frac{1}{\alpha}\left(\int
f_{1}(z,x)\,\mu(x)+W\ast\rho(z)\right)\leq 0$, we have

$\displaystyle\int_{K}|\nabla\rho_{*}|^{2}\mathrm{d}z=c_{1}[\rho_{*},\mu_{*}]^{2}\int_{K}|\nabla\tilde{\rho}+\tilde{\rho}\nabla
T[\gamma_{*}]|^{2}\exp(2T[\gamma_{*}])\mathrm{d}z$ $\displaystyle\quad\leq
2c_{1}[\rho_{*},\mu_{*}]^{2}\int_{K}|\nabla\tilde{\rho}|^{2}\exp(2T[\gamma_{*}])\mathrm{d}z+2c_{1}[\rho_{*},\mu_{*}]^{2}\int_{K}|\nabla
T[\gamma_{*}]|^{2}\tilde{\rho}^{2}\exp(2T[\gamma_{*}])\mathrm{d}z\,,$

which is bounded noting that $\exp(2T[\gamma_{*}])\leq 1$ and that
$T[\gamma_{*}](\cdot)$, $\nabla T[\gamma_{*}](\cdot)$ and $\nabla\tilde{\rho}$
are in $L^{\infty}_{loc}$, where we used that
$f_{1},(\cdot,x),W(\cdot),\tilde{\rho}(\cdot)\in C^{1}(\mathbb{R}^{d})$ by
Assumptions 1-3. We conclude that $\rho_{*}\in W^{1,2}_{loc}$, and indeed
$(\rho_{*},\mu_{*})$ solves (12) in the sense of distributions as a
consequence of (14).

Next, we show that
$\operatorname{supp}(\rho_{*})=\operatorname{supp}(\tilde{\rho})$ using again
the relation (18). Firstly, note that
$\operatorname{supp}(\rho_{*})\subset\operatorname{supp}(\tilde{\rho})$ since
$\rho_{*}$ is absolutely continuous with respect to $\tilde{\rho}$. Secondly,
we claim that $\exp{\left[-\frac{1}{\alpha}\left(\int
f_{1}(z,x)\,\mu_{*}(x)+W\ast\rho_{*}(z)\right)\right]}>0$ for all
$z\in\mathbb{R}^{d}$. In other words, we claim that $\int
f_{1}(z,x)\,\mu_{*}(x)<\infty$ and $W\ast\rho_{*}(z)<\infty$ for all
$z\in\mathbb{R}^{d}$. Indeed, for the first term, fix any $z\in\mathbb{R}^{d}$
and choose $R>0$ large enough such that $z\in B_{R}(0)$. Then, thanks to
continuity of $f_{1}$ according to Assumption 1, we have

$\int f_{1}(z,x)\,\mu_{*}(x)\leq\sup_{z\in B_{R}(0)}\int
f_{1}(z,x)\,\mu_{*}(x)<\infty\,.$

For the second term, note that by Assumption 3, we have for any
$z\in\mathbb{R}^{d}$ and $\epsilon>0$,

$\displaystyle W(z)$ $\displaystyle\leq W(0)+\nabla W(z)\cdot z\leq
W(0)+\frac{1}{2\epsilon}\|\nabla W(z)\|^{2}+\frac{\epsilon}{2}\|z\|^{2}$
$\displaystyle\leq
W(0)+\frac{D^{2}}{2\epsilon}\left(1+\|z\|\right)^{2}+\frac{\epsilon}{2}\|z\|^{2}\leq
W(0)+\frac{D^{2}}{\epsilon}+\left(\frac{D^{2}}{\epsilon}+\frac{\epsilon}{2}\right)\|z\|^{2}$
$\displaystyle=W(0)+\frac{D}{\sqrt{2}}+\sqrt{2}D\|z\|^{2}\,,$

where the last equality follows by choosing the optimal $\epsilon=\sqrt{2}D$.
We conclude that

$\displaystyle W\ast\rho_{*}(z)$ $\displaystyle\leq
W(0)+\frac{D}{\sqrt{2}}+\sqrt{2}D\int\|z-\tilde{z}\|^{2}\,\rho_{*}(\tilde{z})$
$\displaystyle\leq
W(0)+\frac{D}{\sqrt{2}}+2\sqrt{2}D\|z\|^{2}+2\sqrt{2}D\int\|\tilde{z}\|^{2}\,\rho_{*}(\tilde{z})\,,$
(19)

which is finite for any fixed $z\in\mathbb{R}^{d}$ thanks to the fact that
$\rho_{*}\in\mathcal{P}_{2}(\mathbb{R}^{d})$. Hence,
$\operatorname{supp}(\rho_{*})=\operatorname{supp}(\tilde{\rho})$. ∎

######  Remark 7.

If we have in addition that $\tilde{\rho}\in L^{\infty}(\mathbb{R}^{d})$, then
the minimizer $\rho_{*}$ of $G_{a}$ is in $L^{\infty}(\mathbb{R}^{d})$ as
well. This follows directly by bounding the right-hand side of (18).

The following inequality is referred to as HWI inequality and represents the
key result to obtain convergence to equilibrium.

######  Proposition 13 (HWI inequality).

Define the dissipation functional

$\displaystyle
D_{a}(\gamma):=\iint|\nabla_{x,z}\delta_{\gamma}G_{a}(z,x)|^{2}\mathrm{d}\gamma(z,x)\,.$

Assume $\alpha,\beta>0$ and let $\lambda_{a}$ as defined in (15). Let
$\gamma_{0},\gamma_{1}\in\mathcal{P}_{2}\times\mathcal{P}_{2}$ such that
$G_{a}(\gamma_{0}),G_{a}(\gamma_{1}),D_{a}(\gamma_{0})<\infty$. Then

$\displaystyle
G_{a}(\gamma_{0})-G_{a}(\gamma_{1})\leq\overline{W}(\gamma_{0},\gamma_{1})\sqrt{D_{a}(\gamma_{0})}-\frac{\lambda_{a}}{2}\,\overline{W}(\gamma_{0},\gamma_{1})^{2}$
(20)

###### Proof.

For simplicity, consider $\gamma_{0},\gamma_{1}$ that have smooth Lebesgue
densities of compact support. The general case can be recovered using
approximation arguments. Let $(\gamma_{s})_{s\in[0,1]}$ denote a
$\overline{W}$-geodesic between $\gamma_{0},\gamma_{1}$. Following similar
arguments as in [CMV03] and [OV00, Section 5] and making use of the
calculations in the proof of Proposition 10, we have

$\displaystyle\left.\frac{\mathrm{d}}{\mathrm{d}s}G_{a}(\gamma_{s})\right|_{s=0}$
$\displaystyle\geq\iint\begin{bmatrix}\xi_{1}(z)\\\
\xi_{2}(x)\end{bmatrix}\cdot\begin{bmatrix}(\nabla\phi(z)-z)\\\
(\nabla\psi(x)-x)\end{bmatrix}\,\mathrm{d}\gamma_{0}(z,x)\,,$

where

$\displaystyle\xi_{1}[\gamma_{0}](z)$
$\displaystyle:=\alpha\nabla_{z}\log\left(\frac{\rho_{0}(z)}{\tilde{\rho}(z)}\right)+\int\nabla_{z}f_{1}(z,x)\mathrm{d}\mu_{0}(x)+\int\nabla_{z}W(z-z^{\prime})\mathrm{d}\rho_{0}(z^{\prime})\,,$
$\displaystyle\xi_{2}[\gamma_{0}](x)$
$\displaystyle:=\int\nabla_{x}f_{1}(z,x)\mathrm{d}\rho_{0}(z)+\nabla_{x}V(x)\,.$

Note that the dissipation functional can then be written as

$D_{a}(\gamma_{0})=\iint\left(|\xi_{1}(z)|^{2}+|\xi_{2}(x)|^{2}\right)\mathrm{d}\gamma_{0}(z,x)\,.$

Using the double integral Cauchy-Schwarz inequality [Ste04], we obtain

$\displaystyle\left.\frac{\mathrm{d}}{\mathrm{d}s}G_{a}(\gamma_{s})\right|_{s=0}$
$\displaystyle\geq-\left(\sqrt{\iint\left\lVert\begin{bmatrix}\xi_{1}\\\
\xi_{2}\end{bmatrix}\right\rVert_{2}^{2}\mathrm{d}\gamma_{0}}\right)\left(\sqrt{\iint\left\lVert\begin{bmatrix}\nabla\phi(z)-z\\\
\nabla\psi(x)-x\end{bmatrix}\right\rVert_{2}^{2}\mathrm{d}\gamma_{0}}\right)$
$\displaystyle=-\sqrt{D_{a}(\gamma_{0})}\,\sqrt{\int\|\nabla\phi(z)-z\|^{2}\mathrm{d}\rho_{0}+\int\|\nabla\psi(x)-x\|^{2}\mathrm{d}\mu_{0}}$
$\displaystyle=-\sqrt{D_{a}(\gamma_{0})}\,\overline{W}(\gamma_{0},\gamma_{1})\,.$

Next, we compute a Taylor expansion of $G_{a}(\gamma_{s})$ when considered as
a function in $s$ and use the bound on
$\frac{\mathrm{d}^{2}}{\mathrm{d}s^{2}}G_{a}$ from (15):

$\displaystyle G_{a}(\gamma_{1})$
$\displaystyle=G_{a}(\gamma_{0})+\left.\frac{\mathrm{d}}{\mathrm{d}s}G_{a}(\gamma_{s})\right|_{s=0}+\int_{0}^{1}(1-t)\left.\left(\frac{\mathrm{d}^{2}}{\mathrm{d}s^{2}}G_{a}(\gamma_{s})\right)\right|_{s=t}\,\mathrm{d}t$
$\displaystyle\geq
G_{a}(\gamma_{0})-\sqrt{D_{a}(\gamma_{0})}\,\overline{W}(\gamma_{0},\gamma_{1})+\frac{\lambda_{a}}{2}\overline{W}(\gamma_{0},\gamma_{1})^{2}\,.$

∎

######  Remark 8.

The HWI inequality in Proposition 13 immediately implies uniqueness of
minimizers for $G_{a}$ in the set
$\left\\{\gamma\in\mathcal{P}\times\mathcal{P}\,:\,D_{a}(\gamma)<+\infty\right\\}$.
Indeed, if $\gamma_{0}$ is such that $D_{a}(\gamma_{0})=0$, then for any other
$\gamma_{1}$ in the above set we have $G_{a}(\gamma_{0})\leq
G_{a}(\gamma_{1})$ with equality if and only if
$\overline{W}(\gamma_{0},\gamma_{1})=0$.

######  Corollary 14 (Generalized Log-Sobolev inequality).

Denote by $\gamma_{*}$ the unique minimizer of $G_{a}$. With $\lambda_{a}$ as
defined in (15), any product measure
$\gamma\in\mathcal{P}_{2}\times\mathcal{P}_{2}$ such that
$G(\gamma),D_{a}(\gamma)<\infty$ satisfies

$D_{a}(\gamma)\geq 2\lambda_{a}\,G_{a}(\gamma|\gamma_{*})\,.$ (21)

###### Proof.

This statement follows immediately from Proposition 13. Indeed, let
$\gamma_{1}=\gamma_{*}$ and $\gamma_{0}=\gamma$ in (20). Then

$\displaystyle G_{a}(\gamma\,|\,\gamma_{*})$
$\displaystyle\leq\overline{W}(\gamma,\gamma_{*})\sqrt{D_{a}(\gamma)}-\frac{\lambda_{a}}{2}\,\overline{W}(\gamma,\gamma_{*})^{2}$
$\displaystyle\leq\max_{t\geq
0}\left(\sqrt{D_{a}(\gamma)}t-\frac{\lambda_{a}}{2}\,t^{2}\right)=\frac{D_{a}(\gamma)}{2\lambda_{a}}\,.$

∎

######  Corollary 15 (Talagrand inequality).

Denote by $\gamma_{*}$ the unique minimizer of $G_{a}$. With $\lambda_{a}$ as
defined in (15), it holds

$\overline{W}(\gamma,\gamma_{*})^{2}\leq\frac{2}{\lambda_{a}}G_{a}(\gamma\,|\,\gamma_{*})$

for any $\gamma\in\mathcal{P}_{2}\times\mathcal{P}_{2}$ such that
$G_{a}(\gamma)<\infty$.

###### Proof.

This is also a direct consequence of Proposition 13 by setting
$\gamma_{0}=\gamma_{*}$ and $\gamma_{1}=\gamma$. Then
$G_{a}(\gamma_{*})<\infty$ and $D_{a}(\gamma_{*})=0$, and the result follows.
∎

###### Proof of Theorem 2.

The entropy term $\int\rho\log\rho$ produces diffusion in $\rho$ for the
corresponding PDE in (4). As a consequence, solutions $\rho_{t}$ to (4) and
minimizers $\rho^{*}$ for $G_{a}$ have to be $L^{1}$ functions. As there is no
diffusion for the evolution of $\mu_{t}$, solutions may have a singular part.
In fact, for initial condition $\mu_{0}=\delta_{x_{0}}$, the corresponding
solution will be of the form $\mu_{t}=\delta_{x(t)}$, where $x(t)$ solves the
ODE (3) with initial condition $x_{0}$. This follows from the fact that the
evolution for $\mu_{t}$ is a transport equation (also see Section A.1 for more
details). Results (a) and (b) are the statements in Proposition 11, Corollary
12 and Corollary 15. To obtain (c), we differentiate the energy $G_{a}$ along
solutions $\gamma_{t}$ to the equation (4):

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}G_{a}(\gamma_{t})$
$\displaystyle=\int\delta_{\rho}G_{a}[\gamma_{t}](z)\partial_{t}\rho_{t}\mathrm{d}z+\int\delta_{\mu}G_{a}[\gamma_{t}](x)\partial_{t}\mu_{t}\mathrm{d}x$
$\displaystyle=-\int\left\|\nabla_{z}\delta_{\rho}G_{a}[\gamma_{t}](z)\right\|^{2}\mathrm{d}\rho_{t}(z)-\int\left\|\nabla_{x}\delta_{\mu}G_{a}[\gamma_{t}](x)\right\|^{2}\mathrm{d}\mu_{t}(x)$
$\displaystyle=-D_{a}(\gamma_{t})\leq-2\lambda_{a}G_{a}(\gamma_{t}\,|\,\gamma_{*})\,,$

where the last bound follows from Corollary 14. Applying Gronwall’s
inequality, we immediately obtain decay in energy,

$G_{a}(\gamma_{t}\,|\,\gamma_{*})\leq
e^{-2\lambda_{a}t}G_{a}(\gamma_{0}\,|\,\gamma_{*})\,.$

Finally, applying Talagrand’s inequality (Corollary 15), the decay in energy
implies decay in the product Wasserstein metric,

$\overline{W}(\gamma_{t},\gamma_{*})\leq ce^{-\lambda_{a}t}$

where $c>0$ is a constant only depending on $\gamma_{0}$, $\gamma_{*}$ and the
parameter $\lambda_{a}$. ∎

##  Appendix C Proof of Theorem 3

In the case of competing objectives, we rewrite the energy
$G_{c}(\rho,x):\mathcal{P}(\mathbb{R}^{d})\times\mathbb{R}^{d}\mapsto[-\infty,\infty]$
as follows:

$\displaystyle G_{c}(\rho,x)$ $\displaystyle=\int
f_{1}(z,x)\mathrm{d}\rho(z)+\int
f_{2}(z,x)\mathrm{d}\bar{\rho}(z)+\frac{\beta}{2}\left\lVert
x-x_{0}\right\rVert^{2}-P(\rho)\,,$

where

$\displaystyle P(\rho)$ $\displaystyle:=\alpha
KL(\rho|\tilde{\rho})+\frac{1}{2}\int\rho W\ast\rho\,.$

Note that for any fixed $\rho\in\mathcal{P}$, the energy $G_{c}(\rho,\cdot)$
is strictly convex in $x$, and therefore has a unique minimizer. Define the
best response by

$b(\rho)\coloneqq\operatorname*{argmin}_{\bar{x}}G_{c}(\rho,\bar{x})$

and denote $G_{b}(\rho)\coloneqq G_{c}(\rho,b(\rho))$. We begin with auxiliary
results computing the first variations of the best response $b$ and then the
different terms in $G_{b}(\rho)$ using Definition 1.

######  Lemma 16 (First variation of the best response).

The first variation of the best response of the classifier at $\rho$ (if it
exists) is

$\displaystyle\delta_{\rho}b[\rho](z)=-Q(\rho)^{-1}\nabla_{x}f_{1}(z,b(\rho))\quad\text{
for almost every }z\in\mathbb{R}^{d}\,,$

where $Q(\rho)\succeq(\beta+\lambda_{1}+\lambda_{2})\operatorname{I}_{d}$ is a
symmetric matrix, constant in $z$ and $x$, defined as

$Q(\rho)\coloneqq\beta\operatorname{I}_{d}+\int\nabla_{x}^{2}f_{1}(z,b(\rho))\mathrm{d}\rho(z)+\int\nabla_{x}^{2}f_{2}(z,b(\rho))\mathrm{d}\bar{\rho}(z)\,.$

In particular, we then have for any $\psi\in C_{c}^{\infty}(\mathbb{R}^{d})$
with $\int\psi\,\mathrm{d}z=0$ that

$\lim_{\epsilon\to
0}\frac{1}{\epsilon}\left\|b[\rho+\epsilon\psi]-b[\rho]-\epsilon\int\delta_{\rho}b[\rho](z)\psi(z)\mathrm{d}z\right\|=0\,.$

###### Proof.

Let $\psi\in C_{c}^{\infty}(\mathbb{R}^{d})$ with $\int\psi\,\mathrm{d}z=0$
and fix $\epsilon>0$. Any minimizer of $G_{c}(\rho+\epsilon\psi,x)$ for fixed
$\rho$ must satisfy

$\displaystyle\nabla_{x}G_{c}(\rho+\epsilon\psi,b(\rho+\epsilon\psi))=0\,.$

Differentiating in $\epsilon$, we obtain

$\displaystyle\int\delta_{\rho}\nabla_{x}G_{c}[\rho+\epsilon\psi,b(\rho+\epsilon\psi)]\psi(z)\,\mathrm{d}z+\nabla_{x}^{2}G_{c}(\rho+\epsilon\psi,b(\rho+\epsilon\psi))\int\delta_{\rho}b[\rho+\epsilon\psi](z)\psi(z)\,\mathrm{d}z=0\,.$
(22)

Next, we explicitly compute all terms involved in (22). Computing the
derivatives yields

$\displaystyle\nabla_{x}G_{c}(\rho,x)$
$\displaystyle=\int\nabla_{x}f_{1}(z,x)\mathrm{d}\rho(z)+\int\nabla_{x}f_{2}(z,x)\mathrm{d}\bar{\rho}(z)+\beta(x-x_{0})$
$\displaystyle\delta_{\rho}\nabla_{x}G_{c}[\rho,x](z)$
$\displaystyle=\nabla_{x}f_{1}(z,x)$
$\displaystyle\nabla_{x}^{2}G_{c}(\rho,x)$
$\displaystyle=\int\nabla_{x}^{2}f_{1}(z,x)\mathrm{d}\rho(z)+\int\nabla_{x}^{2}f_{2}(z,x)\mathrm{d}\bar{\rho}(z)+\beta\operatorname{I}_{d}.$

Note that $\nabla_{x}^{2}G_{c}$ is invertible by Assumption 1, which states
that $f_{1}$ and $f_{2}$ have positive-definite Hessians. Inverting this term
and substituting these expressions into (22) for $\epsilon=0$ gives

$\displaystyle\int\delta_{\rho}b[\rho](z)\psi(z)\,\mathrm{d}z$
$\displaystyle=-\left[\beta\operatorname{I}_{d}+\int\nabla_{x}^{2}f_{1}(z,b(\rho))\mathrm{d}\rho(z)+\int\nabla_{x}^{2}f_{2}(z,b(\rho))\mathrm{d}\bar{\rho}(z)\right]^{-1}\int\nabla_{x}f_{1}(z,b(\rho))\psi(z)\,\mathrm{d}z$
$\displaystyle=-\int
Q(\rho)^{-1}\nabla_{x}f_{1}(z,b(\rho))\psi(z)\,\mathrm{d}z\,.$

Finally, the lower bound on $Q(\rho)$ follows thanks to Assumption 1. ∎

######  Remark 9.

If we include the additional assumption that $f_{i}\in
C^{3}(\mathbb{R}^{d}\times\mathbb{R}^{d};[0,\infty))$ for $i=1,2$, then the
Hessian of $b[\rho]$ is well-defined. More precisely, the Hessian is given by

$\displaystyle\frac{\mathrm{d}^{2}}{\mathrm{d}\epsilon^{2}}b[\rho+\epsilon\psi]|_{\epsilon=0}=Q(\rho)^{-1}\left(\frac{\mathrm{d}}{\mathrm{d}\epsilon}Q(\rho+\epsilon\psi)|_{\epsilon=0}+\int\nabla_{x}^{2}f_{1}(z,b[\rho])\psi(z)\mathrm{d}z\right)Q(\rho)^{-1}u[\rho,\psi]$

where $u[\rho,\psi]=\int\nabla_{x}f_{1}(z,b[\rho])\psi(z)\mathrm{d}z$ and

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}\epsilon}Q_{ij}(\rho+\epsilon\psi)|_{\epsilon=0}$
$\displaystyle=\int\partial_{x_{i}}\partial_{x_{j}}f_{1}(z,b[\rho])\psi(z)\mathrm{d}z-\int\partial_{x_{i}}\partial_{x_{j}}\nabla_{x}f_{1}(z,b[\rho])\psi(z)\rho(z)\mathrm{d}z\,Q(\rho)^{-1}u[\rho,\psi]$
$\displaystyle\quad-\int\partial_{x_{i}}\partial_{x_{j}}\nabla_{x}f_{2}(z,b[\rho])\psi(z)\bar{\rho}(z)\mathrm{d}z\,Q(\rho)^{-1}u[\rho,\psi].$

Therefore, we can Taylor expand $b[\rho]$ up to second order and control the
remainder term of order $\epsilon^{2}$.

######  Lemma 17 (First variation of $G_{b}$).

The first variation of $G_{b}$ is given by

$\delta_{\rho}G_{b}[\rho](z)=h_{1}(z)+h_{2}(z)+\beta
h_{3}(z)-\delta_{\rho}P[\rho](z)\,,$

where

$\displaystyle h_{1}(z):=\frac{\delta}{\delta\rho}\left(\int
f_{1}(\tilde{z},b(\rho))\mathrm{d}\rho(\tilde{z})\right)(z)=\left\langle\int\nabla_{x}f_{1}(\tilde{z},\,b(\rho))\mathrm{d}\rho(\tilde{z}),\frac{\delta
b}{\delta\rho}[\rho](z)\right\rangle+f_{1}(z,b(\rho))\,,$ $\displaystyle
h_{2}(z):=\frac{\delta}{\delta\rho}\left(\int
f_{2}(\tilde{z},b(\rho))\mathrm{d}\bar{\rho}(\tilde{z})\right)(z)=\left\langle\int\nabla_{x}f_{2}(\tilde{z},\,b(\rho))\mathrm{d}\bar{\rho}(\tilde{z}),\frac{\delta
b}{\delta\rho}[\rho](z)\right\rangle\,,$ $\displaystyle
h_{3}(z):=\frac{1}{2}\frac{\delta}{\delta\rho}\|b(\rho)-x_{0}\|^{2}=\left\langle
b(\rho)-x_{0},\,\frac{\delta b}{\delta\rho}[\rho](z)\right\rangle\,,$

and

$\displaystyle\delta_{\rho}P[\rho](z)=\alpha\log(\rho(z)/\tilde{\rho}(z))+(W\ast\rho)(z)\,.$

###### Proof.

We begin with general expressions for Taylor expansions of
$b:\mathcal{P}(\mathbb{R}^{d})\to\mathbb{R}^{d}$ and
$f_{i}(z,b(\cdot)):\mathcal{P}(\mathbb{R}^{d})\to\mathbb{R}$ for $i=1,2$
around $\rho$. Let $\psi\in\mathcal{T}$ with
$\mathcal{T}=\\{\psi:\int\psi(z)\mathrm{d}z=0\\}$. Then

$\displaystyle b(\rho+\epsilon\psi)=b(\rho)+\epsilon\int\frac{\delta
b}{\delta\rho}[\rho](z^{\prime})\psi(z^{\prime})\mathrm{d}z^{\prime}+O(\epsilon^{2})$
(23)

and

$\displaystyle
f_{i}(z,b(\rho+\epsilon\psi))=f_{i}(z,b(\rho))+\epsilon\left\langle\nabla_{x}f_{i}(z,\,b(\rho)),\int\frac{\delta
b}{\delta\rho}[\rho](z^{\prime})\psi(z^{\prime})\mathrm{d}z^{\prime}\right\rangle+O(\epsilon^{2})\,.$
(24)

We compute explicitly each of the first variations:

  1. (i)

Using (24), we have

$\displaystyle\int\psi(z)h_{1}(z)\mathrm{d}z$
$\displaystyle=\lim_{\epsilon\rightarrow 0}\frac{1}{\epsilon}\bigg{[}\int
f_{1}(z,b(\rho+\epsilon\psi))(\rho(z)+\epsilon\psi(z))\mathrm{d}z-\int
f_{1}(z,b(\rho))\rho(z)\mathrm{d}z\bigg{]}$
$\displaystyle=\left\langle\int\nabla_{x}f_{1}(z,\,b(\rho))\mathrm{d}\rho(z),\int\frac{\delta
b(\rho)}{\delta\rho}[\rho](z^{\prime})\psi(z^{\prime})\mathrm{d}z^{\prime}\right\rangle+\int
f_{1}(z,b(\rho))\psi(z)\mathrm{d}z$
$\displaystyle=\int\left\langle\int\nabla_{x}f_{1}(z,\,b(\rho))\mathrm{d}\rho(z),\frac{\delta
b(\rho)}{\delta\rho}[\rho](z^{\prime})\right\rangle\psi(z^{\prime})\mathrm{d}z^{\prime}+\int
f_{1}(z,b(\rho))\psi(z)\mathrm{d}z$ $\displaystyle\Rightarrow h_{1}(z)$
$\displaystyle=\left\langle\int\nabla_{x}f_{1}(\tilde{z},\,b(\rho))\mathrm{d}\rho(\tilde{z}),\frac{\delta
b}{\delta\rho}[\rho](z)\right\rangle+f_{1}(z,b(\rho))\,.$

  2. (ii)

Similarly, using again (24),

$\displaystyle\int\psi(z)h_{2}(z)\mathrm{d}z$
$\displaystyle=\lim_{\epsilon\rightarrow 0}\frac{1}{\epsilon}\bigg{[}\int
f_{2}(z,b(\rho+\epsilon\psi))\mathrm{d}\bar{\rho}(z)-\int
f_{2}(z,b(\rho))\bar{\rho}(z)\mathrm{d}z\bigg{]}$
$\displaystyle=\int\left\langle\int\nabla_{x}f_{2}(\tilde{z},\,b(\rho))\mathrm{d}\bar{\rho}(\tilde{z}),\frac{\delta
b}{\delta\rho}[\rho](z)\right\rangle\psi(z)\mathrm{d}z$
$\displaystyle\Rightarrow h_{2}(z)$
$\displaystyle=\left\langle\int\nabla_{x}f_{2}(\tilde{z},\,b(\rho))\mathrm{d}\bar{\rho}(\tilde{z}),\frac{\delta
b}{\delta\rho}[\rho](z)\right\rangle\,.$

  3. (iii)

Finally, from (23) it follows that

$\displaystyle\int\psi(z)h_{3}(z)\mathrm{d}z$
$\displaystyle=\lim_{\epsilon\rightarrow
0}\frac{1}{2\epsilon}\bigg{[}\left\langle
b(\rho+\epsilon\psi)-x_{0},\,b(\rho+\epsilon\psi)-x_{0}\right\rangle-\left\langle
b(\rho)-x_{0},\,b(\rho)-x_{0}\right\rangle\bigg{]}$
$\displaystyle=\int\left\langle b(\rho)-x_{0},\,\frac{\delta
b}{\delta\rho}[\rho](z)\right\rangle\psi(z)\mathrm{d}z$
$\displaystyle\Rightarrow h_{3}(z)$ $\displaystyle=\left\langle
b(\rho)-x_{0},\,\frac{\delta b}{\delta\rho}[\rho](z)\right\rangle\,.$

Finally, the expression for $\delta_{\rho}P[\rho]$ follows by direct
computation ∎

######  Lemma 18.

Denote $G_{b}(\rho)\coloneqq G_{c}(\rho,b(\rho))$ with $b(\rho)$ given by (5).
Then
$\delta_{\rho}G_{b}[\rho]=\left.\delta_{\rho}G_{c}[\rho]\right|_{x=b(\rho)}$.

###### Proof.

We start by computing $\delta_{\rho}G_{c}(\cdot,x)[\rho](z)$ for any
$z,x\in\mathbb{R}^{d}$:

$\displaystyle\delta_{\rho}G_{c}(\cdot,x)[\rho](z)$
$\displaystyle=f_{1}(z,x)-\delta_{\rho}P[\rho](z).$ (25)

Next, we compute $\delta_{\rho}G_{b}$. Using Lemma 17, the first variation of
$G_{b}$ is given by

$\displaystyle\delta_{\rho}G_{b}[\rho](z)$
$\displaystyle=h_{1}(z)+h_{2}(z)+\beta h_{3}(z)-\delta_{\rho}P[\rho](z)$
$\displaystyle=-\left\langle\left[\int\nabla_{x}f_{1}(\tilde{z},\,b(\rho))\mathrm{d}\rho(\tilde{z})+\int\nabla_{x}f_{2}(\tilde{z},b(\rho))\mathrm{d}\bar{\rho}(\tilde{z})+\beta(b(\rho)-x_{0})\right],\delta_{\rho}b[\rho](z)\right\rangle$
$\displaystyle\quad+f_{1}(z,b(\rho))-\delta_{\rho}P[\rho](z)\,.$

Note that

$\displaystyle\nabla_{x}G_{c}(\rho,x)=\int\nabla_{x}f_{1}(\tilde{z},x)\mathrm{d}\rho(\tilde{z})+\int\nabla_{x}f_{2}(\tilde{z},x)\mathrm{d}\bar{\rho}(\tilde{z})+\beta(x-x_{0})\,,$
(26)

and by the definition of the best response $b(\rho)$, we have
$\nabla_{x}G_{x}(\rho,x)|_{x=b(\rho)}=0$. Substituting into the expression for
$\delta_{\rho}G_{b}$ and using (25), we obtain

$\displaystyle\delta_{\rho}G_{b}[\rho](z)=f_{1}(z,b(\rho))-\delta_{\rho}P[\rho](z)=\delta_{\rho}G_{c}(\cdot,x)[\rho](z)\bigg{|}_{x=b(\rho)}\,.$

This concludes the proof. ∎

######  Lemma 19 (Uniform boundedness of the best response).

Let Assumption 1 hold. Then for any $\rho\in\mathcal{P}(\mathbb{R}^{d})$, we
have

$\left\lVert b(\rho)\right\rVert^{2}\leq\left\lVert
x_{0}\right\rVert^{2}+\frac{2(a_{1}+a_{2})}{\beta}\,.$

###### Proof.

By definition of the best response $b(\rho)$, we have

$\displaystyle\int\nabla_{x}f_{1}(z,b(\rho))\mathrm{d}\rho_{t}+\int\nabla_{x}f_{2}(z,b(\rho))\mathrm{d}\bar{\rho}(z)+\beta(b(\rho)-x_{0})=0\,.$

To show that that $b(\rho)$ is uniformly bounded, we take the inner product of
the above expression with $b(\rho)$ itself

$\displaystyle\beta\|b(\rho)\|^{2}=\beta x_{0}\cdot
b(\rho)-\int\nabla_{x}f_{1}(z,b(\rho))\cdot
b(\rho)\mathrm{d}\rho(z)-\int\nabla_{x}f_{2}(z,b(\rho))\cdot
b(\rho)\mathrm{d}\bar{\rho}(z)\,.$

Using Assumption 1 to bound the two integrals, together with using Young’s
inequality to bound the first term on the right-hand side, we obtain

$\displaystyle\beta\|b(\rho)\|^{2}\leq\frac{\beta}{2}\left\lVert
x_{0}\right\rVert^{2}+\frac{\beta}{2}\left\lVert
b(\rho)\right\rVert+a_{1}+a_{2}\,,$

which concludes the proof after rearranging terms. ∎

######  Lemma 20 (Upper semi-continuity).

Let Assumptions 1-3 hold. The functional
$G_{c}:\mathcal{P}(\mathbb{R}^{d})\times\mathbb{R}^{d}\to[-\infty,+\infty]$ is
upper semi-continuous when $\mathcal{P}(\mathbb{R}^{d})\times\mathbb{R}^{d}$
is endowed with the product topology of the weak-$*$ topology and the
Euclidean topology. Moreover, the functional
$G_{b}:\mathcal{P}(\mathbb{R}^{d})\to[-\infty,+\infty]$ is upper semi-
continuous with respect to the weak-$*$ topology.

###### Proof.

The functional
$G_{c}:\mathcal{P}(\mathbb{R}^{d})\times\mathbb{R}^{d}\to[-\infty,+\infty]$ is
continuous in the second variable thanks to Assumption 1. Similarly, $\int
f_{1}(z,x)\mathrm{d}\rho(z)+\int f_{2}(z,x)\mathrm{d}\bar{\rho}(z)$ is
continuous in $\rho$ thanks to [San15, Proposition 7.1] using the continuity
of $f_{1}$ and $f_{2}$. Further, $-P$ is upper semi-continuous using [Pos75]
and [San15, Proposition 7.2] thanks to Assumptions 2 and 3. This concludes the
continuity properties for $G_{c}$.

The upper semi-continuity of $G_{b}$ then follows from a direct application of
a version of Berge’s maximum theorem [AB06, Lemma 16.30]. Let $R:=\left\lVert
x_{0}\right\rVert^{2}+\frac{2(a_{1}+a_{2})}{\beta}>0$. We define
$\varphi:(\mathcal{P}(\mathbb{R}^{d}),W_{2})\twoheadrightarrow\mathbb{R}^{d}$
as the correspondence that maps any $\rho\in\mathcal{P}(\mathbb{R}^{d})$ to
the closed ball $\overline{B_{R}(0)}\subset\mathbb{R}^{d}$. Then the graph of
$\varphi$ is
$\text{Gr}\,\varphi=\mathcal{P}(\mathbb{R}^{d})\times\\{\overline{B_{R}(0)}\\}$.
With this definition of $\varphi$, the range of $\varphi$ is compact and
$\varphi$ is continuous with respect to weak-$*$ convergence, and so it is in
particular upper hemicontinuous. Thanks to Lemma 19, the best response
function $b(\rho)$ is always contained in $\overline{B_{R}(0)}$ for any choice
of $\rho\in\mathcal{P}(\mathbb{R}^{d})$. As a result, maximizing
$-G_{c}(\rho,x)$ in $x$ over $\mathbb{R}^{d}$ for a fixed
$\rho\in\mathcal{P}(\mathbb{R}^{d})$ reduces to maximizing it over
$\overline{B_{R}(0)}$. Using the notation introduced above, we can restrict
$G_{c}$ to $G_{c}:\text{Gr}\,\varphi\rightarrow\mathbb{R}$ and write

$\displaystyle
G_{b}(\rho)\coloneqq\max_{\hat{x}\in\varphi(\rho)}-G_{c}(\rho,\hat{x}).$

Because $G_{c}(\rho,x)$ is upper semi-continuous when
$\mathcal{P}(\mathbb{R}^{2})\times\mathbb{R}^{d}$ is endowed with the product
topology of the weak-$*$ topology and the Euclidean topology, [AB06, Lemma
16.30] guarantees that $G_{b}(\cdot)$ is upper semi-continuous in the weak-$*$
topology. ∎

######  Proposition 21.

Let $\alpha,\beta>0$ and assume Assumptions 1-4 hold with the parameters
satisfying $\alpha\tilde{\lambda}>\Lambda_{1}$. Fix
$\rho_{0},\rho_{1}\in\mathcal{P}(\mathbb{R}^{d})$. Along any geodesic
$(\rho_{s})_{s\in[0,1]}\in\mathcal{P}_{2}(\mathbb{R}^{d})$ connecting
$\rho_{0}$ to $\rho_{1}$, we have for all $s\in[0,1]$

$\displaystyle\frac{\mathrm{d}^{2}}{\mathrm{d}s^{2}}G_{b}(\rho_{s})\leq-\lambda_{b}W_{1}(\rho_{0},\rho_{1})^{2}\,,\qquad\lambda_{b}:=\alpha\tilde{\lambda}-\Lambda_{1},.$
(27)

As a result, the functional
$G_{b}:\mathcal{P}_{2}(\mathbb{R}^{d})\to[-\infty,+\infty]$ is uniformly
displacement concave with constant $\lambda_{b}>0$.

###### Proof.

Consider any $\rho_{0},\rho_{1}\in\mathcal{P}_{2}(\mathbb{R}^{d})$. Then any
$W_{2}$-geodesic $(\rho_{s})_{s\in[0,1]}$ connecting $\rho_{0}$ with
$\rho_{1}$ solves the following system of geodesic equations:

$\begin{cases}\partial_{s}\rho_{s}+\mathrm{div}\left(\rho_{s}v_{s}\right)=0\,,\\\
\partial_{s}(\rho_{s}v_{s})+\mathrm{div}\left(\rho_{s}v_{s}\otimes
v_{s}\right)=0\,,\end{cases}$ (28)

where $\rho_{s}:\mathbb{R}^{d}\to\mathbb{R}$ and
$v_{s}:\mathbb{R}^{d}\mapsto\mathbb{R}^{d}$ . The first derivative of $G_{b}$
along geodesics can be computed explicitly as

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}s}G_{b}(\rho_{s})$
$\displaystyle=\int\nabla_{z}f_{1}(z,b(\rho_{s}))\cdot
v_{s}(z)\rho_{s}(z)\mathrm{d}z-\frac{\mathrm{d}}{\mathrm{d}s}P(\rho_{s})$
$\displaystyle+\left\langle\bigg{[}\int\nabla_{x}f_{1}(z,\,x)\mathrm{d}\rho_{s}(z)+\int\nabla_{x}f_{2}(z,x)\mathrm{d}\bar{\rho}(z)+\beta(x-x_{0})\bigg{]}\bigg{|}_{x=b(\rho_{s})},\frac{\mathrm{d}}{\mathrm{d}s}b(\rho_{s})\right\rangle.$

The left-hand side of the inner product is zero by definition of the best
response $b(\rho_{s})$ to $\rho_{s}$, see (26). Therefore

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}s}G_{b}(\rho_{s})=\int\nabla_{z}f_{1}(z,b(\rho_{s}))\cdot
v_{s}(z)\rho_{s}(z)\mathrm{d}z-\frac{\mathrm{d}}{\mathrm{d}s}P(\rho_{s})\,.$

Differentiating a second time, using (28) and integration by parts, we obtain

$\displaystyle\frac{\mathrm{d}^{2}}{\mathrm{d}s^{2}}G_{b}(\rho_{s})=L_{1}(\rho_{s})+L_{2}(\rho_{s})-\frac{\mathrm{d}^{2}}{\mathrm{d}s^{2}}P(\rho_{s})\,,$

where

$\displaystyle L_{1}(\rho_{s})$
$\displaystyle:=\int\nabla_{z}^{2}f_{1}(z,b(\rho_{s}))\cdot(v_{s}\otimes
v_{s})\,\rho_{s}\mathrm{d}z=\int\left\langle
v_{s},\,\nabla_{z}^{2}f_{1}(z,b(\rho_{s}))\cdot
v_{s}\right\rangle\,\rho_{s}\mathrm{d}z\,,$ $\displaystyle L_{2}(\rho_{s})$
$\displaystyle:=\int\frac{\mathrm{d}}{\mathrm{d}s}b(\rho_{s})\cdot\nabla_{x}\nabla_{z}f_{1}(z,b(\rho_{s}))\cdot
v_{s}(z)\,\rho_{s}(z)\mathrm{d}z\,.$

From (16), we have that

$\frac{\mathrm{d}^{2}}{\mathrm{d}s^{2}}\tilde{P}(\rho_{s})\geq\alpha\tilde{\lambda}\,W_{2}(\rho_{0},\rho_{1})^{2}\,,$

and thanks to Assumption 4 it follows that

$\displaystyle L_{1}(s)$
$\displaystyle\leq\Lambda_{1}W_{2}(\rho_{0},\rho_{1})^{2}.$

This leaves $L_{2}$ to bound; we first consider the term
$\frac{\mathrm{d}}{\mathrm{d}s}b(\rho_{s})$:

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}s}b(\rho_{s})$
$\displaystyle=\int\delta_{\rho}b[\rho_{s}](\tilde{z})\partial_{s}\rho_{s}(\mathrm{d}\tilde{z})=-\int\delta_{\rho}b[\rho_{s}](\tilde{z})\mathrm{div}\left(\rho_{s}v_{s}\right)\mathrm{d}\tilde{z}$
$\displaystyle=\int\nabla_{z}\delta_{\rho}b[\rho_{s}](\tilde{z})\cdot
v_{s}(\tilde{z})\mathrm{d}\rho_{s}(\tilde{z}).$

Defining $u(\rho_{s})\in\mathbb{R}^{d}$ by

$\displaystyle
u(\rho_{s})\coloneqq\int\nabla_{x}\nabla_{z}f_{1}(z,b(\rho_{s}))\cdot
v_{s}(z)\mathrm{d}\rho_{s}(z)\,,$

using the results from Lemma 16 for $\nabla_{z}\delta_{\rho}b[\rho_{s}]$,
Assumption 1 and the fact that $Q(\rho)$ is constant in $z$ and $x$, we have

$\displaystyle L_{2}(\rho_{s})$
$\displaystyle=-\iint\left[Q(\rho_{s})^{-1}\nabla_{x}\nabla_{z}f_{1}(\tilde{z},b(\rho_{s}))\cdot
v_{s}(\tilde{z})\right]\cdot\nabla_{x}\nabla_{z}f_{1}(z,b(\rho_{s}))\cdot
v_{s}(z)\,\mathrm{d}\rho_{s}(z)\mathrm{d}\rho_{s}(\tilde{z})$
$\displaystyle=-\left\langle
u(\rho_{s}),\,Q(\rho_{s})^{-1}u(\rho_{s})\right\rangle\leq 0$

Combining all terms together, we obtain

$\displaystyle\frac{\mathrm{d}^{2}}{\mathrm{d}s^{2}}G_{b}(\rho_{s})\leq-\left(\alpha\tilde{\lambda}-\Lambda_{1}\right)W_{2}(\rho_{0},\rho_{1})^{2}\,.$

∎

######  Remark 10.

Under some additional assumptions on the functions $f_{1}$ and $f_{2}$, we can
obtain an improved convergence rate. In particular, assume that for all
$z,x\in\mathbb{R}^{d}$,

  * •

there exists a constant $\Lambda_{2}\geq\lambda_{2}\geq 0$ such that
$\nabla^{2}_{x}f_{2}(z,x)\preceq\Lambda_{2}\operatorname{I}_{d}$;

  * •

there exists a constant $\sigma\geq 0$ such that
$\left\lVert\nabla_{x}\nabla_{z}f_{1}(z,x)\right\rVert\geq\sigma$.

Then we have
$-Q(\rho_{s})^{-1}\preceq-1/(\beta+\Lambda_{1}+\Lambda_{2})\operatorname{I}_{d}$.
Using Lemma 16, we then obtain a stronger bound on $L_{2}$ as follows:

$\displaystyle L_{2}(\rho_{s})$
$\displaystyle\leq-\frac{1}{\beta+\Lambda_{1}+\Lambda_{2}}\left\lVert
u(\rho_{s})\right\rVert^{2}\leq-\frac{1}{\beta+\Lambda_{1}+\Lambda_{2}}\int\left\lVert\nabla_{x}\nabla_{z}f_{1}(z,b(\rho_{s}))\right\rVert^{2}\mathrm{d}\rho_{s}(z)\int\left\lVert
v_{s}(z)\right\rVert^{2}\mathrm{d}\rho_{s}(z)$
$\displaystyle\leq-\frac{\sigma^{2}}{\beta+\Lambda_{1}+\Lambda_{2}}W_{2}(\rho_{0},\rho_{1})^{2}.$

This means we can improve the convergence rate in (27) to
$\lambda_{b}:=\alpha\tilde{\lambda}+\frac{\sigma^{2}}{\beta+\Lambda_{1}+\Lambda_{2}}-\Lambda_{1}$.

######  Proposition 22 (Ground state).

Let Assumptions 1-4 hold for $\alpha\tilde{\lambda}>\Lambda_{1}\geq 0$ and
$\beta>0$. Then there exists a unique maximizer $\rho_{*}$ for the functional
$G_{b}$ over $\mathcal{P}(\mathbb{R}^{d})$, and it satisfies
$\rho_{*}\in\mathcal{P}_{2}(\mathbb{R}^{d})\cap L^{1}(\mathbb{R}^{d})$ and
$\rho_{*}$ is absolutely continuous with respect to $\tilde{\rho}$.

###### Proof.

Uniqueness of the maximizer (if it exists) is guaranteed by the uniform
concavity provided by Lemma 21. To show existence of a maximizer, we use the
direct method in the calculus of variations, requiring the following key
properties for $G_{b}$: (1) boundedness from above, (2) upper semi-continuity,
and (3) tightness of any minimizing sequence. To show (1), note that
$\nabla_{z}^{2}(f_{1}(z,x)+\alpha\log\tilde{\rho}(z))\preceq-(\alpha\tilde{\lambda}-\Lambda_{1})\operatorname{I}_{d}$
for all $z,x\in\mathbb{R}^{d}\times\mathbb{R}^{d}$ by Assumptions 2 and 4, and
so

$f_{1}(z,x)+\alpha\log\tilde{\rho}(z)\leq
c_{0}(x)-\frac{(\alpha\tilde{\lambda}-\Lambda_{1})}{4}|z|^{2}\qquad\forall(z,x)\in\mathbb{R}^{d}\times\mathbb{R}^{d}$
(29)

with
$c_{0}(x):=f_{1}(0,x)+\alpha\log\tilde{\rho}(0)+\frac{1}{\alpha\tilde{\lambda}-\Lambda_{1}}\|\nabla_{z}\left[f_{1}(0,x)+\alpha\log\tilde{\rho}(0)\right]\|^{2}$.
Therefore,

$\displaystyle G_{b}(\rho)$
$\displaystyle=\int\left[f_{1}(z,b(\rho))+\alpha\log\tilde{\rho}(z)\right]\mathrm{d}\rho(z)+\int
f_{2}(z,b(\rho))\mathrm{d}\bar{\rho}(z)+\frac{\beta}{2}\|b(\rho)-x_{0}\|^{2}$
$\displaystyle\qquad-\alpha\int\rho\log\rho-\int\rho W\ast\rho$
$\displaystyle\leq c_{0}(b(\rho))+\int
f_{2}(z,b(\rho))\mathrm{d}\bar{\rho}(z)+\frac{\beta}{2}\|b(\rho)-x_{0}\|^{2}\,.$

To estimate each of the remaining terms on the right-hand side, denote
$R:=\left\lVert x_{0}\right\rVert^{2}+\frac{2(a_{1}+a_{2})}{\beta}$ and recall
that $\|b(\rho)\|\leq R$ for any $\rho\in\mathcal{P}(\mathbb{R}^{d})$ thanks
to Lemma 19. By continuity of $f_{1}$ and $\log\tilde{\rho}$, there exists a
constant $c_{1}\in\mathbb{R}$ such that

$\displaystyle\sup_{x\in B_{R}(0)}c_{0}(x)=\sup_{x\in
B_{R}(0)}\left[f_{1}(0,x)+\alpha\log\tilde{\rho}(0)+\frac{1}{\alpha\tilde{\lambda}-\Lambda_{1}}\left\|\nabla_{z}\left(f_{1}(0,x)+\alpha\log\tilde{\rho}(0)\right)\right\|^{2}\right]\leq
c_{1}\,.$ (30)

The second term is controlled by $c_{2}$ thanks to Assumption 4. And the third
term can be bounded directly to obtain

$\displaystyle G_{b}(\rho)$ $\displaystyle\leq
c_{1}+c_{2}+\beta(R^{2}+\|x_{0}\|^{2})\,.$

This concludes the proof of (1). Statement (2) was shown in Lemma 20. Then we
obtain a minimizing sequence $(\rho_{n})\in\mathcal{P}(\mathbb{R}^{d})$ which
is in the closed unit ball of $C_{0}(\mathbb{R}^{d})^{*}$ and so the Banach-
Anaoglu theorem [Rud91, Theorem 3.15] there exists a limit $\rho_{*}$ in the
Radon measures and a subsequence (not relabeled) such that
$\rho_{n}\overset{\ast}{\rightharpoonup}\rho_{*}$. In fact, $\rho_{*}$ is
absolutely continuous with respect to $\tilde{\rho}$ as otherwise
$G_{b}(\rho_{*})=-\infty$, which contradicts that $G_{b}(\cdot)>-\infty$
somewhere. We conclude that $\rho_{*}\in L^{1}(\mathbb{R}^{d})$ since
$\tilde{\rho}\in L^{1}(\mathbb{R}^{d})$ by Assumption 2. To ensure
$\rho_{*}\in\mathcal{P}(\mathbb{R}^{d})$, we require (3) tightness of the
minimizing sequence $(\rho_{n})$. By Markov’s inequality [Gho02] it is
sufficient to establish a uniform bound on the second moments:

$\int\|z\|^{2}\mathrm{d}\rho_{n}(z)<C\qquad\forall n\in\mathrm{N}\,.$ (31)

To see this we proceed in a similar way as in the proof of Proposition 10.
Defining

$\displaystyle
K(\rho):=-\int\left[f_{1}(z,b(\rho))+\alpha\log\tilde{\rho}(z)\right]\mathrm{d}\rho(z)+\alpha\int\rho\log\rho\,\mathrm{d}z+\frac{1}{2}\int\rho
W\ast\rho\,\mathrm{d}z\,,$

we have $K(\rho)=-G_{b}(\rho)+\int
f_{2}(z,b(\rho))\mathrm{d}\bar{\rho}(z)+\frac{\beta}{2}\left\lVert
b(\rho)-x_{0}\right\rVert^{2}$. Then using again the bound on $b(\rho)$ from
Lemma 19,

$\displaystyle K(\rho)$ $\displaystyle\leq-G_{b}(\rho)+\sup_{x\in
B_{R}(0)}\int f_{2}(z,x)\mathrm{d}\bar{\rho}(z)+\beta\left(R^{2}+\left\lVert
x_{0}\right\rVert^{2}\right)$ $\displaystyle\leq-
G_{b}(\rho)+c_{2}+\beta\left(R^{2}+\left\lVert
x_{0}\right\rVert^{2}\right)\,,$

where the last inequality is thanks to Assumption 4. Hence, using the
estimates (29) and (30) from above, and noting that the sequence $(\rho_{n})$
is minimizing $(-G_{b})$, we have

$\displaystyle\frac{(\alpha\tilde{\lambda}-\Lambda_{1})}{4}\int\left\lVert
z\right\rVert^{2}\mathrm{d}\rho_{n}(z)$ $\displaystyle\leq
c_{0}(b(\rho_{n}))-\int\left[f_{1}(z,b(\rho_{n}))+\alpha\log\tilde{\rho}(z)\right]\mathrm{d}\rho_{n}(z)$
$\displaystyle\leq c_{1}+K(\rho_{n})\leq
c_{1}-G_{b}(\rho_{n})+c_{2}+\beta\left(R^{2}+\left\lVert
x_{0}\right\rVert^{2}\right)$ $\displaystyle\leq
c_{1}-G_{b}(\rho_{1})+c_{2}+\beta\left(R^{2}+\left\lVert
x_{0}\right\rVert^{2}\right)<\infty\,.$

which uniformly bounds the second moments of $(\rho_{n})$. This concludes the
proof for the estimate (31) and also ensures that
$\rho_{*}\in\mathcal{P}_{2}(\mathbb{R}^{d})$.

∎

######  Corollary 23.

Any maximizer $\rho_{*}$ of $G_{b}$ is a steady state for equation (5)
according to Definition 9, and satisfies
$\operatorname{supp}(\rho_{*})=\operatorname{supp}(\tilde{\rho})$.

###### Proof.

To show that $\rho_{*}$ is a steady state we can follow exactly the same
argument as in the proof of Corollary 12, just replacing
$-\frac{1}{\alpha}\int f_{1}(z,x)\,\mu_{*}(x)$ with $+\frac{1}{\alpha}\int
f_{1}(z,b(\rho_{*})$. It remains to show that
$\operatorname{supp}(\rho_{*})=\operatorname{supp}(\tilde{\rho})$. As
$\rho^{*}$ is a maximizer, it is in particular a critical point, and therefore
satisfies that $\delta_{\rho}G_{b}[\rho_{*}](z)$ is constant on all connected
components of $\operatorname{supp}(\rho_{*})$. Thanks to Lemma 18, this means
there exists a constant $c[\rho_{*}]$ (which may be different on different
components of $\operatorname{supp}(\rho_{*})$) such that

$\displaystyle
f_{1}(z,b(\rho_{*}))-\alpha\log\left(\frac{\rho_{*}(z)}{\tilde{\rho}(z)}\right)-W\ast\rho_{*}(z)=c[\rho_{*}]\qquad\text{
on }\operatorname{supp}(\rho_{*})\,.$

Rearranging, we obtain (for a possible different constant $c[\rho_{*}]\neq 0$)

$\displaystyle\rho_{*}(z)=c[\rho_{*}]\tilde{\rho}(z)\exp{\left[\frac{1}{\alpha}\left(f_{1}(z,b(\rho_{*}))-W\ast\rho_{*}(z)\right)\right]}\qquad\text{
on }\operatorname{supp}(\rho_{*})\,.$ (32)

Firstly, note that