The Deep Minimizing Movement SchemeH. J. Hwang, C. Kim, M. S. Park, H. Son

# The Deep Minimizing Movement Scheme

Hyung Ju Hwang222Department of Mathematics, Pohang University of Science and
Technology, Pohang 790-784, Republic of Korea (,,). 444Graduate School of
Artificial Intelligence, Pohang University of Science and Technology, Pohang
790-784, Republic of Korea (). 555Corresponding Author Cheolhyeong
Kim222Department of Mathematics, Pohang University of Science and Technology,
Pohang 790-784, Republic of Korea (,,). Min Sue Park222Department of
Mathematics, Pohang University of Science and Technology, Pohang 790-784,
Republic of Korea (,,). Hwijae Son333Stochastic Analysis and Application
Research Center, Korea Advanced Institute of Science and Technology, Daejeon,
Republic of Korea<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Solutions of certain partial differential equations (PDEs) are often
represented by the steepest descent curves of corresponding functionals.
Minimizing movement scheme was developed in order to study such curves in
metric spaces. Especially, Jordan-Kinderlehrer-Otto studied the Fokker-Planck
equation in this way with respect to the Wasserstein metric space. In this
paper, we propose a deep learning-based minimizing movement scheme for
approximating the solutions of PDEs. The proposed method is highly scalable
for high-dimensional problems as it is free of mesh generation. We demonstrate
through various kinds of numerical examples that the proposed method
accurately approximates the solutions of PDEs by finding the steepest descent
direction of a functional even in high dimensions.

###### keywords:

Minimizing Movement Scheme, JKO Scheme, Neural Networks

68T07

##  1 Introduction

Starting from the pioneering concerns given in [12], the theory of gradient
flow evolution in general metric spaces has been intensively studied to unify
many problems in calculus of variations, differential equations, mean
curvature flow, etc. In this trend, minimizing movement scheme [11] was
devised as a general method for the study of steepest descent curves of a
functional $F$ in a metric space $(X,d_{X})$, and this method was extensively
applied to various differential equations as a time-implicit semi-discrete
scheme.

Historically, the application of minimizing movement scheme to differential
equations was exclusively based on gradient flows arising from a separable
Hilbert space, particularly $L^{2}(\mathbb{R}^{n})$. Later, in [18], Jordan,
Kinderlehrer, and Otto observed that the Fokker-Planck equation can be
understood as a steepest descent for the free energy with respect to the
Wasserstein distance and utilized the minimizing movement scheme for
Wasserstein spaces, and that is why the minimizing movement scheme is now
widely known as the JKO-scheme (Jordan, Kinderlehrer and Otto). In recent
years, it has been known that a large class of diffusion equations can be
variationally reformulated as gradient flows in $L^{2}$-Wasserstein spaces,
and the minimizing movement scheme has now become a popular tool for deriving
numerical schemes for these equations.

Neural networks have been considered as candidates for approximators of the
solutions of PDEs in past decades. The earliest work dates back to the 1990s,
for example, [21, 13]. However, due to the limitation of computing resources
and the absence of an efficient training algorithm, the proposed neural
network solvers were often undervalued at that time. Along with the recent
advances in deep learning theory such as the back-propagation and the
stochastic gradient descent algorithm, several approaches for solving PDEs via
neural networks have been revisited. One major line of research is Physics
Informed Neural Networks (PINNs) [39], which proposed to train a neural
network to minimize the sum of PDE residual functions. Plenty of works are
then reported regarding the convergence property of continuous loss
formulation [17, 16, 22], dealing with the exact imposition of the boundary
condition [4, 35], works on improving training efficiency [42, 46, 30, 44],
and applications to high-dimensional problems [41] to name a few. Another
approach, called the deep Ritz method, reformulates a PDE problem into an
optimization problem of a variational functional [47], and a corresponding
convergence result is reported in [34]. Neural network solutions of
variational problems with essential boundary conditions are studied in [25],
and a pre-training strategy for the deep Ritz method condition is proposed in
[9].

Traditional numerical schemes such as Finite Difference Method (FDM), Finite
Element Method (FEM), Finite Volume Method (FVM) require fine mesh generation
for an accurate approximation. Since mesh generation and computation become
intractable as the computational dimension increases, they often fail to solve
PDEs in high dimensions [27]. On the other hand, the aforementioned neural
network methods have successful approximations in high-dimensional settings.
For example, [41] proposed an iterative sampling technique together with the
stochastic gradient descent to tackle the high-dimensional PDEs, and [41, 42]
showed that neural networks could accurately approximate the solution in high
dimensions. Additionally, [47] considers variational problems in high
dimensions for elliptic PDEs. In this paper, we consider several high-
dimensional PDEs as example problems to demonstrate the empirical success of
our method.

Optimizing a functional and finding a minimizer have been widely studied both
in the machine learning community and in the PDE literature. Recently,
optimization of a functional defined on a space of probability measures
equipped with $L^{2}$-Wasserstein distance arises as an important problem.
Therefore, Wasserstein gradient flow plays an important role in studying the
steepest descent of a target functional in the space of probability measures.
[1] and [32] proposed to applying JKO scheme with neural networks with
tractable algorithms. However, their usage of JKO scheme mainly focused on a
fast approximation of the steepest descent of the target functional, not on
the accurate prediction of a whole dynamic system. For example, the experiment
results of [1] on the Wasserstein gradient flow type PDEs were mainly
discussed on its asymptotic convergence of the objective functions, but the
accurate prediction of the solution $u(t,x)$ was not mentioned. Meanwhile,
[26] proposed to solve Fokker-Planck equation with neural networks by emerging
JKO scheme, Wasserstein Gradient Flow, and its geometric interpretation.
However, its usage is limited to Fokker-Planck equation, since [26] projected
the Fokker-Planck equation onto a parametric manifold with a new metric tensor
and viewed the equation as a new ODE system on the manifold.

In this paper, we propose a neural network model which reflects the minimizing
movement scheme for solving PDEs. Our proposed method has several advantages
in the following aspects. Firstly, our method is more suitable for solving
high-dimensional problems than the traditional mesh-based schemes [26] because
the mesh generation is not necessary for the proposed method, and neural
networks are capable of handling high-dimensional objects in general.
Secondly, while specifying boundary condition often causes a problem in both
residual minimization [35] and the deep Ritz method [9], it often follows
naturally from the minimizing movement scheme in our method. Therefore, we do
not need to explicitly specify boundary conditions for those PDEs derivable
from the minimizing movement scheme. Lastly, our model easily enables time
extrapolation of the numerical solution in contrast to residual minimization;
we only need to solve one more optimization problem for one-step time
extrapolation, while the residual minimization requires additional training in
the whole time-space domain. Furthermore, an increment of the error at each
time step is bounded by $\mathcal{O}(\tau^{1/2})$, where $\tau$ denotes the
step size. See, theorem 3.11.

##  2 Preliminaries

###  2.1 Gradient flow

In general, given a Riemannian manifold $(M,g)$, a curve
$u:(0,\infty)\rightarrow M$ is said to be the gradient flow of a smooth
functional $\mathcal{F}:M\rightarrow\mathbb{R}$ if it follows the direction in
which $\mathcal{F}$ decreases at most. Formally, this definition can be
interpreted by the evolution equation

(1) $\frac{d}{dt}u(t)=-\nabla_{M}\mathcal{F}(u(t)).$

If $M$ is a Hilbert space, the gradient of $\mathcal{F}$ is given by its
functional derivative, which is derived from the Riesz representation of the
Fréchet derivative $D\mathcal{F}$ of $\mathcal{F}$:
$D\mathcal{F}(x)(y)=\langle\nabla_{M}\mathcal{F}(x),y\rangle$ for every
$x,y\in M$. Therefore, for $\Omega\subset\mathbb{R}^{n}$,
$\nabla_{L^{2}(\Omega)}\mathcal{F}(u)$ is merely the functional derivative
$\frac{\delta\mathcal{F}}{\delta u}(u)$ since $L^{2}(\Omega)$ is itself a
Hilbert space.

For example, the heat equation $\frac{d}{dt}u(t)=\Delta u$ is the gradient
flow of the Dirichlet energy $\mathcal{F}(u):=\frac{1}{2}\int||\nabla
u||^{2}dx$ in $L^{2}(\mathbb{R}^{n})$ because
$\nabla_{L^{2}(\mathbb{R}^{n})}\mathcal{F}(u)=\frac{\delta\mathcal{F}}{\delta
u}(u)=-\Delta u$.

As illustrated in the above example, to reformulate an evolution equation into
its variational form as in (1), finding an appropriate Riemannian manifold
$(M,g)$ would be helpful. An important and interesting example of such
Riemannian manifold was observed in [37] that there is a metric tensor $g$
which induces the $L^{2}$-Wasserstein metric $\mathbb{W}_{2}$ on
$\mathcal{P}_{2}(\Omega)$, the space of Borel probability distributions on
$\Omega\subset\mathbb{R}^{n}$ with finite second moments. Here, $\Omega$ is
assumed to be either an open bounded subset of $\mathbb{R}^{n}$ with smooth
boundary or $\mathbb{R}^{n}$. Abusing notation and identifying the absolutely
continuous measure $u$ with its density as $du(x)=u(x)dx$, it was shown in
[37] that, under some assumptions on $\mathcal{E}$, the gradient
$\nabla_{\mathbb{W}_{2}}\mathcal{E}(u)$ can be explicitly expressed by

(2)
$\nabla_{\mathbb{W}_{2}}\mathcal{F}(u)=-\nabla\cdot(u\nabla\frac{\delta\mathcal{F}}{\delta
u}(u)).$

In general, many interesting linear functionals $\mathcal{F}$ on
$(\mathcal{P}_{2}(\mathbb{R}^{n}),\mathbb{W}_{2})$ are given as linear
combinations of three types of basic functionals; potential energy
$\mathcal{V}(u):=\int Vu$, interaction energy $\mathcal{W}(u):=\int(W\ast
u)u$, and internal energy $\mathcal{E}(u):=\int u\phi(u)$, where these
energies are defined for absolutely continuous $u$’s and defined by $\infty$
otherwise. For example, the heat equation $\frac{d}{dt}u(t)=\nabla u$ is also
the gradient flow of the entropy
$\mathcal{E}:(\mathcal{P}_{2}(\mathbb{R}^{n}),\mathbb{W}_{2})\rightarrow[0,+\infty]$
defined by $\int u(x)\log u(x)dx$ for absolutely continuous $u$ and $\infty$
otherwise, because $\nabla\cdot(u\nabla\frac{\delta\mathcal{E}}{\delta
u}(u))=\Delta u$.

Examples of basic functionals and their corresponding gradient flows on
$L^{2}(\mathbb{R}^{n})$ and $(\mathcal{P}_{2}(\Omega),\mathbb{W}_{2})$ are
given in Table 1.

| Functional | Gradient Flow
---|---|---
$L^{2}(\mathbb{R}^{n})$ |  $\mathcal{F}(u)=\int\phi(u)$ |  $\frac{d}{dt}u(t)=-\phi^{\prime}(u)$
$\mathcal{F}(u)=\int||\nabla u||^{2}$ |  $\frac{d}{dt}u(t)=\Delta u$
$\mathcal{F}(u)=\int u(x)f(x)$ |  $\frac{d}{dt}u(t)=-f(x)$
$\mathbb{W}_{2}$ |  $\mathcal{F}(u)=\int u\log u$ |  $\frac{d}{dt}u(t)=\Delta u$
$\mathcal{F}(u)=\frac{1}{m-1}\int u^{m}$ |  $\frac{d}{dt}u(t)=\Delta u^{m}$
$\mathcal{F}(u)=\int Vu$ |  $\frac{d}{dt}u(t)=\nabla\cdot(\nabla Vu)$
$\mathcal{F}(u)=\int(W\ast u)u$ |  $\frac{d}{dt}u(t)=\nabla\cdot(\nabla(W\ast u)u)$
Table 1: Examples of functionals and their corresponding gradient flows.

###  2.2 Minimizing movement scheme (JKO scheme)

A metric space $(X,d_{X})$ is said to be a geodesic space if $d_{X}(x,y)$ is a
geodesic between $x$ and $y\in X$ for every $x$ and $y\in X$. It is well-known
that if a metric space $X$ is a geodesic Polish space, then
$(\mathcal{P}_{2}(X),\mathbb{W}_{2})$ is a geodesic Polish space as well, and
any Hilbert space is a geodesic space [2]. Therefore, for any convex subset
$\Omega\subset\mathbb{R}^{n}$, $L^{2}(\Omega)$ and
$(\mathcal{P}_{2}(\Omega),\mathbb{W}_{2})$ are geodesic Polish spaces.

Let $(X,d_{X})$ be a geodesic Polish space with Riemannian structure and
$F:X\rightarrow\mathbb{R}\cup\\{+\infty\\}$ be a lower semi-continuous
$\lambda$-geodesically convex function (under suitable compactness assumptions
to guarantee the existence of a minimum), and iteratively define

(3) $x_{k+1}^{\tau}=\operatorname*{arg\,min}_{x\in
X}\\{\mathcal{F}(x)+\frac{d_{X}(x,x_{k}^{\tau})^{2}}{2\tau}\\},\quad
x_{0}^{\tau}=x_{0},$

for $\tau>0$ and some initial point $x_{0}\in X$ such that $E(x_{0})<+\infty$.

This scheme is said to be the minimizing movement scheme. If we define a curve
such that $x^{\tau}(k\tau)=x_{k}^{\tau}$ and such that $x^{\tau}$ restricted
on any interval $[k\tau,(k+1)\tau]$ is a constant-speed geodesic with speed
$\frac{d(x_{k}^{\tau},d_{k+1}^{\tau})}{\tau}$, up to a subsequence
$\tau_{j}\to 0$, then $x^{\tau}$ converges uniformly to the solution of (1)
when some regularity conditions on $\mathcal{F}$ and $(X,d_{X})$ are assumed.

