<figure><img src="image_1.png" /></figure>

HAL
open science

# Wasserstein gradient flow of the Fisher information from a non-smooth convex minimization viewpoint

Guillaume Carlier, Jean-David Benamou, Daniel Matthes

► To cite this version:

Guillaume Carlier, Jean-David Benamou, Daniel Matthes. Wasserstein gradient flow of the Fisher information from a non-smooth convex minimization viewpoint. 2023. hal-04178416

HAL Id: hal-04178416

https://inria.hal.science/hal-04178416

Preprint submitted on 8 Aug 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# WASSERSTEIN GRADIENT FLOW OF THE FISHER INFORMATION FROM A NON-SMOOTH CONVEX MINIMIZATION VIEWPOINT

GUILLAUME CARLIER, JEAN-DAVID BENAMOU, AND DANIEL MATTHES

ABSTRACT. Motivated by the Derrida-Lebowitz-Speer-Spohn (DLSS) quantum drift equation, which is the Wasserstein gradient flow of the Fisher information, we study in details solutions of the corresponding implicit Euler scheme. We also take advantage of the convex (but non-smooth) nature of the corresponding variational problem to propose a numerical method based on the Chambolle-Pock primal-dual algorithm.

*Dedicated to Giuseppe Buttazzo, a master of the calculus of variations, on the occasion of his 70th birthday.*

## 1. INTRODUCTION

The non-linear fourth order parabolic equation

$$ (1.1) \qquad \partial_t \rho + \operatorname{div} \left( \rho \nabla \frac{\Delta \sqrt{\rho}}{\sqrt{\rho}} \right) = 0 $$

appears in the description of interface fluctuations in spin systems [22, 23], as quantum correction to semi-conductor models [2], and in the diffusive limit of a quantum mechanical Boltzmann equation [21]. The mathematical analysis of (1.1) has been initiated in [8]. By now, well-posedness and long-time asymptotics are fairly well understood.

In the course of the analysis, a remarkable number of structural properties for (1.1) has been discovered. There exists a whole zoo of integral-type Lyapunov functionals [8, 34], these have been the basis for proving global existence of non-negative weak solutions [36, 35, 28]. Further, the flow induced by (1.1) is contractive in the Hellinger distance [36], leading to uniqueness of weak solutions [26]. Significant for the paper at hand is the observation by Gianazza, Savaré and Toscani that (1.1) is the (non-convex) gradient flow of the Fisher information in the $L^2$-Wasserstein distance [28]. Exploiting that gradient flow representation and in particular the intimate relation to the contractive heat flow [17], self-similar long-time asymptotics towards the Gaussian could be proven [41]. Most importantly, the gradient

<footnote n="2020">2020 Mathematics Subject Classification. 49Q22, 65K10.</footnote>

<footnote>Key words and phrases. DLSS equation, Wasserstein distance, Fisher information, convex duality, Chambolle-Pock algorithm.</footnote>
flow formulation of (1.1) paves the way to the construction of weak solutions by means of the variational JKO scheme [32], see (1.3) below, both analytically and, after an additional spatial discretization, also numerically.

In the paper at hand, we derive a novel discrete-in-time variational formulation of (1.1), prove its unique solvability, and obtain a practical numerical scheme. We thus complement the existing literature both on semi-discretizations [37, 10, 9, 38] and on full discretizations [24, 40, 43, 39, 15, 31] of (1.1). In comparison to the very sophisticated numerical methods that have been proposed before, our scheme sticks out by simplicity yet robustness. In a nutshell, a time step for (1.1) is obtained by iterating the following two simple operations: first, the solution of the linear elliptic equation ($D_h$ and $-D_h^*$ denote discrete gradient and divergence, respectively),

$$ \phi + \delta D_h^*(\rho_0 D_h \phi) = f $$

with data $\rho_0$, an input $f$ from the previous iterate, and a parameter $\delta > 0$; and second, the projection — pointwise at each vertex of the spatial grid — of a certain pair $(\alpha, \beta)$ of functions defined from the previous iterate to the convex set

$$ (1.2) \qquad K = \left\{ (a, b) \in \mathbb{R} \times \mathbb{R}^d : a + \frac{1}{2} |b|^2 \le 0 \right\}. $$

The parabola $K$ plays a pivotal role in our considerations, since both the Fisher information and the action integral in the dynamical formulation of the Wasserstein metric can be written with the help of $K$'s support function; see Section 2 for further details.

We obtain our time-discrete variational formulation of (1.1) as dual to the JKO time stepping. We recall the latter: given a time step $\tau > 0$ and an initial datum $\rho_0$, the semi-discrete approximation $\rho_k$ of the solution to (1.1) at time $t = k\tau$ is obtained by inductive solution of

$$ (1.3) \qquad \rho_k \in \underset{\rho \in \mathcal{P}}{\operatorname{argmin}} \left\{ \frac{1}{2\tau} W_2^2(\rho, \rho_{k-1}) + F(\rho) \right\}, \qquad F(\rho) = \int |\nabla \sqrt{\rho}|^2 dx, $$

where $W_2$ is the $L^2$-Wasserstein distance. The latter involves the solution of an inner minimization problem, see Section 2 below. The novel dual formulation amounts to choosing $\rho_k \in \mathcal{P}$ such that $\nabla \log \rho_k = \beta_k$, where

$$ (1.4) \qquad (\alpha_k, \beta_k) \in \underset{(\alpha, \beta) \in K}{\operatorname{argmax}} \int \left( \frac{\tau}{2} (\operatorname{div}\beta - \alpha) \right)^c \rho_{k-1}. $$

Also this problem involves an inner minimization via the $c$-transform $(\cdot)^c$, recalled in Section 2.

Our results concern solutions $\rho_k \in \mathcal{P}(\mathbb{T}^d)$ on the flat $d$-dimensional torus $\mathbb{T}^d := \mathbb{R}^d/\mathbb{Z}^d$. We show that the minimization problem (1.3) possesses a unique solution
$\rho_k$, and that this solution is positive and classical. Then we prove that

$$ (1.5) \qquad \alpha_k = -\frac{1}{2}|\nabla \log \rho_k|^2 \quad \text{and} \quad \beta_k = \nabla \log \rho_k $$

form the unique solution pair for the maximization problem (1.4). Finally, we perform a spatial discretization of (1.4), apply the Chambolle-Pock algorithm [19] for numerical solution, and use the resulting scheme for approximation of solutions to (1.1).

Convex minimization algorithms have been used extensively for practical implementation of the JKO time stepping for (non-)linear Fokker-Planck and similar second order parabolic equations. The most common algorithm is the augmented Lagrangian method, also known as ADMM for Alternating Direction Method of Multipliers. The first application of ADMM to an optimal transport problem has been proposed in [4] for numerical solution of the dynamic Benamou-Brenier formulation. This approach has then been extended to transport problems with congestion by Buttazzo, Jimenez and Oudet [12], see also [5, 46]. By now, there are plenty of applications of ADMM to JKO discretizations of gradient flows [6]. For gradient flows of functionals involving derivatives of the density, ADMM is unpractical since it involves the solution of auxiliary nonlinear PDEs at each time step. In the specific case of the Fisher information F from above, the Chambolle-Pock algorithm provides a much better alternative: it is easy to implement and requires no regularization of the datum $\rho_{k-1}$. In particular, densities that vanish on sets of positive measure are admissible data. The price of robustness, however, is a relatively slow speed of convergence.

We briefly review some of the various alternatives to the aforementioned augmented Lagrangian methods for numerical solution of the JKO stepping for Wasserstein gradient flows. On the one hand, there are fully Lagrangian methods for second [16, 30, 42, 24, 33, 27, 7] and also fourth order [18, 45] parabolic equations, that take advantage of the geometry behind mass transport by approximating the Wasserstein distance directly from mass displacement. Specifically for discretization of (1.1) in one space dimension, a convergent Lagrangian scheme is given in [43], which also preserves several structural properties like the long-time asymptotics. Concerning variational methods for (1.1) using finite differences, we mention [39], where minimization in the original JKO formulation (1.3) is carried out by means of Newton's method. The spatial discretization of the Fisher information is performed such that it blows up for vanishing densities, which keeps the Newton iterates in the region where densities are positive and the functional is strictly convex. A spatial discretization of the JKO scheme by finite elements has been performed in [11, 13], based on an approximation of the Wasserstein metric by an $H^{-1}$-distance with explicitly [11] or implicitly [13] defined weight. The resulting method is convergent and numerically robust for Fokker-Planck equations, but is apparently difficult to adapt to fourth order equations.
The paper is organized as follows. Preliminaries from the theory of optimal mass transport are gathered in Section 2 below. In Section 3, we prove existence, uniqueness, positivity and smoothness of the minimizer to the JKO scheme (1.3). Then, in Section 4, the correspondence (1.5) between (1.3) and the dual formulation (1.4) is established. Finally, in Section 5, the dual problem (1.4) is spatially discretized and numerically solved.

