# Mean-Field Langevin Dynamics :
Exponential Convergence and Annealing

Lénaïc Chizat<EMAIL_ADDRESS>
EPFL

###### Abstract

_Noisy particle gradient descent_ (NPGD) is an algorithm to minimize convex
functions over the space of measures that include an entropy term. In the
many-particle limit, this algorithm is described by a _Mean-Field Langevin_
dynamics—a generalization of the Langevin dynamic with a non-linear
drift—which is our main object of study. Previous work have shown its
convergence to the unique minimizer via non-quantitative arguments. We prove
that this dynamics converges at an exponential rate, under the assumption that
a certain family of Log-Sobolev inequalities holds. This assumption holds for
instance for the minimization of the risk of certain two-layer neural
networks, where NPGD is equivalent to standard noisy gradient descent. We also
study the annealed dynamics, and show that for a noise decaying at a
logarithmic rate, the dynamics converges in value to the global minimizer of
the unregularized objective function.

##  1 Introduction

Let $\mathcal{P}_{2}(\mathbb{R}^{d})$ (resp.
$\mathcal{P}^{a}_{2}(\mathbb{R}^{d})$) be the set of probability measures
(resp. absolutely continuous probability measures) with finite second moment
on $\mathbb{R}^{d}$ and let $G:\mathcal{P}_{2}(\mathbb{R}^{d})\to\mathbb{R}$
be a convex function which is “smooth” in the sense of Assumption 1 below. Our
goal is to solve problems of the form

$\min_{\mu\in\mathcal{P}_{2}^{\text{a}}(\mathbb{R}^{d})}F_{\tau}(\mu)\quad\text{where}\quad
F_{\tau}(\mu):=G(\mu)+\tau H(\mu)$ (1)

with $H(\mu)\coloneqq\int\log(\frac{\mathrm{d}\mu}{\mathrm{d}x})\mathrm{d}\mu$
the entropy of $\mu$ and $\tau>0$ the regularization/temperature parameter.
See Section 5 for examples of problems of this form (typically with $\tau=0$)
that arise in machine learning such as the regularized risk functional of wide
two-layer neural networks, or Maximum Mean Discrepancy (MMD) minimization.

##### Noisy Particle Gradient Descent (NPGD)

The starting idea of NPGD is to parameterize the measure $\mu$ as a mixture of
$m$ particles $\mu=\frac{1}{m}\sum_{i=1}^{m}\delta_{X_{i}}$. Let
$\mathbf{X}=(X_{1},\dots,X_{m})\in(\mathbb{R}^{d})^{m}$ encode the position of
all particles and consider the function

$G_{m}(\mathbf{X})\coloneqq
G\Big{(}\frac{1}{m}\sum_{i=1}^{m}\delta_{X_{i}}\Big{)}.$ (2)

Then, NPGD is just noisy gradient descent on $G_{m}$ with initialization
sampled from $\mu_{0}\in\mathcal{P}_{2}(\mathbb{R}^{d})$. It is defined, for
$k\geq 0$, as

$\mathbf{X}[k+1]=\mathbf{X}[k]-m\eta\nabla
G_{m}(\mathbf{X}[k])+\sqrt{2\eta\tau}\mathbf{Z}[k],\quad\mathbf{X}[0]\sim\mu_{0}^{\otimes
m}$ (3)

where $\eta>0$ is the step-size and $\mathbf{Z}[1],\mathbf{Z}[2],\dots$ are
i.i.d. standard Gaussian vectors (see Eq. (10) for an equivalent definition of
NPGD directly in terms of $G$ and its first-variation).

When $G$ is linear, i.e. $G(\mu)=\int V\mathrm{d}\mu$ for some smooth
$V:\mathbb{R}^{d}\to\mathbb{R}$, the particles $X_{i}$ are independent and
each follows the (unadjusted) Langevin algorithm [Ermak, 1975, Roberts and
Tweedie, 1996, Durmus and Moulines, 2017] given by the stochastic recursion

$\displaystyle X[k+1]=X[k]-\eta\nabla V(X[k])+\sqrt{2\eta\tau}Z[k],\quad
X[0]\sim\mu_{0}$ (4)

and it is thus sufficient to choose $m=1$ in that case. In the general case of
a convex and non-linear $G$, the particles will interact in non-trivial ways
and $m$ should be taken large, so that a mean-field behavior emerges.

##### Mean-Field Langevin

The dynamics obtained in the many-particle $m\to\infty$ and vanishing step-
size $\eta\to 0$ limit was called the _Mean-Field Langevin_ dynamics in Hu et
al. [2021] and is our object of interest. In this limit, the distribution
$\mu_{t}$ of particles at time $t=k\eta$ solves the following drift-diffusion
partial differential equation (PDE) of _McKean-Vlasov_ type:

$\partial_{t}\mu_{t}=\nabla\cdot\big{(}\mu_{t}\nabla
V[\mu_{t}]\big{)}+\tau\Delta\mu_{t}$ (5)

where $\nabla\cdot$ stands for the divergence operator and
$V[\mu]\in\mathcal{C}^{1}(\mathbb{R}^{d})$ is the _first-variation_ of $G$ at
$\mu$ (see Definition 2.1). This dynamics, which can be interpreted as the
gradient flow of $F_{\tau}$ under the $W_{2}$ Wasserstein metric [Ambrosio and
Savaré, 2007], is a generalization of the Langevin dynamics to a specific form
of non-linear drift term.

There is a long line of work around mean-field dynamics [Dobrushin, 1979,
Sznitman, 1991] (see Lacker [2018] for an introduction and references) which
guarantee that NPGD (3) indeed converges to the Mean-Field Langevin dynamics,
sometimes with fine quantitative bounds [Lacker, 2021, Mei et al., 2019]. As
for the behavior of the Mean-Field Langevin dynamics (5) itself, it is shown
in [Mei et al., 2018, Hu et al., 2021] that, under suitable coercivity
assumptions, $(\mu_{t})$ weakly converges to the unique minimizer of
$F_{\tau}$ as $t\to\infty$. Moreover, Hu et al. [2021] remarks that the
results from Eberle et al. [2019] to obtain quantitative rates apply here, but
this argument is restricted to the large noise regime and does not exploit the
convexity of $G$. These works leave open the question of quantitative
guarantees without a strong noise assumption.

###  1.1 Contributions and related work

Our contributions are the following:

  * –

We prove that, under a certain uniform log-Sobolev inequality assumption
(which is in particular satisfied in the settings of Mei et al. [2018], Hu et
al. [2021]), solutions to (5) converge at a global exponential rate to the
minimizer of $F_{\tau}$ (Theorem 3.2). The known convergence rate of the
Langevin dynamics under a log-Sobolev inequality is recovered as a particular
case when $G$ is linear.

  * –

We study the annealed dynamics where the noise $\tau=\tau_{t}$ is time-
dependent and decays as $\alpha/\log(t)$ and prove that for $\alpha>0$ large
enough, $G(\mu_{t})$ converges towards the minimum of the _unregularized_
functional $F_{0}=G$ (Theorem 4.1).

  * –

In Section 5, we show that our results apply to noisy gradient descent on
infinitely wide two-layer neural networks and we provide numerical experiments
for $G$ being a kernel Maximum Mean Discrepancy (MMD).

Let us mention that other algorithms to solve problems of the form (1) are
possible. Nitanda et al. [2021] proposed a dual averaging scheme which
involves a sequence of Langevin diffusions and enjoys a $O(1/t)$ convergence
rate in the mean-field limit. For low-dimensional problems, one can resort to
discretizing the measure on a fixed grid, which leads to a convex problem
amenable to standard (Bregman) gradient descent algorithms [Tseng, 2010].

The long-time behavior of drift-diffusion PDEs of the form Eq. (5) has been
studied in the mathematical physics literature. General convergence rates
(under assumptions that imply a large noise in our context) are proved in
Eberle et al. [2019]. For interacting particule systems with an interaction
kernel $k$ (discussed in Section 5.2), the case where $k(x,y)=h(x-y)$ with $h$
convex can be dealt with using the notion of _displacement convexity_ , see
[Villani, 2021, Chap. 9.6]. In contrast, we rely on standard convexity, which
corresponds to a positive semi-definite interaction kernel $k$.

Upon completion of this work, we became aware of the paper [Nitanda et al.,
2022] which also proves the exponential convergence of the Mean-Field Langevin
dynamics with the same proof technique. The main differences between these two
works is that they perform a discrete time analysis while we study the
annealed dynamics. These works were conducted independently and
simultaneously.

###  1.2 Notations

We use $\|\cdot\|$ for the Euclidean norm on $\mathbb{R}^{d}$. For
$\mu,\nu\in\mathcal{P}_{2}(\mathbb{R}^{d})$, $\Pi(\mu,\nu)$ is the set of
transport plans, that is, probability measures on
$\mathbb{R}^{d}\times\mathbb{R}^{d}$ with marginals $\mu$ and $\nu$
respectively. The Wasserstein distance
$W_{2}:\mathcal{P}_{2}(\mathbb{R}^{d})^{2}\to\mathbb{R}_{+}$ is defined as the
square-root of

$\displaystyle
W_{2}(\mu,\nu)^{2}\coloneqq\min_{\gamma\in\Pi(\mu,\nu)}\int_{(\mathbb{R}^{d})^{2}}\|y-x\|^{2}\mathrm{d}\gamma(x,y).$
(6)

Relevant background on the Wasserstein distance can be found in Ambrosio and
Savaré [2007]. We often identify absolutely continuous probability measures
with their density with respect to the Lebesgue measure. A function
$f:\mathbb{R}^{d}\to\mathbb{R}$ is said $L$-smooth if its gradient is a
$L$-Lipschitz continuous function.

##  2 Assumptions and preliminaries

###  2.1 First-variation and smoothness of $G$

The Mean-Field Langevin dynamics in Eq. (5) involves the _first-variation_ $V$
of $G$, defined as follows.

######  Definition 2.1 (First-variation).

We say that $G:\mathcal{P}_{2}(\mathbb{R}^{d})\to\mathbb{R}$ admits a _first-
variation_ at $\mu\in\mathcal{P}_{2}(\mathbb{R}^{d})$ if there exists a
continuous function $V[\mu]:\mathbb{R}^{d}\to\mathbb{R}$ such that

$\forall\nu\in\mathcal{P}_{2}(\mathbb{R}^{d}),\
\lim_{\epsilon\,\downarrow\,0}\frac{1}{\epsilon}\big{(}G((1-\epsilon)\mu+\epsilon\nu)-G(\mu)\big{)}=\int_{\mathbb{R}^{d}}V[\mu](x)\mathrm{d}(\nu-\mu)(x).$
(7)

If it exists, the first-variation $V[\mu]$ is unique up to an additive
constant.

The notion of first-variation appears naturally when studying variational
problems over $\mathcal{P}_{2}(\mathbb{R}^{d})$ and its precise definition
varies across references, see e.g. [Santambrogio, 2015, Def. 7.12]. Throughout
our work, we make the following regularity assumptions on $G$.

######  Assumption 1 (Smoothness of $G$).

For all $\mu\in\mathcal{P}_{2}(\mathbb{R}^{d})$, $G$ admits a first-variation
$V[\mu]\in\mathcal{C}^{1}(\mathbb{R}^{d})$ and $(\mu,x)\to\nabla V[\mu](x)$ is
Lipschitz continuous in the following sense: there exists $L>0$ such that
$\forall\mu,\nu\in\mathcal{P}_{2}(\mathbb{R}^{d}),\quad\forall
x,y\in\mathbb{R}^{d},\quad\|\nabla V[\mu](x)-\nabla V[\nu](y)\|_{2}\leq
L\big{(}\|x-y\|_{2}+W_{2}(\mu,\nu)\big{)}.$

Let us now state a lemma that is useful in our proofs, that gives the
evolution of $G$ and $H$ along dynamics $(\mu_{t})_{t\in(a,b)}$ in
$\mathcal{P}_{2}^{a}(\mathbb{R}^{d})$ that solve the _continuity equation_ (in
the sense of distributions):

$\partial_{t}\mu_{t}=-\nabla\cdot(\mu_{t}v_{t})$ (8)

for some time-dependent velocity field $v\in L^{2}((a,b),L^{2}(\mu_{t}))$.
Observe that Eq. (5) is an equation of this form with $v_{t}=-\nabla
V[\mu_{t}]-\tau\nabla\log(\mu_{t})$.

######  Lemma 2.2 (Chain rule).

Let $(\mu_{t})_{t\in(a,b)}$ be a weakly continuous solution to Eq. (8) such
that $\nabla\log(\mu_{t})\in L^{2}((a,b),L^{2}(\mu_{t}))$. Then $G(\mu_{t})$
and $H(\mu_{t})$ are absolutely continuous functions of $t$ and it holds for
a.e. $t\in(a,b)$,

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}G(\mu_{t})=\int_{\mathbb{R}^{d}}\nabla
V[\mu_{t}]^{\top}v_{t}\mathrm{d}\mu_{t}$ and
$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}H(\mu_{t})=\int_{\mathbb{R}^{d}}(\nabla\log(\mu_{t}))^{\top}v_{t}\mathrm{d}\mu_{t}.$

###### Proof.