Moreover, a gradient flow derived through the minimizing movement scheme in
$L^{2}(\Omega)$ is always accompanied by the boundary condition
$\frac{\delta\phi}{\delta(\nabla u)}\cdot\hat{\boldsymbol{n}}=0$ on
$\partial\Omega$ when the functional $\mathcal{F}$ is given in the form
$\mathcal{F}(u):=\int\phi(u,\nabla u)$, if $\Omega$ is not $\mathbb{R}^{n}$
itself [5]. Likewise, a gradient flow derived through the minimizing movement
scheme in $(\mathcal{P}_{2}(\Omega),\mathbb{W}_{2})$ is always accompanied by
no-flux boundary condition on $\partial\Omega$, i.e.
$u\nabla\frac{\delta\mathcal{F}}{\delta u}(u)\cdot\hat{\boldsymbol{n}}=0$, if
$\Omega$ is not $\mathbb{R}^{n}$ itself [40].

In particular, the gradient flows of the heat equation and the Allen-Cahn
equation derived through the minimizing movement scheme in the space
$L^{2}(\Omega)$ satisfy the Neumann boundary condition [5, 31]. Similarly, the
gradient flow of the heat equation derived through the minimizing movement
scheme in the space $(\mathcal{P}_{2}(\Omega),\mathbb{W}_{2})$ satisfies the
Neumann boundary condition [2]. In this paper, we consider Neumann boundary
conditions when $\Omega\neq\mathbb{R}^{n}$.

###  2.3 Optimal transport and $\boldsymbol{L^{2}}$-Wasserstein distance

Let $\mu,\nu$ be Borel probability measures on Polish spaces $X,Y$
respectively. Given a cost function $c:X\times Y\rightarrow[0,+\infty]$, where
$c(x,y)$ measures the cost of transporting one unit of mass from $x\in X$ to
$y\in Y$, the optimal transport problem is how to tranposrt $\mu$ to $\nu$
while minimizing the cost $c$. This was firstly formulated by Monge [33] as
follows:

(4) $(MP)\quad\inf_{T_{\ast}\mu=\nu}\int_{X}c(x,T(x))d\mu(x),$

where $T:X\rightarrow Y$ is a measurable map (transport map).

However, this problem is ill-posed because of the constraint
$T_{\ast}\mu=\nu$. For instance, if $\mu$ is a Dirac measure and $\nu$ is not,
there is no admissible $T$. Later, Kantorovich proposed the following way to
relax the Monge problem:

(5) $(KP)\quad\min_{\gamma\in\Pi(\mu,\nu)}\int_{X\times Y}c(x,y)d\gamma(x,y),$

where $\Pi(\mu,\nu)$ is the collection of joint distributions (transport
plans) of $\mu$ and $\nu$.

Moreover, if $\mu$ is non-atomic and $c$ is continuous, then there always
exist optimal transport plans satisfying (5) and it coincides with the infimum
of (4) [38]. In particular, this is the case when $X=Y=\mathbb{R}^{n}$ and
$c(x,y):=||x-y||_{2}^{2}$. The $L^{2}$-Wasserstein distance between $\mu$ and
$\nu$ is defined to be the square root of the minimum of (5) when
$c(x,y):=||x-y||_{2}^{2}$, and is denoted by $\mathbb{W}_{2}(\mu,\nu)$.

Estimation of $L^{2}$-Wasserstein distance is essential for solving PDEs with
JKO scheme and applying optimal transport in a machine learning area. Previous
studies [43, 28] proposed to estimate this $L^{2}$-Wasserstein distance by
applying Kantorovich-Rubinstein duality theorem [45] and approximating
primal/conjugate convex functions. However, their computations were based on
minimax optimization, which leads to a heavy computational cost and a slow
estimation speed.

Meanwhile, the authors in [20] proposed Wasserstein-2 generative networks
(W2GN), an end-to-end non-minimax algorithm for training optimal transport
mappings for $L^{2}$-Wasserstein distance. They avoided the minimax
optimization by introducing a regularization term on the convex conjugate
potential. Instead of optimizing minimax problem, W2GN simply minimizes the
cost and controls the conjugate function. Compared with the previous minimax-
based computation algorithms in [43, 28], W2GN significantly improves the
computational cost for estimating $L^{2}$-Wasserstein distance.

In the rest of this subsection, we briefly describe the method for estimating
$\mathbb{W}_{2}(\mu,\nu)$, that was proposed in [20].

When absolutely continuous distributions $u,\rho$ with finite second moments
are given, by re-arranging the Monge’s formulation (4),
$\mathbb{W}_{2}(u,\rho)$ can be expressed by

(6)
$\mathbb{W}_{2}^{2}(u,\rho)=\int_{\mathbb{R}^{n}}||x||^{2}du(x)+\int_{\mathbb{R}^{n}}||y||^{2}d\rho(y)-2\cdot\max_{T_{\ast}u=v}\int_{\mathbb{R}^{n}}\langle
x,T(x)\rangle du(x).$

We denote the maximum term on the right hand side of (6) by
$\text{Corr}(u,\rho)$, and $-\text{Corr}(u,\rho)$ can be regarded as (4) for
$c(x,y):=-\langle x,y\rangle$ [29]. Therefore, by applying the Kantorovich-
Rubinstein duality theorem, we have

(7)
$\text{Corr}(u,\rho)=\min_{\psi\in\text{Convex}}\bigg{[}\int_{\mathbb{R}^{n}}\psi(x)du(x)+\int_{\mathbb{R}^{n}}\overline{\psi}(y)d\rho(y)\bigg{]},$

where the minimum is taken over all the convex functions
$\psi:\mathbb{R}^{n}\rightarrow\mathbb{R}\cup\\{+\infty\\}$, and
$\overline{\psi}(y):=\max_{x\in\mathbb{R}^{n}}(\langle x,y\rangle-\psi(x))$ is
the Fenchel convex conjugate to $\psi$.

According to [45, 29], the gradient $T^{*}:=\nabla\psi^{*}$ of the optimal
$\psi^{*}$ readily gives the maximizer of (6), and we have

(8) $(T^{*})^{-1}(y)=(\nabla\psi^{*})^{-1}(y)=\nabla\overline{\psi^{*}}(y).$

Therefore, it suffices to find the minimizer $\psi$ of (7) to evaluate the
$L^{2}$-Wasserstein distance $\mathbb{W}_{2}(u,\rho)$. However, to do this, we
need to solve the optimization sub-problem
$\overline{\psi}(y)=\max_{x\in\mathbb{R}^{n}}(\langle x,y\rangle-\psi(x))$,
and this would lead to a minimax problem

$\displaystyle\min_{\psi\in\text{Convex}}\bigg{[}\int_{\mathbb{R}^{n}}\psi_{\theta}(x)du(x)+\int_{\mathbb{R}^{n}}\bar{\psi}(y)d\rho(y)\bigg{]}$
$\displaystyle=\min_{\psi\in\text{Convex}}\bigg{[}\int_{\mathbb{R}^{n}}\psi_{\theta}(x)du(x)+\int_{\mathbb{R}^{n}}[\max_{x\in\mathbb{R}^{n}}\langle
x,y\rangle-\psi(x)]d\rho(y)\bigg{]}$ (9)
$\displaystyle=\min_{\psi\in\text{Convex}}\bigg{[}\int_{\mathbb{R}^{n}}\psi_{\theta}(x)du(x)+\max_{S:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}}\int_{\mathbb{R}^{n}}[\langle
S(y),y\rangle-\psi(S(y))]d\rho(y)\bigg{]}$

where the maximum is taken over arbitrary measurable functions $S$.

Nevertheless, [20] proposed a trick to convert this problem into an end-to-end
non-minimax problem by considering the optimal $\psi$ and $\bar{\psi}$. For
the variational approximation in (2.3), we obtain a variational lower bound
which matches the entire value for
$S=(\nabla\psi)^{-1}(y)=\nabla\overline{\psi}(y)$.

Therefore, the primal potential $\psi$ and its conjugate $\overline{\psi}$ can
be approximated by two parametrized convex functions $\psi_{\theta}$ and
$\overline{\psi_{\omega}}$ respectively, if we minimize the following
objective:

(10)
$\min_{\theta,\omega}\bigg{[}\int_{\mathbb{R}^{n}}\psi_{\theta}(x)du(x)+\int_{\mathbb{R}^{n}}[\langle\nabla\overline{\psi_{\omega}}(y),y\rangle-\psi_{\theta}(\nabla\overline{\psi_{\omega}}(y))]d\rho(y)\bigg{]}.$

Moreover, the relation (8) must be imposed on the optimized pair
$(\psi_{\theta},\overline{\psi_{\omega}})$, and this can be done by
additionally minimizing the following regularization:

(11)
$R(\theta,\omega):=\int_{\mathbb{R}^{n}}||\nabla\psi_{\theta}\circ\nabla\overline{\psi_{\omega}}(y)-y||^{2}d\rho(y).$

In conclusion, our final objective would be the following:

(12)
$\min_{\theta,\omega}\bigg{[}\int_{\mathbb{R}^{n}}\psi_{\theta}(x)du(x)+\int_{\mathbb{R}^{n}}[\langle\nabla\overline{\psi_{\omega}}(y),y\rangle-\psi_{\theta}(\nabla\overline{\psi_{\omega}}(y))]d\rho(y)+\lambda\cdot
R(\theta,\omega)\bigg{]},$

where $\lambda>0$ is a hyperparameter.

Using the identities (6) and (7), once the optimal
$\psi_{\theta^{*}},\overline{\psi_{\omega^{*}}}$ are found,
$\mathbb{W}_{2}(u,\rho)$ can now be approximated by the following:

(13)
$\int_{\mathbb{R}^{n}}||x||^{2}du(x)+\int_{\mathbb{R}^{n}}||y||^{2}d\rho(y)-2\cdot\bigg{[}\int_{\mathbb{R}^{n}}\psi_{\theta^{*}}(x)du(x)+\int_{\mathbb{R}^{n}}\overline{\psi_{\omega^{*}}}(y)d\rho(y)\bigg{]}$

The remaining obstacle is to find the optimal $\psi_{\theta^{*}}$ and
$\overline{\psi_{\omega^{*}}}$ with a tractable computation. In order to find
them and to use them as approximated solutions of $\psi^{*}$ and
$\overline{\psi^{*}}$, a large enough but tractable space of parametrized
convex functions should be introduced. Here, we used input convex neural
network (ICNN) architectures to approximate the optimal convex functions,
which were firstly presented in [3]. Recently, ICNN architectures are widely
chosen to approximate the optimal convex potential and its conjugate in the
previous studies [43, 28, 20, 1, 32]. Among them, a fully-convex, $k$-layer,
fully-connected ICNN (FICNN) architecture was proposed

$\displaystyle z_{i+1}=g_{i}(W_{i}^{(z)}z_{i}+W_{i}^{(y)}y+b_{i}),\quad 0\leq
i\leq k-1,\quad W_{0}^{(z)}=0,$

where $y$ denotes an input, $g_{i}$ denotes the activation function, and
$z_{k}$ denotes the final output. Let
$\theta:=\\{W_{0:k-1}^{(y)},W_{1:k-1}^{(z)},b_{0:k-1}\\}$ be the parameters of
$k$-layer FICNN $f(y;\theta)=z_{k}$. The authors in [3] claimed that the
function $f$ is convex with respect to $y$ given $W_{1:k-1}^{(z)}$ are non-
negative and $g_{i}$ are non-decreasing convex activation functions.

###  2.4 Universal approximation theorem

The approximation theory of the neural network has been widely studied in past
decades after a seminal work [10] which states the universal approximation
property of one hidden layer neural network with sigmoidal activations.

The original result had been generalized in many directions. For example, [15]
showed the approximation property for more general activation functions.

######  Theorem 2.1 (Theorem 2 in [15]).

Whenever $\sigma$ is continuous, bounded and non-constant, the the finite sums
of the form

(14) $G(x)=\sum_{j=1}^{N}a_{j}\sigma(y_{j}^{T}x+\theta_{j})$

are dense in $C(X)$ for all compact sets $X\subset\mathbb{R}^{n}$.

Later, the authors in [23] showed that being a universal approximator is
equivalent to having a non-polynomial activation function $\sigma$.
Approximation theorems in the space of differentiable functions are
established in [24] through the following theorem.

######  Theorem 2.2 (Theorem 2.1 in [24]).

Let X be a compact subset of $\mathbb{R}^{n}$ and
$f\in\hat{C}^{m_{1}}(X)\cap\hat{C}^{m_{2}}(X)\cap\cdots\cap\hat{C}^{m_{q}}(X)$,
where $m_{i}\in\mathbb{Z}^{n}_{+}$. Also, let $\sigma$ be any non-polynomial
function in $C^{n}(\mathbb{R})$, where $n=\max\\{|m_{i}|,1\leq i\leq q\\}$.
Then for any $\epsilon>0$, there is a neural network defined as in (14) such
that

$\|D^{k}f-D^{k}G\|_{L^{\infty}(X)}<\epsilon,\quad k\leq m_{i},\quad\text{ for
some }i\quad,1\leq i\leq q.$

The density argument on simultaneous approximations of multivariate functions
and their partial derivatives have been applied to show the existence of a
neural network that arbitrarily minimizes the PDE residual loss function, see
[17, 41].

##  3 The Deep Minimizing Movement Method

We consider the following sequential minimization problem with time step
$\tau$ :

(15) $u_{k+1}=\operatorname*{arg\,min}_{u\in
X}\\{\mathcal{F}(u)+\frac{d_{X}^{2}(u,u_{k})}{2\tau}\\},\quad u_{0}=u(0,x),$

where $X$ denotes either $L^{2}(\Omega)$, or $\mathcal{P}_{2}(\Omega)$ with
$L^{2}$-Wasserstein distance.