## 2. SETTING AND PRELIMINARIES

**Spaces.** Throughout the paper, we will work in the periodic setting, i.e., on the flat $d$-dimensional torus $\mathbb{T}^d := \mathbb{R}^d/\mathbb{Z}^d$, equipped with the distance

$$ \mathrm{dist}(x, y) = \mathrm{dist}(x - y, 0) := \min_{k \in \mathbb{Z}^d} |x - y + k|. $$

We denote by $\mathcal{P}(\mathbb{T}^d)$ and $\mathcal{M}(\mathbb{T}^d)$, respectively, the sets of Borel probability measures and of Borel measures on $\mathbb{T}^d$. By abuse of notation, we shall not distinguish between an absolutely continuous measure $\mu \in \mathcal{P}(\mathbb{T}^d)$ and its probability density $\rho \in L^1(\mathbb{T}^d)$.

**Wasserstein metric: primal formulation.** The $L^2$-Wasserstein distance $W_2(\mu, \nu)$ between probability measures $\mu, \nu \in \mathcal{P}(\mathbb{T}^d)$ is given by

$$ W_2^2(\mu, \nu) := \inf_{\gamma \in \Pi(\mu, \nu)} \int_{\mathbb{T}^d \times \mathbb{T}^d} \mathrm{dist}^2(x, y) \mathrm{d}\gamma(x, y), \quad \forall (\mu, \nu) \in \mathcal{P}(\mathbb{T}^d)^2 $$

where $\Pi(\mu, \nu)$ denotes the set of transport plans between $\mu$ and $\nu$, i.e. the set of probability measures on $\mathbb{T}^d \times \mathbb{T}^d$ having $\mu$ and $\nu$ as marginals.

**Wasserstein metric: dual formulation.** The well-known Kantorovich duality formula, c.f. [48], [49], expresses the Wasserstein distance as a maximization problem,

$$ \frac{1}{2} W_2^2(\mu, \nu) := \sup \left\{ \int_{\mathbb{T}^d} \phi \mathrm{d}\mu + \int_{\mathbb{T}^d} \psi \mathrm{d}\nu : \phi \oplus \psi \le \frac{\mathrm{dist}^2}{2}, \right\} $$

where $\phi \oplus \psi$ denotes the separable function $(x, y) \mapsto \phi(x) + \psi(y)$. It is well-known that this dual formulation admits optimal potentials which can be chosen to be $c$-conjugates of each other, i.e.

$$ \phi(x) = \psi^c(x) := \inf_{y \in \mathbb{T}^d} \left\{ \frac{\mathrm{dist}^2(x, y)}{2} - \psi(y) \right\}, \quad \forall x \in \mathbb{T}^d, $$

$$ \psi(y) = \phi^c(y) := \inf_{x \in \mathbb{T}^d} \left\{ \frac{\mathrm{dist}^2(x, y)}{2} - \phi(x) \right\}, \quad \forall y \in \mathbb{T}^d. $$

In particular $\phi$ and $\psi$ are $\sqrt{d}$-Lipschitz and 1/2-semi-concave, i.e., the (non-periodic) map $x \in \mathbb{R}^d \mapsto \frac{1}{2}|x|^2 - \phi(x)$ is convex.

If the pair $(\phi, \psi) = (\phi, \phi^c)$ solves the Kantorovich dual above, one says that $\phi$ is a Kantorovich potential between $\mu$ and $\nu$. This notion is useful when one wants
to differentiate $W_2^2$ with respect to one of its marginals. The following result is well-known in the euclidean case, for the sake of completeness, we give a proof in the periodic case (which is covered by the analysis of McCann for the general Riemannian case [44], see also Cordero-Erausquin [20] for more on the periodic case):

**Lemma 2.1.** Let $\mu$ and $\nu$ be in $\mathcal{P}(\mathbb{T}^d)$ with $\mu$ absolutely continuous with respect to the Lebesgue measure and with an a.e. positive density. Then there exists a unique (up to an additive constant) Kantorovich potential $\phi$ between $\mu$ and $\nu$ and

$$ \frac{1}{2} \frac{d}{dt} W_2^2(\mu + t(\alpha - \mu), \nu)|_{t=0^+} = \int_{\mathbb{T}^d} \phi d(\alpha - \mu). $$

*Proof.* Let $\gamma \in \Pi(\mu, \nu)$ and $\phi$ be respectively an optimal plan, and a Kantorovich potential between $\mu$ and $\nu$, then we have:

$$ \phi(x) + \phi^c(y) \le \frac{\mathrm{dist}^2(x, y)}{2}, \quad \forall (x, y) \in \mathbb{T}^d \times \mathbb{T}^d \quad \text{with equality on } \mathrm{spt}(\gamma) $$

since $\phi$ is Lipschitz it is differentiable a.e. hence $\mu$ a.e. and $\mathrm{dist}^2$ is superdifferen- tiable everywhere. Now if $(x, y) \in \mathrm{spt}(\gamma)$ is such that $\phi$ is differentiable at $x$, and that $q$ is in the superdifferential of $\frac{1}{2} \mathrm{dist}^2(\cdot, y)$ at $x$, then

$$ \begin{aligned} \phi(x+h) &\le \frac{\mathrm{dist}^2(x+h, y)}{2} - \phi^c(y) \\ &\le \frac{\mathrm{dist}^2(x, y)}{2} - \phi^c(y) + q \cdot h + o(h) = \phi(x) + q \cdot h + o(h). \end{aligned} $$

Hence $q = \nabla\phi(x)$, so the superdifferential of $\frac{1}{2} \mathrm{dist}^2(\cdot, y)$ is reduced to the singleton $\nabla\phi(x)$. This means that $\frac{1}{2} \mathrm{dist}^2(\cdot, y)$ is in fact differentiable at $x$, with $\nabla\phi(x) = \frac{1}{2} \nabla_x \mathrm{dist}^2(x, y)$. Since this holds at any $y$ such that $(x, y) \in \mathrm{spt}\gamma$, disintegrating $\gamma$ with respect to $\mu$ as $\gamma = \mu \otimes \gamma^x$, we get

$$ (2.1) \qquad \nabla\phi(x) = \int_{\mathbb{T}^d} \frac{1}{2} \nabla_x \mathrm{dist}^2(x, y) d\gamma^x(y) $$

for $\mu$-a.e. $x \in \mathbb{T}^d$. But since the right-hand side of (2.1) only depends on $\gamma$, any other Kantorovich potential should satisfy (2.1) as well, and then should have $\mu$-a.e. and thus also Lebesgue a.e. the same gradient as $\phi$. This shows uniqueness of the Kantorovich potential $\phi$ up to a constant. Now the differentiability claim is easy to deduce, since denoting by $\phi_t$ ($t \in (0, 1]$) a Kantorovich potential between $\mu + t(\alpha - \mu)$ and $\nu$, one has

$$ \int_{\mathbb{T}^d} \phi d(\alpha - \mu) \le \frac{1}{2t} [W_2^2(\mu + t(\alpha - \mu), \nu) - W_2^2(\mu, \nu)] \le \int_{\mathbb{T}^d} \phi_t d(\alpha - \mu) $$

we can normalize $\phi_t$ and $\phi$ to have zero mean and then, since cluster points as $t \to 0^+$ of the (uniformly Lipschitz and bounded) family $\phi_t$ should be Kantorovich
potentials between $\mu$ and $\nu$, by the uniqueness shown previously, we deduce that $\phi_t$ in fact converges uniformly to $\phi$ which shows the desired result. $\square$

**Wasserstein metric: dynamic formulation.** For later use in the numerical scheme, we also need the dynamical formulation from [4] of optimal transport. For $\rho_0, \rho_1 \in \mathcal{P}(\mathbb{T}^d)$, we have

$$ (2.2) \qquad \frac{1}{2} W_2^2(\rho_0, \rho_1) = \inf \left\{ \int_0^1 \int_{\mathbb{T}^d} \frac{|\mathbf{m}_s|^2}{2\rho_s} \, dx \, ds : \partial_s \rho_s + \mathrm{div} \mathbf{m}_s = 0 \right\}, $$