Using the vocabulary of analysis in Wasserstein space, the function $H$ is
displacement convex with subdifferential $\nabla\log\mu$ [Ambrosio and Savaré,
2007, Thm. 4.16]. Also we prove in Lemma A.2 that $G$ is $(-2L)$-displacement
convex with subdifferential $\nabla V[\mu_{t}]$. Then the claim is a
consequence of [Ambrosio and Savaré, 2007, Sec. 4.4.E]. ∎

###  2.2 Characterization of the minimizer

We recall the optimality conditions for $F_{\tau}$ which have been proved in
several works (see e.g. [Mei et al., 2018, Lem. 10.4] or [Hu et al., 2021,
Prop. 2.5]) and require the following assumptions.

######  Assumption 2.

The function $G$ is convex and $F_{\tau}=G+\tau H$ admits a minimizer
$\mu^{*}_{\tau}$.

We stress that by convexity we mean standard convexity for the linear
structure in $\mathcal{P}_{2}(\mathbb{R}^{d})$, i.e.

$\forall\mu,\nu\in\mathcal{P}_{2}(\mathbb{R}^{d}),\forall\alpha\in[0,1],\quad
G(\alpha\mu+(1-\alpha)\nu)\leq\alpha G(\mu)+(1-\alpha)G(\nu).$

######  Proposition 2.3.

Under Assumption 1 and 2, the minimizer $\mu^{*}_{\tau}$ of $F$ is unique and
satisfies

$\mu^{*}_{\tau}\propto e^{-V[\mu^{*}_{\tau}]/\tau}.$ (9)

The uniqueness comes from the strict convexity of $H$. For Eq. (9), one first
derives the first order optimality condition, which require that
$V[\mu_{\tau}^{*}]+\tau\log(\mu_{\tau}^{*})$ must be a constant
$\mu_{\tau}^{*}$-almost everywhere. Then one shows that $\mu_{\tau}^{*}$ has
positive density everywhere due to the entropy term, and concludes. We refer
to [Mei et al., 2018, Lem. 10.4] for details.

###  2.3 Noisy Particle Gradient Descent (NPGD)

The NPGD algorithm has been defined in Section 1 via the function $G_{m}$. We
now give an alternative definition of this algorithm involving the first-
variation $V$ of $G$.

Assume that $G$ satisfies Assumption 1, let $V$ be its first-variation and fix
$m\in\mathbb{N}^{*}$. For $i\in[m]$, initialize randomly
$X_{i,0}\overset{iid}{\sim}\mu_{0}\in\mathcal{P}_{2}(\mathbb{R}^{d})$ and
define recursively for $k\geq 0$

$\left\\{\begin{aligned} X_{i,k+1}&=X_{i,k}-\eta\nabla
V[\hat{\mu}_{k}](X_{i,k})+\sqrt{2\eta\tau}Z_{i,k}\\\
\hat{\mu}_{k}&=\frac{1}{m}\sum_{i=1}^{m}\delta_{X_{i,k}}\end{aligned}\right.$
(10)

where $\eta>0$ is the step-size and $Z_{i,k}\sim\mathcal{N}(0,I)$ are iid
standard Gaussian random variables.

######  Proposition 2.4.

Under Assumption (1), the two definitions of NPGD in Eq. (3) and in Eq. (10)
are equivalent.

###### Proof.

For $\mathbf{X}\in(\mathbb{R}^{d})^{m}$ and
$\mathbf{Y}\in(\mathbb{R}^{d})^{m}$ define
$\mu_{t}=\frac{1}{m}\sum_{i=1}^{m}\delta_{X_{i}+tY_{i}}$. It satisfies the
continuity equation (8) with velocity field $v_{t}(X_{i}+tY_{i})=Y_{i}$. By
Lemma 2.2, it holds

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}G_{m}(\mathbf{X}+t\mathbf{Y})|_{t=0}$
$\displaystyle=\frac{\mathrm{d}}{\mathrm{d}t}G(\mu_{t})|_{t=0}=\int_{\mathbb{R}^{d}}\nabla
V[\mu_{0}](x)^{\top}v_{0}(x)\mathrm{d}\mu_{0}(x)=\frac{1}{m}\sum_{i=1}^{m}\nabla
V[\mu_{0}](X_{i})^{\top}Y_{i}.$

This proves that $\forall i\in[m]$, $m\nabla_{X_{i}}G_{m}(\mathbf{X})=\nabla
V[\mu_{0}](X_{i})$ and thus the update equations in Eq. (3) and Eq. (10) are
the same. ∎

###  2.4 Mean-Field Langevin dynamics

Given Eq. (10) standard results about mean-field systems tell us that as
$m\to\infty$, the random measure $\hat{\mu}_{k}$ becomes deterministic, so
that in the limit (and taking also the small-step size limit $\eta\to 0$) the
particles trajectories are given by i.i.d. samples from the following
stochastic differential equation (SDE)