The main idea of this paper is to optimize (15) by parametrizing $u$ by the
neural network $u_{\theta}$. We present some theoretical results to support
our method in the following sequel.

First of all, we show that there exists a neural network $u_{\theta}$ which
minimizes the loss
$\mathcal{F}(u_{\theta})+\frac{d_{L^{2}(\Omega)}^{2}(u_{\theta},u_{k})}{2\tau}$
in (14) as much as we want in most cases.

######  Theorem 3.1.

Let $\Omega$ be a compact subset of $\mathbb{R}^{n}$ and
$\phi:\Omega\times\mathbb{R}\times\mathbb{R}^{n}\rightarrow\mathbb{R}$ be a
continuous function and $\tau>0$. Let $g,h$ be $C^{k}$-functions on $\Omega$.
Define $\mathcal{I}(u):=\int_{\Omega}\phi(x,u(x),\nabla
u(x))dx+\frac{d_{L^{2}(\Omega)}^{2}(u,h)}{2\tau}$ for every $C^{k}$-function
$u$ on $\Omega$. Then, for $\epsilon>0$, there is a neural network
$u_{\theta}$ defined as in (14) such that
$||u_{\theta}-g||_{L^{\infty}(\Omega)}\leq\epsilon$ and
$\mathcal{I}(u_{\theta})-\mathcal{I}(g)\leq\epsilon$.

######  Proof 3.2.

For notational convenience, let us define
$\mathcal{J}(u):=\int_{\Omega}\phi(x,u(x),\nabla u(x))dx$ for every
$C^{k}$-function $u$ on $\Omega$, so that
$\mathcal{I}(u)=\mathcal{J}(u)+\frac{d_{L^{2}(\Omega)}^{2}(u,h)}{2\tau}$ for
all $u$.

Let $\epsilon>0$ be given and pick $\epsilon_{1}>0$ such that
$(1+\frac{\epsilon_{1}}{2\tau})\cdot
vol(\Omega)\cdot\epsilon_{1}+\frac{1}{\tau}\cdot\sqrt{vol(\Omega)}\cdot\epsilon_{1}\cdot
d_{L^{2}(\Omega)}(g,h)<\epsilon$.

Since $g$ is $C^{k}$ on $\Omega$, both $g$ and $\nabla g$ are bounded on
$\Omega$. Thus, there exists $R>0$ such that the range of $g$ is contained in
the interval $(-R,R)$ and the range of $\nabla g$ is contained in the open
ball $B(0,R)$. Since $\phi$ is continuous, it is uniformly continuous on
$\Omega\times[-R-\epsilon_{1},R+\epsilon_{1}]\times\overline{B(0,R+\epsilon_{1})}$.
Thus, there exists $\delta$ such that $|x-y|+||z-w||_{2}<\delta$ implies
$|\phi(a,x,z)-\phi(a,y,w)|<\epsilon_{1}$ for any
$(a,x,z),(a,y,w)\in\Omega\times[-R-\epsilon_{1},R+\epsilon_{1}]\times\overline{B(0,R+\epsilon_{1})}$.
Then, by Theorem 2.2, there exists a neural network $u_{\theta}$ satisfying
$||u_{\theta}-g||_{L^{\infty}(\Omega)}+||\nabla u_{\theta}-\nabla
g||_{L^{\infty}(\Omega)}<\min\\{\delta,\epsilon_{1},\epsilon\\}$, so that
$||u_{\theta}-g||_{L^{\infty}(\Omega)}\leq\epsilon$.

Next, we pick $x\in\Omega$. Since $||\nabla u_{\theta}(x)||_{2}-||\nabla
g(x)||_{2}\leq||\nabla u_{\theta}(x)-\nabla g(x)||_{2}\leq||\nabla
u_{\theta}-\nabla g||_{L^{\infty}(\Omega)}<\epsilon_{1}$, we know that $\nabla
u_{\theta}(x)\in\overline{B(0,R+\epsilon_{1})}$. Analogously,
$u_{\theta}(x)\in[-R-\epsilon_{1},R+\epsilon_{1}]$. Therefore, both
$(u_{\theta}(x),\nabla u_{\theta}(x))$ and $(g(x),\nabla g(x))$ are elements
of $[-R-\epsilon_{1},R+\epsilon_{1}]\times\overline{B(0,R+\epsilon_{1})}$.
Since $|u_{\theta}(x)-g(x)|+||\nabla u_{\theta}(x)-\nabla g(x)||_{2}<\delta$,
we conclude that $|\phi(x,u_{\theta}(x),\nabla
u_{\theta}(x))-\phi(x,g(x),\nabla g(x))|<\epsilon_{1}$. Since this holds for
every $x$, we have

$\displaystyle|\mathcal{J}(u_{\theta})-\mathcal{J}(g)|$
$\displaystyle=\bigg{|}\int_{\Omega}\phi(x,u_{\theta}(x),\nabla
u_{\theta}(x))-\phi(x,g(x),\nabla g(x))dx\bigg{|}$
$\displaystyle\leq\int_{\Omega}|\phi(x,u_{\theta}(x),\nabla
u_{\theta}(x))-\phi(x,g(x),\nabla g(x))|dx$
$\displaystyle\leq\epsilon_{1}\cdot vol(\Omega)$

Now, note that
$d_{L^{2}(\Omega)}^{2}(u_{\theta},g)=\int_{\Omega}|u_{\theta}-g|^{2}dx\leq
vol(\Omega)\cdot||u_{\theta}-g||_{L^{\infty}(\Omega)}^{2}\leq
vol(\Omega)\cdot\epsilon_{1}^{2}$. Thus, we have

$\displaystyle d_{L^{2}(\Omega)}^{2}(u_{\theta},h)-d_{L^{2}(\Omega)}^{2}(g,h)$
$\displaystyle\leq\bigg{(}d_{L^{2}(\Omega)}(u_{\theta},g)+d_{L^{2}(\Omega)}(g,h)\bigg{)}^{2}-d_{L^{2}(\Omega)}^{2}(g,h)$
$\displaystyle=d_{L^{2}(\Omega)}^{2}(u_{\theta},g)+2\cdot
d_{L^{2}(\Omega)}(u_{\theta},g)\cdot d_{L^{2}(\Omega)}(g,h)$
$\displaystyle\leq
vol(\Omega)\cdot\epsilon_{1}^{2}+2\cdot\sqrt{vol(\Omega)}\cdot\epsilon_{1}\cdot
d_{L^{2}(\Omega)}(g,h)$

Hence,

$\displaystyle\mathcal{I}(u_{\theta})-\mathcal{I}(g)$
$\displaystyle\leq(1+\frac{\epsilon_{1}}{2\tau})\cdot
vol(\Omega)\cdot\epsilon_{1}+\frac{1}{\tau}\cdot\sqrt{vol(\Omega)}\cdot\epsilon_{1}\cdot
d_{L^{2}(\Omega)}(g,h)$ $\displaystyle\leq\epsilon.$

Theorem 3.1 shows that if we take $g$ as the global minimizer $u_{k+1}$ of
(14) with $h=u_{k}$, there always exist a neural network $u_{\theta}$ such
that $\mathcal{I}(u_{\theta})\leq\mathcal{I}(u_{k+1})+\epsilon$ with
$||u_{\theta}-u_{k+1}||_{L^{\infty}(\Omega)}\leq\epsilon$. Therefore, we can
expect that we can find a neural network solution approximating $u_{k+1}$ by
optimizing the parametric form of (14).

Next, we will justify our method for $\mathbb{W}_{2}$-gradient flows. Yet, we
need a variation of the standard universal approximation theorem for
probability density functions before we prove the analogous result of Theorem
3.1 for $\mathbb{W}_{2}$ spaces. This is because we will approximate a given
probability density function by neural networks which behave as if probability
density functions. We start by proving that any non-negative $C^{k}$-function
can be approximated by a non-negative neural network.

######  Theorem 3.3.

Let $X,f,\sigma$ and $m_{i}$’s be given as in Theorem 2.2 with an additional
assumption that $f$ is non-negative. Then, for any $\epsilon>0$, there is a
non-negative neural network defined as in $\eqref{nn_dense}$ such that

$\|D^{k}f-D^{k}G\|_{L^{\infty}(X)}<\epsilon,\quad k\leq m_{i},\quad\text{ for
some }i,1\leq i\leq q.$

######  Proof 3.4.

Fix $\epsilon>0$ and define $\tilde{f}:=f+\epsilon$. Then, by Theorem 2.2,
there exists a neural network $\tilde{G}$ of the form
$\sum_{j=1}^{N}a_{j}\sigma(y_{j}^{T}x+\theta_{j})$ such that

$\|D^{k}\tilde{f}-D^{k}\tilde{G}\|_{L^{\infty}(X)}<\epsilon/2,\quad k\leq
m_{i},\quad\text{ for some }i,1\leq i\leq q.$

Since
$\tilde{f}(x)-\tilde{G}(x)\leq||\tilde{f}-\tilde{G}||_{L^{\infty}(\Omega)}<\epsilon/2$
for any $x\in X$, we have
$\tilde{G}(x)>\tilde{f}-\epsilon/2=f(x)+\epsilon/2\geq\epsilon/2$ for any
$x\in X$.

Since $\sigma$ is non-polynomial, there exists $\theta_{N+1}\in\mathbb{R}$
such that $\sigma(\theta_{N+1})\neq 0$. Now, we define $y_{N+1}:=0$ and
$a_{N+1}:=\frac{-\epsilon}{2\cdot\sigma(\theta_{N+1})}$ and
$G(x):=\sum_{j=1}^{N+1}a_{j}\sigma(y_{j}^{T}x+\theta_{j})$ for every $x\in X$,
so that $G$ is a neural network defined as in (14).

Now, we pick $x\in X$. Then,

$\displaystyle G(x)$
$\displaystyle=\tilde{G}(x)+a_{N+1}\sigma(y_{N+1}^{T}x+\theta_{N+1})$
$\displaystyle=\tilde{G}(x)+a_{N+1}\sigma(\theta_{N+1})$ (16)
$\displaystyle=\tilde{G}(x)-\epsilon/2\geq 0.$

Therefore, we have

$\displaystyle f(x)-G(x)$ $\displaystyle=f(x)-\tilde{G}+\epsilon/2$
$\displaystyle=f(x)+\epsilon-\tilde{G}(x)-\epsilon/2$
$\displaystyle=\tilde{f}(x)-\tilde{G}(x)-\epsilon/2.$

Thus, $|f(x)-G(x)|\leq|\tilde{f}(x)-\tilde{G}(x)|+\epsilon/2<\epsilon$ and
$G(x)\geq 0$ for every $x\in X$ by (3.4), meaning that
$||f-G||_{L^{\infty}(\Omega)}<\epsilon$ and $G$ is a non-negative neural
network as in (14). Since $\tilde{f}=f+\epsilon$ and $G=\tilde{G}-\epsilon/2$,
$D^{k}f=D^{k}\tilde{f}$ and $D^{k}G=D^{k}\tilde{G}$ for every $k\neq 0$.
Consequently, $G$ is the desired non-negative neural network.

Now, we prove that there exists a probability density function
$\widetilde{u_{\theta}}$ which can be produced by manipulating a neural
network, which approximates a given density $f$ as much as we want.

######  Theorem 3.5.

Let $X,f,\sigma$ and $m_{i}$’s be given as in Theorem 2.2 with an additional
assumption that $f$ is a probability density function whose domain is $X$.
Then, for any $\epsilon>0$, there is a non-negative neural network $G$ defined
as in (14) such that

$\|f-\tilde{G}\|_{L^{\infty}(X)}<\epsilon,$

where $\tilde{G}=\frac{1}{\int_{X}G}G$.

######  Proof 3.6.

Fix $\epsilon>0$. Then, take $\epsilon_{0}>0$ and $0<\epsilon_{1}<1$ such that
$\epsilon_{0}\cdot||f||_{L^{\infty}(\Omega)}+\frac{\epsilon_{0}^{2}}{vol(X)}<\epsilon/2$
and
$\max\\{\frac{\epsilon_{1}}{1-\epsilon_{1}},\frac{\epsilon_{1}}{1+\epsilon_{1}},\epsilon_{1}\\}<\epsilon_{0}$.

Then, by Theorem 3.3, there exists a non-negative neural network $G$ defined
as in (14) such that
$||f-G||_{L^{\infty}(\Omega)}<\min\\{\frac{\epsilon_{1}}{vol(X)},\epsilon/2\\}$.
Thus,
$||G||_{L^{\infty}(\Omega)}\leq||f||_{L^{\infty}(\Omega)}+\frac{\epsilon_{1}}{vol(X)}<||f||_{L^{\infty}(\Omega)}+\frac{\epsilon_{0}}{vol(X)}$.
Since

$\displaystyle\bigg{|}1-\int_{X}G\bigg{|}$
$\displaystyle=\bigg{|}\int_{X}f-\int_{X}G\bigg{|}=\bigg{|}\int_{X}(f-G)\bigg{|}\leq\int_{X}|f-G|$
$\displaystyle\leq
vol(X)\cdot||f-G||_{L^{\infty}(\Omega)}<vol(X)\cdot\frac{\epsilon_{1}}{vol(X)}=\epsilon_{1},$

we have $1-\epsilon_{1}<\int_{X}G<1+\epsilon_{1}$. Therefore,
$-\epsilon_{0}<-\frac{\epsilon_{1}}{1-\epsilon_{1}}<1-\frac{1}{\int_{X}G}<\frac{\epsilon_{1}}{1+\epsilon_{1}}<\epsilon_{0}$,
leading to $\bigg{|}1-\frac{1}{\int_{X}G}\bigg{|}\leq\epsilon_{0}$.

Hence,