where the infimum runs over all pairs $(\rho_s, \mathbf{m}_s)_{s \in [0,1]}$ of weakly continuous curves $\rho(\cdot) : [0,1] \to \mathcal{P}(\mathbb{T}^d)$ with prescribed end points $\rho_0$ and $\rho_1$ and $\mathbb{R}^d$-valued Radon measures $\mathbf{m}$ on $(0,1) \times \mathbb{T}^d$ such that the continuity equation $\partial_s \rho_s + \mathrm{div} \mathbf{m}_s = 0$ holds in the sense of distribution. The quotient

$$ (2.3) \qquad \mathrm{BB}(r, m) := \frac{|m|^2}{2r} $$

inside the integral in (2.2) is understood as the convex and lower semi-continuous envelope of the respective expression for $r > 0$, which is zero where $r \ge 0$ and $m = 0$, and infinity where either $r < 0$, or $r = 0$ but $m \ne 0$. This ensures in particular non-negativity of the curve $\rho_s$. For the fully rigorous formulation of (2.2), we refer to [1, Chapter 8].

Of particular relevance to us is the fact that the quotient $\mathrm{BB}$ is the support function of the convex set $K$ from (1.2) that is to say that

$$ (2.4) \qquad \mathrm{BB}(r, m) = \sup_{(a,b) \in K} \{ar + b \cdot m\} \quad \text{for all } (r, m) \in \mathbb{R} \times \mathbb{R}^d. $$

This confirms the claimed convexity and lower semi-continuity of $\mathrm{BB}$. The Legendre transform of $\mathrm{BB}$ is thus the indicator function

$$ \mathrm{BB}^*(a, b) = \chi_K(a, b) = \begin{cases} 0 & \text{if } (a, b) \in K \\ +\infty & \text{otherwise.} \end{cases} $$

**Fisher information.** In (1.3), we have (informally) introduced the Fisher information functional $F : \mathcal{M}(\mathbb{T}^d) \to \mathbb{R} \cup \{+\infty\}$ as

$$ (2.5) \qquad F(\rho) := \begin{cases} \int_{\mathbb{T}^d} |\nabla \sqrt{\rho}|^2 & \text{if } \rho \ge 0 \text{ and } \sqrt{\rho} \in H^1(\mathbb{T}^d) \\ +\infty & \text{otherwise.} \end{cases} $$

This functional is convex and lower semi-continuous. Note that if $\rho \ge 0$ and $\sqrt{\rho} \in H^1(\mathbb{T}^d)$, then $\nabla \rho = 2\sqrt{\rho}\nabla \sqrt{\rho} \in L^1$, and we may write

$$ F(\rho) = \frac{1}{4} \int_{\mathbb{T}^d} \frac{|\nabla \rho|^2}{\rho}. $$

The integrand above is to be understood as zero for $\nabla \rho = 0$ and $\rho = 0$, and as $+\infty$ for $\rho = 0$ but $\nabla \rho \ne 0$. In other words, the quotient is again of the type (2.3),
that is, the support function of the parabola K from (1.2). This gives rise to the following alternative representation (which can be obtained by similar arguments as in Proposition 5.18 in [48]):

$$ (2.6) \qquad F(\rho) = \sup \left\{ \int_{\mathbb{T}^d} \frac{(\alpha - \operatorname{div}\beta)}{2} d\rho \ : \ (\alpha, \beta) \in C(\mathbb{T}^d) \times C^1(\mathbb{T}^d, \mathbb{R}^d), \ (\alpha, \beta) \in K \right\}. $$

**JKO functional.** Finally, let $\tau > 0$ and $\rho_0 \in \mathcal{P}(\mathbb{T}^d)$ be given. We are interested in the (convex) minimization

$$ (2.7) \qquad \inf_{\rho \in \mathcal{P}(\mathbb{T}^d)} J(\rho) := \frac{1}{2\tau} W_2^2(\rho_0, \rho) + F(\rho). $$

Existence of a minimizer for (2.7) is easily shown by the direct method. Uniqueness requires more work, see next section, since neither $W^2(\rho_0, \cdot)$ (which is linear if $\rho_0$ is a Dirac mass) nor $F$ (which is positively homogeneous) is strictly convex.

### 3. UNIQUENESS, POSITIVITY AND REGULARITY OF THE MINIMIZER

In this section, we prove:

**Theorem 3.1.** The minimization problem (2.7) admits a unique solution $\rho$. Moreover, $\rho \in W^{3,p}(\mathbb{T}^d)$ for every $p \in [1, +\infty)$ $\rho$ is strictly positive, and $\rho$ solves the Euler-Lagrange PDE

$$ (3.1) \qquad \frac{\phi}{\tau} = \frac{1}{2} \Delta \log \rho + \frac{1}{4} |\nabla \log \rho|^2 = \frac{\Delta \sqrt{\rho}}{\sqrt{\rho}}. $$

where $\phi$ is a Kantorovich potential between $\rho$ and $\rho_0$.

Notice that $\rho \in W^{3,p}(\mathbb{T}^d)$ for all $p \ge 1$ implies that $\rho \in C^{2,\alpha}(\mathbb{T}^d)$ for every $\alpha \in [0, 1]$. Hence all derivatives in (3.1) are actually classical.

**3.1. Approximation.** To enforce positivity, given $\delta > 0$ we consider (as in [14] in the context of the TV-JKO scheme) the entropic approximation of (2.7):

$$ (3.2) \qquad \inf_{\rho \in \mathcal{P}(\mathbb{T}^d)} J_{\delta}(\rho) := \frac{1}{2\tau} W_2^2(\rho_0, \rho) + F(\rho) + \delta \int_{\mathbb{T}^d} \rho \log \rho. $$

**Lemma 3.2.** Given $\delta > 0$, (3.2) admits a unique solution $\rho_\delta$, and $\rho_\delta$ satisfies $\rho_\delta + \frac{1}{\rho_\delta} \in L^\infty(\mathbb{T}^d)$.

*Proof.* Existence follows from weak * compactness and lower semicontinuity of $J_\delta$ and uniqueness by strict convexity. If $\rho_\delta$ is constant, there is nothing to prove, setting $m := \operatorname{essinf}\rho_\delta$ and $M := \operatorname{essup}\rho_\delta$ we may therefore assume that $m < M$
and our aim is to show $m > 0$ and $M < +\infty$. To do so, consider the two functions $\underline{G}$ and $\overline{G}$ defined by

$$ \underline{G}(\alpha) := \int_{\mathbb{T}^d} (\alpha - \rho_{\delta})_{+}, \quad \overline{G}(\alpha) := \int_{\mathbb{T}^d} (\rho_{\delta} - \alpha)_{+}, \quad \forall \alpha \in \mathbb{R}_{+}. $$

Both $\underline{G}$ and $\overline{G}$ are 1-Lipschitz, $\underline{G}$ is nondecreasing with $= \underline{G}(m) = 0$, $\underline{G}(\alpha) > 0$ for $\alpha > m$, $\overline{G}$ is nonincreasing with $\overline{G}(M) = 0$ and $\overline{G}(\alpha) > 0$ for $\alpha < M$. We can therefore choose sequences $(m_n)_n$ and $(M_n)_n$ such that

$$ (3.3) \qquad m < m_n < M_n < M, \quad M_n \nearrow M, \quad m_n \searrow m, \quad \underline{G}(m_n) = \overline{G}(M_n), \quad \forall n. $$

With this choice, we define the sequence of competitors to $\rho_{\delta}$ by suitably truncating $\rho_{\delta}$:

$$ \mu_n := \min(M_n, \max(m_n, \rho_{\delta})) $$

by construction $\mu_n$ is a probability measure and

$$ \|\mu_n - \rho_{\delta}\|_{L^1} = \underline{G}(m_n) + \overline{G}(M_n) =: 2\varepsilon_n \to 0 \text{ as } n \to \infty. $$

Since by construction $\sqrt{\mu_n}$ is obtained by truncating $\sqrt{\rho_{\delta}}$, we have $F(\mu_n) \le F(\rho_{\delta})$, the fact that $\rho_{\delta}$ minimizes $J_{\delta}$ thus gives

$$ 0 \le \frac{1}{2\tau} [W_2^2(\mu_n, \rho_0) - W_2^2(\rho_{\delta}, \rho_0)] + \delta \int_{\mathbb{T}^d} \mu_n \log \mu_n - \delta \int_{\mathbb{T}^d} \rho_{\delta} \log \rho_{\delta} $$