$\displaystyle\left\\{\begin{aligned} \mathrm{d}X_{t}&=-\nabla
V[\mu_{t}](X_{t})\mathrm{d}t+\sqrt{2\tau}\mathrm{d}B_{t},\quad
X_{0}\sim\mu_{0}\\\ \mu_{t}&=\mathrm{Law}(X_{t})\end{aligned}\right.$ (11)

where $(B_{t})_{t\geq 0}$ is a Brownian motion. As mentioned in the
introduction, the law $(\mu_{t})$ of a solution to this SDE solves the
following PDE which is our main object of study:

$\displaystyle\partial_{t}\mu_{t}=\nabla\cdot\big{(}\mu_{t}\nabla
V[\mu_{t}]\big{)}+\tau\Delta\mu_{t}$ (12)

where $\nabla\cdot$ stands for the divergence operator. Standard results about
this class of PDEs guarantee its well-posedness, i.e. the existence of a
unique solution, under Assumption 1 (see e.g. [Huang et al., 2021, Thm. 3.3]
or Ambrosio and Savaré [2007] for an approach based on the gradient flow
structure which applies here thanks to Lemma A.2 which states that $G$ is
$(-2L)$-displacement convex).

Let us now study the convergence of $(\mu_{t})_{t\geq 0}$ to the global
minimum of $F_{\tau}$.

##  3 Exponential convergence of Mean-Field Langevin dynamics

For $\mu,\nu\in\mathcal{P}_{2}^{a}(\mathbb{R}^{d})$ with $\mu$ absolutely
continuous w.r.t. $\nu$ we define the _relative entropy_ (a.k.a. Kullback-
Leibler divergence) by

$H(\mu|\nu)\coloneqq\int_{\mathbb{R}^{d}}\log\Big{(}\frac{\mathrm{d}\mu}{\mathrm{d}\nu}\Big{)}\mathrm{d}\mu,$

and the _relative Fisher information_ by

$I(\mu|\nu)\coloneqq\int_{\mathbb{R}^{d}}\left\|\nabla\log\frac{\mathrm{d}\mu}{\mathrm{d}\nu}\right\|^{2}\mathrm{d}\mu.$

######  Definition 3.1 (Log-Sobolev inequality).

We say that $\nu\in\mathcal{P}_{2}^{a}(\mathbb{R}^{d})$ satisfies a
logarithmic Sobolev inequality with constant $\rho>0$ (in short
$\mathrm{LSI}(\rho)$) if for all $\mu\in\mathcal{P}_{2}^{a}(\mathbb{R}^{d})$
absolutely continuous w.r.t. $\nu$, it holds

$\displaystyle H(\mu|\nu)\leq\frac{1}{2\rho}I(\mu|\nu).$ (13)

This inequality can be interpreted as a $2$-Łojasiewicz gradient inequality
for the functional $\mu\mapsto H(\mu|\nu)=\int V\mathrm{d}\mu+H(\mu)$ (where
we have posed $V=-\log\nu$) in the Wasserstein geometry [Otto and Villani,
2000] and thus directly implies the exponential convergence of its Wasserstein
gradient flow. This corresponds to our objective function in the linear case
$G(\mu)=\int V\mathrm{d}\mu$, and in this case exponential convergence towards
minimizers is thus guaranteed when $\nu=e^{-V}$ satisfies a Log-Sobolev
inequality.

In the general case, we make an analogous assumption that such an inequality
holds _uniformly_ for $e^{-V[\mu_{t}]/\tau}$ throughout the dynamics.

######  Assumption 3 (Uniform log-Sobolev).

There exists $\rho_{\tau}>0$ such that
$\forall\mu\in\mathcal{P}_{2}(\mathbb{R}^{d})$ it holds $e^{-V[\mu]/\tau}\in
L^{1}(\mathbb{R}^{d})$ and the probability measure $\nu\propto
e^{-V[\mu]/\tau}$ satisfies $\mathrm{LSI}(\rho_{\tau})$.

Remember that $V[\mu]$ is defined up to a constant term, and in this section,
we fix this constant so that
$e^{-V[\mu]/\tau}\in\mathcal{P}_{2}(\mathbb{R}^{d})$. Let us recall two
criteria for a probability measure to satisfy a Log-Sobolev inequality:

  * –

If $\nabla^{2}V\succeq\rho I_{d}$ then $e^{-V}\in\mathcal{P}(\mathbb{R}^{d})$
satisfies $\mathrm{LSI}(\rho)$ [Bakry and Émery, 1985];

  * –

if $\nu$ satisfies $\mathrm{LSI}(\rho)$ and
$\tilde{\nu}=e^{-\psi}\nu\in\mathcal{P}(\mathbb{R}^{d})$ is a perturbation of
$\nu$ with $\psi\in L^{\infty}(\mathbb{R}^{d})$ then $\tilde{\nu}$ satisfies
$\mathrm{LSI}(\tilde{\rho})$ with $\tilde{\rho}=\rho e^{\inf\psi-\sup\psi}$
[Holley and Stroock, 1987].

These two criteria are standard, but many finer criteria are known, such as
integral conditions [Wang, 2001], Lyapunov conditions [Cattiaux et al., 2010]
or criteria for mixture distributions [Chen et al., 2021] (see also Section
5.1). Our main result regarding the Mean-Field Langevin dynamics (11) is the
following.

######  Theorem 3.2.

Under Assumptions 1, 2 and 3, let $\mu_{0}\in\mathcal{P}_{2}(\mathbb{R}^{d})$
be such that $F_{\tau}(\mu_{0})<\infty$. For $t\geq 0$, it holds
$\displaystyle F_{\tau}(\mu_{t})-F_{\tau}(\mu^{*}_{\tau})\leq
e^{-2\tau\rho_{\tau}t}(F_{\tau}(\mu_{0})-F_{\tau}(\mu^{*}_{\tau})).$ (14)

###### Proof.

Let $\nu_{t}=e^{-V[\mu_{t}]/\tau}$. By Lemma 2.2 applied to $v_{t}=-\nabla
V[\mu_{t}]-\tau\nabla\log(\mu_{t})$, we have

$\displaystyle\frac{d}{dt}F_{\tau}(\mu_{t})=-\int_{\mathbb{R}^{d}}\left\|\nabla
V[\mu_{t}]+\tau\nabla\log(\mu_{t})\right\|^{2}\mathrm{d}\mu_{t}=-\tau^{2}I(\mu_{t}|\nu_{t}).$
(15)

Note that although Lemma 2.2 requires some regularity estimates, they can be
bypassed here thanks to general results about Wasserstein gradient flows
[Ambrosio and Savaré, 2007, Thm. 5.3 (v)]. Combining this energy identity with
the log-Sobolev inequality and Lemma 3.4, it follows

$\displaystyle\frac{d}{dt}\big{(}F_{\tau}(\mu_{t})-F_{\tau}(\mu^{*}_{\tau})\big{)}$
$\displaystyle=-\tau^{2}I(\mu_{t}|\nu_{t})\leq-2\rho_{\tau}\tau^{2}H(\mu_{t}|\nu_{t})\leq-2\rho_{\tau}\tau(F(\mu_{t})-F(\mu^{*}))$

which is a $2$-Łojasiewicz gradient inequality for $F_{\tau}$. By integrating
in time we get Eq. (14). ∎

In the proof, we see that we could relax Assumption 3 and require the Log-
Sobolev inequality to hold only for all
$\mu\in\mathcal{P}_{2}(\mathbb{R}^{d})$ such that $F_{\tau}(\mu)\leq
F_{\tau}(\mu_{0})$. Also, Assumption 1 is only a general assumption that
guarantees well-posedness of the dynamics and the energy decay formula Eq.
(15); this regularity assumption can be relaxed on a case by case basis.
Convergence guarantees in parameter space directly follow from the previous
theorem.

######  Corollary 3.3.

Under the assumptions of Theorem 3.2, for $t\geq 0$ we have

$\displaystyle
H(\mu_{t}|\mu^{*}_{\tau})\leq\frac{1}{\tau}e^{-2\tau\rho_{\tau}t}(F_{\tau}(\mu_{0})-F_{\tau}(\mu^{*}))$
and $\displaystyle
W_{2}^{2}(\mu_{t},\mu^{*}_{\tau})\leq\frac{2e^{-2\tau\rho_{\tau}t}}{\tau\rho_{\tau}}\big{(}F_{\tau}(\mu_{0})-F_{\tau}(\mu^{*})\big{)}.$

###### Proof.

The first inequality follows from Theorem 3.2 and Lemma 3.4. For the second
one, it follows from the fact that if $\nu$ satisfies $\mathrm{LSI}(\rho)$,
then it satisfies the _Talagrand inequality_ , which states that
$\forall\mu\in\mathcal{P}_{2}(\mathbb{R}^{d})$,
$W_{2}^{2}(\mu,\nu)\leq\frac{2}{\rho}H(\mu|\nu)$, as proved in Otto and
Villani [2000]. ∎

The following lemma establishes inequalities which are key to handle the non-
linear aspect of the dynamics (when $G$ is linear, they become trivial
equalities).

######  Lemma 3.4 (Entropy Sandwich).

Under Assumption 1, 2 and 3, let $\mu^{*}_{\tau}$ be the unique minimizer of
$F_{\tau}$. For all $\mu\in\mathcal{P}_{2}(\mathbb{R}^{d})$, letting
$\nu\coloneqq e^{-V[\mu]/\tau}\in\mathcal{P}_{2}(\mathbb{R}^{d})$, it holds

$\tau H(\mu|\mu^{*}_{\tau})\leq F_{\tau}(\mu)-F_{\tau}(\mu^{*}_{\tau})\leq\tau
H(\mu|\nu).$

###### Proof.

The convexity of $G$ implies that,
$\forall\mu,\nu\in\mathcal{P}_{2}(\mathbb{R}^{d})$,
$\frac{1}{\epsilon}(G((1-\epsilon)\mu+\epsilon\nu)-G(\mu))\leq G(\nu)-G(\mu).$
So, passing to the limit in the definition of the first-variation (Definition
2.1), we recover the usual convexity inequality (interpreting $V[\mu]$ as the
gradient of $G$ at $\mu$):

$G(\nu)\geq G(\mu)+\int V[\mu]\mathrm{d}(\nu-\mu).$ (16)

Invoking this inequality twice with the role of $\mu$ and $\mu^{*}$ exchanged,
it holds

$\displaystyle\int V[\mu^{*}]\mathrm{d}(\mu-\mu^{*})\leq
G(\mu)-G(\mu^{*})\leq\int V[\mu]\mathrm{d}(\mu-\mu^{*}).$

Recalling $F_{\tau}(\mu)=G(\mu)+\tau H(\mu)$, it holds, on the one hand,

$\displaystyle F_{\tau}(\mu)-F_{\tau}(\mu^{*})$ $\displaystyle\leq\int
V[\mu]\mathrm{d}\mu+\tau H(\mu)-\int V[\mu]\mathrm{d}\mu^{*}-\tau H(\mu^{*})$
$\displaystyle=\tau H(\mu|\nu)-\tau H(\mu^{*}|\nu)\leq\tau H(\mu|\nu).$

On the other hand, using the fact that $\mu^{*}_{\tau}=e^{-V[\mu^{*}]/\tau}$
(Proposition 2.3), it holds

$\displaystyle F_{\tau}(\mu)-F_{\tau}(\mu^{*}_{\tau})$ $\displaystyle\geq\int
V[\mu^{*}]\mathrm{d}\mu+\tau H(\mu)-\int V[\mu^{*}]\mathrm{d}\mu^{*}-\tau
H(\mu^{*})$ $\displaystyle=\tau H(\mu|\mu^{*})-\tau H(\mu^{*}|\mu^{*})=\tau
H(\mu|\mu^{*}).\qed$

##  4 Convergence of the annealed dynamics

We now turn our attention to the “annealed” Mean-Field Langevin dynamics

$\displaystyle\partial\mu_{t}=\nabla\cdot\big{(}\mu_{t}\nabla
V[\mu_{t}]\big{)}+\tau_{t}\Delta\mu_{t}$ (17)

with a time-dependent _temperature_ parameter $\tau_{t}$ that converges to
$0$. The existence of a unique solution from any
$\mu_{0}\in\mathcal{P}_{2}(\mathbb{R}^{d})$ follows again from the theory of
McKean-Vlasov equations, now with time inhomogeneous coefficients (see e.g.
[Huang et al., 2021, Thm. 3.3]). As a side note, notice that (17) cannot
strictly be interpreted as a Wasserstein gradient flow anymore, but some
aspects of the theory of Wasserstein gradient flows have been extended to
cover the case of time-dependent diffusion coefficients [Ferreira and
Valencia-Guevara, 2018, Sec. 6.2].

The linear case when $G(\mu)=\int V\mathrm{d}\mu$ has been considered in
numerous works (e.g. Holley et al. [1989], Geman and Hwang [1986], Miclo
[1992], Raginsky et al. [2017], Tang and Zhou [2021]). It is known in
particular [Miclo, 1992] that under suitable coercivity assumptions for $V$
and if $\tau_{t}=C/\log(t)$ for some $C>0$ large enough, then $G(\mu_{t})$
converges to
$\min_{\mu\in\mathcal{P}_{2}(\mathbb{R}^{d})}G(\mu)=\min_{x\in\mathbb{R}^{d}}V(x)$.

Here we show that a similar guarantee holds in our more general context.

######  Theorem 4.1 (Convergence of annealed dynamics).

Suppose Assumptions 1, 2 and 3 hold for all $\tau>0$, and moreover assume
that: – the Log-Sobolev constants satisfy $\rho_{\tau}\geq
C_{0}e^{-\alpha^{*}/\tau}$ for some $\alpha^{*},C_{0}>0$, – $G$ is lower-
bounded, – $(\tau_{t})_{t}$ is smooth, decreases, and for $t$ large it holds
$\tau_{t}=\alpha/\log(t)$ for some $\alpha>\alpha^{*}$. Let
$\mu_{0}\in\mathcal{P}_{2}(\mathbb{R}^{d})$ be such that
$F_{\tau_{0}}(\mu_{0})<\infty$. Then for each $\epsilon>0$, there exists
$C,C^{\prime}>0$ such that $\displaystyle
F_{\tau_{t}}(\mu_{t})-F_{\tau_{t}}(\mu^{*}_{\tau_{t}})$ $\displaystyle\leq
Ct^{-\big{(}1-\frac{\alpha^{*}}{\alpha}-\epsilon\big{)}},$ (18) and
$\displaystyle G(\mu_{t})-\inf G$ $\displaystyle\leq C^{\prime}\frac{\log\log
t}{\log t}.$ (19)

We can make the following comments:

  * –

The lower-bound assumed on $\rho_{\tau}$ is natural when one has in mind the
Holley and Stroock criterion given in Section 3. In Section 5, we show a lower
bound of this form on a concrete example related to two-layer neural networks.

  * –

The bounds of Theorem 4.1 exhibit a two time-scales phenomenon: the dynamics
$(\mu_{t})$ converges at a polynomial rate to the regularization path
$(\mu^{*}_{\tau_{t}})$ (in relative entropy or $W_{2}^{2}$ distance, thanks to
the “entropy sandwich” Lemma 3.4 or the Talagrand inequality) but the
regularization path only converges at a logarithmic rate to the optimal value
$\inf G$, because of the slow decay of $\tau_{t}$.

  * –

The slow decay of $\tau_{t}$ is an inconvenience but it cannot be improved. It
is known that in the linear case $G(\mu)=\int V\mathrm{d}\mu$, convergence is
lost if $\tau_{t}$ decays faster [Holley et al., 1989, Sec. 3] (in fact,
taking $\tau_{t}=\alpha/\log(t)$ with $\alpha>0$ too small already breaks
convergence).

###### Proof.

Our proof is partly inspired by Miclo [1992], as revisited by Tang and Zhou
[2021].

Step 1. Consider the function that returns the values of the regularization
path

$h(\tau)\coloneqq
F_{\tau}(\mu_{\tau}^{*})=\min_{\mu\in\mathcal{P}_{2}(\mathbb{R}^{d})}G(\mu)+\tau
H(\mu).$

As an infimum of affine functions, $h$ is concave and since the minimizer
$\mu^{*}_{\tau}$ is unique, $h$ is differentiable for $\tau>0$ and its
derivative is $h^{\prime}(\tau)=H(\mu^{*}_{\tau})$. We focus on $t\geq t_{0}$
so that $\tau_{t}=\alpha/\log(t)$. By Lemma 2.2 applied to $v_{t}=-\nabla
V[\mu_{t}]+\tau_{t}\nabla\log(\mu_{t})$ (here again, the regularity
assumptions of Lemma 2.2 can be bypassed using the gradient flow-like
structure, see [Ferreira and Valencia-Guevara, 2018, Thm. 6.9]), we have

$\displaystyle\frac{d}{dt}\big{(}F_{\tau_{t}}(\mu_{t})-F_{\tau_{t}}(\mu_{\tau_{t}}^{*})\big{)}$
$\displaystyle=-\int\|\nabla
V[\mu_{t}]+\tau_{t}\nabla\log(\mu_{t})\|^{2}\mathrm{d}\mu_{t}+\tau^{\prime}_{t}H(\mu_{t})-\tau^{\prime}_{t}h^{\prime}(\tau_{t})$
$\displaystyle\leq-\tau_{t}^{2}I(\mu_{t}|\nu_{t})+\tau^{\prime}_{t}\big{(}H(\mu_{t})-H(\mu^{*}_{\tau_{t}})\big{)}$

where we introduced the probability measure $\nu_{t}\propto
e^{-V[\mu_{t}]/\tau_{t}}$. On the one hand, we have by the Log-Sobolev
inequality and the “entropy sandwich” Lemma 3.4,

$\tau^{2}_{t}I(\mu_{t}|\nu_{t})\geq
2\tau^{2}_{t}\rho_{\tau_{t}}H(\mu_{t}|\nu_{t})\geq
2\rho_{\tau_{t}}\tau_{t}\big{(}F_{\tau_{t}}(\mu_{t})-F_{\tau_{t}}(\mu_{\tau_{t}}^{*})\big{)}.$

On the other hand, by Lemma 4.2 below, it holds for some $C_{2},C_{3}>0$
independent from $\mu$ and $t$,

$\displaystyle-\tau_{t}(H(\mu_{t})-H(\mu^{*}_{\tau_{t}}))$ $\displaystyle\leq
C_{2}\tau_{t}F_{\tau_{t}}(\mu_{t})+\tau_{t}C_{3}+F_{\tau_{t}}(\mu^{*}_{\tau_{t}})-G(\mu^{*}_{\tau_{t}})$
$\displaystyle\leq
C_{2}\tau_{t}(F_{\tau_{t}}(\mu_{t})-F_{\tau_{t}}(\mu^{*}_{\tau_{t}}))+(1+C_{2}\tau_{t})F_{\tau_{t}}(\mu^{*}_{\tau_{t}})-G(\mu^{*}_{\tau_{t}})+\tau_{t}C_{3}$
$\displaystyle\leq
C^{\prime}_{2}(F_{\tau_{t}}(\mu_{t})-F_{\tau_{t}}(\mu^{*}_{\tau_{t}}))+C^{\prime}_{3}$

where in the last step, we used that $G$ is lower bounded and
$h(\tau)=F_{\tau}(\mu^{*}_{\tau})$ is bounded for $\tau\in[0,\tau_{0}]$.
Combining the previous estimates, we get that for any $\epsilon>0$, there
exists $C_{1},C_{2},C_{3}>0$ such that

$\displaystyle\frac{d}{dt}\big{(}F_{\tau_{t}}(\mu_{t})-F_{\tau_{t}}(\mu_{\tau_{t}}^{*})\big{)}$
$\displaystyle\leq-2\rho_{\tau_{t}}\tau_{t}\big{(}F_{\tau_{t}}(\mu_{t})-F_{\tau_{t}}(\mu_{\tau_{t}}^{*})\big{)}-C_{2}\frac{\tau_{t}^{\prime}}{\tau_{t}}\big{(}F_{\tau_{t}}(\mu_{t})-F_{\tau_{t}}(\mu_{\tau_{t}}^{*})\big{)}-\frac{\tau_{t}^{\prime}}{\tau_{t}}C_{3}$
$\displaystyle\leq\frac{-2C_{1}\alpha
t^{-\frac{\alpha^{*}}{\alpha}}+C_{2}t^{-1}}{\log
t}\big{(}F_{\tau_{t}}(\mu_{t})-F_{\tau_{t}}(\mu^{*}_{\tau_{t}})\big{)}+C_{3}\frac{t^{-1}}{\log
t}$

where we used $\tau_{t}=\alpha/\log(t)$, $\tau^{\prime}_{t}=-\alpha/(t(\log
t)^{2})$ and $\rho_{\tau_{t}}\geq C_{0}t^{-\alpha^{*}/\alpha}$. In passing,
the first inequality in the above display guarantees that
$F_{\tau_{t}}(\mu_{t})-F_{\tau_{t}}(\mu_{\tau_{t}}^{*})$ remains finite at all
time because $\log\tau_{t}\in\mathcal{C}^{1}$, which justifies the fact that
we can consider only $t$ large enough in the rest of the proof.

It follows that for any $\epsilon>0$ such that $\epsilon<1-\alpha^{*}/\alpha$,
for $t$ large enough and some $C,C^{\prime}>0$,

$\displaystyle\frac{d}{dt}\big{(}F_{\tau_{t}}(\mu_{t})-F_{\tau_{t}}(\mu_{\tau_{t}}^{*})\big{)}$
$\displaystyle\leq-
Ct^{-\frac{\alpha^{*}}{\alpha}-\epsilon}\big{(}F_{\tau_{t}}(\mu_{t})-F_{\tau_{t}}(\mu_{\tau_{t}}^{*})\big{)}+C^{\prime}t^{-1}{/\log(t)}.$

Now define

$Q(t)\coloneqq\big{(}F_{\tau_{t}}(\mu_{t})-F_{\tau_{t}}(\mu_{\tau_{t}}^{*})\big{)}-\frac{C^{\prime}}{C}t^{-1+\frac{\alpha^{*}}{\alpha}+\epsilon}$

which satisfies

$\displaystyle\frac{d}{dt}Q(t)\leq-
Ct^{-\frac{\alpha^{*}}{\alpha}-\epsilon}Q(t)-C^{\prime}t^{-1}+Ct^{-1}{/\log(t)}+\frac{C^{\prime}(1-\frac{\alpha^{*}}{\alpha}-\epsilon)}{C}t^{-2+\frac{\alpha^{*}}{\alpha}+\epsilon}.$
(20)

Observe that the term $-C^{\prime}t^{-1}$ dominates the two last terms for $t$
large enough. Thus for $t\geq t_{*}$ large enough, $\frac{d}{dt}Q(t)\leq-
Ct^{-\frac{\alpha^{*}}{\alpha}-\epsilon}Q(t)$ which implies $Q(t)\leq
Q(t_{*})\exp(-C\int_{t_{*}}^{t}s^{-\frac{\alpha^{*}}{\alpha}-\epsilon}\mathrm{d}s)$.
As a consequence

$F_{\tau_{t}}(\mu_{t})-F_{\tau_{t}}(\mu_{\tau_{t}}^{*})\leq\frac{C^{\prime}}{C}t^{-1+\frac{\alpha^{*}}{\alpha}+\epsilon}+Q(t_{*})\exp\Big{(}-\frac{C}{\kappa}(t^{\kappa}-t_{*}^{\kappa})\Big{)}$

and thus $F_{\tau_{t}}(\mu_{t})-F_{\tau_{t}}(\mu_{\tau_{t}}^{*})\leq
C^{\prime\prime}t^{-\kappa}$ because $\kappa\coloneqq
1-\frac{\alpha^{*}}{\alpha}-\epsilon>0$ and $Q(t^{*})$ is finite. This proves
Eq. (18).

Step 2. Let us now prove Eq. (19), under the assumption that $G$ admits a
minimizer $\mu_{0}^{*}\in\mathcal{P}_{2}(\mathbb{R}^{d})$. The proof can be
easily adapted to the general case by choosing $\mu_{0}^{*}$ as a quasi-
minimizer such that $G(\mu_{0}^{*})\leq\inf G+\epsilon$ and taking $\epsilon$
arbitrarily small. Remember that $h(0)=G(\mu_{0}^{*})=F_{0}(\mu^{*}_{0})$, so

$\displaystyle G(\mu_{t})-G(\mu_{0}^{*})$
$\displaystyle=F_{\tau_{t}}(\mu_{t})-F_{\tau_{t}}(\mu_{\tau_{t}}^{*})+F_{\tau_{t}}(\mu_{\tau_{t}}^{*})-F_{0}(\mu_{0}^{*})-\tau_{t}H(\mu_{t})$
$\displaystyle\leq Ct^{-\kappa}+(h(\tau_{t})-h(0))+C^{\prime}\tau_{t}$

where we have used the bound $-H(\mu_{t})\leq
C_{1}F_{\tau_{t}}(\mu_{t})+C_{2}$ from Lemma 4.2, which is uniformly bounded
for $t\geq 0$ by some $C^{\prime}$ thanks to Step 1.

The rest of the proof consists in bounding $h(\tau)-h(0)$ via an approximation
argument. Let
$g_{\sigma}(x)=(2\pi\sigma^{2})^{-d/2}\exp(-\|x\|^{2}/(2\sigma^{2}))$ be the
standard Gaussian kernel and let $\tilde{\mu}_{\sigma}(x)=\int
g_{\sigma}(x-y)\mathrm{d}\mu_{0}^{*}(y)$. We consider the transport plan
$\gamma\in\Pi(\tilde{\mu}_{0},\tilde{\mu}_{\sigma})$ given by the joint law of
$(X,X+Z)$ for $\mathrm{Law}(X)=\tilde{\mu}_{0}=\mu_{0}^{*}$ and
$\mathrm{Law}(Z)=g_{\sigma}$. On the one hand, it holds by convexity of $G$

$\displaystyle 0\geq G(\tilde{\mu}_{0})-G(\tilde{\mu}_{\sigma})$
$\displaystyle\geq\int
V[\tilde{\mu}_{\sigma}]\mathrm{d}[\tilde{\mu}_{0}-\tilde{\mu}_{\sigma}]$
$\displaystyle=\int(V[\tilde{\mu}_{\sigma}](y)-V[\tilde{\mu}_{\sigma}](x))\mathrm{d}\gamma(x,y).$

It follows, using the smoothness bound
$V[\tilde{\mu}_{\sigma}](x)-V[\tilde{\mu}_{\sigma}](y)\leq\nabla
V[\tilde{\mu}_{\sigma}](y)^{\top}(x-y)+\frac{L}{2}\|y-x\|^{2}$ and the fact
that the Gaussian kernel is centered, that

$\displaystyle|G(\tilde{\mu}_{0})-G(\tilde{\mu}_{\sigma})|$
$\displaystyle\leq\int(V[\tilde{\mu}_{\sigma}](x)-V[\tilde{\mu}_{\sigma}](y))\mathrm{d}\gamma(x,y)$
$\displaystyle\leq\int\nabla
V[\tilde{\mu}_{\sigma}](x)^{\top}(y-x)\mathrm{d}\gamma(x,y)+\frac{L}{2}\int\|y-x\|^{2}\mathrm{d}\gamma(x,y)$
$\displaystyle=0+\frac{L}{2}\sigma^{2}.$

On the other hand, we have by Jensen’s inequality for the convex function
$\varphi:s\mapsto s\log(s)$ and Fubini’s theorem:

$\displaystyle H(\tilde{\mu}_{\sigma})$ $\displaystyle=\int\varphi\Big{(}\int
g_{\sigma}(x-y)\mathrm{d}\mu_{0}(y)\Big{)}\mathrm{d}x$
$\displaystyle\leq\int\Big{(}\int\varphi(g_{\sigma}(x-y))\mathrm{d}x\Big{)}\mathrm{d}\mu_{0}(y)=-\frac{1}{2}\big{(}1+\log(2\pi\sigma^{2})\big{)}$

which is the entropy of the Gaussian distribution $g_{\sigma}$. Thus we have

$h(\tau)-h(0)\leq\inf_{\sigma>0}\frac{L}{2}\sigma^{2}-\frac{\tau}{2}\big{(}1+\log(2\pi\sigma^{2})\big{)}\leq-\frac{\tau}{2}\log(\pi\tau)$

by choosing $\sigma^{2}=\tau/L$. Plugging the value of
$\tau_{t}=\alpha/\log(t)$ we get, for some $C,C^{\prime}>0$,

$G(\mu_{t})-G(\mu_{0}^{*})\leq\frac{\alpha}{2}\frac{(\log\log
t-\log(\pi\alpha)}{\log t}+Ct^{-\kappa}+C^{\prime}\frac{\alpha}{\log(t)}\leq
C^{\prime\prime}\frac{\log\log t}{\log t}.\qed$

In the proof of Theorem 4.1, we used a lower bound on the value of $H(\mu)$ in
terms of the functional value that is provided in the following lemma.

######  Lemma 4.2.

Under the assumptions of Theorem 4.1, there exists $C_{1},C_{2},C_{3}>0$ such
that for all $0<\tau\leq\tau_{0}$ and
$\mu\in\mathcal{P}_{2}^{a}(\mathbb{R}^{d})$, it holds $F_{\tau}(\mu)\geq-
C_{3}$ and

$\displaystyle-H(\mu)\leq C_{1}F_{\tau}(\mu)+C_{2}.$

###### Proof.

In the following proof, $C_{i},C_{i}^{\prime},C_{i}^{\prime\prime}>0$ are
constants independent from $\mu$ which value may change from line to line.
Since by assumption the probability measure $\nu$ proportional to
$e^{-V[\mu_{0}]}$ satisfies a logarithmic Sobolev inequality, there exists
$C_{1},C_{2}>0$ such that $\forall x\in\mathbb{R}^{d}$, $V[\mu_{0}](x)\geq
C_{1}\|x\|^{2}-C_{2}$. Indeed, by Herbst argument [Bakry et al., 2014, Prop.
5.4.1], there exists $C_{1}>0$ such that $\int
e^{C_{1}\|x\|^{2}-V[\mu_{0}](x)}\mathrm{d}x<\infty$ and we conclude using the
fact that if $f\in\mathcal{C}^{1}(\mathbb{R}^{d})$ has a Lipschitz gradient
and $\int e^{f}\mathrm{d}x<\infty$ then $f$ must be upper-bounded.

Letting $M_{2}(\mu)\coloneqq\int\|x\|^{2}\mathrm{d}\mu(x)$, it follows, using
convexity of $G$, that

$G(\mu)\geq G(\mu_{0})+\int V[\mu_{0}]\mathrm{d}(\mu-\mu_{0})\geq
2C_{1}M_{2}(\mu)-C_{2}.$

Invoking Lemma 4.3 with $\sigma^{2}=\tau/C_{1}$ we have

$\tau H(\mu)\geq-C_{1}M_{2}(\mu)-\tau-\tau d\log(2\tau\pi/C_{1}).$

Summing the two previous equations (with the same value of $C_{1}$), we get
that for $\tau\leq\tau_{0}$,

$F_{\tau}(\mu)\geq C_{1}M_{2}(\mu)-C_{2}^{\prime}.$

Combined with the fact that $-H(\mu)\leq
C_{1}^{\prime}M_{2}(\mu)+C_{2}^{\prime}$, we get $-H(\mu)\leq
C_{1}^{\prime\prime}F_{\tau}(\mu)+C_{2}^{\prime\prime}$.∎

See e.g. [Mei et al., 2018, Lem. 10.1] for a proof of the following lemma.

######  Lemma 4.3.

For $\mu\in\mathcal{P}^{a}_{2}(\mathbb{R}^{d})$, let
$M_{2}(\mu)\coloneqq\int\|x\|^{2}\mathrm{d}\mu(x)$. For any $\sigma^{2}>0$, it
holds

$-H(\mu)\leq\frac{1}{\sigma^{2}}M_{2}(\mu)+1+d\log(2\pi\sigma^{2}).$

##  5 Applications and experiments

###  5.1 Noisy GD on a wide two-layer neural network

We now show that our results apply to the training dynamics of certain wide
$2$-layer neural networks trained with noisy gradient descent.

Let us introduce the formulation of two-neural networks of arbitrary width
parameterized by a probability measure, which is at the heart of the mean-
field analysis of the training dynamics [Nitanda and Suzuki, 2017, Mei et al.,
2018, Sirignano and Spiliopoulos, 2020, Rotskoff and Vanden-Eijnden, 2018,
Chizat and Bach, 2018]. Consider a input/output data distribution
$(z,y)\sim\mathcal{D}\in\mathcal{P}(\mathbb{R}^{n}\times\mathbb{R})$, a loss
function $\ell:\mathbb{R}^{2}\to\mathbb{R}_{+}$, a “feature function”
$\Phi(z,x)\in\mathcal{C}(\mathbb{R}^{n}\times\mathbb{R}^{d})$ and let

$G(\mu)\coloneqq\mathbf{E}_{(z,y)\sim\mathcal{D}}\,\ell\Big{(}y,\int\Phi(z,x)\mathrm{d}\mu(x)\Big{)}+\frac{\lambda}{2}\int\|x\|^{2}\mathrm{d}\mu(x)$
(21)

where $\lambda>0$ is regularization parameter. Typical choices for the loss
are the logistic loss $\ell(y,y^{\prime})=\log(1+\exp(-yy^{\prime}))$ and the
square loss $\ell(y,y^{\prime})=\frac{1}{2}|y-y^{\prime}|^{2}$ and in what
follows, $\ell^{\prime}$ denotes the derivative of $\ell$ with respect to
$y^{\prime}$.

When $\mu=\frac{1}{m}\sum_{i=1}^{m}\delta_{x_{i}}$ is an empirical
distribution with $m$ atoms/particles, the function $G_{m}$ derived from $G$
as in Eq. (2) is exactly the risk with weight decay regularization for a two-
layer neural network of width $m$. Thus noisy gradient descent for two-layer
neural networks is equivalent to NPGD with $G$ defined in Eq (21).

Let us give simple conditions under which our convergence theorems apply in
this case.

######  Proposition 5.1.

Assume that $\ell$ is the square or the logistic loss, that $|\Phi|$ is
bounded by $K>0$ and that $\Phi$ smooth in $x$, uniformly in $z$. Then
Assumptions 1 and 2 are satisfied and the first variation of $G$ is given, for
$\mu\in\mathcal{P}_{2}(\mathbb{R}^{d})$ and $x\in\mathbb{R}^{d}$, by

$V[\mu](x)=\mathbf{E}_{(z,y)\sim\mathcal{D}}\,\ell^{\prime}\Big{(}y,\int\Phi(z,x^{\prime})\mathrm{d}\mu(x^{\prime})\Big{)}\Phi(z,x)+\frac{\lambda}{2}\|x\|^{2}.$

Moreover Assumption 3 is satisfied when:

  * –

$\ell$ is the logistic loss. Then we have
$\rho_{\tau}\geq\frac{\lambda}{\tau}e^{-2K/\tau}$, or

  * –

$\ell$ is the square loss and $|\mathbf{E}[y|z]|\leq K^{\prime}$ a.s. Then we
have $\rho_{\tau}\geq\frac{\lambda}{\tau}e^{-2K(K+K^{\prime})/\tau}$.

###### Proof.

For the computation of the first variation and Assumptions 1, we refer e.g. to
Hu et al. [2021]. For Assumption 2, $G$ is convex as a composition of a linear
operator and a convex function. To see that $F_{\tau}$ admits a minimizer,
notice that thanks to the regularization term, the sublevel sets of $G$ are
tight and thus weakly-precompact by Prokorov’s theorem. Moreover, the loss
term in $G$ is weakly continuous, the regularization term is weakly lower-
semicontinuous (lsc) and $H$ is weakly lsc [Ambrosio and Savaré, 2007, Sec.
3.2] so, overall, $F_{\tau}$ is lsc. Thus a minimizer $\mu^{*}_{\tau}$ exists
for all $\tau\geq 0$ by the Direct Method in the calculus of variations.

Let us derive the lower-bound on the log-Sobolev constant $\rho_{\tau}$ using
the criteria given below Assumption 3. First, by the Bakry-Émery criterion,
the probability measure $\propto e^{-\frac{\lambda}{2\tau}\|x\|^{2}}$
satisfies $\mathrm{LSI(\lambda/\tau)}$. Also, our assumptions guarantee that
the first term in $V$ is uniformly bounded by $K$ – in case of the logistic
loss because $|\ell^{\prime}(y,y^{\prime})|\leq 1$ – or by $K(K+K^{\prime})$ –
in case of the square loss. We conclude by applying the Holley-Stroock
criterion with a perturbation $\psi$ that satisfies $\sup\psi-\inf\psi\leq
2K/\tau$ (for the logistic loss) or $\sup\psi-\inf\psi\leq
2K(K+K^{\prime})/\tau$ (for the square loss). ∎

##### Limitations of this approach

While the previous proposition, combined with our theorems, gives new
convergence guarantees for noisy gradient descent on neural networks (in a
certain limit), let us stress on the limitations of these results. First, the
convergence proof fundamentally relies on the existence of noise and our
analysis misses the fact that, in certain contexts, $G$ has a specific
structure (of a different nature than the uniform LSI) that in practice seems
to ease convergence and even make the noiseless dynamics converge, see e.g.
Bach and Chizat [2021]. Second, this approach introduces, in addition to
weight decay, an entropic regularization which might be detrimental to the
statistical performance. Finally, the risk for a vanilla two-layer neural
network with non-linearity $\phi:\mathbb{R}\to\mathbb{R}$ is obtained from Eq.
(21) by taking

$\displaystyle\Phi(z,x)=a\phi(b^{\top}z)\quad\text{where}\quad
x=(a,b)\in\mathbb{R}\times\mathbb{R}^{d-1}$ (22)

which is not covered by our assumptions (in particular because it is not
bounded). In the case of the ReLU non-linearity $\phi(s)=\max\\{0,s\\}$, there
is in addition a lack of smoothness issue. An interesting direction for future
research would be to adapt this algorithm and analysis in order to cover the
case of ReLU non-linearities.

##### Relaxing the boundedness assumption

In Proposition 5.1, we assumed that $\Phi$ is bounded in order to obtain easy
quantitative bounds on the LSI constant, but this assumption is not necessary.
For instance, the Lyapunov condition in [Cattiaux et al., 2010, Cor. 2.1]
implies that if there exists $C_{1},C_{2},C_{3}>0$ such that for all
$\mu\in\mathcal{P}_{2}(\mathbb{R}^{d})$,
$V[\mu]\in\mathcal{C}^{2}(\mathbb{R}^{d})$ and

$\displaystyle\forall\|x\|_{2}\geq C_{1},\quad x^{\top}\nabla V[\mu](x)\geq
C_{2}\|x\|^{2}_{2}\quad\text{and}\quad\forall
x\in\mathbb{R}^{d},\quad\nabla^{2}V[\mu](x)\succeq-C_{3}\mathrm{Id}$ (23)

then a uniform LSI holds (i.e. Assumption 3 holds).

As an illustration, consider a positively $2$-homogeneous and
$\mathcal{C}^{2}$ (uniformly in $z$) function $x\mapsto\Phi(x,z)$ (this does
not cover functions of the form Eq. (22) except in the simple case
$\phi=\mathrm{id}$; a valid example is the square non-linearity
$\Phi(x,z)=(x^{\top}z)^{2}$). Given the expression of $V$ in Proposition 5.1
and invoking Euler’s identity $x^{\top}\nabla\Phi(x,z)=2\Phi(x,z)$, we have

$x^{\top}\nabla
V[\mu](x)=2\mathbf{E}_{(z,y)\sim\mathcal{D}}\,\ell^{\prime}\Big{(}y,\int\Phi(z,x^{\prime})\mathrm{d}\mu(x^{\prime})\Big{)}\Phi(z,x)+\lambda\|x\|^{2}.$

So, assuming $\ell^{\prime}$ is absolutely bounded by some $M>0$, the first
part of Eq. (23) holds for $\lambda$ large enough, namely
$\lambda>2M\mathbf{E}|\Phi(z,u)|$ for all $u\in\SS_{d-1}$. Moreover, the
Hessian lower-bound is also satisfied since $\nabla^{2}\Phi$ is
$0$-homogeneous. Thus, in this simple example, Assumption 3 holds although
$\Phi$ is unbounded with a quadratic growth (but this is at the price of
requiring a strong weight decay regularization).

###  5.2 Numerical illustration: kernel Maximum Mean Discrepancy

We conclude this paper with numerical experiments exploring the behavior of
NPGD111Link to Julia code to reproduce the experiments:
https://github.com/lchizat/2022-mean-field-langevin-rate. defined in (3). Let
us stress that our theoretical guarantees only apply to the mean-field
Langevin dynamics – recovered in the many-particle and continuous time limit –
so there remains a gap between the theory and the NPGD algorithm.

We consider the torus $\mathcal{X}\coloneqq(\mathbb{R}/(2\pi\mathbb{Z}))^{d}$
and the convex function defined on $\mathcal{P}(\mathcal{X})$ by

$\displaystyle G(\mu)\coloneqq\frac{1}{2}\int
k(x,y)\mathrm{d}\mu(x)\mathrm{d}\mu(y)-\int
k(x,y)\mathrm{d}\mu(x)\mathrm{d}\nu(y)+\frac{1}{2}\int
k(x,y)\mathrm{d}\nu(x)\mathrm{d}\nu(y)$ (24)

where $k\in\mathcal{C}^{2}(\mathcal{X}\times\mathcal{X})$ is a smooth positive
semi-definite kernel and $\nu\in\mathcal{P}(\mathcal{X})$ a fixed probability
measure. This function $G$ can be interpreted as the square kernel Maximum
Mean Discrepancy (kMMD) [Gretton et al., 2008] between $\mu$ and $\nu$. This
choice of function is convenient for numerical experiments because its minimum
value is known and is $0$, attained in particular for $\mu=\nu$. Although our
theory was developed for $\mathcal{X}=\mathbb{R}^{d}$, it is straightforward
to adapt it to the torus, and our main convergence results apply, as shown
below.

######  Proposition 5.2.

The first variation of $G$ is given, for $\mu\in\mathcal{P}(\mathcal{X})$ and
$x\in\mathcal{X}$ by

$V[\mu](x)=\int_{\mathcal{X}}k(x,y)\mathrm{d}(\mu-\nu)(y).$

Moreover, Assumptions 1,2 and 3 are satisfied with
$\rho_{\tau}\geq((1+2\pi)d)^{-1}e^{(\inf k-\sup k)/\tau}$.

###### Proof.

The properties of $G$ and $V$ are obtained by standard arguments. For the log-
Sobolev inequality, we note that the normalized volume measure on
$\mathcal{X}$ satifies LSI with a constant larger than
$\frac{\pi^{2}}{(1+2\pi)\mathrm{diam}(\mathcal{X})^{2}}$ [Ledoux, 1999, Thm.
7.3] and here $\mathrm{diam}(\mathcal{X})=\pi\sqrt{d}$. The lower bound on
$\rho_{\tau}$ follows by the Holley-Stroock criterion. ∎

In our experiments, we consider $d=2$ and the translation invariant kernel
$k(x,y)=\prod_{i=1}^{d}\big{(}1+2\sum_{k=1}^{n}(1+k)^{-1}\cos(k(x_{i}-y_{i}))\big{)}$
with $n=5$ frequency components. Because of the frequency cut-off, this kernel
is not strictly positive definite, so $G$ admits minimizers other than $\nu$
[Sriperumbudur et al., 2010, Simon-Gabriel et al., 2020] although in practice
we observed that NPGD was in general attracted towards $\nu$. We take $\nu$ as
a random empirical distribution of $m^{*}=10$ samples from the uniform
distribution on $\mathcal{X}$. We run NPGD with $m=50$ particles, a step-size
$\eta=0.08$ and $\mu_{0}$ being the uniform distribution on $\mathcal{X}$.

Figure 1(a) shows an example of a large-time particle configuration, with the
atoms of $\nu$ is red and the atoms of $\hat{\mu}_{t}$ in black (with $t$
large), with a noise temperature $\tau=0.1$. Here the measure $\hat{\mu}_{t}$
is a noisy version of $\nu^{*}$.

Figure 1(b) shows the evolution of the objective $F_{\tau}=G+\tau H$ (up to a
constant, adjusted for ease of comparison) along the iterations, where the
entropy $H$ is estimated using the $1$-nearest-neighbor estimator [Kozachenko
and Leonenko, 1987, Singh et al., 2003]. We observe the exponential decay of
$F_{\tau}$ towards a plateau which we expect to be the global minimum of
$F_{\tau}$, up to discretization errors. For small values of $\tau$, it is not
excluded that the plateau corresponds instead to a suboptimal metastable
state.

Finally, Figure 1(c) shows the advantage of NPGD with simulated annealing vs.
PGD to minimize the unregularized function $G$. We used a noise temperature
that decays polynomially as $\tau_{t}=20(t+1)^{-1}$ where $t$ is the iteration
count, which is a faster decay than what the theory suggests. At iteration
$800$, we stopped the noise in order to observe the “quality” of the
configuration of particles. We see that the NPGD with simulated annealing
consistently outperforms PGD, which gets stuck in poorer local minima.

(a) Example of a large-time configuration of NPGD for $\tau=0.1$. (red) atoms
of $\nu$ (black) atoms of $\hat{\mu}_{t}$.

(b) Evolution of $F_{\tau}(\mu_{t})$ (averaged over $10$ random experiments).
Curves’ height adjusted to end at $1$.

(c) Evolution of $G(\mu_{t})$ for NPGD with simulated annealing vs. PGD
(averaged over $10$ experiments).

##  6 Conclusion

We have proved the convergence of the Mean-Field Langevin dynamics to the
global minimizer at an exponential rate, under natural assumptions that
include all settings where (non-quantitative) convergence was previously
shown. We have also proved the convergence of the annealed dynamics for a
suitable noise decay.

From a higher perspective, our analysis—in particular the simple “entropy
sandwich” Lemma 3.4—suggests that often, the guarantees about Langevin
dynamics obtained via log-Sobolev inequalities can be generalized to _mean-
field_ Langevin dynamics. In this paper, we focused on exponential convergence
and on simulated annealing, but other aspects could be considered, such as a
direct analysis of the discrete dynamics, which could lead to computational
bounds, as done in e.g. [Vempala and Wibisono, 2019, Ma et al., 2019] for the
Langevin algorithm.

Another interesting direction for future work is to develop and study more
applications of Mean-Field Langevin dynamics, since many problems can be cast
as optimization problems of the form Eq. (1). This includes sparse
deconvolution problems, mixture models fitting [Boyd et al., 2017] or problems
involving optimal transport [Peyré and Cuturi, 2019, Chap. 9].

#### Acknowledgments

I would like to thank Loucas Pillaud-Vivien for orienting me through the
literature of simulated annealing and Mo Zhou for noticing a gap in a previous
version of the proof of Theorem 4.1. I would also like to thank the anonymous
reviewers for their insightful comments and for suggesting the discussion on
LSI in Section 5.1.

## References

  * Ambrosio and Savaré [2007] Luigi Ambrosio and Giuseppe Savaré.  Gradient flows of probability measures.  _Handbook of differential equations: evolutionary equations_ , 3:1–136, 2007.
  * Bach and Chizat [2021] Francis Bach and Lénaïc Chizat.  Gradient descent on infinitely wide neural networks: Global convergence and generalization.  _arXiv preprint arXiv:2110.08084_ , 2021.
  * Bakry and Émery [1985] Dominique Bakry and Michel Émery.  Diffusions hypercontractives.  In _Seminaire de probabilités XIX 1983/84_ , pages 177–206. Springer, 1985.
  * Bakry et al. [2014] Dominique Bakry, Ivan Gentil, and Michel Ledoux.  _Analysis and geometry of Markov diffusion operators_ , volume 103\.  Springer, 2014.
  * Boyd et al. [2017] Nicholas Boyd, Geoffrey Schiebinger, and Benjamin Recht.  The alternating descent conditional gradient method for sparse inverse problems.  _SIAM Journal on Optimization_ , 27(2):616–639, 2017.
  * Cattiaux et al. [2010] Patrick Cattiaux, Arnaud Guillin, and Li-Ming Wu.  A note on Talagrand’s transportation inequality and logarithmic Sobolev inequality.  _Probability theory and related fields_ , 148(1):285–304, 2010.
  * Chen et al. [2021] Hong-Bin Chen, Sinho Chewi, and Jonathan Niles-Weed.  Dimension-free log-Sobolev inequalities for mixture distributions.  _Journal of Functional Analysis_ , 281(11):109236, 2021.
  * Chizat and Bach [2018] Lénaïc Chizat and Francis Bach.  On the global convergence of gradient descent for over-parameterized models using optimal transport.  _Advances in Neural Information Processing Systems_ , 31:3036–3046, 2018.
  * Dobrushin [1979] R. L. Dobrushin.  Vlasov equations.  _Functional Analysis and its Applications_ , 13(2):115–123, 1979.
  * Durmus and Moulines [2017] Alain Durmus and Eric Moulines.  Nonasymptotic convergence analysis for the unadjusted Langevin algorithm.  _The Annals of Applied Probability_ , 27(3):1551–1587, 2017.
  * Eberle et al. [2019] Andreas Eberle, Arnaud Guillin, and Raphael Zimmer.  Quantitative Harris-type theorems for diffusions and McKean–Vlasov processes.  _Transactions of the American Mathematical Society_ , 371(10):7135–7173, 2019.
  * Ermak [1975] Donald L. Ermak.  A computer simulation of charged particles in solution. I. Technique and equilibrium properties.  _The Journal of Chemical Physics_ , 62(10):4189–4196, 1975.
  * Ferreira and Valencia-Guevara [2018] Lucas C.F. Ferreira and Julio C. Valencia-Guevara.  Gradient flows of time-dependent functionals in metric spaces and applications to PDEs.  _Monatshefte für Mathematik_ , 185(2):231–268, 2018.
  * Geman and Hwang [1986] Stuart Geman and Chii-Ruey Hwang.  Diffusions for global optimization.  _SIAM Journal on Control and Optimization_ , 24(5):1031–1043, 1986.
  * Gretton et al. [2008] Arthur Gretton, Karsten M Borgwardt, Malte J. Rasch, Bernhard Schölkopf, and Alexander Smola.  A kernel method for the two-sample problem.  _Journal of Machine Learning Research_ , 1:1–10, 2008.
  * Holley and Stroock [1987] Richard Holley and Daniel Stroock.  Logarithmic Sobolev inequalities and stochastic Ising models.  _Journal of Statistical Physics_ , 46(5-6):1159–1194, 1987.
  * Holley et al. [1989] Richard A. Holley, Shigeo Kusuoka, and Daniel W. Stroock.  Asymptotics of the spectral gap with applications to the theory of simulated annealing.  _Journal of functional analysis_ , 83(2):333–347, 1989.
  * Hu et al. [2021] Kaitong Hu, Zhenjie Ren, David Šiška, and Łukasz Szpruch.  Mean-field Langevin dynamics and energy landscape of neural networks.  _Annales de l’Institut Henri Poincaré, Probabilités et Statistiques_ , 57(4):2043 – 2065, 2021.  doi: 10.1214/20-AIHP1140.  URL https://doi.org/10.1214/20-AIHP1140.
  * Huang et al. [2021] Xing Huang, Panpan Ren, and Feng-Yu Wang.  Distribution dependent stochastic differential equations.  _Frontiers of Mathematics in China_ , 16(2):257–301, 2021.
  * Kozachenko and Leonenko [1987] Lyudmyla F. Kozachenko and Nikolai N. Leonenko.  Sample estimate of the entropy of a random vector.  _Problemy Peredachi Informatsii_ , 23(2):9–16, 1987.
  * Lacker [2018] Daniel Lacker.  Mean field games and interacting particle systems.  _Preprint_ , 2018.
  * Lacker [2021] Daniel Lacker.  Hierarchies, entropy, and quantitative propagation of chaos for mean field diffusions.  _arXiv preprint arXiv:2105.02983_ , 2021.
  * Ledoux [1999] Michel Ledoux.  Concentration of measure and logarithmic Sobolev inequalities.  In _Seminaire de probabilites XXXIII_ , pages 120–216. Springer, 1999\.
  * Ma et al. [2019] Yi-An Ma, Yuansi Chen, Chi Jin, Nicolas Flammarion, and Michael Jordan.  Sampling can be faster than optimization.  _Proceedings of the National Academy of Sciences_ , 116(42):20881–20885, 2019.
  * Mei et al. [2018] Song Mei, Andrea Montanari, and Phan-Minh Nguyen.  A mean field view of the landscape of two-layer neural networks.  _Proceedings of the National Academy of Sciences_ , 115(33):E7665–E7671, 2018.
  * Mei et al. [2019] Song Mei, Theodor Misiakiewicz, and Andrea Montanari.  Mean-field theory of two-layers neural networks: dimension-free bounds and kernel limit.  In _Conference on Learning Theory_ , pages 2388–2464. PMLR, 2019\.
  * Miclo [1992] Laurent Miclo.  Recuit simulé sur $\mathbb{R}^{n}$. étude de l’évolution de l’énergie libre.  In _Annales de l’IHP Probabilités et statistiques_ , volume 28, pages 235–266, 1992.
  * Nitanda and Suzuki [2017] Atsushi Nitanda and Taiji Suzuki.  Stochastic particle gradient descent for infinite ensembles.  _arXiv preprint arXiv:1712.05438_ , 2017.
  * Nitanda et al. [2021] Atsushi Nitanda, Denny Wu, and Taiji Suzuki.  Particle dual averaging: Optimization of mean field neural networks with global convergence rate analysis.  In _Neural Information Processing Systems_ , 2021.
  * Nitanda et al. [2022] Atsushi Nitanda, Denny Wu, and Taiji Suzuki.  Convex analysis of the mean field langevin dynamics.  In _International Conference on Artificial Intelligence and Statistics_ , pages 9741–9757. PMLR, 2022.
  * Otto and Villani [2000] Felix Otto and Cédric Villani.  Generalization of an inequality by Talagrand and links with the logarithmic Sobolev inequality.  _Journal of Functional Analysis_ , 173(2):361–400, 2000.
  * Peyré and Cuturi [2019] Gabriel Peyré and Marco Cuturi.  Computational optimal transport: With applications to data science.  _Foundations and Trends® in Machine Learning_ , 11(5-6):355–607, 2019.
  * Raginsky et al. [2017] Maxim Raginsky, Alexander Rakhlin, and Matus Telgarsky.  Non-convex learning via stochastic gradient Langevin dynamics: a nonasymptotic analysis.  In _Conference on Learning Theory_ , pages 1674–1703. PMLR, 2017\.
  * Roberts and Tweedie [1996] Gareth O. Roberts and Richard L. Tweedie.  Exponential convergence of Langevin distributions and their discrete approximations.  _Bernoulli_ , pages 341–363, 1996.
  * Rotskoff and Vanden-Eijnden [2018] Grant M. Rotskoff and Eric Vanden-Eijnden.  Parameters as interacting particles: long time convergence and asymptotic error scaling of neural networks.  In _Proceedings of the 32nd International Conference on Neural Information Processing Systems_ , pages 7146–7155, 2018.
  * Santambrogio [2015] Filippo Santambrogio.  Optimal transport for applied mathematicians.  _Birkäuser, NY_ , 55(58-63):94, 2015.
  * Simon-Gabriel et al. [2020] Carl-Johann Simon-Gabriel, Alessandro Barp, Bernhard Schölkopf, and Lester Mackey.  Metrizing weak convergence with Maximum Mean Discrepancies.  _arXiv preprint arXiv:2006.09268_ , 2020.
  * Singh et al. [2003] Harshinder Singh, Neeraj Misra, Vladimir Hnizdo, Adam Fedorowicz, and Eugene Demchuk.  Nearest neighbor estimates of entropy.  _American journal of mathematical and management sciences_ , 23(3-4):301–321, 2003.
  * Sirignano and Spiliopoulos [2020] Justin Sirignano and Konstantinos Spiliopoulos.  Mean field analysis of neural networks: A law of large numbers.  _SIAM Journal on Applied Mathematics_ , 80(2):725–752, 2020.
  * Sriperumbudur et al. [2010] Bharath K. Sriperumbudur, Arthur Gretton, Kenji Fukumizu, Bernhard Schölkopf, and Gert R.G. Lanckriet.  Hilbert space embeddings and metrics on probability measures.  _The Journal of Machine Learning Research_ , 11:1517–1561, 2010.
  * Sznitman [1991] Alain-Sol Sznitman.  Topics in propagation of chaos.  In _Ecole d’été de probabilités de Saint-Flour XIX—1989_ , pages 165–251. Springer, 1991.
  * Tang and Zhou [2021] Wenpin Tang and Xun Yu Zhou.  Simulated annealing from continuum to discretization: a convergence analysis via the Eyring–Kramers law.  _arXiv preprint arXiv:2102.02339_ , 2021.
  * Tseng [2010] Paul Tseng.  Approximation accuracy, gradient methods, and error bound for structured convex optimization.  _Mathematical Programming_ , 125(2):263–295, 2010\.
  * Vempala and Wibisono [2019] Santosh Vempala and Andre Wibisono.  Rapid convergence of the Unadjusted Langevin Algorithm: Isoperimetry suffices.  _Advances in Neural Information Processing Systems_ , 32:8094–8106, 2019.
  * Villani [2021] Cédric Villani.  _Topics in optimal transportation_ , volume 58.  American Mathematical Soc., 2021.
  * Wang [2001] Feng-Yu Wang.  Logarithmic Sobolev inequalities: conditions and counterexamples.  _Journal of Operator Theory_ , pages 183–197, 2001.

##  Appendix A Additional proofs

Let us start with a relation between $G$ and its first-variation $V$ that is
more convenient for proofs.

######  Lemma A.1 (Integral formula).

Under Assumption (1), for $\mu_{0},\mu_{1}\in\mathcal{P}_{2}(\mathbb{R}^{d})$,
one has

$G(\mu_{1})-G(\mu_{0})=\int_{0}^{1}\int_{\mathbb{R}^{d}}V[\mu_{t}]\mathrm{d}(\mu_{1}-\mu_{0})\mathrm{d}t$

where $\mu_{t}=(1-t)\mu_{0}+t\mu_{1}$ for $t\in[0,1]$.

###### Proof.

Let $h(t)=G(\mu_{t})$. By definition of the first-variation, $h$ is right
(resp. left) continuous at $t=0$ (resp. $t=1$). We just need to prove that $h$
is differentiable on $]0,1[$ with $h^{\prime}(t)=\int
V[\mu_{t}]\mathrm{d}(\mu_{1}-\mu_{0})$. Then, because this expression is
continuous in $t$ under Assumption 1, the fundamental theorem of calculus
would imply $h(1)-h(0)=\int_{0}^{1}h^{\prime}(t)\mathrm{d}t$, which is our
claim. For $t,\epsilon\in{]0,1[}$ one has
$(1-\epsilon)\mu_{t}+\epsilon\mu_{0}=\mu_{t-t\epsilon}$ and thus

$-th^{\prime}_{-}(t)=\lim_{\epsilon\to
0_{+}}\frac{h(t-t\epsilon)-h(t)}{\epsilon}=\lim_{\epsilon\to
0+}\frac{G((1-\epsilon)\mu_{t}+\epsilon\mu_{0})-G(\mu_{t})}{\epsilon}\\\ =\int
V[\mu_{t}]\mathrm{d}(\mu_{0}-\mu_{t})=-t\int
V[\mu_{t}]\mathrm{d}(\mu_{1}-\mu_{0})$

where $h^{\prime}_{-}(t)$ stands for the left-derivative of $h$ at $t$. This
shows that $h^{\prime}_{-}(t)=\int V[\mu_{t}]\mathrm{d}(\mu_{1}-\mu_{0})$ for
$t\in{]0,1[}$. A similar computation using
$(1-\epsilon)\mu_{t}+\epsilon\mu_{1}=\mu_{t+(1-t)\epsilon}$ shows that the
right derivative $h^{\prime}_{+}(t)$ has the same value, and thus
$h^{\prime}(t)=\int V[\mu_{t}]\mathrm{d}(\mu_{1}-\mu_{0})$ for $t\in[0,1]$
which concludes the proof of the formula. ∎

In the following lemma, we verify that $G$ is well-behaved (in fact smooth) as
function in the Wasserstein space $\mathcal{P}_{2}(\mathbb{R}^{d})$, using the
vocabulary and results from Ambrosio and Savaré [2007].

######  Lemma A.2.

Let $\mu\in\mathcal{P}_{2}(\mathbb{R}^{d})$, let $v\in L^{2}(\mu)$ and let
$\mu_{t}=(\mathrm{id}+tv)_{\\#}\mu$. Then

$\frac{d}{dt}G(\mu_{t})=\int_{\mathbb{R}^{d}}\nabla
V[\mu_{t}](x+tv(x))^{\top}v(x)\mathrm{d}\mu(x).$

Moreover, $G$ is $(-2L)$-semiconvex along any interpolating curve in
$\mathcal{P}_{2}^{a}(\mathbb{R}^{d})$ and the $W_{2}$-derivative of $G$ at
$\mu$ is $V[\mu]$.

Since the same holds true for $-G$, we could say that $G$ is $2L$-smooth in
the Wasserstein geometry, in the sense that it is both $2L$-semiconvex and
$2L$-semiconvex.

###### Proof.

For $\epsilon>0$ and $s\in[0,1]$, let $\mu_{s}=(1-s)\mu+s\mu_{t+\epsilon}$. It
holds by Lemma A.1,

$\displaystyle\frac{1}{\epsilon}(G(\mu_{t+\epsilon})-G(\mu_{t}))$
$\displaystyle=\frac{1}{\epsilon}\int_{0}^{1}\int_{\mathbb{R}^{d}}V[\mu_{s}]\mathrm{d}(\mu_{t+\epsilon}-\mu_{t})$
$\displaystyle=\frac{1}{\epsilon}\int_{0}^{1}\int_{\mathbb{R}^{d}}(V[\mu_{s}](x+(t+\epsilon)v(x))-V[\mu_{s}](x+tv(x)))\mathrm{d}\mu(x)$
$\displaystyle=\int_{0}^{1}\int_{\mathbb{R}^{d}}\nabla
V[\mu_{s}](x+tv(x))^{\top}v(x)\mathrm{d}\mu(x)+O(L\epsilon\|v\|_{L^{2}(\mu)})$
$\displaystyle=\int_{\mathbb{R}^{d}}\nabla
V[\mu](x+tv(x))^{\top}v(x)\mathrm{d}\mu(x)+O(L\epsilon\|v\|_{L^{2}(\mu)})$

where we used successively the Lipschitz continuity of $x\mapsto\nabla
V[\mu](x)$ and of $\mu\mapsto\nabla V[\mu](x)$ in the last two lines. The
first claim follows by taking the limit $\epsilon\to 0$. This also shows that
$V[\mu]$ is the unique (strong) $W_{2}$-differential of $W_{2}$ at $\mu$, in
the sense of [Ambrosio and Savaré, 2007, Def. 4.1].

For the semi-convexity claim, let $h(t)\coloneqq G(\mu_{t})$. For
$s,t\in[0,1]$, it holds by Cauchy-Schwarz

$\displaystyle|h^{\prime}(t)-h^{\prime}(s)|^{2}$
$\displaystyle\leq\|v\|^{2}_{L^{2}(\mu)}\int_{\mathbb{R}^{d}}\|\nabla
V[\mu_{t}](x+tv(x))-\nabla V[\mu_{s}](x+sv(x))\|^{2}\mathrm{d}\mu(x)$
$\displaystyle\leq\|v\|^{2}_{L^{2}(\mu)}L^{2}\int_{\mathbb{R}^{d}}\big{(}W_{2}(\mu_{s},\mu_{t})+|t-s|\|v(x)\|\big{)}^{2}\mathrm{d}\mu(x)$
$\displaystyle\leq\|v\|^{2}_{L^{2}(\mu)}L^{2}(2W_{2}^{2}(\mu_{s},\mu_{t})+2|t-s|^{2}\|v\|^{2}_{L^{2}(\mu)}).$

Since $W_{2}(\mu_{s},\mu_{t})\leq|t-s|\|v\|_{L^{2}(\mu)}$, it follows

$\displaystyle|h^{\prime}(t)-h^{\prime}(s)|\leq 2L|t-s|\|v\|^{2}_{L^{2}(\mu)}$

which proves that $G$ is $(-2L)$-convex in the sense of [Ambrosio and Savaré,
2007, Remark 3.2]. ∎


11institutetext: Institute of Theoretical Astrophysics, Universitetet i Oslo,
P.O. Box 1029 Blindern, N-0315 Oslo, Norway 22institutetext: Department of
Physics, McGill University, 3600 rue University, Montreal QC H3A 2T8, Canada
33institutetext: Institute for Computational Science, Universität Zürich,
Winterthurerstrasse 190, 8057 Zürich, Switzerland 44institutetext: Département
de Physique Théorique, Université de Genève, 24 quai Ernest-Ansermet, 1211
Genève 4, Switzerland

# Biased tracers as a probe of beyond-$\Lambda$CDM cosmologies

Farbod Hassani , 1122<EMAIL_ADDRESS>Julian Adamek 33 Ruth Durrer
44 Martin Kunz 44

Cosmological models beyond $\Lambda$CDM, such as those featuring massive
neutrinos or modifications of gravity, often display a characteristic change
(scale-dependent suppression or enhancement) in the matter power spectrum when
compared to a six-parameter $\Lambda$CDM baseline. It is therefore a widely
held view that constraints on those models can be obtained by searching for
such features in the clustering statistics of large-scale structure. However,
when using biased tracers of matter in the analysis, the situation is
complicated by the fact that the bias also depends on cosmology. Here we
investigate how the selection of tracers affects the observed signatures for
two examples of beyond-$\Lambda$CDM cosmologies: massive neutrinos and
clustering dark energy ($k$-essence). We study the signatures in the monopole,
quadrupole, and hexadecapole of the redshift-space power spectra for halo
catalogues from large $N$-body simulations and argue that a fixed selection
criterion based on local attributes, such as tracer mass, leads to a near loss
of signal in most cases. Instead, the full signal is recovered only if the
selection of tracers is done at fixed bias. This emphasises the need to model
or measure the bias parameters accurately in order to get meaningful
constraints on the cosmological model.

###### Key Words.:

Cosmology: large-scale structure of Universe – Methods: $N$-body simulations,
statistical

##  1 Introduction

The six-parameter $\Lambda$-cold-dark-matter ($\Lambda$CDM) model of cosmology
fits most current cosmological observations (Riess et al. 1998; Perlmutter et
al. 1999; Alam et al. 2017; Aghanim et al. 2020) but it relies on two
unexplained ingredients, dark matter (DM) and dark energy ($\Lambda$), which
make up about $95\%$ of the energy budget of the Universe. For this reason,
cosmologists are investigating alternative models of dark energy and DM and
also theories of gravity which deviate from general relativity on large
scales. In the coming years, the next generation of large-scale structure
(LSS) surveys (Amendola et al. 2018; Laureijs et al. 2011; Abell et al. 2009)
will reach an extraordinary statistical power that will allow us to
discriminate between different theories and possibly reveal the nature of the
dark components. These surveys will map the three-dimensional galaxy
distribution up to redshift $z\simeq 2$. Under the assumption that galaxies
and their host halos trace the underlying DM distribution, we can measure the
three-dimensional mass distribution in the Universe.

To keep up with these future observations we need sufficiently accurate
theoretical predictions for the models we wish to study. While perturbation
theory is sufficient to model the anisotropies of the cosmic microwave
background (CMB), the galaxy distribution is more complicated and must be
modelled with cosmological $N$-body simulations that also provide access to a
large number of modes in the non-linear regime (Angulo & Hahn 2021). As a
result, $N$-body simulations of so-called non-standard cosmologies have been
developed over recent years (Baldi et al. 2010; Barreira et al. 2013; Llinares
et al. 2014; Adamek et al. 2017; Li 2018; Hassani et al. 2019; Hassani &
Lombriser 2020). While far from exhausting all possibilities, these help us to
place stringent constraints on modified theories of gravity or models of dark
energy and DM by comparing the details of structure formation in the different
scenarios.

At the most basic level, such a comparison often starts with the matter power
spectrum. The different scenarios typically lead to a scale-dependent and
redshift-dependent modification of the power spectrum when compared to a
$\Lambda$CDM baseline, and this signature can then be used to constrain the
models. When applying this reasoning to observed clustering in LSS, which is
always measured from biased tracers of matter such as galaxies, one needs to
keep in mind that the bias of a fixed type of tracer is not independent of
cosmology either. In the simplest halo bias model (Kaiser 1984), the linear
bias parameter for a fixed halo mass is approximately inversely proportional
to the amplitude of fluctuations. This is because it becomes proportionally
more difficult to form objects of a given mass when the initial fluctuations
are smaller. Consequently, for a fixed type of halo defined through a mass
threshold, a change in matter power is expected to be partially compensated by
a change in bias, such that the clustering of the tracer is less affected by a
change in cosmology.

Here, we present a study of this issue in the context of two
beyond-$\Lambda$CDM scenarios: a cosmology with massive neutrinos and a
cosmology with $k$-essence-type clustering dark energy. In the first case, the
matter power gets suppressed on small scales due to the free streaming of
neutrinos that make up a certain fraction of the matter as determined by their
mass. In the second case, the growth rate gets modified at all scales due to a
change in the expansion history of the late Universe. Additionally, the
clustering of dark energy leaves a certain scale-dependent imprint. We use
large $N$-body simulations to produce halo catalogues for the different
scenarios and show that the model-specific signatures get almost completely
wiped out in the halo power spectra if the comparison is carried out using
halo populations with a fixed mass threshold. We further show, rather
trivially, that the expected signal is recovered if the comparison is instead
carried out at fixed bias value. What is perhaps less trivial is that these
statements also apply for the multipoles of the power spectra in redshift
space, and we demonstrate this for the monopole, quadrupole, and hexadecapole.

##  2 Theory

We define the halo bias at a given wave number as

$b(k)=b_{1}+b_{\mathrm{NL}}(k)=\sqrt{\frac{P_{hh}(k)}{P_{mm}(k)}}\,,$ (1)

where $P_{hh}(k)$ and $P_{mm}(k)$ are, respectively, the halo and matter
density power spectra without redshift-space distortions (RSDs), and $b_{1}$
is the linear bias parameter that captures the bias at large scales where
matter and halo distributions are maximally correlated. At these scales, the
non-linear bias $b_{\mathrm{NL}}(k)\to 0$. The bias $b(k)$ is a measure of how
well the halo density field traces the underlying matter density at a given
scale; see Desjacques et al. (2018) for a review. Alternative definitions
would be $b(k)=P_{mh}(k)/P_{mm}(k)$ or $b(k)=P_{hh}(k)/P_{mh}(k)$, where
$P_{mh}(k)$ is the matter-halo cross power spectrum. We note that in order to
suppress the shot noise in our halo power spectra we employ a simple jack-
knife method as described for example in Inman et al. (2015). In Fig. 1 we
show the bias as a function of wavenumber obtained from these different
definitions for a $\Lambda$CDM simulation. While the non-linear bias depends
on the definition because of the non-trivial correlation between matter and
halos on small scales, the linear bias is almost the same within statistical
error bars, for all definitions. In our simulations, we measure the linear
bias by taking averages over $b(k)$ on modes in the linear and quasi-linear
regime.

Figure 1:  Halo bias obtained at redshift $z=0$ from three different
definitions for a $\Lambda$CDM simulation with box size $4032$ Mpc/$h$ and
$N_{p}=4608^{3}$ particles. The bias is computed for halos with mass
$M_{h}\geq 2\times 10^{13}$${M}_{\odot}$/$h$. Bands show the confidence
intervals for the linear bias $b_{1}$.

We consider two different selection procedures for the halo populations. In
the first case, the population is selected according to a fixed mass
threshold. We should note that our particle-mesh $N$-body simulations have a
fixed effective force resolution (no adaptive mesh refinement) and the
measured mass of the halos is not expected to be fully converged numerically.
The quoted value of the mass threshold should therefore not be taken too
literally, but rather considered a mass proxy. The important point is that the
force resolution was not changed between runs of different cosmologies, and
hence the mass proxies can be considered comparable. This selection method
therefore selects, in some specific sense, the same types of objects across
different models.

In the second case, we want to select the halo populations at a fixed value of
the linear bias. This is achieved by measuring the bias for several choices of
mass threshold, and then fitting for the threshold value that provides the
desired bias. This selection procedure guarantees that any modifications of
the matter power spectrum will appear in the corresponding halo power spectrum
on all scales where the linear bias model is accurate. We note that carrying
out a similar selection in a catalogue from an actual LSS survey would require
some procedure for measuring the bias from observations.

##  3 Non-standard cosmologies

The two non-standard scenarios we consider in this work are massive neutrino
cosmology and clustering dark energy ($k$-essence). In both cases, we use
existing simulation suites that include suitable $\Lambda$CDM reference runs
such that the impact of the non-standard scenarios on the clustering
statistics can be quantified. To minimise the contamination from cosmic
variance, the simulations are initialised on the same realisation of the
Gaussian random field that sets the initial perturbations.

### Massive neutrinos

Laboratory measurements of neutrino flavour oscillations show that neutrinos
have mass (Esteban et al. 2019). However, the absolute neutrino mass scale has
not yet been determined, but only mass differences are known. It is expected
that we will obtain the strongest constraints on the sum of neutrino masses
through cosmological probes (Lesgourgues & Pastor 2006). This motivates
cosmologists to study the impact of massive neutrinos on LSS through $N$-body
simulations. In this work we use a suite of large $N$-body simulations from
Adamek et al. (2017) that include neutrinos with the sum of masses given by
$0$, $0.06$, $0.2,$ and $0.3$ eV. These simulations have a box size of $2048$
Mpc/$h$ with $4096^{3}$ particles, corresponding to a spatial resolution of
$0.5$ Mpc/$h$ and mass resolution of about $10^{10}$${M}_{\odot}$/$h$. The
$\Lambda$CDM baseline model in this case is the one with vanishing neutrino
masses.

###  $k$-essence dark energy

The $k$-essence model is a viable candidate for the late-time accelerated
expansion of the Universe. This model was first introduced in Armendariz-Picon
et al. (2001) to avoid problems such as fine tuning, coincidence, or anthropic
reasoning. The expansion history is slightly different from $\Lambda$CDM at
low redshift due to the fact that $k$-essence has an effective equation of
state of $w\neq-1$. Moreover, the $k$-essence field clusters around matter
over-densities depending on its speed of sound. This model was studied
thoroughly by Hassani et al. (2019, 2020b, 2020a) using cosmological $N$-body
simulations based on the code $k-$evolution. Here we use two of these
simulations with a fixed equation of state parameter $w=-0.9$ and two choices
of the speed of sound, $c_{s}^{2}=1$ and $c_{s}^{2}=10^{-4}$. The simulations
have a box size of $4032$ Mpc/$h$ with $4608^{3}$ particles, corresponding to
a spatial resolution of $0.875$ Mpc/$h$ and mass resolution of about $6\times
10^{10}$${M}_{\odot}$/$h$.

##  4 Results

In this section, we discuss the signatures of the $k$-essence and massive
neutrino models in the matter and halo power spectra. We study the even
moments of the spectra in redshift space, that is, the monopole, quadrupole,
and hexadecapole, which contain the most information about clustering
(including bias) and standard RSDs on linear and quasi-linear scales (Kaiser
1987). The spectra are computed using the code
Pylians3111https://pylians3.readthedocs.io/en/master/ on the snapshots at
$z=0$, and Doppler RSD are included using the distant-observer approximation
where the velocity-induced shifts in redshift are parallel to one of the
coordinate axes. Our halo catalogues are produced with the Rockstar halo
finder (Behroozi et al. 2013). From each halo catalogue, we select two samples
as explained in Section 2: one using a fixed mass threshold, and one using a
fixed bias value by adjusting the mass threshold accordingly. For the matter
power spectra, as well as for the halo power spectra from each selection
method, we compute the ratios with respect to the one found in the
corresponding $\Lambda$CDM simulation. As the initial conditions were chosen
to be perfectly correlated across simulations, cosmic variance largely drops
out in these ratios.

In Fig. 2 we show the results for the monopole power spectra from the
$k$-essence simulations. Blue colour is used to show the case where
$c_{s}^{2}=1$ and orange colour is used for the case where
$c_{s}^{2}=10^{-4}$. The stars correspond to the samples selected at fixed
bias while the circles are for halos with fixed mass threshold. The figure
shows clearly that the first selection can recover the signature in the
monopole of the power spectrum very well in the linear and quasi-linear
regime, $k\lesssim 0.1h$/Mpc. On the other hand, in the power spectra of the
halos selected according to fixed mass threshold, the signal is completely
removed and the ratio of the spectra is close to unity. We observe a similar
behaviour for the other moments of the power spectrum for the $k-$essence
models, which, on the other hand, is less visible due to the weakness of the
effect and to the large noise in the higher multipoles.

Figure 2:  Ratio of $k$-essence to $\Lambda$CDM halo monopole power spectrum
(symbols) compared to the ratio of the corresponding matter power spectra
(solid lines). The stars represent the result for the fixed bias selection,
while the circles correspond to the selection at fixed mass threshold.

In Fig. 3 we show the results for the monopole spectra of the massive neutrino
models with the sum of neutrino masses $\sum m_{\nu}=0.06$ eV, $0.2$ eV, and
$0.3$ eV in blue, orange, and green, respectively. Similar to the $k$-essence
model, comparing halos selected by the same mass threshold within different
models obscures the difference to $\Lambda$CDM at the intermediate scales,
$0.03\,h/\mathrm{Mpc}\lesssim k\lesssim 0.3\,h/\mathrm{Mpc}$, where clustering
measurements are most sensitive. While some differences are noticeable at
larger and smaller scales, they bear no resemblance to the signatures seen in
the matter, apparently even reversing some of the usual trends. These
differences are somewhat difficult to interpret due to the relatively large
contribution from RSD. On the other hand, a comparison within the fixed bias
halo samples recovers the signature of the matter power spectrum at the
intermediate scales which is also where the bias is measured. We also find the
same behaviour in the higher moments of the power spectra as shown in Figs. 4
and 5. We see that comparing quadrupole and hexadecapole for the fixed bias
halo selection recovers the signatures found in matter, although less clearly
than in the monopole. This is mainly due to the larger uncertainties in the
measurements of these higher moments which have significantly smaller
amplitudes. On the other hand, we do not see any significant effect from
different neutrino masses when we consider a fixed mass selection in different
theories. It is also interesting to note that the ratios of the quadrupoles
and hexadecapoles at fixed bias value agree relatively well with the ones of
matter even in the nonlinear regime, $k>0.2h/$Mpc, where the corresponding
monopole ratios deviate.

Figure 3: Same as Fig. 2, but comparing models with different neutrino masses
to $\Lambda$CDM (where neutrinos are massless). Only a halo selection at fixed
bias (stars) recovers the monopole power ratio of matter (solid lines) on
intermediate scales. Bias values and mass thresholds are chosen in a different
way to in the $k$-essence comparison because of the better mass resolution.
Figure 4:  Ratio with respect to the massless case ($\Lambda$CDM) of the
quadrupole of the halo power spectra (symbols with error bars) compared to the
corresponding ratios for matter (solid lines) for different neutrino masses.
The stars represent the result for the fixed bias selection, while the circles
correspond to the selection at fixed mass threshold. The divergence at
$k\simeq 0.3h^{-1}$Mpc is caused by a zero-crossing of the quadrupole. Figure
5:  Same as Fig. 4, but for the hexadecapole of the power spectra.

From the bias definition in Eq. (1) we can understand why comparing fixed bias
halo samples in different models can recover the behaviour of matter at least
for the scales on which the bias is linear. By selecting halo samples in
$\Lambda$CDM and in model $X$ at the same bias value, we can write

$\frac{P^{\Lambda{\rm{CDM}}}_{hh}(k)}{P^{\Lambda{\rm{CDM}}}_{mm}(k)}=\frac{P^{X}_{hh}(k)}{P^{X}_{mm}(k)}\quad\Rightarrow\quad\frac{P^{X}_{hh}(k)}{P^{\Lambda{\rm{CDM}}}_{hh}(k)}=\frac{P^{X}_{mm}(k)}{P^{\Lambda{\rm{CDM}}}_{mm}(k)}.$
(2)

This relation is valid on scales where the linear bias model is a good
description in both cosmologies, and it can explain the results seen in Figs.
2 and 3 at linear and mildly non-linear scales, except that they show the
monopole in redshift space which includes RSD. In linear theory, the RSD can
be modelled following Kaiser (1987), leading to the well-known expressions for
the multipoles of the power spectrum

$\displaystyle P^{(0)}(k)$ $\displaystyle=$
$\displaystyle\left(1+\frac{2}{3}\beta+\frac{1}{5}\beta^{2}\right)P(k)\,,$ (3)
$\displaystyle P^{(2)}(k)$ $\displaystyle=$
$\displaystyle\left(\frac{4}{3}\beta+\frac{4}{7}\beta^{2}\right)P(k)\,,$ (4)
$\displaystyle P^{(4)}(k)$ $\displaystyle=$
$\displaystyle\frac{8}{35}\beta^{2}P(k)\,,$ (5)

where $\beta=f/b_{1}$ and $f$ is the logarithmic growth rate in linear theory.
Having fixed $b_{1}$ we can see from Eq. (2) that the halo power spectrum
ratio would deviate from the matter power spectrum ratio by less than
$\frac{2}{3}\Delta f$ for the monopole, by less than $\frac{3}{7}\Delta f$ for
the quadrupole, and the hexadecapole ratio remains unaffected. Here, $\Delta
f$ is the change in the growth rate between the two cosmologies. Large changes
in the growth rate might therefore affect the monopole and quadrupole ratios,
but for the models studied here we have $\Delta f$ of approximately a few
percent at most, and therefore smaller than the typical change of bias when
considering a fixed type of tracer. It should also be noted that the linear
RSD model is not very accurate, and becomes completely useless at short
distance scales where the so-called fingers-of-God effect becomes relevant. It
is therefore interesting to see our results play out on those non-linear
scales.

##  5 Conclusions

Near-future LSS surveys will put very tight bounds on cosmological parameters.
They will measure the sum of the neutrino masses and will help us to
understand the nature of dark energy. To exploit the full potential of these
surveys in better constraining our theories using biased tracers, we need to
understand how a signal in the matter distribution is translated into the halo
and the galaxy clustering statistics.

In this work, we show for the first time that a naive selection of the halo
sample based on a fixed mass threshold within different theories leads to the
almost complete removal of the signature in the halo power spectrum. For this
purpose, we used suites of large-scale cosmological simulations for two
different models, namely $k$-essence and massive neutrinos. We demonstrate
that, even in the higher moments of the power spectra, one sees a removal of
signal for a fixed mass selection. On the other hand, when considering a
selection at fixed value of the linear bias, the signal is recovered in the
halo power spectra. The simple argument leading to this conclusion is
summarised in Eq. (2) and we expect it to also hold for other model
comparisons. On non-linear scales, the relation between the matter power
spectrum and the halo power spectrum is difficult to model and is beyond the
scope of this paper.

The results of this paper also show the importance of modelling or measuring
the bias accurately. The bias can in principle be measured by correlating weak
lensing shear (which depends on total matter) and clustering statistics of the
tracers, but the current sensitivity is still rather poor. Using a measured
bias, we can construct the tracer catalogue accordingly such that we see the
expected effect in the tracer clustering.

###### Acknowledgements.

We thank Vincent Desjacques for helpful discussions. This work is supported by
a grant from the Swiss National Supercomputing Centre (CSCS) under project ID
s1051. We have benefited from the Baobab and Yggdrasil clusters at the
University of Geneva. Part of the numerical calculations were performed on the
computational facilities at University of Oslo. FH is supported by the
Overseas Research Fellowship from the Research Council of Norway. JA, RD and
MK acknowledge financial support from the Swiss National Science Foundation.

## References

  * Abell et al. (2009) Abell, P. A. et al. 2009 [arXiv:0912.0201]
  * Adamek et al. (2017) Adamek, J., Durrer, R., & Kunz, M. 2017, JCAP, 11, 004
  * Aghanim et al. (2020) Aghanim, N. et al. 2020, Astron. Astrophys., 641, A6, [Erratum: Astron.Astrophys. 652, C4 (2021)]
  * Alam et al. (2017) Alam, S. et al. 2017, Mon. Not. Roy. Astron. Soc., 470, 2617
  * Amendola et al. (2018) Amendola, L. et al. 2018, Living Rev. Rel., 21, 2
  * Angulo & Hahn (2021) Angulo, R. E. & Hahn, O. 2021, Living Rev. Comput. Astrophys., 8, 1
  * Armendariz-Picon et al. (2001) Armendariz-Picon, C., Mukhanov, V. F., & Steinhardt, P. J. 2001, Phys. Rev. D, 63, 103510
  * Baldi et al. (2010) Baldi, M., Pettorino, V., Robbers, G., & Springel, V. 2010, Mon. Not. Roy. Astron. Soc., 403, 1684
  * Barreira et al. (2013) Barreira, A., Li, B., Hellwing, W. A., Baugh, C. M., & Pascoli, S. 2013, JCAP, 10, 027
  * Behroozi et al. (2013) Behroozi, P. S., Wechsler, R. H., & Wu, H.-Y. 2013, Astrophys. J., 762, 109
  * Desjacques et al. (2018) Desjacques, V., Jeong, D., & Schmidt, F. 2018, Phys. Rept., 733, 1
  * Esteban et al. (2019) Esteban, I., Gonzalez-Garcia, M. C., Hernandez-Cabezudo, A., Maltoni, M., & Schwetz, T. 2019, JHEP, 01, 106
  * Hassani et al. (2020a) Hassani, F., Adamek, J., & Kunz, M. 2020a, Mon. Not. Roy. Astron. Soc., 500, 4514
  * Hassani et al. (2019) Hassani, F., Adamek, J., Kunz, M., & Vernizzi, F. 2019, JCAP, 12, 011
  * Hassani et al. (2020b) Hassani, F., L’Huillier, B., Shafieloo, A., Kunz, M., & Adamek, J. 2020b, JCAP, 04, 039
  * Hassani & Lombriser (2020) Hassani, F. & Lombriser, L. 2020, Mon. Not. Roy. Astron. Soc., 497, 1885
  * Inman et al. (2015) Inman, D., Emberson, J. D., Pen, U.-L., et al. 2015, Phys. Rev. D, 92, 023502
  * Kaiser (1984) Kaiser, N. 1984, Astrophys. J. Lett., 284, L9
  * Kaiser (1987) Kaiser, N. 1987, Mon. Not. Roy. Astron. Soc., 227, 1
  * Laureijs et al. (2011) Laureijs, R. et al. 2011 [arXiv:1110.3193]
  * Lesgourgues & Pastor (2006) Lesgourgues, J. & Pastor, S. 2006, Phys. Rept., 429, 307
  * Li (2018) Li, B. 2018, Simulating Large-Scale Structure for Models of Cosmic Acceleration, 2514-3433 (IOP Publishing)
  * Llinares et al. (2014) Llinares, C., Mota, D. F., & Winther, H. A. 2014, Astron. Astrophys., 562, A78
  * Perlmutter et al. (1999) Perlmutter, S. et al. 1999, Astrophys. J., 517, 565
  * Riess et al. (1998) Riess, A. G. et al. 1998, Astron. J., 116, 1009