$\displaystyle\bigg{|}\bigg{|}G-\frac{G}{\int_{X}G}\bigg{|}\bigg{|}_{L^{\infty}(\Omega)}$
$\displaystyle\leq||G||_{L^{\infty}(\Omega)}\cdot\bigg{|}1-\frac{1}{\int_{X}G}\bigg{|}$
$\displaystyle\leq||G||_{L^{\infty}(\Omega)}\cdot\epsilon_{0}$
$\displaystyle\leq(||f||_{L^{\infty}(\Omega)}+\frac{\epsilon_{0}}{vol(X)})\cdot\epsilon_{0}$
$\displaystyle=\epsilon_{0}\cdot||f||_{L^{\infty}(\Omega)}+\frac{\epsilon_{0}^{2}}{vol(X)}$
$\displaystyle<\epsilon/2.$

Therefore, we have

$\displaystyle\bigg{|}\bigg{|}f-\frac{G}{\int_{X}G}\bigg{|}\bigg{|}_{L^{\infty}(\Omega)}\leq||f-G||_{L^{\infty}(\Omega)}+\bigg{|}\bigg{|}G-\frac{G}{\int_{X}G}\bigg{|}\bigg{|}_{L^{\infty}(\Omega)}<\epsilon/2+\epsilon/2=\epsilon.$

Throughout, we identify absolutely continuous measures with their densities,
$d\rho(x)=\rho(x)dx$. We will now show that there exists a non-negative neural
network $u_{\theta}$ which minimizes the loss
$\mathcal{F}(\widetilde{u_{\theta}})+\frac{\mathbb{W}_{2}^{2}(\widetilde{u_{\theta}},u_{k})}{2\tau}$
in (15) as much as we want in most cases, where
$\widetilde{u_{\theta}}:=\frac{1}{\int_{\Omega}u_{\theta}}u_{\theta}$.

######  Lemma 3.7.

Assume that $\Omega$ is a compact subset of $\mathbb{R}^{n}$, and let $u,\rho$
be absolutely continuous probability distributions on $\Omega$ with finite
second moments such that $||u-\rho||_{L^{\infty}(\Omega)}\leq\epsilon$. Then,
$\mathbb{W}_{2}(u,\rho)\leq C\cdot\epsilon^{1/2}$ for some constant $C$ which
only depends on $\Omega$.

######  Proof 3.8.

Let $\gamma$ be a joint distribution of $u$ and $\rho$. Then, we have

$\displaystyle\mathbb{W}_{2}(u,\rho)$
$\displaystyle=\bigg{(}\int_{\Omega}||x-y||_{2}^{2}d\gamma(x,y)\bigg{)}^{1/2}$
$\displaystyle=\bigg{(}\int_{\Omega}||x-y||_{2}\cdot||x-y||_{2}d\gamma(x,y)\bigg{)}^{1/2}$
$\displaystyle\leq
diam(\Omega)^{\frac{1}{2}}\cdot\bigg{(}\int_{\Omega}||x-y||_{2}d\gamma(x,y)\bigg{)}^{1/2}.$

Since $\gamma$ was arbitrary, we obtain

(17) $\mathbb{W}_{2}(u,\rho)\leq
diam(\Omega)^{\frac{1}{2}}\cdot\bigg{(}\inf_{\gamma\in\Pi(u,\rho)}\int_{\Omega}||x-y||_{2}d\gamma(x,y)\bigg{)}^{\frac{1}{2}},$

where $\Pi(u,\rho)$ is the collection of joint distributions of $u$ and
$\rho$.

The term $\inf_{\gamma\in\Pi(u,\rho)}\int_{\Omega}||x-y||_{2}d\gamma(x,y)$ is
called the Wasserstein-1 metric and is denoted by $\mathbb{W}_{1}(u,\rho)$.
Now, by the Kantorovich-Rubinstein duality theorem, we have

(18) $\mathbb{W}_{1}(u,\rho)=\sup\biggl{\\{}\int_{\Omega}f(x)d(u-\rho)(x):f\in
C(\Omega)\text{ with }||f||_{Lip}\,\leq 1\biggr{\\}},$

where $||\cdot||_{Lip}$ denotes the minimal Lipschitz constant.

Define the set $\mathcal{E}:=\biggl{\\{}\int_{\Omega}g(x)d(u-\rho)(x):g\in
C(\Omega)\text{ with }||g||_{L^{\infty}(\Omega)}\leq 1\biggr{\\}}$. Let $f\in
C(\Omega)$ with $||f||_{Lipz}\leq 1$, and pick a point $x_{0}\in\Omega$. Now,
define $g:=\frac{1}{diam(\Omega)}\cdot(f-f(x_{0}))$. Since
$|f(x)-f(x_{0})|\leq||x-x_{0}||_{2}\leq diam(\Omega)$ for every $x$, we have
$||g||_{L^{\infty}(\Omega)}\leq 1$. Since

$\displaystyle\int_{\Omega}g(x)d(u-\rho)(x)$
$\displaystyle=\frac{1}{diam(\Omega)}\int_{\Omega}(f(x)-f(x_{0}))d(u-\rho)(x)$
$\displaystyle=\frac{1}{diam(\Omega)}\int_{\Omega}fd(u-\rho)(x),$

we have
$\frac{1}{diam(\Omega)}\int_{\Omega}f(x)d(u-\rho)(x)=\int_{\Omega}g(x)d(u-\rho)(x)\in\mathcal{E}$.
Since $f$ was arbitrary, by (18), we have

(19) $\mathbb{W}_{1}(u,\rho)\leq diam(\Omega)\cdot\sup\mathcal{E}.$

Combining with (17), we get

(20) $\mathbb{W}_{2}(u,\rho)\leq
diam(\Omega)\cdot(\sup\mathcal{E})^{\frac{1}{2}},$

Next, consider $g\in C(\Omega)$ with $||g||_{L^{\infty}(\Omega)}\leq 1$. Then,

$\displaystyle\int_{\Omega}g(x)d(u-\rho)(x)$
$\displaystyle=\int_{\Omega}g(x)(u(x)-\rho(x))dx$
$\displaystyle\leq\int_{\Omega}|g(x)|\cdot|u(x)-\rho(x)|dx$
$\displaystyle\leq\int_{\Omega}|u(x)-\rho(x)|dx$ $\displaystyle\leq
vol(\Omega)\cdot||u-\rho||_{L^{\infty}(\Omega)}$ $\displaystyle\leq
vol(\Omega)\cdot\epsilon.$

Since this holds for all such $g$, we have $\sup\mathcal{E}\leq
vol(\Omega)\cdot\epsilon$.

Consequently, we have $\mathbb{W}_{2}^{2}(u,\rho)\leq C\cdot\epsilon$ where
$C$ only depends on $\Omega$.

######  Theorem 3.9.

Let $\Omega$ be a compact subset of $\mathbb{R}^{n}$ and
$\phi:\Omega\times[0,1]\rightarrow\mathbb{R}$ be a continuous function and
$\tau>0$. Let $\rho,\pi$ be absolutely continous $C^{k}$ probability
distributions on $\Omega$ with finite second moments. Define
$\mathcal{I}(u):=\int_{\Omega}\phi(x,u(x))dx+\frac{\mathbb{W}_{2}^{2}(u,\pi)}{2\tau}$
for every $C^{k}$ probability distribution $u$ on $\Omega$ with finite second
moment. Then, for $\epsilon>0$, there is a non-negative neural network
$u_{\theta}$ defined as in (14) such that
$||\widetilde{u_{\theta}}-\rho||_{L^{\infty}(\Omega)}\leq\epsilon$ and
$\mathcal{I}(\widetilde{u_{\theta}})-\mathcal{I}(\rho)\leq\epsilon$, where
$\widetilde{u_{\theta}}=\frac{1}{\int_{\Omega}u_{\theta}}u_{\theta}$.

######  Proof 3.10.

Let $\epsilon>0$ be given, and $C$ be the constant given in Lemma 3.7 which
only depends on $\Omega$. Let $\epsilon_{0}>0$ be a constant such that
$\epsilon_{0}+\frac{1}{2\tau}(\epsilon_{0}+2\cdot\sqrt{\epsilon_{0}}\cdot\mathbb{W}_{2}(\rho,\pi))<\epsilon$.
Since $\Omega\times[0,1]$ is compact, $\phi$ is uniformly continuous.
Therefore, there exists $\delta$ such that $|x-y|\leq\delta$ implies
$|\phi(z,x)-\phi(z,y)|\leq\frac{1}{vol(\Omega)}\cdot\epsilon_{0}$, for every
$z\in\Omega$ and $x,y\in[0,1]$. Now, by Theorem 3.5, there is a non-negative
neural network $u_{\theta}$ satisfying
$||\widetilde{u_{\theta}}-\rho||_{L^{\infty}(\Omega)}<\min\\{\delta,\frac{1}{C^{2}}\cdot\epsilon_{0},\epsilon\\}$.
Therefore, we obtain

$\displaystyle\bigg{|}\int_{\Omega}\phi(x,\widetilde{u_{\theta}}(x))dx-\int_{\Omega}\phi(x,\rho(x))dx\bigg{|}$
$\displaystyle\leq\int_{\Omega}|\phi(x,\widetilde{u_{\theta}}(x))-\phi(x,\rho(x))|dx$
$\displaystyle\leq vol(\Omega)\cdot\frac{1}{vol(\Omega)}\cdot\epsilon_{0}$
$\displaystyle=\epsilon_{0}$

Since
$||\widetilde{u_{\theta}}-\rho||_{L^{\infty}(\Omega)}<\frac{1}{C^{2}}\cdot\epsilon_{0}$,
by Lemma 3.7, we have
$\mathbb{W}_{2}(\widetilde{u_{\theta}},\rho)\leq\sqrt{\epsilon_{0}}$.
Therefore,

$\displaystyle\mathbb{W}_{2}^{2}(\widetilde{u_{\theta}},\pi)-\mathbb{W}_{2}^{2}(\rho,\pi)$
$\displaystyle\leq\bigg{(}\mathbb{W}_{2}(\widetilde{u_{\theta}},\rho)+\mathbb{W}_{2}(\rho,\pi)\bigg{)}^{2}-\mathbb{W}_{2}^{2}(\rho,\pi)$
$\displaystyle=\mathbb{W}_{2}^{2}(\widetilde{u_{\theta}},\rho)+2\cdot\mathbb{W}_{2}(\widetilde{u_{\theta}},\rho)\cdot\mathbb{W}_{2}(\rho,\pi)$
$\displaystyle\leq\epsilon_{0}+2\cdot\sqrt{\epsilon_{0}}\cdot\mathbb{W}_{2}(\rho,\pi)$
$\displaystyle\leq 2\tau\cdot(\epsilon-\epsilon_{0})$

Consequently, we have

$\displaystyle\mathcal{I}(\widetilde{u_{\theta}})-\mathcal{I}(\rho)$
$\displaystyle=\int_{\Omega}\phi(\widetilde{u_{\theta}}(x))-\phi(\rho(x))dx+\frac{1}{2\tau}\cdot(\mathbb{W}_{2}^{2}(\widetilde{u_{\theta}},\pi)-\mathbb{W}_{2}^{2}(\rho,\pi))$
$\displaystyle\leq\epsilon_{0}+(\epsilon-\epsilon_{0})$
$\displaystyle=\epsilon$

Theorem 3.9 shows that if we take $\rho$ as the global minimizer $u_{k+1}$ of
(14) with $\pi=u_{k}$, there always exist a neural network $u_{\theta}$ such
that $\mathcal{I}(\widetilde{u_{\theta}})\leq\mathcal{I}(u_{k+1})+\epsilon$
with $||\widetilde{u_{\theta}}-u_{k+1}||_{L^{\infty}(\Omega)}\leq\epsilon$,
where $\widetilde{u_{\theta}}=\frac{1}{\int u_{\theta}}u_{\theta}$. Therefore,
we can expect that we can find a neural network solution approximating
$u_{k+1}$ by optimizing the parametric form of (14).

Lastly, we present a result which shows that the increment of the error at
each update is bounded above.

######  Theorem 3.11.

Let $X$ be a metric space and
$\mathcal{F}:X\rightarrow\mathbb{R}\cup\\{+\infty\\}$ be a lower semi-
continuous function with some lower bounds to guarantee the existence of
minimizers of (15) for small $\tau$. Given $u_{0}\in X$ with
$\mathcal{F}(u_{0})<\infty$, let $\\{u_{k}\\}$ be the true minimizers of (15)
for each $k$, and let $\\{u_{\theta,k}\\}$ be a sequence of neural networks
obtained by solving the optimization problem (15) iteratively with a natrual
assumption that
$\mathcal{F}(u_{\theta,k+1})+\frac{d_{X}^{2}(u_{\theta,k+1},u_{\theta,k})}{2\tau}\leq\mathcal{F}(u_{\theta,k})+\frac{d_{X}^{2}(u_{\theta,k},u_{\theta,k})}{2\tau}$
for each $k$. Then, we obtain

$d_{X}(u_{\theta,k+1},u_{k+1})\leq d_{X}(u_{\theta,k},u_{k})+C\tau^{1/2}$

for some constant $C$ depending on $\mathcal{F}$ and initial $u_{0}$.

######  Proof 3.12.

Firstly, since $u_{k+1}$ is a minimizer of (15), we get

$\mathcal{F}(u_{k+1})+\frac{d_{X}^{2}(u_{k+1},u_{k})}{2\tau}\leq\mathcal{F}(u_{k})+\frac{d_{X}^{2}(u_{k},u_{k})}{2\tau}=\mathcal{F}(u_{k})$

Therefore, we get $\mathcal{F}(u_{k+1})\leq\mathcal{F}(u_{k})$ for each $k$.
Since this is valid for every $k$, we have
$\mathcal{F}(u_{k})\leq\mathcal{F}(u_{0})<+\infty$ for each $k$. Analogous
arguments apply to $u_{\theta,k}$, and we conclude that
$\mathcal{F}(u_{\theta,k})<+\infty$ for each $k$.