By Kantorovich duality formula, denoting by $\phi_n$ a Kantorovich potential ($\sqrt{d}$-Lipschitz and with zero mean, say) between $\mu_n$ and $\rho_0$, the first term in the right-hand side can be bounded from above by $\tau^{-1} \int_{\mathbb{T}^d} \phi_n(\mu_n - \rho_{\delta}) \le \tau^{-1} \|\phi_n\|_{L^{\infty}} \|\mu_n - \rho_{\delta}\|_{L^1} \le C_d \tau^{-1} \varepsilon_n$ where $C_d$ is a just a dimensional constant, using this with the convexity of $t \mapsto t \log t$ and the construction of $\mu^n$, we arrive at

$$ \begin{aligned} 0 &\le \frac{C_d}{\tau} \varepsilon_n + \delta \int_{\rho_{\delta} \le m_n} (m_n \log(m_n) - \rho_{\delta} \log(\rho_{\delta})) + \delta \int_{\rho_{\delta} \ge M_n} (M_n \log(M_n) - \rho_{\delta} \log(\rho_{\delta})) \\ &\le \frac{C_d}{\tau} \varepsilon_n + \delta \left( (1 + \log(m_n)) \int_{\mathbb{T}^d} (m_n - \rho_{\delta})_+ - (1 + \log(M_n)) \int_{\mathbb{T}^d} (\rho_{\delta} - M_n)_+ \right) \\ &= \varepsilon_n \left[ \frac{C_d}{\tau} + \delta \log(\frac{m_n}{M_n}) \right] \end{aligned} $$

which since $\varepsilon_n = \underline{G}(m_n) = \overline{G}(M_n) > 0$ implies that $\frac{m_n}{M_n} \ge e^{-C_d \tau \delta}$ preventing $m_n$ to approach 0 and $M_n$ to approach $+\infty$ and thus enabling us to conclude that $m > 0$ and $M < +\infty$ (with $\frac{m}{M} \ge e^{-C_d \tau \delta}$). $\square$

Since $\sqrt{\rho_{\delta}} \in H^1(\mathbb{T}^d)$, the fact that $\rho_{\delta}$ is bounded and bounded away from 0 implies that $\rho_{\delta}$, $\frac{1}{\sqrt{\rho}}$ and $\log \rho_{\delta}$ belong to $H^1(\mathbb{T}^d) \cap L^{\infty}(\mathbb{T}^d)$ as well. More importantly, it enables us to justify the Euler-Lagrange satisfied by $\rho_{\delta}$ as follows. Thanks to Lemma 2.1, there exists a unique (up to an additive constant) Kantorovich potential $\phi_{\delta}$ between $\rho_{\delta}$ and $\rho_0$. Let $\theta \in C^1(\mathbb{T}^d)$ such that $\int_{\mathbb{T}^d} \theta = 0$, for $t \in \mathbb{R}$
with $|t|$ small enough $\rho_\delta + t\theta$ is admissible for (3.2) so that $J_\delta(\rho_\delta + t\theta)$ is minimal at $t=0$. Using Lemma 2.1, we first have

$$ (3.4) \qquad \frac{1}{2} \frac{d}{dt} W_2^2(\rho_\delta + t\theta, \rho_0)|_{t=0} = \int_{\mathbb{T}^d} \phi_\delta \theta. $$

Moreover, thanks to Lemma 3.2, we can write

$$ \begin{aligned} \nabla \sqrt{\rho_{\delta} + t\theta} &= \nabla \sqrt{\rho_{\delta}} + \frac{1}{2} \int_{0}^{t} \nabla \left( \frac{\theta}{\sqrt{\rho_{\delta} + s\theta}} \right) ds \\ &= \nabla \sqrt{\rho_{\delta}} + \frac{t}{2} \nabla \left( \frac{\theta}{\sqrt{\rho_{\delta}}} \right) + t\epsilon_t \quad \text{with } \epsilon_t \to 0 \text{ as } t \to 0 \text{ in } L^2(\mathbb{T}^d). \end{aligned} $$

whence

$$ (3.5) \qquad \frac{d}{dt} F(\rho_{\delta} + t\theta)|_{t=0} = \int_{\mathbb{T}^d} \nabla \sqrt{\rho_{\delta}} \cdot \nabla \left( \frac{\theta}{\sqrt{\rho_{\delta}}} \right) $$

and since

$$ (3.6) \qquad \frac{d}{dt} \left[ \int_{\mathbb{T}^d} (\rho_\delta + t\theta) \log(\rho_\delta + t\theta) \right]_{t=0} = \int_{\mathbb{T}^d} \log(\rho_\delta) \theta $$

we deduce

$$ (3.7) \qquad \int_{\mathbb{T}^d} \left( \frac{\phi_\delta}{\tau} + \delta \log \rho_\delta \right) \theta + \int_{\mathbb{T}^d} \nabla \sqrt{\rho_\delta} \cdot \nabla \left( \frac{\theta}{\sqrt{\rho_\delta}} \right) = 0, \quad \forall \theta \in C^1(\mathbb{T}^d) \text{ such that } \int_{\mathbb{T}^d} \theta = 0. $$

In other words there is some constant $C_\delta$ such that $\sqrt{\rho_\delta}$ is a weak solution of

$$ \frac{\phi_\delta}{\tau} + \delta \log \rho_\delta - \frac{\Delta \sqrt{\rho_\delta}}{\sqrt{\rho_\delta}} = C_\delta $$

up to absorbing this constant in the Kantorovich potential $\phi_\delta$ we thus have

$$ (3.8) \qquad \left( \frac{\phi_\delta}{\tau} + \delta \log \rho_\delta \right) \sqrt{\rho_\delta} = \Delta \sqrt{\rho_\delta}. $$

By standard elliptic regularity, the fact that $\phi_\delta \in W^{1,\infty}(\mathbb{T}^d)$ implies that $\sqrt{\rho_\delta}$ belongs to $W^{3,p}(\mathbb{T}^d)$ for every $p \in [1, +\infty)$ hence $\sqrt{\rho_\delta}$ (and $\rho_\delta$) belong to $C^{2,\alpha}(\mathbb{T}^d)$ for every $\alpha \in [0, 1)$, in particular $\sqrt{\rho_\delta}$ is a classical solution of (3.8).

## 3.2. Proof of Theorem 3.1. We are now in position to prove Theorem 3.1.

**Step 1:** $\delta \to 0$. Note that if $\mu \in \mathcal{P}(\mathbb{T}^d)$ and $\sqrt{\mu} \in H^1(\mathbb{T}^d)$ then by the Sobolev imbedding $\mu \in L^{q(d)}$ for some $q(d) > 1$ (for $d = 1, 2$ one can take any exponent in $(1, +\infty)$) and for $d \ge 3$ one can take $q(d) = \frac{d}{d-2}$) and in particular $\mu$ has finite entropy. From the inequality $J(\rho_\delta) + \delta \int_{\mathbb{T}^d} \rho_\delta \log \rho_\delta \le J(\mu) + \delta \int_{\mathbb{T}^d} \mu \log \mu$ and the fact that $\int_{\mathbb{T}^d} \rho_\delta \log \rho_\delta$ is bounded from below, one readily deduces that $\sqrt{\rho_\delta}$ is bounded in $H^1(\mathbb{T}^d)$ and therefore (up to extracting a vanishing sequence) converges weakly in $H^1(\mathbb{T}^d)$ as $\delta \to 0^+$ (hence strongly in $L^{p(d)}$ for some $p(d) > 2$)
to some $\sqrt{\rho}$, and then $\rho_\delta$ converges strongly to $\rho$ in $L^{p(d)/2}$ and $\rho$ obviously solves (2.7). Note also that both $\rho_\delta \log \rho_\delta$ and $\sqrt{\rho_\delta} \log \rho_\delta$ are bounded in $L^1$ so that

$$ (3.9) \qquad \delta \rho_\delta \log \rho_\delta \text{ and } \delta \sqrt{\rho_\delta} \log \rho_\delta \text{ tend to 0 strongly in } L^1. $$

As for the family of Kantorovich potentials $\phi_\delta$, we already know that $\|\nabla\phi_\delta\|_{L^\infty} \le \sqrt{d}$ and integrating (3.8) and using the periodicity of $\rho_\delta$, we get

$$ \int_{\mathbb{T}^d} \left( \frac{\phi_\delta}{\tau} + \delta \log \rho_\delta \right) \sqrt{\rho_\delta} = 0 $$