Again, by the definition of $u_{k}$’s, we get

$\displaystyle d_{X}^{2}(u_{k},u_{k+1})\leq
2\tau(\mathcal{F}(u_{k})-\mathcal{F}(u_{k+1}))\leq
2\tau(\mathcal{F}(u_{0})-\inf\mathcal{F})=2A\tau,$

where $A=A(u_{0},\mathcal{F})\coloneqq(\mathcal{F}(u_{0})-\inf\mathcal{F})$
only depends on the functional $\mathcal{F}$ and the initial $u_{0}$.
Analogously, by the assumption, we get
$d_{X}^{2}(u_{\theta,k},u_{\theta,k+1})<2A\tau$. Then,

$\displaystyle d_{X}(u_{\theta,k+1},u_{k+1})$ $\displaystyle\leq
d_{X}(u_{\theta,k+1},u_{\theta,k})+d_{X}(u_{\theta,k},u_{k})+d_{X}(u_{k},u_{k+1})$
$\displaystyle<\sqrt{2A\tau}+d_{X}(u_{\theta,k},u_{k})+\sqrt{2A\tau}=d_{X}(u_{\theta,k},u_{k})+C\tau^{1/2},$

where $C=2\sqrt{2A}$.

In the remaining of this section, we give a detailed description of the
framework of our method for the spaces $L^{2}(\Omega)$ and
$\mathcal{P}_{2}(\Omega)$ with $L^{2}$-Wasserstein distance.

###  3.1 Case when $\boldsymbol{X=L^{2}(\Omega)}$

If $X=L^{2}(\Omega)$ and the functional $\mathcal{F}(u)$ is given in the form
$\int\phi(u,\nabla u)$, this integration can be evaluated by Monte Carlo
estimate over mini-batches from the uniform distribution on $\Omega$ if
$\Omega$ is bounded. If $\Omega=\mathbb{R}^{n}$, we can still approximate this
integration by restricting the domain to a sufficiently large cropped bounded
subset, i.e. $\int_{\mathbb{R}^{n}}\phi(u,\nabla
u)\approx\int_{[-M,M]^{n}}\phi(u,\nabla u)$. Likewise, the distance
$d_{L^{2}(\Omega)}^{2}(u,v)=\int_{\Omega}||u-v||^{2}$ can also be evaluated by
Monte Carlo estimate in the same manner. Now, by Theorem 3.1, it would suffice
to minimize

(21) $\operatorname*{arg\,min}_{\theta}\\{\int\phi(u_{\theta},\nabla
u_{\theta})+\frac{d_{L^{2}(\Omega)}^{2}(u_{\theta},u_{k})}{2\tau}\\},$

for a given $u_{k}$ to obtain the minimizer $u_{\theta^{*(k+1)}}$, which is an
approximation of the minimizer of (15). Then, we take the optimized
$u_{\theta^{*(k+1)}}$ in place of $u_{k}$ and repeat this process iteratively
to derive the corresponding time-discretized gradient flow.

###  3.2 Case when $\boldsymbol{X=\mathcal{P}_{2}(\Omega)}$

If $X=\mathcal{P}_{2}(\Omega)$ and the functional $\mathcal{F}(u)$ is given as
a linear combination of potential energy $\mathcal{V}$ and interaction energy
$\mathcal{W}$ and internal energy $\mathcal{E}$, there is no problem
approximating $\mathcal{F}(u)$ by Monte Carlo estimate as described in Section
3.1. Moreover, by Theorem 3.5, any probability density function can be
approximated by manipulating a neural network $u_{\theta}$ having non-negative
outputs. Therefore, if we define $\widetilde{u_{\theta}}:=\frac{1}{\int
u_{\theta}}u_{\theta}$ for non-negative neural networks $u_{\theta}$ so that
$\int\widetilde{u_{\theta}}=1$, it suffices to minimize

(22)
$\operatorname*{arg\,min}_{\theta}\\{\mathcal{F}(\widetilde{u_{\theta}})+\frac{\mathbb{W}_{2}^{2}(\widetilde{u_{\theta}},u_{k})}{2\tau}\\},$

for a given $u_{k}$ to obtain the minimizer $\widetilde{u_{\theta^{*(k+1)}}}$,
which is an approximation of the minimizer of (15). As in the case
$X=L^{2}(\Omega)$, we take the optimized $\widetilde{u_{\theta^{*(k+1)}}}$ in
place of $u_{k}$ and repeat this process iteratively to derive the
corresponding time-discretized gradient flow.

However, unlike the case $X=L^{2}(\Omega)$, there is no direct easy way to
evaluate the squared metric, i.e. $\mathbb{W}_{2}^{2}(u,\rho)$, when
minimizing the objective (22). There are many algorithms to approximate the
$L^{2}$-Wasserstein distance, but we take an idea proposed in [20] for the
reasons described in Section 2.3. In other words, we approximate it by solving
the optimization problem (12) at first and then by evaluating (13) using the
optimal $\psi_{\theta^{*}}$ and $\overline{\psi_{w^{*}}}$. All the
integrations in formulae (13) and (12) are again evaluated by Monte Carlo as
described in Section 3.1. Moreover, the non-negativity assumption on neural
networks can be easily satisfied by choosing non-negative activation functions
such as Softplus [14], ReLU [36], ELU [8], etc.

##  4 Experiments

In this section, we show through various kinds of examples, that one can
accurately approximate the solutions of PDEs which have gradient flow
structures via neural networks. We divide this section into two subsections,
each of which is devoted to gradients flows for the case $L^{2}$ and the case
$\mathbb{W}_{2}$ respectively. For the $L^{2}$ gradient flow case, we
demonstrate the scalability of the proposed method to the high-dimensional
problems through the heat equation. Throughout this section, we apply the
domain truncation technique for the unbounded domains (in our cases
$\mathbb{R}^{n}$) as applied in the literature (See, for example, [16, 7]). We
compute the $L^{2}$ error of a neural network solution $u_{NN}$ by
$\|u-u_{NN}\|_{L^{2}}$ and the relative $L^{2}$ error by
$\frac{\|u-u_{NN}\|_{L^{2}}}{\|u\|_{L^{2}}}$, where $u$ denotes either an
analytic solution or a numerical solution.

##### Motivating Example

We first consider a simple motivating example with the $L^{2}$ energy
functional $\frac{\kappa}{2}\int_{\Omega}u^{2}dx$. The corresponding
minimization problem reads as :

(23)
$\operatorname*{arg\,min}_{\theta}\\{\frac{\kappa}{2}\int_{\Omega}u_{\theta}^{2}dx+\frac{d^{2}_{L^{2}(\Omega)}(u_{\theta},u_{k})}{2\tau}\\},$

where $u_{\theta}$ is a neural network solution and $\theta$ denotes the set
of parameters. We set an initial condition $u(0,x)=\sin(x)$,
$\Omega=[-\pi,\pi]$, then the corresponding solution of the PDE is
$u(t,x)=e^{-\kappa t}\sin(x)$. In Figure 1, we demonstrate the minimizer
$u_{\theta,k+1}$ approximately satisfies the Euler-Lagrange equation of (23),
i.e.,

$\frac{u_{\theta,k+1}-u_{\theta,k}}{\tau}=-\nabla F(u_{k+1}),$

during $40$ timesteps with $\tau=0.01$.

Figure 1: Left: Scatter plot for the values of
$\frac{u_{\theta,k+1}-u_{\theta,k}}{\tau}$ on the uniform samples. Right:
Scatter plot for the values of $-\nabla F(u_{\theta,k+1})$ on the same sample
points.

###  4.1 $\boldsymbol{L^{2}}$-gradient flows

In this section, we consider two examples of $L^{2}$-gradient flows, the heat
equation and the Allen-Cahn equation.

####  4.1.1 Heat equation: 2-dimensional, 8-dimensional

The heat equation

(24) $\begin{cases}\partial_{t}u=\kappa\Delta u\quad\text{in }\Omega,\\\
\nabla u\cdot\hat{\boldsymbol{n}}=0\quad\text{on }\partial\Omega,\end{cases}$

where $\boldsymbol{n}$ denotes the unit normal vector, is the $L^{2}$-gradient
flow of the Dirichlet energy
$\mathcal{F}(u):=\frac{\kappa}{2}\int_{\Omega}||\nabla u||^{2}dx$.

For the special case $\Omega=[-a,a]^{n}$, where the domain is given by an
$n$-dimensional cube, the ground-truth solution of (24) can be computed. Let
the initial condition be given by
$u(0,x)=\sum_{[i]}A_{[i]}\prod_{k}\cos(\frac{i_{k}\pi x_{k}}{a})$, where
$[i]=(i_{1},i_{2},...,i_{n})$ denotes the multi-index and $A_{[i]}$ denotes
the coefficient of $\prod_{k}\cos(\frac{i_{k}\pi x_{k}}{a})$. Assuming the
convergence of the initial condition $u(0,x)$, then

(25) $\displaystyle u(t,x)$
$\displaystyle=\sum_{[i]}A_{[i]}\biggr{(}\prod_{k}\cos(\frac{i_{k}\pi
x_{k}}{a})\biggr{)}e^{-\kappa\mu_{[i]}^{2}t},\quad\mu_{[i]}^{2}=\pi^{2}\sum_{k}(\frac{i_{k}}{a})^{2}$

solves (24) and satisfies the boundary condition.

##### 2-dimensional case

In the experiment for 2-dimensional heat equation, we set $\kappa=0.1$ and
$a=\pi/4$. The domain $\Omega$ is given by $[-\pi/4,\pi/4]^{2}$. For an
initial condition
$u_{0}(x)=4-0.5\cos{4x_{1}}+0.6\cos{4x_{2}}+\cos{8x_{1}}\cos{4x_{2}}+\cos{8x_{1}}\cos{8x_{2}}$,
its ground-truth solution $u(t,x)$ is computed as in (25). A size of the time-
step for the minimizing movement scheme is given by $\tau=0.005$.

##### Model training details for 2d case

We use a 2-layer fully connected neural network with 256 hidden units for the
model that approximates the solution $u(x)$. A tangent hyperbolic (tanh)
activation function is applied. We choose the learning rate as $10^{-3}$ with
an Adam optimizer [19]. While training the neural network model, we uniformly
picked 10,000 samples from $\Omega=[-\pi/4,\pi/4]^{2}$ and computed the
minimization target (21) by the Monte-Carlo approximation.

##### Results - 2d case

The approximated results are illustrated in Figure 2. We visualized the
ground-truth solution and the approximated neural network solution at the time
$t=2\tau,20\tau,40\tau$. Overall, the deep minimizing movement scheme achieved
a good result on the 2d heat equation. Our approximated neural network
solution achieved the relative $L^{2}$ error around $3.15\times 10^{-3}$ at
$t=0.2=40\tau$.

Figure 2: Ground-truth solution (top row) and our neural network solution
(middle row) with deep minimizing movement scheme at time $t=0.01,0.1,0.2$.
(At $t=2\tau,20\tau,40\tau$) The relative $L^{2}$ error for each timestep is
plotted in the graph at the bottom row.

##### 8-dimensional case

In the experiment for 8-dimensional heat equation, we set $\kappa=0.25$, and
$a=\pi/4$. The domain $\Omega$ is given by $[-\pi/4,\pi/4]^{8}$. For an
initial condition
$u_{0}(x)=1-0.2\cos{8x_{2}}+2\cos{4x_{4}}-\cos{4x_{6}}+\cos{4x_{8}}+2\cos{4x_{1}}\cos{4x_{3}}\cos{4x_{8}}+2\cos{4x_{2}}\cos{4x_{5}}\cos{4x_{7}}$,
its ground-truth solution $u(t,x)$ is computed by (25), the same as the
2-dimensional case. A size of the time-step for the minimizing movement scheme
is given by $\tau=0.005$.

##### Model training details for 8d case

For the neural network model that approximates the solution $u(x)$, a 3-layer
fully connected neural network with 256 hidden units is used. A tangent
hyperbolic (Tanh) activation function is applied. We choose the learning rate
as $10^{-5}$ with an Adam optimizer. While training the neural network model,
we uniformly picked 50,000 samples from $\Omega=[-\pi/4,\pi/4]^{8}$ and
computed the minimization target (21) by the Monte-Carlo approximation.

Figure 3: Ground-truth solution (top row) and our neural network solution
(middle row) with deep minimizing movement scheme at time $t=0.2=40\tau$. The
relative $L^{2}$ error for each timestep is plotted in the graph at the bottom
row. In the relative $L^{2}$ error plot, we estimated the error with 500
trials, each evaluated on 50,000 points that are uniformly sampled from
$\Omega$. Blue dotted line indicates the maximum error and the red dotted line
indicates the minimum error among the 500 trials for each timestep. As the
plots of the maximum error and the minimum error are almost coincident, it can
be deduced that the evaluated result is surely robust.

##### Results - 8d case

For the high dimensional heat equation, we illustrated our approximated neural
network solution and its relative $L^{2}$ error result in Figure 3, by
visualizing the solution on the 4 hyperplanes $x_{1}x_{2}$, $x_{3}x_{4}$,
$x_{5}x_{6}$, $x_{7}x_{8}$-planes. Each solution $u_{sliced}(x_{i},x_{i+1})$
was evaluated and visualized by fixing the other 6 variables to zero and
varying $(x_{i},x_{i+1})\in[-\pi/4,\pi/4]^{2}$ for $i=1,3,5,7$. We plotted the
ground-truth solution and the approximated neural network solution at the time
$t=40\tau$. Its relative $L^{2}$ error is estimated 500 times, each with
50,000 points that are uniformly sampled from the domain $\Omega$. The deep
minimizing movement scheme also achieved a good result on the 8d heat
equation. Our approximated neural network solution achieved the relative
$L^{2}$ error around $3.67\times 10^{-2}$ at $t=0.2=40\tau$. These results
clearly showed the applicability of our proposed method as a cornerstone of
dealing with high dimensional PDEs.

####  4.1.2 Allen-Cahn equation: 2-dimensional

The Allen-Cahn equation

(26) $\begin{cases}\partial_{t}u=\Delta
u-\epsilon^{-2}F^{\prime}(u)\quad\text{in }\Omega,\\\ \nabla
u\cdot\hat{\boldsymbol{n}}=0\quad\text{on }\partial\Omega,\end{cases}$

where $\boldsymbol{n}$ denotes the unit normal vector, is the $L^{2}$-gradient
flow of the functional $I_{\epsilon}(u):=\frac{1}{2}\int||\nabla
u||^{2}dx+\epsilon^{-2}\int F(u)dx$, where $\epsilon>0$ and $F$ is a double
well potential. In our experiment, we consider the case where the double well
potential $F(s)$ is given by $\frac{(s^{2}-1)^{2}}{4}$. We set $\epsilon=0.25$
and the domain $\Omega=[-2,2]\times[-2,2]$. The ground-truth solution of (26)
is computed with a forward Euler scheme, where the initial condition in Fig. 4
is examined. A size of the time-step for the minimizing movement scheme is
given by $\tau=0.005$.

Figure 4: A hard torus-shape initial condition (left) and a 4-point initial
condition (right) are used for the experiment on the Allen-Cahn equation.
Black color indicates -1 and white color indicates 1.

##### Model training details

For the neural network model that approximates the solution $u(x)$, a 2-layer
fully connected neural network with 256 hidden units is used. An exponential
linear unit (ELU) activation function is applied. We choose the learning rate
as $10^{-3}$ with an Adam optimizer. As in the 2d heat equation case, we
uniformly picked 10,000 samples from $\Omega=[-2,2]^{2}$ and computed the
minimization target.

##### Results

The results for the hard torus-shape initial condition and the 4-point initial
condition in Fig. 4 are visualized in Fig. 5 and Fig. 6 respectively. The
ground-truth solution is obtained with the forward Euler scheme. We visualized
the ground-truth solution and the approximated neural network solution at the
time $t=2\tau,20\tau,40\tau$ and plotted the relative $L^{2}$ error for each
timestep. The trained models for the two initial conditions at $t=0.2=40\tau$
achieved its relative $L^{2}$ error $3.76\times 10^{-2}$ and $6.69\times
10^{-2}$ respectively. Although the borderline area between -1 and +1 evolves
rapidly and has a steep slope, the neural network approximator and our
proposed deep minimizing movement scheme achieved good results for the given
two initial conditions.

Figure 5: Ground-truth numerical solution (top row) and our neural network
solution (middle row) with deep minimizing movement scheme of the Allen-Cahn
equation at time $t=0.01,0.1,0.2$. (At $t=2\tau,20\tau,40\tau$), based on the
hard torus-shape initial condition. Its relative $L^{2}$ error for each
timestep is plotted in the graph at the bottom row.

Figure 6: Ground-truth numerical solution (top row) and our neural network
solution (middle row) with deep minimizing movement scheme of the Allen-Cahn
equation at time $t=0.01,0.1,0.2$. (At $t=2\tau,20\tau,40\tau$), based on the
4-point initial condition. Its relative $L^{2}$ error for each timestep is
plotted in the graph at the bottom row.

###  4.2 $\boldsymbol{\mathbb{W}_{2}}$-gradient flows

In this section, we consider several examples of $\mathbb{W}_{2}$-gradient
flows. In the following experiments, we used the Softplus activation function
for the last layer to guarantee the positivity of neural networks. Activation
functions used in the remaining layers except for the last layer are specified
for each experiment.

####  4.2.1 Heat equation: 2-dimensional

We consider

$\begin{cases}\partial_{t}u=\kappa\Delta u\quad\text{in
}[0,T]\times[-\pi,\pi]^{2},\\\ \nabla
u\cdot\hat{\boldsymbol{n}}=0\quad\text{on
}[0,T]\times\partial[-\pi,\pi]^{2},\\\
u(0,x_{1},x_{2})=(1+\frac{1}{5}\cos(x_{1})+\frac{1}{5}\cos(x_{2}))/(4\pi^{2}),\end{cases}$

which corresponds to a $\mathbb{W}_{2}$-gradient flow of the entropy
functional $\mathcal{F}(u):=\int_{\Omega}u(x)\log u(x)dx$.

##### Model training details

In this numerical experiment, we employ a 4-layer fully connected neural
network with 512 hidden units with cosine activation function. We use Adam
optimizer with a learning rate $2\times 10^{-4}$. We set $\tau=0.005$. For
each training epoch, we sample 20,000 points from $[-\pi,\pi]\times[-\pi,\pi]$
uniformly. Numerical integrations are computed by using the Monte-Carlo
method.

Figure 7: Left: Analytic solutions and approximated solutions of (4.2.1) at
$t=\tau,51\tau$. Right: $L^{2}$ error and relative $L^{2}$ error of an
approximated solution.

####  4.2.2 Porous Medium: 1-dimensional

The porous medium equation

(27) $\partial_{t}u(t,x)=\Delta u^{m},\quad\text{in
}[0,T]\times\mathbb{R}^{n}$

is the $\mathbb{W}_{2}$-gradient flow of the energy
$\mathcal{F}(u):=\frac{1}{m-1}\int_{\mathbb{R}^{n}}u^{m}(x)dx$.

A fundamental example of exact solution of this equation was obtained
independently by Barenblatt and Pattle [6], which are densities of the form

(28)
$u(t,x)=(t+t_{0})^{-kn}\bigg{(}C-\frac{(m-1)k}{2m}||x||^{2}t^{-2k}\bigg{)}_{+}^{\frac{1}{m-1}},$

where $k=(n(m-1)+2)^{-1}$ and the positive constant $C$ is defined by the
identity $\int_{\mathbb{R}^{n}}u(t,x)dx=1$.

In this experiment, we set $t_{0}=10^{-1}$, $m=2$, and
$C=(\frac{3}{64})^{\frac{1}{3}}$. Moreover, we set $u(x,0)$ in the equation
(28) to be the initial condition.

##### Model training details

In this numerical experiment, we employ a 5-layer fully connected neural
network with 512 hidden units with ELU activation function. We use Adam
optimizer with a learning rate $2\times 10^{-4}$. We set $\tau=0.005$ and
truncated the domain $\mathbb{R}$ by $[-3,3]$ so that the numerical
integration of the initial condition has a negligible error. For each training
epoch, we sample $3,000$ points from $[-3,3]$ uniformly. Numerical
integrations in the algorithm are computed by using the Monte-Carlo method.

##### Results

Both the neural network solution and the ground-truth solution over time until
$t=60\tau$ are plotted in Figure 8 and error estimation is given in Figure 9.
As shown in Figure 8, despite the accumulation of errors over time, it fits
well for a fairly long time.

Figure 8: Neural network solution and the ground-truth solution for the
1-dimensional porous medium equation over timesteps $t=\tau,\cdots,60\tau$.
Figure 9: Error estimation of the neural network solution for the porous
medium equation.

####  4.2.3 Porous medium equation: 2-dimensional

In this experiment, we set $t_{0}=10^{-1},m=2$, and $C=\frac{1}{\sqrt{8\pi}}$.
As in the $1$-dimensional case, we take the $u(x,0)$ in the equation
$\eqref{porous_sol}$ to be the initial condition.

##### Model training details

In this numerical experiment, we employ a 5-layer fully connected neural
network with 512 hidden units with ELU activation function. We use Adam
optimizer with a learning rate $2\times 10^{-4}$. We set $\tau=0.005$ and
truncated the domain $\mathbb{R}^{2}$ by $[-3,3]\times[-3,3]$ so that the
numerical integration of the initial condition has a negligible error. For
each training epoch, we sample $10,000$ points from $[-3,3]\times[-3,3]$
uniformly. Numerical integrations in the algorithm are computed by using the
Monte-Carlo method.

##### Results

We summarize the results in Figure 10 and Figure 11. In Figure 10, the neural
network solution and the ground-truth solution are both plotted for timesteps
$t=\tau$ and $t=60\tau$. As in the 1-dimensional case, neural network solution
fits well for a fairly long time. It worths noticing that since we sampled
more points in this experiment than in the 1-dimensional experiment, the
errors are smaller.

Figure 10: Neural network solution and the ground-truth solution for the
2-dimensional porous medium equation for timesteps $t=\tau$ and $t=60\tau$.
Figure 11: Error estimation of the neural network solution for the
2-dimensional porous medium equation.

####  4.2.4 Fokker-Planck equation: 2-dimensional

The Fokker-Planck equation

(29) $\partial_{t}u(t,x)=\nabla\cdot(u\nabla V)+\Delta u,\quad\text{in
}[0,T]\times\mathbb{R}^{n}$

is the $\mathbb{W}_{2}$-gradient flow of the energy
$\mathcal{F}(u):=\int_{\mathbb{R}^{2}}V(x)u(x)+u(x)\log u(x)dx$. In this
numerical example, we set the initial condition to be the standard Gaussian
distribution and set $V(x)=\frac{1}{2}(x-\mu)^{T}\Sigma^{-1}(x-\mu)$, with
$\mu=[\frac{1}{3},\frac{1}{3}]$, and
$\Sigma=\begin{bmatrix}\frac{5}{8}&-\frac{3}{8}\\\
-\frac{3}{8}&\frac{5}{8}\end{bmatrix}$. Then, an analytic solution of (29)
reads

(30)
$u(t,x)=\mathcal{N}(\mu(t),\Sigma(t)):=\frac{1}{2\pi\sqrt{|\Sigma(t)|}}\exp{-\frac{1}{2}(x-\mu(t))^{T}\Sigma(t)^{-1}(x-\mu(t))},$

where $\mu(t)=(1-e^{-4t})\mu$,
$\Sigma(t)=\begin{bmatrix}\frac{5}{8}+\frac{3}{8}\times
e^{-8t}&-\frac{3}{8}+\frac{3}{8}\times e^{-8t}\\\
-\frac{3}{8}+\frac{3}{8}\times e^{-8t}&\frac{5}{8}+\frac{3}{8}\times
e^{-8t}\end{bmatrix}$.

##### Model training details

In this numerical experiment, we employ a 4-layer fully connected neural
network with 512 hidden units with ReLU activation function. We use Adam
optimizer with a learning rate $2\times 10^{-4}$. We set $\tau=0.005$ and
truncated the domain $\mathbb{R}^{2}$ by $[-5,5]\times[-5,5]$ so that the
numerical integration of the initial condition has a negligible error. For
each training epoch, we sample 10,000 points from $[-5,5]\times[-5,5]$
uniformly. Numerical integrations in the algorithm are computed by using the
Monte-Carlo method.

##### Results

We summarize the results in Figure 12. The left panel of Figure 12 shows both
an analytic solution (30) and a neural network solution at $t=\tau$, and
$t=21\tau$. The right panel of the figure shows trajectories of mean vector of
an analytic solution and a neural network solution. We achieved a relative
error less than 0.01 (averaged in time) in the truncated domain, where
$u_{NN}$ denotes a neural network solution.

Figure 12: Left: Analytic solutions and approximated solutions of (29) at
$t=\tau,\cdots,21\tau$. Right: Mean trajectories of an analytic solution and
an approximated solution.

####  4.2.5 Fokker-Planck equation: 4-dimensional

In this experiment, we set the initial condition to be the standard Gaussian
distribution and $V(x)=\frac{1}{2}(x-\mu)^{T}\Sigma^{-1}(x-\mu)$, where
$\mu=[\frac{1}{3},\frac{1}{3},0,0]$ and
$\Sigma=\begin{bmatrix}\frac{5}{8}&-\frac{3}{8}\\\
-\frac{3}{8}&\frac{5}{8}\end{bmatrix}\bigoplus\begin{bmatrix}1&0\\\
0&1\end{bmatrix}$. In this case, the analytic solution of (29) is given by
(30), where $\mu(t)=[\frac{1}{3}(1-e^{-4t}),\frac{1}{3}(1-e^{-4t}),0,0]$ and
$\Sigma(t)=\begin{bmatrix}\frac{5}{8}+\frac{3}{8}\times
e^{-8t}&-\frac{3}{8}+\frac{3}{8}\times e^{-8t}\\\
-\frac{3}{8}+\frac{3}{8}\times e^{-8t}&\frac{5}{8}+\frac{3}{8}\times
e^{-8t}\end{bmatrix}\bigoplus\begin{bmatrix}1&0\\\ 0&1\end{bmatrix}$.

##### Model training details

In this numerical experiment, we employ a 5-layer fully connected neural
network with 256 hidden units with ELU activation function. We use Adam
optimizer with a learning rate $2\times 10^{-4}$. We set $\tau=0.005$ and
truncated the domain $\mathbb{R}^{4}$ by $[-3.5,3.5]^{4}$, and in this case,
the error of the numerical integration of the initial condition on this domain
and the whole domain is less than $2\%$. For each training epoch, we sampled
160,000 points from $[-3.5,3.5]^{4}$ uniformly. Numerical integrations in the
algorithm are computed by using the Monte-Carlo method.

##### Results

The mean trajectories of the neural network solution and the ground-truth
solution along timesteps $t=\tau,\cdots,60\tau$ are given in Figure 13. The
first and second coordinates of the mean trajectories of both neural network
solution and the ground-truth solution are in the relation $y=x$ as
illustrated in the left top of Figure 13. The right panel of Figure 13 is the
mean trajectories of both neural network solution and the ground-truth
solution over $t=0,\cdots,20\tau$, and it can be checked that they are
overlapping.

Figure 13: Mean trajectories of our optimized neural network and the ground-
truth solution. The trajectories taken from $\tau$ to $60\tau$ are plotted
with respect to the first and the second coordinates on the top left, and with
respect to the third and fourth coordinates on the bottom left. On the right,
trajectories of neural network and the ground-truth solution taken for $t=0$
$20\tau$ are plotted in the same plot together for detailed comparison.