which, with (3.9), the fact that $\sqrt{\rho_\delta}$ converges strongly in $L^2$ to $\sqrt{\rho}$ with $\int_{\mathbb{T}^d} \rho = 1$ and $\max_{[0,1]^d} \phi_\delta - \min_{[0,1]^d} \phi_\delta \le \sqrt{d} \operatorname{diam}([0,1]^d) = d$ yields a uniform $L^\infty$ bound on $\phi_\delta$. Thanks to the Arzelà-Ascoli's Theorem, we may assume (after an extraction again) that $\phi_\delta$ converges uniformly to some $\phi$ which is a Kantorovich potential between $\rho$ and $\rho_0$. Since $\Delta\sqrt{\rho_\delta}$ converges to $\Delta\sqrt{\rho}$ in $H^{-1}$, $\phi_\delta\sqrt{\rho_\delta}$ converges to $\phi\sqrt{\rho}$ in $H^1$ and recalling (3.9), we may pass to the limit in (3.8) and get

$$ (3.10) \qquad \frac{\phi}{\tau} \sqrt{\rho} = \Delta \sqrt{\rho}. $$

**Step 2: positivity and regularity of $\rho$.** The announced $W^{3,p}$ regularity of $\rho$ then follows from (3.10) and elliptic regularity. Setting $u = \sqrt{\rho} \in C^{2,\alpha}(\mathbb{T}^d)$, $f = \tau^{-1}\phi \in W^{1,\infty}$, we have $Lu := \Delta u - f_+u = -f_-u \le 0$, hence it follows from the strong maximum principle (Theorem 3.5 in [29]) that if $u$ takes the (global minimum) value 0 one should have $u = 0$ contradicting $\int_{\mathbb{T}^d} u^2 = 1$. We thus have found a smooth and everywhere strictly positive solution $\rho$ of (2.7) and shown that it is a classical solution of (3.1).

**Step 3: uniqueness.** If $\mu$ solves (2.7), so does $\frac{\mu+\rho}{2}$ by convexity of $F$ and $W^2(\rho_0, \cdot)$ yielding that $F(\rho+\mu) = F(\rho) + F(\mu)$ which in particular implies that a.e. on $\{\mu > 0\}$, one has

$$ |\nabla \mu|^2 \frac{\rho}{\mu} + |\nabla \rho|^2 \frac{\mu}{\rho} = 2 \nabla \mu \cdot \nabla \rho $$

so that

$$ (3.11) \qquad \sqrt{\rho} \nabla \sqrt{\mu} = \sqrt{\mu} \nabla \sqrt{\rho} $$

and since $\sqrt{\mu} \in H^1$, $\nabla \sqrt{\mu} = 0$ a.e. on $\{\mu = 0\}$ hence (3.11) actually holds a.e.; taking the divergence of (3.11), we get

$$ \Delta \sqrt{\mu} = \frac{\sqrt{\mu}}{\sqrt{\rho}} \Delta \sqrt{\rho} = \frac{\phi}{\tau} \sqrt{\mu} $$

and arguing exactly as in step 2, we deduce that $\mu$ is $W^{3,p}$ for every $p \in [1, +\infty)$ and everywhere strictly positive. Getting back, to (3.11) we deduce that the smooth functions $\log \mu$ and $\log \rho$ differ by an additive constant which should be 0 since $\mu$
and $\rho$ are probability measures. This proves that $\mu = \rho$ and so that the minimizer of (2.7) is unique.

Remark 3.3. Note that the full sequence $\rho_\delta$ in fact converges to the solution $\rho$ of (2.7). Let us also remark that the Kantorovich potential $\phi$ in (3.1) is defined uniquely (not only up to an additive constant) by the condition

$$ \int_{\mathbb{T}^d} \phi \sqrt{\rho} = 0 $$

obtained by multiplying (3.1) by $\sqrt{\rho}$ and integrating.

### 4. DUALITY

For $\varphi \in C(\mathbb{T}^d)$, set

$$ \varphi^{c_{\tau}}(y) := \inf_{y \in \mathbb{T}^d} \left\{ \frac{1}{2\tau} \operatorname{dist}^2(x, y) - \varphi(y) \right\}, x \in \mathbb{T}^d $$

and

$$ V(\varphi) := - \int_{\mathbb{T}^d} \varphi^{c_{\tau}} d\rho_0 $$

note that $V$ is convex and everywhere finite and continuous on $C(\mathbb{T}^d)$. Moreover, the Kantorovich duality formula can be conveniently expressed in terms of the Legendre transform of $V$, $V^*: \mathcal{M}(\mathbb{T}^d) \to \mathbb{R} \cup \{+\infty\}$:

$$ V^*(\mu) = \begin{cases} \frac{1}{2\tau} W_2^2(\mu, \rho_0) & \text{if } \mu \in \mathcal{P}(\mathbb{T}^d), \\ +\infty & \text{otherwise.} \end{cases} $$

Consider now the concave maximization problem

$$ (4.1) \sup \left\{ \int_{\mathbb{T}^d} \left( \frac{\operatorname{div}\beta - \alpha}{2} \right)^{c_{\tau}} d\rho_0, (\alpha, \beta) \in C(\mathbb{T}^d) \times C^1(\mathbb{T}^d, \mathbb{R}^d) : \alpha + \frac{1}{2} |\beta|^2 \le 0 \right\} $$

which can also be written as

$$ - \inf_{(\alpha, \beta) \in C(\mathbb{T}^d) \times C^1(\mathbb{T}^d, \mathbb{R}^d)} \left\{ U(\alpha, \beta) + V\left(\frac{\operatorname{div}\beta - \alpha}{2}\right) \right\} $$

where

$$ U(\alpha, \beta) = \begin{cases} 0 & \text{if } \alpha + \frac{1}{2} |\beta|^2 \le 0 \text{ on } \mathbb{T}^d \\ +\infty & \text{otherwise} \end{cases} $$

whose dual in the sense of Fenchel-Rockafellar see [25] reads:

$$ \inf_{\mu \in \mathcal{M}(\mathbb{T}^d)} U^* \left( \frac{1}{2} (\mu, \nabla \mu) \right) + V^*(\mu) $$

and since, by (2.6), we have

$$ U^* \left( \frac{1}{2} (\mu, \nabla \mu) \right) = \begin{cases} F(\mu), & \text{if } \mu \ge 0 \text{ and } \sqrt{\mu} \in H^1(\mathbb{T}^d) \\ +\infty & \text{otherwise} \end{cases} $$
it follows from the Fenchel-Rockafellar theorem (Theorem 4.1 in [25]) that:

$$ \min(2.7) = \sup(4.1). $$

In other words, (2.7) is the dual of (4.1). Note that (4.1) is a slightly unusual optimization problem and it is a priori unclear whether it admits a solution $(\alpha, \beta) \in C(\mathbb{T}^d) \times C^1(\mathbb{T}^d, \mathbb{R}^d)$. But thanks to Theorem 3.1, we can easily find a (smooth) solution of (4.1):

**Proposition 4.1.** Let $\rho$ be the solution of (2.7) then the pair (of $C^{1,\alpha}$ functions)

$$ (4.2) \qquad \alpha := -\frac{1}{2}|\nabla \log \rho|^2, \quad \beta := \nabla \log \rho $$

is the only solution of (4.1).

*Proof.* Note that by construction $(\alpha, \beta)$ is admissible for (4.1) and from (3.1), we have

$$ \frac{\phi}{\tau} = \frac{1}{2}\Delta \log \rho + \frac{1}{4}|\nabla \log \rho|^2 = \frac{\mathrm{div}\beta - \alpha}{2} $$

where $\phi$ is a Kantorovich potential between $\rho$ and $\rho_0$ so that, on the one hand

$$ \frac{1}{2\tau}W_2^2(\rho, \rho_0) = \int_{\mathbb{T}^d} \frac{\phi}{\tau}\rho + \int_{\mathbb{T}^d} \left(\frac{\phi}{\tau}\right)^{c_{\tau}} d\rho_0 = \int_{\mathbb{T}^d} \frac{\phi}{\tau}\rho + \int_{\mathbb{T}^d} \left(\frac{\mathrm{div}\beta - \alpha}{2}\right)^{c_{\tau}} d\rho_0. $$

On the other hand

$$ F(\rho) = \frac{1}{4} \int_{\mathbb{T}^d} \frac{|\nabla \rho|^2}{\rho} = \int_{\mathbb{T}^d} \frac{\alpha - \mathrm{div}\beta}{2} \rho = - \int_{\mathbb{T}^d} \frac{\phi}{\tau} \rho. $$

Summing the two equalities above, duality and the fact that $\rho$ solves (2.7) thus imply