The heat maps of cross-sectional $u_{sliced}$’s are given in Figure 14. The
ground-truth solution of $u_{sliced}(x_{3},x_{4})$ converges to a Gaussian
distribution with the variance $\begin{bmatrix}\frac{5}{8}&-\frac{3}{8}\\\
-\frac{3}{8}&\frac{5}{8}\end{bmatrix}$ while the ground-truth solution of
$u_{sliced}(x_{1},x_{2})$ remains unchanged. As time goes by, the shape of the
variance of the neural network solution of $u_{sliced}(x_{3},x_{4})$ follows
$\begin{bmatrix}\frac{5}{8}&-\frac{3}{8}\\\
-\frac{3}{8}&\frac{5}{8}\end{bmatrix}$. Meanwhile, even with the passage of
time, the neural network solution of $u_{sliced}(x_{3},x_{4})$ remains the
same as the standard Gaussian distribution.

Figure 14: Heat maps of cross-sectional $u_{sliced}$ for the neural network
solution and the ground-truth solution at the origin for timesteps $t=\tau$
and $t=60\tau$ respectively. Our neural network solution for
$u_{sliced}(x_{1},x_{2})$ becomes anisotropic, leaving
$u_{sliced}(x_{3},u_{4})$ isotropic, as time passes. Since the ground-truth
solution has this property, this pattern is a desiderata for approximated
solutions.

Error estimation is given in Figure 15. The right panel of Figure 15 is the
plot of the mean of the relative $L^{2}$ error over time, i.e.
$\frac{1}{K}\sum_{k=1}^{K}[\text{relative }L^{2}\text{ error at }t=k\tau]$.

Figure 15: Errors of the neural network solution for $4$-dimensional Fokker-
Planck equation.

##  5 Discussion

In the experiments on several gradient flow equations above, the proposed
mesh-free deep minimizing movement scheme achieved outstanding results for
both low and high dimensions. Despite lack of mature error analysis, it has
been known that neural networks are more fluent than the classical numerical
methods because they are free from mesh-generation. By the series of
experiments in this work, we verified that neural networks perform great for
both low and high dimensions, and this indicates the possibility of applying
our method to solving various PDEs related to gradient flows.

As pointed out in [4], mesh generation often fails to scale when the domain
has a complex geometry. Our approach for both $L^{2}$ and $\mathbb{W}_{2}$
cases, which is completely free of mesh generation, is therefore directly
applicable for the gradient-flow type equation in a domain with complex
geometry. Our method also takes advantage in that it can be directly applied
not only to $L^{2}$ and $\mathbb{W}_{2}$ spaces, but also to any general
space. Therefore, if some classes of PDEs are realized as gradient flows in
some spaces other than $L^{2}$ and $\mathbb{W}_{2}$ in future, our method can
be a breakthrough for numerically solving those equations.

## Acknowledgement

H. J. Hwang is supported by the National Research Foundation of Korea (NRF)
grant funded by the Korea government (MSIT) (No. NRF-2017R1E1A1A03070105 and
NRF-2019R1A5A1028324), Institute for Information & Communications Technology
Promotion (IITP) grant funded by the Korea government (MSIT) (No.
2019-0-01906, Artificial Intelligence Graduate School Program (POSTECH)), and
the Information Technology Research Center (ITRC) support program (No.
IITP-2018-0-01441). H. Son is supported by National Research Foundation of
Korea (NRF) grants funded by the Korean government (MSIT) (No.
NRF-2019R1A5A1028324).

## References

  * [1] D. Alvarez-Melis, Y. Schiff, and Y. Mroueh, Optimizing functionals on the space of probabilities with input convex neural networks, arXiv preprint arXiv:2106.00774, (2021).
  * [2] L. Ambrosio and N. Gigli, A user’s guide to optimal transport, in CIME summer school, Italy, 2009, https://hal.archives-ouvertes.fr/hal-00769391.
  * [3] B. Amos, L. Xu, and J. Z. Kolter, Input convex neural networks, in International Conference on Machine Learning, PMLR, 2017, pp. 146–155.
  * [4] J. Berg and K. Nyström, A unified deep artificial neural network approach to partial differential equations in complex geometries, Neurocomputing, 317 (2018), pp. 28–41.
  * [5] A. Braides, Local minimization, variational evolution and $\Gamma$-convergence, vol. 2094, Springer, 2014.
  * [6] J. A. Carrillo and G. Toscani, Asymptotic l 1-decay of solutions of the porous medium equation to self-similarity, Indiana University Mathematics Journal, (2000), pp. 113–142.
  * [7] S. W. Cho, H. J. Hwang, and H. Son, Traveling wave solutions of partial differential equations via neural networks, Journal of Scientific Computing, 89 (2021).
  * [8] D.-A. Clevert, T. Unterthiner, and S. Hochreiter, Fast and accurate deep network learning by exponential linear units (elus), arXiv preprint arXiv:1511.07289, (2015).
  * [9] L. Courte and M. Zeinhofer, Robin pre-training for the deep ritz method, arXiv preprint arXiv:2106.06219, (2021).
  * [10] G. Cybenko, Approximation by superpositions of a sigmoidal function, Mathematics of control, signals and systems, 2 (1989), pp. 303–314.
  * [11] E. De Giorgi, Movimenti minimizzanti, in Aspetti e problemi della Matematica oggi, Proc. of Conference held in Lecce, 1992.
  * [12] E. De Giorgi, A. Marino, and M. Tosques, Problems of evolution in metric spaces and maximal decreasing curve, Atti Accad. Naz. Lincei Rend. Cl. Sci. Fis. Mat. Natur.(8), 68 (1980), pp. 180–187.
  * [13] M. Dissanayake and N. Phan-Thien, Neural-network-based approximations for solving partial differential equations, communications in Numerical Methods in Engineering, 10 (1994), pp. 195–201.
  * [14] C. Dugas, Y. Bengio, F. Bélisle, C. Nadeau, and R. Garcia, Incorporating second-order functional knowledge for better option pricing, Advances in neural information processing systems, (2001), pp. 472–478.
  * [15] K. Hornik, Approximation capabilities of multilayer feedforward networks, Neural networks, 4 (1991), pp. 251–257.
  * [16] H. J. Hwang, J. W. Jang, H. Jo, and J. Y. Lee, Trend to equilibrium for the kinetic fokker-planck equation via the neural network approach, Journal of Computational Physics, (2020), p. 109665.
  * [17] H. Jo, H. Son, H. J. Hwang, and E. H. Kim, Deep neural network approach to forward-inverse problems, Networks & Heterogeneous Media, 15 (2020), pp. 247–259.
  * [18] R. Jordan, D. Kinderlehrer, and F. Otto, The variational formulation of the fokker–planck equation, SIAM journal on mathematical analysis, 29 (1998), pp. 1–17.
  * [19] D. P. Kingma and J. Ba, Adam: A method for stochastic optimization, arXiv preprint arXiv:1412.6980, (2014).
  * [20] A. Korotin, V. Egiazarian, A. Asadulaev, A. Safin, and E. Burnaev, Wasserstein-2 generative networks, arXiv preprint arXiv:1909.13082, (2019).
  * [21] I. E. Lagaris, A. Likas, and D. I. Fotiadis, Artificial neural networks for solving ordinary and partial differential equations, IEEE transactions on neural networks, 9 (1998), pp. 987–1000.
  * [22] J. Y. Lee, J. W. Jang, and H. J. Hwang, The model reduction of the vlasov-poisson-fokker-planck system to the poisson-nernst-planck system via the deep neural network approach, arXiv preprint arXiv:2009.13280, (2020).
  * [23] M. Leshno, V. Y. Lin, A. Pinkus, and S. Schocken, Multilayer feedforward networks with a nonpolynomial activation function can approximate any function, Neural networks, 6 (1993), pp. 861–867.
  * [24] X. Li, Simultaneous approximations of multivariate functions and their derivatives by neural networks with one hidden layer, Neurocomputing, 12 (1996), pp. 327–343.
  * [25] Y. Liao and P. Ming, Deep nitsche method: Deep ritz method with essential boundary conditions, arXiv preprint arXiv:1912.01309, (2019).
  * [26] S. Liu, W. Li, H. Zha, and H. Zhou, Neural parametric fokker-planck equations, arXiv preprint arXiv:2002.11309, (2020).
  * [27] L. Lu, X. Meng, Z. Mao, and G. E. Karniadakis, Deepxde: A deep learning library for solving differential equations, SIAM Review, 63 (2021), pp. 208–228.
  * [28] A. Makkuva, A. Taghvaei, S. Oh, and J. Lee, Optimal transport mapping via input convex neural networks, in International Conference on Machine Learning, PMLR, 2020, pp. 6672–6681.
  * [29] R. J. McCann, Existence and uniqueness of monotone measure-preserving maps, Duke Mathematical Journal, 80 (1995), pp. 309–323.
  * [30] L. McClenny and U. Braga-Neto, Self-adaptive physics-informed neural networks using a soft attention mechanism, arXiv preprint arXiv:2009.04544, (2020).
  * [31] M. Mizuno and Y. Tonegawa, Convergence of the allen–cahn equation with neumann boundary conditions, SIAM Journal on Mathematical Analysis, 47 (2015), pp. 1906–1932.
  * [32] P. Mokrov, A. Korotin, L. Li, A. Genevay, J. Solomon, and E. Burnaev, Large-scale wasserstein gradient flows, arXiv preprint arXiv:2106.00736, (2021).
  * [33] G. Monge, Mémoire sur la théorie des déblais et des remblais, Histoire de l’Académie Royale des Sciences de Paris, (1781).
  * [34] J. Müller and M. Zeinhofer, Deep ritz revisited, arXiv preprint arXiv:1912.03937, (2019).
  * [35] J. Müller and M. Zeinhofer, Notes on exact boundary values in residual minimisation, arXiv preprint arXiv:2105.02550, (2021).
  * [36] V. Nair and G. E. Hinton, Rectified linear units improve restricted boltzmann machines, in Icml, 2010.
  * [37] F. Otto, The geometry of dissipative evolution equations: the porous medium equation, Communications in Partial Differential Equations, 26 (2001).
  * [38] A. Pratelli, On the equality between monge’s infimum and kantorovich’s minimum in optimal mass transportation, in Annales de l’Institut Henri Poincare (B) Probability and Statistics, vol. 43, Elsevier, 2007, pp. 1–13.
  * [39] M. Raissi, P. Perdikaris, and G. E. Karniadakis, Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations, Journal of Computational Physics, 378 (2019), pp. 686–707.
  * [40] F. Santambrogio, Optimal transport for applied mathematicians, Birkäuser, NY, 55 (2015), p. 94.
  * [41] J. Sirignano and K. Spiliopoulos, Dgm: A deep learning algorithm for solving partial differential equations, Journal of Computational Physics, 375 (2018), pp. 1339–1364.
  * [42] H. Son, J. W. Jang, W. J. Han, and H. J. Hwang, Sobolev training for the neural network solutions of pdes, arXiv preprint arXiv:2101.08932, (2021).
  * [43] A. Taghvaei and A. Jalali, 2-wasserstein approximation via restricted convex potentials with application to improved training for gans, arXiv preprint arXiv:1902.07197, (2019).
  * [44] R. van der Meer, C. Oosterlee, and A. Borovykh, Optimally weighted loss functions for solving pdes with neural networks, arXiv preprint arXiv:2002.06269, (2020).
  * [45] C. Villani, Topics in optimal transportation, vol. 58, American Mathematical Soc., 2003.
  * [46] S. Wang, X. Yu, and P. Perdikaris, When and why pinns fail to train: A neural tangent kernel perspective, arXiv preprint arXiv:2007.14527, (2020).
  * [47] E. Weinan and B. Yu, The deep ritz method: a deep learning-based numerical algorithm for solving variational problems, Communications in Mathematics and Statistics, 6 (2018), pp. 1–12.


$L_{2}(\Omega,\mathbb{Q})$ as well. Let $f\in\mathcal{F}_{b}(\Omega)$ be
fixed. Because $f$ is bounded, $\mathfrak{R}^{n}f\in\mathcal{F}_{b}(\Omega)$,
and by the Markov property of $W$ it follows that

(B.3) $\mathbb{E}_{w}[f(W_{n})]=\mathfrak{R}^{n}f(w),\quad
w\in\Omega,\,n\in\mathbb{N}.$

Because $\mathfrak{R}$ is self-adjoint, we see that $\mathfrak{R}^{2}$ is a
nonnegative-definite operator on the Hilbert space $L_{2}(\Omega,\mathbb{Q})$
equipped with the natural $L_{2}$ inner product, and this allows us to
conclude from [52, Corollary 3] (see also [48, Theorem 1]) that there exist
$\psi,\widehat{\psi}\in L_{2}(\Omega,\mathbb{Q})$ satisfying

(B.4)
$\psi=\lim\limits_{n\to\infty}\mathfrak{R}^{2n}f,\quad\widehat{\psi}=\lim\limits_{n\to\infty}\mathfrak{R}^{2n+1}f,$

where the convergence is in $L_{2}$-norm and $\mathbb{Q}$ almost everywhere.
It is worth mentioning that the convergence in (B.4), which follows from [52,
Corollary 3], essentially uses the classical Banach principle (see e.g., [1])
along with a maximal ergodic inequality. The convergence, in fact, holds for
any function in $(L\log^{+}L)(\Omega,\mathbb{Q})$. By the almost sure
convergence of $\mathfrak{R}^{2n}f$ (resp. $\mathfrak{R}^{2n+1}f$) and the
$L_{\infty}$ contractivity of $\mathfrak{R}$, we see that
$\psi,\widehat{\psi}\in L_{\infty}(\Omega,\mathbb{Q})$ as well. The $L_{2}$
contractivity of the linear operator $\mathfrak{R}$ also implies that,