$$ \sup(4.1) = \inf(2.7) = \frac{1}{2\tau}W_2^2(\rho, \rho_0) + F(\rho) = \int_{\mathbb{T}^d} \left(\frac{\mathrm{div}\beta - \alpha}{2}\right)^{c_{\tau}} d\rho_0 $$

so that $(\alpha, \beta)$ solves (4.1). Now if $(a, b)$ is another solution of (4.1), one should have

$$ F(\rho) = \int_{\mathbb{T}^d} \frac{a\rho + b \cdot \nabla \rho}{2} $$

implying that $a(x), b(x)$ maximizes $a\rho(x) + b \cdot \nabla\rho(x)$ subject to the constraint $a + \frac{1}{2}|b|^2 \le 0$, but since $\rho(x) > 0$, this constraint should be binding and there is a unique such maximizer which is $(\alpha(x), \beta(x))$ defined in (4.2). $\square$

## 5. NUMERICAL SOLUTIONS BY A PRIMAL-DUAL METHOD

The core idea of our numerical scheme is to discretize both the primal (1.3) and the dual (1.4) JKO step in space, and then apply the Chambolle-Pock algorithm to optimize in both problems alternatingly. It suffices to describe the discretization of the first JKO step, from the datum $\rho_0$ to the minimizer $\rho_1$.
5.1. **The discretized primal problem.** The torus $\mathbb{T}^d$ is modelled by a box $B = [0, 1]^d$, functions on $B$ extend periodically. For simplicity, we explain the discretization in $d = 1$ space dimensions, the translation to $d > 1$ is straight-forward but notationally cumbersome. We use a cartesian grid with $N$ vertices $x_j := (j + 0.5)h$, for $j = 1, \dots, N$, and uniform grid step $h = 1/N$. The centering of the grid points in the cells simplifies the implementation of the periodic boundary conditions. Functions $f: \mathbb{T}^d \to \mathbb{R}$ are discretized by finite sequences $(f_j)_{j=1,\dots,N}$ with $f_j \approx f(x_j)$. Introduce the $L^2$ inner product of functions on the grid by

$$ \langle f, g \rangle = \frac{1}{N} \sum_j f_j g_j. $$

The central finite difference quotient is denoted by $D_h$, i.e.,

$$ (D_h f)_j = \frac{f_{j+1} - f_{j-1}}{2h}. $$

Note that $D_h$ is anti self-adjoint with respect to the inner product, i.e., $-D_h^* = D_h$ is also the divergence.

The primal JKO problem (1.3) is now discretized as follows:

$$ (5.1) \qquad \inf_{\rho, m} \left\{ \frac{1}{N} \sum_{j=1}^{N} \frac{m_j^2}{2\rho_{0,j}} + \frac{\tau}{N} \sum_{j=1}^{N} \frac{(D_h \rho)_j^2}{4\rho_j} \right\} $$

$$ (5.2) \qquad \text{subject to} \quad \rho_j - \rho_{0,j} - (D_h^* m)_j = 0 \quad \text{for } j = 1, \dots, N. $$

In (5.1) above, the discretization of the Fisher information $F$ by finite differences in the second term is straight-forward. The relation of the first term in (5.1) and the constraint (5.2) to the dynamical representation (2.2) of $W_2^2(\rho_0, \rho_1)$ deserves some comment: recall that in (2.2), curves $(\rho_s, m_s)_{s \in [0,1]}$ depending on a fictitious time $s \in [0, 1]$ — unrelated to the physical time $t$ — have been considered to relax the problem. For small $\tau > 0$, we expect the initial and terminal densities $\rho_0$ and $\rho_1$ to be close, and the momentum $m_s \approx \bar{m}$ to be almost $s$-independent. Accordingly, the continuity equation $\partial_s \rho_s + \text{div} m_s = 0$ should be well-approximated by the stationary equation

$$ (5.3) \qquad \rho_1 - \rho_0 + \text{div} \bar{m} = 0. $$

For approximation of the $s$-integral in (2.2), the denominator $\rho_s$ should be replaced by some average of $\rho_0$ and $\rho_1$. The two canonical choices are explicit $(\rho_s \approx \rho_0$ as in [11]) and implicit $(\rho_s \approx \rho_1$ as in [13]). Either choice leads to a linearization of the Wasserstein distance in the form ($i = 0$ or 1)

$$ (5.4) \qquad \frac{1}{2} W_2^2(\rho_0, \rho_1) \approx \inf_{\bar{m}} \left\{ \int_{\mathbb{T}^d} \frac{|\bar{m}|^2}{2\rho_i} : -\text{div} \bar{m} = \rho_1 - \rho_0 \right\}. $$
The minimizer is explicitly given by $\bar{m} = \rho_i \nabla \phi$, where $\phi$ solves the linear elliptic equation $-\mathrm{div}(\rho_i \nabla \phi) = \rho_1 - \rho_0$ with periodic boundary conditions. In summary, this leads to the explicit ($i=0$) or implicit ($i=1$) approximation

$$W_2^2(\rho_0, \rho_1) \approx \|\rho_1 - \rho_0\|_{H_{\rho_i}^{-1}}^2$$

by a $\rho_i$-weighted $H^{-1}$-norm. We refer to [47] for a detailed comparison between the Wasserstein distance and its weighted $H^{-1}$ linearization. The discretization in (5.1)-(5.2) above has been built on the representation (5.4), in which we have chosen the numerically cheaper explicit form ($i=0$).

## 5.2. The discretized dual problem.
Somehow mimicking the continuous duality formulation, using the fact that the (discretized) linearized Wasserstein term can be written in dual form via the support function of the parabola $K$ as in (2.4), we can reformulate (5.1)-(5.2) in inf-sup form

$$(5.5) \quad \inf_{\rho,m} \sum_{j=1}^N \frac{m_j^2}{2\rho_{0,j}} + \sup_{\phi,(\alpha,\beta)\in K} \left\{ \langle \phi, \rho - \rho_0 - D_h^*(m) \rangle + \frac{\tau}{2} \langle \alpha, \rho \rangle + \frac{\tau}{2} \langle \beta, D_h(\rho) \rangle \right\}$$

by standard arguments (see chapter VI in [25]) the above Lagrangian has a saddle point and the dual is obtained by switching the inf and the sup. Doing so, we obtain that the saddle point satisfies

$$m = \rho_0 D_h \phi, \quad \phi + \frac{\tau}{2}(\alpha + D_h^* \beta) = 0$$

from which we can eliminate the variable $\alpha = -\frac{2\phi}{\tau} - D_h^* \beta$ so that $(\phi, \beta)$ solves the dual of (5.1)-(5.2) which can be written as

$$(5.6) \quad \inf_{(\phi,\beta)} \left\{ \frac{1}{2N} \sum_{j=1}^{N} \rho_{0,j} (D_h \phi)_j^2 + \langle \phi, \rho_0 \rangle : \left( -\frac{2\phi}{\tau} - D_h^* \beta, \beta \right) \in K \right\}$$

where of course, we have changed the sign in the dual to have a convex minimization problem instead of a concave maximization one.

As we shall see in the next paragraph, problem (5.6) is amenable to the Chambolle-Pock algorithm, and finding the optimal $\phi$ we can of course recover the optimal $\rho$ (nonnegative by construction since $\sigma_K(\rho, m) = +\infty$ if not) by

$$\rho = \rho_0 + D_h^*(\rho_0 D_h \phi).$$

## 5.3. Solution by the Chambolle-Pock algorithm.
Given two lsc convex functions $\Phi: \mathbb{R}^p \to \mathbb{R} \cup \{+\infty\}$, $\Psi: \mathbb{R}^n \to \mathbb{R} \cup \{+\infty\}$ and $\Lambda$ an $n \times p$ matrix, consider the convex minimization

$$(5.7) \quad \inf_{z \in \mathbb{R}^p} \{\Phi(z) + \Psi(\Lambda z)\} = \inf_{q \in \mathbb{R}^n} \left\{ \Phi(z) + \sup_{q \in \mathbb{R}^n} \{q \cdot \Lambda z - \Psi^*(q)\} \right\}.$$
Recalling that for a convex lsc function $G: \mathbb{R}^p \to \mathbb{R} \cup \{+\infty\}$, the proximal operator of $G$ is defined by

$$ \mathrm{prox}_G(x) = (\mathrm{id} + \partial G)^{-1}(x) := \underset{y \in \mathbb{R}^p}{\operatorname{argmin}} \left\{ \frac{1}{2} |x - y|^2 + G(y) \right\}, \forall x \in \mathbb{R}^p $$