(B.5)
$\psi=\mathfrak{R}\widehat{\psi},\quad\widehat{\psi}=\mathfrak{R}\psi\quad\mathbb{Q}\text{-a.s.,}$

from which we get

(B.6) $\mathfrak{R}^{2}\psi=\psi,\
\mathfrak{R}^{2}\widehat{\psi}=\widehat{\psi},\quad\mathbb{Q}\text{ a.s.}$

We claim that if $-1$ is not an eigenvalue of $\mathfrak{R}$ as an operator on
$L_{\infty}(\Omega,\mathbb{Q})$, then we must have

(B.7) $\psi=\widehat{\psi}=\int_{\Omega}f\,{\rm
d}\mathbb{Q},\quad\mathbb{Q}\text{ a.s.}$

Note that (B.1) will follow once we prove (B.7). Indeed, (B.7) combined with
(B.3)–(B.4) implies that, for $\mathbb{Q}$ a.s. $w\in\Omega$, both the odd and
even subsequence of $(\mathbb{E}_{w}[f(W_{n})])_{n\in\mathbb{N}}$ converge to
the same limit $\mathbb{Q}(f):=\int_{\Omega}f\,{\rm d}\mathbb{Q},$ which
necessarily forces the convergence of $\mathbb{E}_{w}[f(W_{n})]$ to
$\mathbb{Q}(f)$ as $n\to\infty$.

To prove (B.7), it suffices to show that $\psi$ and $\widehat{\psi}$ are
constant $\mathbb{Q}$ a.s., because by the invariance of $\mathfrak{R}$ w.r.t.
$\mathbb{Q}$ and bounded convergence we have

(B.8) $\int_{\Omega}\psi\,{\rm
d}\mathbb{Q}=\lim\limits_{n\to\infty}\int_{\Omega}\mathfrak{R}^{2n}f\,{\rm
d}\mathbb{Q}=\int_{\Omega}f\,{\rm
d}\mathbb{Q}=\lim\limits_{n\to\infty}\int_{\Omega}\mathfrak{R}^{2n+1}f\,{\rm
d}\mathbb{Q}=\int_{\Omega}\widehat{\psi}\,{\rm d}\mathbb{Q}.$

We only prove the claim for $\psi$, as the same argument works for
$\widehat{\psi}$. Let us set $g:=\mathfrak{R}\psi-\psi$. From (B.6), we have

(B.9) $\mathfrak{R}g=-g,\quad\mathbb{Q}\text{ a.s.},$

and also $\|g\|_{\infty}\leq 2\|\psi\|_{\infty}<\infty$. Thus, $g\in
L_{\infty}(\Omega,\mathbb{Q})$ is such that $\mathbb{Q}$-a.s.
$\mathfrak{R}g=-g$, and hence by our assumption we must have $g=0$ a.s. In
other words, $\mathbb{Q}$-a.s. $\mathfrak{R}\psi-\psi=0$ and therefore
ergodicity of $\mathfrak{R}$ in $L_{2}(\Omega,\mathbb{Q})$, which is
equivalent to the ergodicity of $W$ under $\mathbb{Q}$, implies that $\psi$ is
necessarily a constant $\mathbb{Q}$ a.s. ∎

## References

  * [1] A. Bellow and R. L. Jones, A Banach principle for ${L}_{\infty}$, _Adv. Math._ 120 (1996), 155–172. MR1392277
  * [2] S. A. Bethuelsen, M. Birkner, A. Depperschmidt, and T. Schlüter, Local limit theorems for a directed random walk on the backbone of a supercritical oriented percolation cluster, 2021, Preprint. arXiv:2105.09030
  * [3] M. Biskup, Recent progress on the random conductance model, _Probab. Surveys_ 8 (2011), 294–373. MR861133
  * [4] J. Blath, E. Buzzoni, A. G. Casanova, and M. Wilke-Berenguer, Structural properties of the seed bank and the two island diffusion, _J. Math. Biol._ 79 (2019), no. 1, 369–392. MR3975877
  * [5] J. Blath, A. G. Casanova, N. Kurt, and D. Spano, The ancestral process of long-range seed bank models, _J. Appl. Probab._ 50 (2013), no. 3, 741–759. MR3102512
  * [6] J. Blath, A. G. Casanova, N. Kurt, and M. Wilke-Berenguer, A new coalescent for seed-bank models, _Ann. Appl. Probab._ 26 (2016), no. 2, 857–891. MR3476627
  * [7] J. Blath, B. Eldon, A. G. Casanova, and N. Kurt, Genealogy of a Wright-Fisher Model with Strong Seed Bank Component, _XI Symposium on Probability and Stochastic Processes_ , Springer International Publishing, 2015, pp. 81–100. doi: 10.1007/978-3-319-13984-5_4
  * [8] J. Blath, F. Hermann, and M. Slowik, A branching process model for dormancy and seed banks in randomly fluctuating environments, _J. Math. Biol._ 83 (2021), no. 2, 1–40. MR4289473
  * [9] J. Blath and N. Kurt, Population genetic models of dormancy, _Probabilistic Structures in Evolution_ , European Mathematical Society Publishing House, 2021, pp. 247–265. MR4331863
  * [10] E. Bolthausen and I. Goldsheid, Recurrence and transience of random walks in random environments on a strip, _Commun. Math. Phys._ 214 (2000), no. 2, 429–447. MR1796029
  * [11] E. Bolthausen and A. Sznitman, Ten lectures on random media, vol. 32, Springer Science & Business Media, 2002. MR1890289
  * [12] B. L. Burkholder and Y. S. Chow, Iterates of conditional expectation operators, _Proc. Am. Math. Soc._ 12 (1961), no. 3, 490–495. MR0142144
  * [13] G. Carinci, C. Giardinà, C. Giberti, and F. Redig, Dualities in population genetics: A fresh look with new dualities, _Stoch. Proc. Appl._ 125 (2015), 941–969. MR3303963
  * [14] A. G. Casanova and D. Spanò, Duality and fixation in ${\Xi}$-Wright–Fisher processes with frequency-dependent selection, _Ann. Appl. Probab._ 28 (2018), no. 1, 250–284. MR3770877
  * [15] A. G. Casanova, D. Spanò, and M. Wilke-Berenguer, The effective strength of selection in random environment, 2019, Preprint. arXiv:1903.12121
  * [16] D. Cohen, Optimizing reproduction in a randomly varying environment, _J. Theor. Biol._ 12 (1966), no. 1, 119–129. doi: 10.1016/0022-5193(66)90188-3
  * [17] G. Cohen, C. Cuny, and M. Lin, Almost everywhere convergence of powers of some positive ${L}_{p}$ contractions, _J. Math. Anal._ 420 (2014), no. 2, 1129–1153. MR3240070
  * [18] F. den Hollander, Mixing properties for random walk in random scenery, _Ann. Probab._ 16 (1988), no. 4, 1788 – 1802. MR0958216
  * [19] F. den Hollander and Pederzani. G., Multi-colony wright-fisher with seed-bank, _Indag. Math._ 28 (2017), no. 3, 637–669. MR3651340
  * [20] F. den Hollander and S. Nandan, Spatially inhomogeneous populations with seed-banks: I. duality, existence and clustering, _J. Theor. Probab._ (2021), 1–47. doi: 10.1007/s10959-021-01119-z
  * [21] F. den Hollander and S. Nandan, Spatially inhomogeneous populations with seed-banks: II. clustering regime, _Stoc. Proc. Appl._ 150 (2022), 116–146. MR4419571
  * [22] F. den Hollander and J. E. Steif, Random walk in random scenery: A survey of some recent results, Institute of Mathematical Statistics Lecture Notes - Monograph Series, Institute of Mathematical Statistics, 2006, pp. 53–65. MR2306188
  * [23] D. Dolgopyat and I. Goldsheid, Invariant measure for random walks on ergodic environments on a strip, _Ann. Probab._ 47 (2019), no. 4, 2494–2528. MR3980926
  * [24] D. Dolgopyat and I. Goldsheid, Local limit theorems for random walks in a random environment on a strip, 2019, Preprint. arXiv:1910.12961
  * [25] S. Floreani, C. Giardinà, F. den Hollander, S. Nandan, and F. Redig, Switching interacting particle systems: Scaling limits, uphill diffusion and boundary layer, _J. Stat. Phys._ 186 (2022), no. 3, 1–33. MR4372572
  * [26] A. Greven, F. den Hollander, and M. Oomen, Spatial populations with seed-bank: well-posedness, duality and equilibrium, _Electron. J. Probab._ 27 (2022), 1–88. MR4377129
  * [27] A. Greven, F. den Hollander, and M. Oomen, Spatial populations with seed-bank: renormalization on the hierarchical group, 2022, Preprint. arXiv:2110.02714
  * [28] A. Iwanik and R. Shiflett, The root problem for stochastic and doubly stochastic operators, _J. Math. Anal._ 113 (1986), no. 1, 93–112. MR0826661
  * [29] S. Jansen and N. Kurt, On the notion(s) of duality for Markov processes, _Probab. Surveys_ 11 (2014), 59–120. MR3201861
  * [30] I. Kaj, S. M. Krone, and M. Lascoux, Coalescent theory for seed bank models, _J. Appl. Probab._ 38 (2001), no. 2, 285–300. MR1834743
  * [31] O. Kallenberg, Foundations of Modern Probability, vol. 2, Springer, 1997\. MR1464694
  * [32] H. Kesten, A renewal theorem for random walk in a random environment, _Proc. Sympos._ 31 (1977), 67–77. MR0458648
  * [33] I. Kontoyiannis and S. P. Meyn, Geometric ergodicity and the spectral gap of non-reversible Markov chains, _Probab. Theory. Relat. Fields._ 154 (2011), no. 1-2, 327–339. MR2981426
  * [34] S. M. Kozlov, The averaging method and walks in inhomogeneous environments, _Uspekhi Mat. Nauk_ 40 (1985), no. 2(242), 61–120, 238\. MR786087
  * [35] A. Kulik and M. Scheutzow, Generalized couplings and convergence of transition probabilities, _Probab. Theory. Relat. Fields._ 171 (2017), no. 1-2, 333–376. doi: 10.1007/s00440-017-0779-8
  * [36] S. Lalley, An extension of Kesten’s renewal theorem for random walk in a random environment, _Adv. Appl. Math._ 7 (1986), no. 1, 80–100. MR0834222
  * [37] G. F. Lawler and V. Limic, Random walk: A Modern Introduction, Cambridge University Press, 2009. MR2677157
  * [38] J. T. Lennon, F. den Hollander, M. Wilke-Berenguer, and J. Blath, Principles of seed banks and the emergence of complexity from dormancy, _Nat. Comm._ 12 (2021), no. 1, 4807. doi: 10.1038/s41467-021-24733-1
  * [39] J. T. Lennon and S. E. Jones, Microbial seed banks: the ecological and evolutionary implications of dormancy, _Nat. Rev. Microbiol._ 9 (2011), 119–130. doi: 10.1038/nrmicro2504
  * [40] T. M. Liggett, Interacting Particle Systems, Springer, Berlin, 1985\.
  * [41] T. M. Liggett, Continuous Time Markov Processes: An Introduction, vol. 113, American Mathematical Soc., 2010.
  * [42] M. Lin, On the “zero-two” law for conservative markov processes, _Probab. Theory Relat. Fields_ 61 (1982), no. 4, 513–525. MR0682577
  * [43] S. P. Meyn and R. L. Tweedie, Markov Chains and Stochastic Stability, Springer, London, 1993.
  * [44] E. Nummelin, General Irreducible Markov Chains and Non-Negative Operators, Cambridge University Press, 1984.
  * [45] M. Oomen, Spatial Populations with Seed-Bank, Doctoral Thesis, Leiden University, 2022, pp. 1–368.
  * [46] D. Ornstein, On the pointwise behavior of iterates of a self-adjoint operator, _J. Math. Mech._ 18 (1968), no. 5, 473–477. MR0236354
  * [47] G. O. Roberts and J. S. Rosenthal, General state space Markov chains and MCMC algorithms, _Probab. Surv._ 1 (2004), 20–71. MR2095565
  * [48] G. Rota, An “alternierende verfahren” for general positive operators, _Bull. Am. Math. Soc._ 68 (1962), no. 2, 95–102. MR0133847
  * [49] T. Shiga, An interacting system in population genetics, _Kyoto J. Math._ 20 (1980), no. 2, 213–242. MR0582165
  * [50] T. Shiga, An interacting system in population genetics, II, _Kyoto J. Math._ 20 (1980), no. 4, 723–733. MR0592356
  * [51] W. Shoemaker and J. T. Lennon, Evolution with a seed bank: The population genetic consequences of microbial dormancy, _Evol. Appl._ 11 (2018), no. 1, 60–75. doi: 10.1111/eva.12557
  * [52] E. M. Stein, On the maximal ergodic theorem, _Proc. Natl. Acad. Sci. U.S.A._ 47 (1961), no. 12, 1894–1897. MR0131517
  * [53] A. Tellier, S. Laurent, H. Lainer, P. Pavlidis, and W. Stephan, Inference of seed bank parameters in two wild tomato species using ecological and genetic data, _Proc. Natl. Acad. Sci._ 108 (2011), no. 41, 17052–17057. doi: 10.1073/pnas.1111266108
  * [54] N. I. Wisnoski and L. G. Shoemaker, Seed banks alter metacommunity diversity: The interactive effects of competition, dispersal and dormancy, _Ecol. Lett_ 25 (2022), no. 4, 740–753. doi: 10.1111/ele.13944
  * [55] D. Živković and A. Tellier, Germ banks affect the inference of past demographic events, _Mol. Ecol._ 21 (2012), no. 22, 5434–5446. doi: 10.1111/mec.12039