and that $G$ is called simple whenever $\mathrm{prox}_G$ (or $\mathrm{prox}_{\sigma G}$) can be computed in closed form. By Moreau's identity (see [3]):

$$ \mathrm{prox}_{G^*} = \mathrm{id} - \mathrm{prox}_G $$

we see that $G^*$ is simple as soon as $G$ is. Note also that if $G$ is the characteristic function of a closed convex set $C$, $\mathrm{prox}_G$ is the projection onto $C$, hence, if projecting onto $C$ is closed form, both $G$ and its Legendre transform $G^* = \sigma_C$ (the support function of $C$) are simple.

The Chambolle-Pock primal-dual algorithm consists in iterating the proximal steps:

$$ (5.8) \qquad q^{t+1} = \mathrm{prox}_{\sigma\Psi^*}(q^t + \sigma\Lambda x^t) $$

$$ (5.9) \qquad z^{t+1} = \mathrm{prox}_{\delta\Phi}(z^t - \delta\Lambda^* q^{t+1}) $$

and the extrapolation step

$$ (5.10) \qquad x^{t+1} = 2z^{t+1} - z^t. $$

Chambolle and Pock proved in [19], that, provided:

* there exists a saddle point for (5.7) (i.e. a pair $(z, q)$ such that $q \in \partial\Psi(\Lambda z)$ and $-\Lambda^*q \in \partial\Phi(z)$),

* the (positive) steps $\sigma$ and $\delta$ are chosen such that $\sigma\delta\|\Lambda\|_{\mathrm{op}}^2 < 1$,

then, the Chambolle-Pock iterates $(z^t, q^t)$ converge to a saddle point (so that the first component solves (5.7) and the second solves its dual). Chambolle-Pock certainly is a good choice for solving (5.7) if $\Phi$ and $\Psi$ are nonsmooth but simple (note that contrarily to ADMM we never have to project on the range of $\Lambda$ but just to compute $\Lambda x^t$ and $\Lambda^* q^{t+1}$).

Let us now consider (5.6) and write it in the standard form (5.7), setting $z = (\phi, \beta)$, $\Lambda z = (-\frac{2\phi}{\tau} - D_h^*\beta, \beta)$ and remark that (5.6) then corresponds to

$$ \Phi(z) = \frac{1}{2N} \sum_{j=1}^{N} \rho_{0,j} (D_h \phi)_j^2 + \langle \phi, \rho_0 \rangle $$

and

$$ \Psi(\alpha, \beta) = \begin{cases} 0 & \text{if } (\alpha_j, \beta_j) \in K \text{ for } j = \dots, N \\ +\infty & \text{otherwise.} \end{cases} $$

Observe first that $\Psi$ is simple since the projection onto $K$ which can be performed on each grid point is almost explicit (see [4] and [46] for details). As for the weighted discretized Dirichlet energy $\Phi$, it is also simple since computing its
proximal operator amounts to solve a linear system discretizing a linear elliptic equation. More precisely $\text{prox}_{\delta\Phi}(\phi_0, \beta_0) = (\phi, \beta_0)$ where $\phi$ solves

$$ \phi + \delta D_h^*(\rho_0 D_h \phi) = \phi_0 - \delta \rho_0. $$

This shows the feasibility of the Chambolle-Pock algorithm for (5.6) which is illustrated in the next paragraph. Note that the convergence is ensured but can be quite slow, see Figure 1 for convergence curves and a discussion associated with the test cases of section 5.4.

### 5.4. Simulations.
The method described above is used to perform simulations in $d=1$. We use the following set of test cases:

$$ \begin{array}{ll} \text{Test 1:} & \rho_0(x) = \cos^2(\pi x) \\ \text{Test 2:} & \rho_0(x) = 10^{-3} + \cos^{16}(\pi x) \\ \text{Test 3:} & \rho_0(x) = 10^{-2} + \cos^{64}(\pi x) \\ \text{Test 4:} & \rho_0(x) = 10^{-1} + \cos^{16}(\pi x) - \sin^2(2\pi x) \end{array} \qquad (5.11) $$

We have been using $N=128$ grid points in space (i.e. $dx = 1/128$). The evolutions of the densities are shown in figure 2. Test case 1 and 2 were treated in [8] and results are consistent with ours (see figure 1 page 15 in [8]).

As expected the relaxation towards the constant equilibrium is very fast. The Finite Difference approach in [8] uses a time discretization $\tau$ of the order of $dx^4$ (the authors of [8] take $N=100$ and $\tau = 10^{-8}$). For the experiments using the method described above, it seems sufficient to start with $\tau = 10^{-7}$ to reach the first step and then increase it by 10 at each represented snapshot in figure 1. For all presented test cases, Chambolle-Pock Algorithm has been applied to 40 times.

The convergence of Chambolle-Pock for each JKO step is however slow. for the first JKO time step of our simulation, we converge to machine precision in $10^5$ iterations. The convergence curves are presented in figure 2. The situation improves as time $t$ increases and the last JKO steps are computed with “only” $10^3$ Chambolle-Pock iterations. Even though each iteration is of order $O(N)$ the slow convergence represents a difficulty in view of the extension of this method to higher dimensions.

## ACKNOWLEDGMENTS

GC acknowledges the support of the Lagrange Mathematics and Computing Research Center. DM's research was supported by the German Research Foundation (DFG), Collaborative Research Center SFB-TR 109.
<figure><img src="image_3.png" /><figcaption>FIGURE 1. Convergence for one JKO step starting with 5.11). Test 1 to 4, left to right and top to bottom. In $||\cdot||_\infty$ norm we show the iteration increment on $F$, and also the primal and dual gradient error respectively $||\Lambda^*q + \partial\Phi(z)||_\infty$ and $||\Lambda z - \partial\Psi^*(q)||_\infty$</figcaption></figure>

REFERENCES

[1] Luigi Ambrosio, Nicola Gigli, and Giuseppe Savaré. *Gradient flows: in metric spaces and in the space of probability measures*. Springer Science & Business Media, 2005.

[2] MG Ancona. Diffusion-drift modeling of strong inversion layers. *COMPEL-The international journal for computation and mathematics in electrical and electronic engineering*, 6(1):11-18, 1987.

[3] Heinz H. Bauschke and Patrick L. Combettes. *Convex analysis and monotone operator theory in Hilbert spaces*. CMS Books in Mathematics/Ouvrages de Mathématiques de la SMC. Springer, New York, 2011. With a foreword by Hédy Attouch.

[4] Jean-David Benamou and Yann Brenier. A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem. *Numer. Math.*, 84(3):375-393, 2000.

[5] Jean-David Benamou and Guillaume Carlier. Augmented Lagrangian methods for transport optimization, mean field games and degenerate elliptic equations. *J. Optim. Theory Appl.*, 167(1):1-26, 2015.
<figure><img src="image_4.png" /><figcaption>FIGURE 2. DLSS with initialisation (5.11), Test 1 to 4, left to right and top to bottom.</figcaption></figure>

[6] Jean-David Benamou, Guillaume Carlier, and Maxime Laborde. An augmented Lagrangian approach to Wasserstein gradient flows and applications. In *Gradient flows: from theory to application*, volume 54 of *ESAIM Proc. Surveys*, pages 1–17. EDP Sci., Les Ulis, 2016.

[7] Jean-David Benamou, Guillaume Carlier, Quentin Mérigot, and Édouard Oudet. Discretization of functionals involving the Monge-Ampère operator. *Numer. Math.*, 134(3):611–636, 2016.

[8] Pavel M Bleher, Joel L Lebowitz, and Eugene R Speer. Existence and positivity of solutions of a fourth-order nonlinear pde describing interface fluctuations. *Communications on Pure and Applied Mathematics*, 47(7):923–942, 1994.

[9] Mario Bukal. Well-posedness and convergence of a numerical scheme for the corrected derrida-lebowitz-speer-spohn equation using the hellinger distance. *Discrete & Continuous Dynamical Systems*, 41(7):3389, 2021.

[10] Mario Bukal, Etienne Emmrich, and Ansgar Jüngel. Entropy-stable and entropy-dissipative approximations of a fourth-order quantum diffusion equation. *Numerische Mathematik*, 127(2):365–396, October 2013.

[11] Martin Burger, José A Carrillo, and Marie-Therese Wolfram. A mixed finite element method for nonlinear diffusion equations. *Kinetic and Related Models*, 3(1):59–83, 2009.

[12] G. Buttazzo, C. Jimenez, and E. Oudet. An optimization problem for mass transportation with congested dynamics. *SIAM J. Control Optim.*, 48(3):1961–1976, 2009.
[13] Clément Cancès, Thomas Gallouët, and Gabriele Todeschi. A variational finite volume scheme for Wasserstein gradient flows. *Numerische Mathematik*, 146(3):pp 437–480, 2020.

[14] Guillaume Carlier and Clarice Poon. On the total variation Wasserstein gradient flow and the TV-JKO scheme. *ESAIM Control Optim. Calc. Var.*, 25:Paper No. 42, 21, 2019.

[15] José A. Carrillo, Ansgar Jüngel, and Shaoqiang Tang. Positive entropic schemes for a non-linear fourth-order parabolic equation. *Discrete Contin. Dyn. Syst. Ser. B*, 3(1):1–20, 2003.

[16] José A Carrillo, Helene Ranetbauer, and Marie-Therese Wolfram. Numerical simulation of nonlinear continuity equations by evolving diffeomorphisms. *Journal of Computational Physics*, 327:186–202, 2016.

[17] José A Carrillo and Giuseppe Toscani. Long-time asymptotics for strong solutions of the thin film equation. *Communications in mathematical physics*, 225:551–571, 2002.

[18] Fausto Cavalli and Giovanni Naldi. A wasserstein approach to the numerical solution of the one-dimensional cahn-hilliard equation. *Kinet. Relat. Models*, 3(1):123–142, 2010.

[19] Antonin Chambolle and Thomas Pock. A first-order primal-dual algorithm for convex problems with applications to imaging. *J. Math. Imaging Vision*, 40(1):120–145, 2011.

[20] Dario Cordero-Erausquin. Sur le transport de mesures périodiques. *C. R. Acad. Sci. Paris Sér. I Math.*, 329(3):199–202, 1999.

[21] Pierre Degond, Florian Méhats, and Christian Ringhofer. Quantum energy-transport and drift-diffusion models. *Journal of statistical physics*, 118:625–667, 2005.

[22] B. Derrida, J. L. Lebowitz, E. R. Speer, and H. Spohn. Fluctuations of a stationary nonequilibrium interface. *Phys. Rev. Lett.*, 67(2):165–168, 1991.

[23] Bernard Derrida, Joel L Lebowitz, Eugene R Speer, and Herbert Spohn. Dynamics of an anchored toom interface. *Journal of Physics A: Mathematical and General*, 24(20):4805, 1991.

[24] Bertram Düring, Daniel Matthes, and Josipa-Pina Milisić. A gradient flow scheme for nonlinear fourth order equations. *Discrete and Continuous Dynamical Systems-series B*, 14:935–959, 2010.

[25] Ivar Ekeland and Roger Temam. *Convex analysis and variational problems*. Studies in Mathematics and its Applications, Vol. 1. North-Holland Publishing Co., Amsterdam-Oxford; American Elsevier Publishing Co., Inc., New York, 1976. Translated from the French.

[26] Julian Fischer. Uniqueness of solutions of the derrida-lebowitz-speer-spohn equation and quantum drift-diffusion models. *Communications in Partial Differential Equations*, 38(11):2004–2047, 2013.

[27] Thomas O Gallouët, Quentin Merigot, and Andrea Natale. Convergence of a lagrangian discretization for barotropic fluids and porous media flow. *SIAM Journal on Mathematical Analysis*, 54(3):2990–3018, 2022.

[28] Ugo Gianazza, Giuseppe Savaré, and Giuseppe Toscani. The Wasserstein gradient flow of the Fisher information and the quantum drift-diffusion equation. *Arch. Ration. Mech. Anal.*, 194(1):133–220, 2009.

[29] David Gilbarg and Neil S. Trudinger. *Elliptic partial differential equations of second order*. Classics in Mathematics. Springer-Verlag, Berlin, 2001. Reprint of the 1998 edition.

[30] Laurent Gosse and Giuseppe Toscani. Lagrangian numerical approximations to one-dimensional convolution-diffusion equations. *SIAM Journal on Scientific Computing*, 28(4):1203–1227, 2006.

[31] Xiaokai Huo and Hailiang Liu. A positivity-preserving and energy stable scheme for a quantum diffusion equation. *Numer. Methods Partial Differential Equations*, 37(6):2973–2999, 2021.

[32] Richard Jordan, David Kinderlehrer, and Felix Otto. The variational formulation of the Fokker-Planck equation. *SIAM J. Math. Anal.*, 29(1):1–17, 1998.
[33] Oliver Junge, Daniel Matthes, and Horst Osberger. A fully discrete variational scheme for solving nonlinear fokker-planck equations in higher space dimensions. *arXiv preprint arXiv:1509.07721*, 2015.

[34] Ansgar Jüngel and Daniel Matthes. An algorithmic construction of entropies in higher-order nonlinear pdes. *Nonlinearity*, 19(3):633–2006.

[35] Ansgar Jüngel and Daniel Matthes. The derrida–lebowitz–speer–spohn equation: Existence, nonuniqueness, and decay rates of the solutions. *SIAM Journal on Mathematical Analysis*, 39(6):1996–2015, 2008.

[36] Ansgar Jüngel and René Pinnau. Global nonnegative solutions of a nonlinear fourth-order parabolic equation for quantum systems. *SIAM Journal on Mathematical Analysis*, 32(4):760–777, 2000.

[37] Ansgar Jüngel and René Pinnau. A positivity-preserving numerical scheme for a nonlinear fourth order parabolic system. *SIAM J. Numer. Anal.*, 39(2):385–406, 2001.

[38] Ansgar Jüngel and Ingrid Violet. First-order entropies for the Derrida-Lebowitz-Speer-Spohn equation. *Discrete Contin. Dyn. Syst. Ser. B*, 8(4):861–877, 2007.

[39] Wuchen Li, Jianfeng Lu, and Li Wang. Fisher information regularization schemes for wasserstein gradient flows. *Journal of Computational Physics*, 416:109449, sep 2020.

[40] Jan Maas and Daniel Matthes. Long-time behavior of a finite volume discretization for a fourth order diffusion equation. *Nonlinearity*, 29(7):1992, 2016.

[41] Daniel Matthes, Robert J. McCann, and Giuseppe Savaré. A family of nonlinear fourth order equations of gradient flow type. *Comm. Partial Differential Equations*, 34(10-12):1352–1397, 2009.

[42] Daniel Matthes and Horst Osberger. Convergence of a variational lagrangian scheme for a nonlinear drift diffusion equation. *ESAIM: Mathematical Modelling and Numerical Analysis*, 48(3):697–726, 2014.

[43] Daniel Matthes and Horst Osberger. A convergent lagrangian discretization for a nonlinear fourth-order equation. *Foundations of Computational Mathematics*, 17(1):73–126, 2017.

[44] Robert J. McCann. Polar factorization of maps on Riemannian manifolds. *Geom. Funct. Anal.*, 11(3):589–608, 2001.

[45] Horst Osberger and Daniel Matthes. 14 convergence of a fully discrete variational scheme for a thin-film equation. *Topological Optimization and Optimal Transport: In the Applied Sciences*, 17:356, 2017.

[46] Nicolas Papadakis, Gabriel Peyré, and Edouard Oudet. Optimal transport with proximal splitting. *SIAM J. Imaging Sci.*, 7(1):212–238, 2014.

[47] Rémi Peyre. Comparison between $W_2$ distance and $\dot{H}^{-1}$ norm, and localization of Wasserstein distance. *ESAIM Control Optim. Calc. Var.*, 24(4):1489–1501, 2018.

[48] Filippo Santambrogio. *Optimal transport for applied mathematicians*, volume 87 of *Progress in Nonlinear Differential Equations and their Applications*. Birkhäuser/Springer, Cham, 2015. Calculus of variations, PDEs, and modeling.

[49] C. Villani. *Topics in Optimal Transportation*, volume 58 of *Graduate Studies in Mathematics*. American Mathematical Society, Providence, 2003.

CEREMADE, UMR CNRS 7534, UNIVERSITÉ PARIS DAUPHINE, PSL, PL. DE LATTRE DE TASSIGNY, 75775 PARIS CEDEX 16, FRANCE AND INRIA-PARIS, MOKAPLAN

Email address: carlier@ceremade.dauphine.fr
CEREMADE, UMR CNRS 7534, UNIVERSITÉ PARIS DAUPHINE, PSL, INRIA-PARIS, MOKAPLAN

*Email address: jean-david.benamou@inria.fr*

DEPARTMENT OF MATHEMATICS, SCHOOL OF CIT, TECHNISCHE UNIVERSITÄT MÜNCHEN

*Email address: matthes@ma.tum.de*