Entropy-dissipation Informed Neural Network for McKean-Vlasov Type PDEs
Authors: Zebang Shen
Zhenfu Wang
Venue: NeurIPS 2023 poster (2023)
Keywords: Entropy-dissipation
McKean-Vlasov
Navier-Stokes
PDE
Coulomb
singular interaction

Abstract: The McKean-Vlasov equation (MVE) describes the collective behavior of particles subject to drift, diffusion, and mean-field interaction. In physical systems, the interaction term can be singular, i.e. it diverges when two particles collide. Notable examples of such interactions include the Coulomb interaction, fundamental in plasma physics, and the Biot-Savart interaction, present in the vorticity formulation of the 2D Navier-Stokes equation (NSE) in fluid dynamics. Solving MVEs that involve singular interaction kernels presents a significant challenge, especially when aiming to provide rigorous theoretical guarantees. In this work, we propose a novel approach based on the concept of entropy dissipation in the underlying system. We derive a potential function that effectively controls the KL divergence between a hypothesis solution and the ground truth. Building upon this theoretical foundation, we introduce the Entropy-dissipation Informed Neural Network (EINN) framework for solving MVEs. In EINN, we utilize neural networks (NN) to approximate the underlying velocity field and minimize the proposed potential function. By leveraging the expressive power of NNs, our approach offers a promising avenue for tackling the complexities associated with singular interactions. To assess the empirical performance of our method, we compare EINN with SOTA NN-based MVE solvers. The results demonstrate the effectiveness of our approach in solving MVEs across various example problems.

Entropy-dissipation Informed Neural Network for McKean-Vlason Type PDEs
Zebang Shen"
Zhenfu Wang" Peking University zwang0bicmr pku. edu. cn
ETH Zürich zebang. shen@ @n ethz.
Abstract
The McKean-Vlasov equation (MVE) describes the collective behavior parti- cles subject to drift, diffusion, and mean-field interaction. In physical systems, the interaction term can be singular, i.e. diverges when two particles collide. No- table examples of such interactions include the Coulomb interaction, fundamental in plasma physics, and the Biot-Savart interaction, present in the vorticity formu- lation of the 2D Navier-Stokes equation (NSE) in fluid dynamies. Solving MVEs that involve singular interaction kernels presents a significant challenge, especially when aiming to provide rigorous thcoretical guarantees. In this work, we propose a novel approach based on the concept of entropy dissipation in the underlying system. We derive a potential function that effectively controls the KL divergence between hypothesis solution and the ground truth. Building upon this theoret- ical foundation, we introduce the Entropy-diss dissipatinn Informed Neural Network (EINN) framework for solving MVEs. In E EINN, we utilize neural networks (NN) approximate the underlying velocity field and minimize the proposed poten- tial function. By leveraging the expressive power of NNs, our approach offers a promising avenue for tackling the complexities associated with singular interac- tions. To assess the empirical performance of our method, we compare E.NNN with SOTA NN-based MVE solvers. The results demonstrate the effectiveness of our approach in solving MVEs across various example problems.
Introduction
Scientists use Partial Differential Equations (PDEs) to describe natural laws and predict the dynam- ics real-world systems. As PDEs are fundamental importance, a growing area in machine learning is the of neural networks (NN) solve these equations (Han al. 2018 Zhang al.
2018 Raissi et al. 2020 Cai et al. 2021 Karniadakis et al.) 2021 Cuomo et al.) 2022]. An im-
portant category of PDES is the McKean-Vlass equation (MVE), which models the dynamics of
stochastic particle system with mean-field interactions
dx, = vV(x))d
(x(X))d V2ndB; P: =Lw
Here X, X' denotes a random particle' position, X' is either I& or the torus IIdd (a cube [-L, LJd with periodic boundary condition), V : Rd F denotes known potential, K Rd Rd denotes some interaction kernel and the convolution operation is defined as
fx h( y) (y)dy
[Bo is the standard d-dimensional Wiener process with 0 being the diffusion coefficient, and Pe X IR is the law or the probability density function of the random variable X, and the initial data Po is given. Under mild regularity conditions, the density function P, satisfies the MVE
(MVE) (,p(()) div (p.(-VV(x) K px()))) (x),
(2)
Authors are listed in alphahetic order.
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

where div denotes the divergence operator, div h(x)= = 0h./Oz; for velocity field h I R. A denotes the Laplacian operator defined as Aø div(Vó), where Vø denotes the gradient of a scalar function Φ: IR + R. Note that all these operators are applied only on the spatial variable z.
In order to describe dynamics in real-workd phenomena such as electromagnetism [Golse| 2016 and fluid mechanics [Majda et al. 2002], the interaction kernels K the MVE can be highly singular, 0. Two of the most notable examples are the Coulomb interactions
ie.. (x)) oo when
[(d - 2- || (d-2),, d2 3, (-(2T) log |l]||
(Coulomb Kernel) ((x) - V Vg(x), with g(x)
(3)
d= 2,
with (1) denoting the surface area of the unit sphere in TR:", and the vorticity formulation of the 2D Navier-Stokes equation (NSE) where the interaction kernel K is given by the Biot-Savart law
(Biot-Savart Kernel)
K(x)
1 x 2T
(-
22
21
||a|| 2x
||0|| ||0|
where x = (x1, 42) and || denotes the Euclidean norm of vector.
Classical methods for solving MVEs, including finite difference, finite volume, finite element, spec- tral methods, and particle methods, have been developed over time. common drawback of these methods lies in the constraints of their solution representation: Sparse representations, such as less granular grids, cells, meshes, fewer basis functions, or particles, may lead to an inferior solution accuracy; On the other hand, dense representations incur higher computational and memory costs.
As a potent tool for function approximation, NNs are anticipated to overcome these hurdles and handle gigher-dimensionall less regular, and more complex systems efficiently (Weinan et al., 2021]. The most renowned NN-based algorithm is the Physics Informed Neural Network (PINN) Raissi et al. 2019.. The philosophy behind the PINN method is that solving PDE system is equivalent to finding the root of the corresponding differential operators. PINN tackles the latter problem by directly parameterizing the hypothesis solution with NN and training minimize the L func- tional residual of the operators. As a versatile PDE solver, PINN may fail to exploit the underlying dynamics of the PDE, which possibly leads to inferior performance task-specific solvers. For example, on the 2D NSE problem, a recent NN-based development Zhang et al.] [022] surpasses PINN and sets new SOTA empirical performance, which however lacks rigorous theoretical sub- stantiation. Despite the widespread applications of PINN, rigorous етrоr estimation guarantees are scarce in the literature. While we could not find results on the MVE with the Coulomb interaction, only in a very recent paper [De Ryck et al.) 2023 the authors establish for NSE that the PINN loss controls the discrepancy between a candidate solution and the ground truth. We highlight that their result holds average-in-time, meaning that at a particular timestamp e [0. 7). candidate solution with small PINN loss may still significantly differ from the true solution. In contrast, all guarantees in this paper are uniform-in-time. Moreover, there is factor in the aforementioned guarantee that exponentially depends on the total evolving time T, while the factor in our guarantee for the NSE is independent of T. We highlight that these novel improvements are achieved for the proposed EINN framework since we take a completely different route from PINN: Our approach is explicitly designed to exploit the underlying dynamics of the system, as elaborated below.
Our approach Deñine the operator
A[p] def
VV K p- vV log
(5)
By noting Ap; = div(p.V p.), we can rewrite the MVE in the form of a continuity equation
a,p((x) dip(((=))AA) =0.
(6)
For simplicity, we will refer to A[p as the underlying velocity. Consider another time-varying hypothesis velocity field f: R x Rd + IR and let of be the solution to the continuity equation
(hypothesis solution) O.p,( + div(p) (x)) (x) (l,, 0, /6 =
for t [0, T]. where we recall that the initial law Po is known. We will refer to Pf as the hypothesis solution and use the superscript to emphasize its dependence on the hypothesis velocity field f. We propose an entropy-dissipate Informed Neural Network framework (E (EINN), which trains an NN parameterized hypothesis velocity field fo by minimizing the following EINN loss
(EINN loss) R(Jo) def ST" J. (1 A
[x)) Pe (x))ddd)
(8)
2

---

The objective is obtained by studying the stability of carefully constructed Lyapunov functions. These Lyapunov functions draw inspiration from the concept of entropy dissipation in the system, leading to the name of our framework. We highlight that we provide a rigorous error estimation guarantee for our framework for MVEs with singular kernels and e. showing that when R(fo)
is sufficiently small, pt" recovers the ground truth Pt in the KL sense, uniform in time.
Theorem 1 (Informal). Suppose that the initial density function Po is sufficiently regular and the
hypothesis velocity field f((-) = f(t.-) is at least three times continuously differentiable both in t and x. We have for the MVE with a bounded interaction kernel K or with the singular Coulomb O or Biot-Savart ( interaction, the KL divergence between the hypothesis solution p and the ground truth Pt is controlled by the EINN loss for any tE [0, T), i.e. there exists some constant C > 0,
sup KL Pc) S CR(S).
re[o,T]
Having stated our main result, elaborate on the difference between EINN and PINN in terms of information flow over time, which explains why EINN achieves better theoretical guarantees: In PINN, the residuals different time stamps are independent of each other and hence there is no information flow from the residual at time L to the one at time (2(> (1). In contrast, in the EINN loss
80, incorrect estimation made in & will also affect the error at f2 through the hypothesis solution pf- Such an information flow gives a stronger gradient signal when we are trying to minimize the EINN loss, compared to the PINN loss. It partially explains why we can obtain the novel uniform-in-time estimation for PINN and why the constant C. in the NSE case is independent of T for EINN (Theorem 2 but exponential in T for PINN
estimation as opposed to the average-in-tinee
Contributions. In summary, we present novel NN-based framework for solving the MVEs. Our method capitalizes on the entropy dissipation property of the underlying system, ensuring robust theoretical guarantees even when dealing with singular interaction kernels. We elaborate on the contributions of our work from theory, algorithm, and empirical perspectives as follows.
1. (Theory-wise) By studying the stability of the MVEs with bounded interaction kernels or with singular interaction kernels in the Coulomb 3 and the Biot-Savar case 44p (the 2D NSE) via en- tropy dissipation, we establish the error estimation guarantee for the EINN loss on these equations. Specifically, we design potential function R(f) of hypothesis velocity f such that R(f) controls the KL divergence between the hypothesis solution pf (defined in equation 7 and the ground truth solution Pe for any time stamp within a given time interval [0. A direct consequence of this result is that R(f) can be used assess the quality of a generie hypothesis solution the above
MVEs and of exactly recovers ýe in the KL sense given that R(f)=0) 0.
2. (Algorithm-we wiie)) When the hypothesis velocity field is parameterized by an NN, i.e. f with 9 being some inite-dimensionn of the NN parameters 9. We discuss in detail how an estimator of the gradient V ReR(fe fe) can be computed so that stochastic gradient-based optimizers can be utilized to train the NN. In particular, for the 2D NSE (the Biot-Savart case B. we show that the singularity in the gradient computation can be removed by exploiting the derivaeivve of the Biot-Savart kernel.
parameters, the EINN loss R(fe) can be used as the loss function
3. (Empirical-wise) We compare the proposed approach, derived from our novel theoretical guar- antees, with SOTA NN-based algorithms for solving the MVE with the Coulomb interaction and the 2D NSE (the Biot-Savart interaction). We pick specific instances of the initial density Po. under which explicit solutions are known and can be used the ground truth test the quality of the hypothesis ones. Using NNs with the same complexity (depth, width, and structure), we observe that the proposed method significantly outperforms the included baselines.
2 -itropo-dissipation Informed Neural Network
In this section, we present the proposed EINN framework for the MVE. To understand the intuition behind our design, we first write the continuity equation 7 in a similar form as the MVE 6
a.pl (æ) )x) + div pl (x)
[() (ALo
=0,
(10)
where ʃ is the hypothesis velocity (recall that s²(-) = f((,,)) and (Perturbation) (x(x) def fe(z) A[pf((2)
(11)
can be regarded as perturbation to the original MVE system. Consequently, it is natural to study the deviation of the hypothesis solution of from the true solution ग using an appropriate Lyapunov
3

---

function I(p),Pt) The functional relation between this deviation and the perturbation is termed as the stability of the underlying dynamical system, which allows us to derive the EINN loss ( Following this idea, the design of the ETNN loss can be determined by the choice of the Lyapunov function I. used in the stability analysis. In the following, we describe the Lyapunov function used for the MVE with the Coulomb interaction and the 2D NSE (MVE with Biot-Savart interaction). The proof of the following results are the major theoretical contributions of this paper and will be elaborated in the analysis section B
For the MVE with the Coulomb interaction, we choose L. to be the modulated free energy E (defined in equation (23) which is originally proposed in [Bresch et al. [019a] to establish the mean-field limit of a corresponding interacting particle system. We have (setting L E)
1 Elt. p.) s 2 Sxxx (x)) C CEP( P(), dt
(12)
where C' is a universal constant depending on v and (pt)se(o.7).
For the 2D NSE (MVE with the analysis is inspired by (Jabin and Wang 2018) which for the first time establishes the quantitative mean-field limit of the stochastic interacting particle systems where the interaction kernel can be in some negative Sobolev space. We have
interaction), we choose the KL divergence. Our
dt KL p.) V/ f (f (x)
Lee (a)) A.) Sx of (x))]) (13) Pr.
where again C' is universal constant depending on v and (f:)sefo.1]
After applying Grönwall's inequality on the above results, we can see that the EINNN loss 8 is precisely the term derived by stability analysis of the MVE system with an appropriate Lyapunov function. In the next section, we elaborate on how a stochastic approximation of R(f( can be efficiently computed for a parameterized hypothesis velocity field f fo so that stochastic opti- mization methods can be utilized to minimize R( R(fo).
2.1 Stochastic Gradient Computation with Neural Network Parameterization
While the choice of the EINN loss + is theoretically justified through the above stability study, in this section, we show that it admits an estimator which can be efficiently computed. Define the flow map X2 via the ODE (t) (== with x00) 3y such that (+t) = (x0)). From the definition of the push-forward measure, one has p. = Xx Pn- Recall the definitions of the EINN loss R(f) in equation g and the perturbation 8, in equation
II Use the change of variable formula
of the push- push-forward measure (a) and Fubini's 's theorem in (b). We have
R()) [1111122 (2) [." [|ó, Xelle, dt
J."
|1,, (((x))l2dddpp((a)
(14)
Consequently, by defining the trajectory-wiss loss (recall x<t)= X(((x))
R(f: So 18, x.(oro)
S.¹ (x)()))))
(15)
we can write the potential function as an expectation R(f) We (f)] Similarly, when f is parameterized as f fe. we obtain the expectation form VoR(fo) En [ VoRRfoo [();
We show VoR(fe: æ.) can be computed accurately, via the adjoint method (for completeness the derivation of the adjoint method in appendix D). As a recap. suppose that we can write R( fo; 2)) in standard ODE-constrained form R(fe; ((O) fo² g(t, s(t), d) where (s(t))equ,r) is the solution the ODE d s()) (t, s(t): with s0) so. and vb a known transition function. The adjoint method states that the gradient 2 (0)) can be computed as
de
de
S."" a(t)
0g 00
(Adjoint Method)
8 (t, (t))
s():
(16)
d0
where att) solves the final value problems (a(t))+att)* 0. In the following, we focus on how R(fe: can be written in the above ODE-constrained form.
du (t, (t) 0g (t.
Do

---

Write R(fo:oo) in ODE-conttraiee form Expanding the definition of 6, in equation gives
(x(((() [((((()) - ((-)))---VVV)))) K of (a<(1) -V log pt (ze(1)))
(17)
Note that in the above quantity, f and V are known functions. Moreover, is known that v log p (z(t)) admits a closed form dynamics (e.g. see Proposition 2 in Shen et al. 2022
V log (z((t) - -- (divi((( (()) v log pf (x(()),
(J,, (x(():)))
(18)
di
which allows it to be explicitly computed by starting from V log Po (zo) and integrating over time (recall that Po is known). Here J,. denotes the Jacobian matrix of ft. Consequently, all we need to handle is the convolution term K p (z(t)).
A common choice to approximate the convolution operation is via Monte-Carlo integration: Let (1(()) iid of for i 1,...,,, N and denote empirical approximation of pf by n 인 Σ (y)) where (..()) denotes the Dirac measure at (t). We approximate the convolution term in equation T7 in different ways for the Coulomb and the Biot-Savart interactions:
1. For the Coulomb type kernel 3 first approximate K of with Ke pL. where
Ko(x) def ((x)
if |x C;
(19)
if x Se.
If p is bounded in X, we have
sup I(K Ke) * pl (=)) sup 2ex
f ya-y/s se l 1x-20 yl a pl (z) (z)dyll s lle/ (ell(2))) S [yl| dy.
x x-
wex
To compute the integral on the right-hand side, we will switch to polar coordinates (r, b): 1 dr-= J dup J(r.w) S/² dr =c
1 dy S - S Misc by)
[ddr
(20)
Here, J(r.g) denotes the determinant of the Jacobian matrix resulting from the transformation from the Cartesian system to the polar coordinate system. In inequality (20). we utilize the fact that ((r,))) $ pdh1, which allows us cancel out the factor 1/rd-1, Now that Ke bounded by c- we
of can further approximate Ke Pf. using Ke * PN with error of the order O(c (C-d+////) /VN). Altogether,
have sup || pf (zx - K2 16 (x)] O(c c c (+N+//N which be made arbitrarily small for a sufficiently small c and a sufficiently large N.
2. For Biot-Savart interaction (2D Navier-Stoker equation), there are more structures to exploit and we can completely avoid the singularity: As noted by Jabin and Wang 2018 the convolution kernel K can be written in divergence form:
arctan( ms 0,
K = V. U, with U(x)
(21)
2
arctan 82)
where the divergence of matrix function is applied row-wise, L.c. [K(x)]; div U.(z). Using integration by parts and noticing that the boundary integration vanishes on the torus, one has
K of (x) f K(v) (x -yddy f V. (.(())) y)dy= = f U()) (a y)dy U(x (v)V logo! (y)dy
y) log of (z)).
If the score function V V p. is bounded, then the integrand in the expectation is also bounded. Therefore, 1
avoid integrating singular functions and the Monte Carlo-type estimation N U(x U(=-((t))))) log of (1)( is accurate for sufficiently large value of N.
With the above discussion, we can write R(fe: in an ODE-constrained form in a standard way, which due to space limitation is deferred to Appendix D.I
Remark Let n (0) be the function we obtained using the above approximation of the convolution, where N is the number of Monte-Can Carlo samples. The above discussion shows that E (f) and R(Jo; wo) are close in the 6 sense, which is hence sufficient when the EINN loss used as error quantification since only function value malters. When both (0) and R(Jo; xa) are in C², one can translate the closeness in function value the closeness of their gradients. In our experiments, using v () as approximation of V Re (f;; 2o) gives very good empirical performance already.
5

---

3 Analysis
In this section, we focus on the torus case, i.c. X 11 is a box with the periodic boundary condition. This is typical setting considered in the literature as the universal function approximation of NNs only holds over a compact set. Moreover, the boundary integral resulting from integration by parts vanishes in this setting, making it amenable for analysis purposes. For completeness, we provide a discussion on the unbounded case, i.e. X' = Rd in the Appendix [ which requires additional regularity assumptions. Given the MVE 2. if K is bounded, it is sufficient to choose the Lyapunov functional [pp Pt) as the KL divergence (please see Theorem 6 the appendix). But for the singular Coulomb kernel, we need also to consider the modulated energy as in [Serfaty 2020
(Modulated Energy) P(P:P) def
S
())))))p)((()((>>>)))))))))))))) (v))
(22)
where g is the fundamental solution to the Laplacian equation in R. i.e. - -g ōg. and the Coulomb interaction reads K = g (see its closed form expression in equation
If we are only interested
in the deterministic dynamics with Coulomb interactions, i.e.
0 in equation P it suffices to
choose L(p] p) as F(p Pt) (please see Theorem 3 But if we consider the system with Coulomb interactions and diffusions, i.c. O, shall combine the KL divergence and the modulated energy to form the modulated free energy as in [Bresch et al. 2019b], which reads
(Modulated Free Energy) E(p,p) def vKL((, + F(p.
(23)
This definition agrees with the physical meaning that "Free Energy Temperature x Entropy + Energy". and we note that the temperature is proportional to the diffusion coefficient v. We remark also for two probability densities and P. F)) 2 O since by looking in the Fourier domain F(p.p)
dE 2 0 as +{E) 2 O. Moreover, (p..p) can be regarded as a negative Sobolev norm for p p. which metricizes weak convergence.
To obtain our main stability estimate, we first obtain the time evolution of the KL divergence.
Lemma (Time Evolution of the KL divergence). Given the hypothesis velocity field f (1, e CLa Assume that (p') tofo, and (pe) te[o.T] are classical solutions to equation 1 and equation O respectively. It holds that (recall the definition of òe in equation
dt Sx pl log
P. = SSx A..lVlogg log Pr
p dip2 Sx S.ot K (ot pe) v log
PL. Pe
oe Vlog
of Pt.
Pt.
where X' is the tours ITd. All the integrands are evaluated at xc.
We refer the proof of this lemma and all other lemmas and theorems in this section to the appendix We remark that to have the existence of classical solution (o,).. we definitely need the regularity assumptions on V and on K. But the linear term - VV will not contribute to the evolution of the relative entropy. See (Jabin and Wang 2018 for detailed discussions.
Similarly, we have the time evolution of the modulated energy as follows.
Lemma 2 (Time evolution of the modulated energy). Under the same assumptions as in Lemmal 7 given the diffusion coefficient v 2 0, it holds that (recall the definition of ôn in equation (7)
d F(P(.P)) - Sx llKK (ol -2))1 - S- p6 5, K (of p.)+ +v Sx pl K * (ol - -c)) v log di
P
Pt
- SSx K(=- v)
(ATa.)a) ((
)a(pk (x..)
where we recall that the operator A is defined in equation
By Lemma and careful analysis, in particular by rewriting the Biot-Savart law in the divergence of a bounded matrix-valued function ( we obtain the following estimate for the 2D NSE.
Theorem 2 (Stability estimate of the 2D NSE). Notice that when K is the Biot-Savart kernel, div = O. Assume that the initial data pa e (IIII)) and there exists e > such that S Po $ c. Assume further the hypothesis velocity field ft, E Ct Then it holds that
pf dx R()),
e
sup S pl log te[0.T]
Pt
where C= 0000 M(fddt with M(t) def
lo- 2/ pel

---

We remark that given Po is smooth enough and fully supported on X' one can propagate the regu- larity to finally show the finiteness of C. See detailed computations as in Guillim al. [02T]. We give the complete proof in the appendix F. This theorem tells us that as long as RJ) is small, the
KL divergence between p! and P: is small and the control is uniform in time [O. for any T.
Moreover, we highlight that C' is independent of T. and our result on the NSE is significantly better than the average-in-time and exponential-in-12 results from (De Ryck al. 2026
To treat the MVE 2 with Coulomb interactions, we exploit the time evolution of the modulated free energy E(p', .t)) Indeed, combining Lemma | and Lemma
2
we arrive at the following identity.
Lemma (Time evolution of the modulated free energy). Under the same assumptions as in Lemma one has (recall the definitions of ó₁ and A in (77 and
respectively)
d E(p( = - fxx pl I||K<* (ol - -p)) - vV log
2 Pt
- Sxpt 6. (K* (p.- - vV log
Pk Pt.
dt
Pr.
3/2 Jxp K(x v)
1
(1)--)) - A (((-))-(A())()))) )a(of= o.) 2)))((r)))
2
Inspired by the mean-field convergence results in Serfaty 2020 and Bresch al. 2019b], we finally can control the growth of E(p! pe) in the case when v and F(pl Pr) in the case when v= = O. Note also that E(p! pt) can also control the KL divergence when v > 0.
Theorem 3 (Stability estimate of MVE with Coulomb interactions). Assume for [0. TT, the underlying velocity field (x) is Lipschitz in x and supre(o.r] A.pppll))=== = C1 < co. Then there exists C > such that
sup v WLL (p). s sup E(rE Pc) s Exx((CC))))) te[u,]]
De[U.]]
In the deterministic case when under the same assumptions, it holds that sup [e(u.]]
F(P(.p)) s (x)((CC.TR()).
Given that X = II and Po is smooth enough and bounded from below, one can propagate regularity to obtain the Lipschitz condition for A[p.). See the proof and the discussion on the Lipschitz assumptions on (.) in the appendix E
Recall the definition of the operator A in equation 5
Approximation Error of Neural Network Theorems 2 and 3 provide the error estimation guar- antee for the proposed EINN loss (8. Suppose that we parameterize the velocity field f fe with an NN parameterized by 6, as we did in Section 2.1 and let f be the output of an optimization pro- cedure when R(fe) is used as objective. In order the explicitly quantify the mismatch between pf and Pr. we need to quantify two errors: (i) Approximation error, reflecting how well the ground truth solution can be approximated among the NN function class of choice; (ii) Optimization error, in- volving minimization of a highly nonlinear non-convex objective. In the following, we show that for function class F with sufficient capacity, there exists least one element f F that can reduce the loss function R(}) as much as desired. We will not discuss how identify such element in the function class F as it is independent of our research and remains possibly the largest open problem in modern Al research. To establish our result, we make the following assumptions.
Assumption 1. Pu is sufficiently regular, such that v log E L (X) and Jt A W2, VV is Lipschitz continuous. Here W2 (X) stands for the Sobolev norm of order (2, 00) over X.
We here again need to propagate the regularity for Je at least for a time interval [0. It is easy to do so for the torus case, but for the unbounded domain, there are some technical issues to be overcome. Similar assumptions are also needed in some mathematical works for instance in Dabin and Wang (2018. We also make the following assumption on the capacity of the function class F. which is satisfied for example by NNs with tanh activation function (De Ryck et al. 2021].
Assumption 2. The function class is sufficiently large, such that there exists f E F satisfying f e Ca (X) and l Jel (illw==x)) (x) S & for all e [0, T).
Theorem 4. Consider the case where the domain is the torus. Suppose that Assumptions and
hold. For both the Coulomb and the Biot-Savart cases, there exists J e F such that R <
C(T) (e. In 1/ where C(T) is some constant independent of e. Here R is the EINN loss 8
7

---

The major difficulty to overcome is the lack of Lipschitz continuity due to the singular interaction. We successfully address this challenge by establishing that the contribution of the singular region (||=] S to R(f) f) can be bounded by O((E log ))²). Please see the detailed proof in Appendix
4 Related Works on NN-based PDE solvers
Due to space limitations, a detailed discussion about the classical PDE solvers is deferred to Appendix In this section, we focus on the NN-based approaches they are more related to our research.
Solving PDEs is key aspect of scientific research, with wealth of literature [Evans)
As previously mentioned, PINN is possibly the most well-known method of this type. PINN regards the solution to PDE system as the root of the corresponding operators (Di Y- and expresses the time and space boundary conditions as (g)) = O, where g is a candidate solution and D; and B operators acting on g. Parameterizing g ge using NN, PINN optimizes its parameters 0 by
minimizing the residual L()) def Σx= (. + ())))) The hyperparameters X; balance the validity of PDEs and boundary conditions under consideration and must be adjusted for optimal performance. In contrast, E EINN requires no hyperparameter tuning. PINN is versatile and can be applied to a wide range of PDEs, but its performance may not be as good as other NN- based solvers tailored for a particular class of PDEs, as it does not take into account other in-depth properties of the system, phenomenon observed in the literature (Krishnapriyan et al. 2021 Wang et al. 2022]. [Shin and Em Karniadakis 2020 initiates the work of theoretically establishing the consistency of PINN by considering the linear elliptic and parabolic PDEs, for which they prove that a vanishing PINN loss L(0) asymptotically implics 9o recovers the true solution. A similar result is extended to the linear advection equations in (Shin et al. 2020]. . Leveraging the stability of the
operators D; (corresponding to PDEs of interest), non-asymptotic error estimations are established for linear Kolmogorov equations in (De Ryck and Mishra 2022 for semi-linear and quasi-linear parabolic equations and the incompressible Buler in [Mishra and Molinaro] 2022]. and for the NSE in (De Ryck et al. 2023]. We highlight these non-asymptotic results are all average-in-time, meaning that even when the PINN loss is small the deviation of the candidate solution to the true solution may be sigmificant at a particular timestamp E [0, T). In comparison, our results are uniform- in-time, i.c. the supremum of the deviation is strictly bounded by the EINN loss. Moreover, we show in Theorem 2 for the NSE our error estimation holds for any T uniformly, while the results in De Ryck et al. 2023) have an exponential dependence on T.
proposes the Random Deep Vortex Network (RDVN) method for solving the 2D NSE and achieves SOTA performance for this task. Let ue be an estimation of the interaction term K * P in the SDE and use pi to denote the law of the particle driven by the SDE dX,, ui )dt 22dBB To train uf. RDVN minimizes the loss L(0) = fe fx (x) Note that in order simulate the SDE, one needs to discretize the time variable
Recent work from Zhang al. 2022
K* (2) in loss function L. After training RDVN outputs of as solution. However, no error estimation guarantee is provided that controls the discrepancy between Pe and Pt using L(O).
Shen et al. 2022 propose the concept of self-consistency for the FPE. However, unlike our work where the EINN loss is derived via the stability analysis, they construct the potential R(f) for the hypothesis velocity field f by observing that the underlying velocity field f" is the fixed point of some coosisty-nnnsitent transformation A and they construct R(f) to be more complicated Sobolev norm of the residual f-) In their result, they bound the Wasserstein distance between of and p by R(F), which is weaker than our KL type control. The improved KL type control for the Fokker-Plann equation has also been discussed in [offi and Vanden- Eijnden| 2023 A very recent work [Li et al. 2023) extends the celf-consis approach to compute the general Wasserstein gradient flow numerically, without providing further theoretical justification.
5 Experiments
To show the efficacy and efficiency of the proposed approach, we conduct numerical studies on example problems that admit explicit solutions and compare the results with SOTA NN-based PDE solvers. The included baselines are PINN [Raissi et al. 2019 and DRVN Zhang et al. [022]. Note that these baselines only considered the 2D NSE. We extend them to solve the MVE with the Coulomb interaction for comparison, and the details are discussed in Appendix C.T
8

---

PM04 Rebul et 20131 that mestiae
average error average error
rela. 12 error at time
objective value
di
iterations (Lamb Oseen)
itterations (Lamb Oseen)
iterations (Lamb Oseen)
rala. 12 error at time
objective value
2010
25 25
Iterations (Barenblatt)
iterations (Barenblatt)
Iterations (Barenblatt)
Figure 1: The first row contains results for the 2D NSE and the second row contains the results for the 3D MVE with Coulomb interaction. The first column reports the objective losses, while the second and third columns report the average and last-time-stamy relative &g error.
Equations with an explicit solution We consider the following two instances that admit explicit solutions. We verify these solutions in Appendix C.2
Lamb-Oseen Vortex (2D NSE) [Oseen 1912]: Consider the whole domain case where X' IR² and the Biot-Savart kernel Let (H, Σ) be the Gaussian distribution with mean μa and covariance
Σ. po =N 22) for some to O, then we have (r(ar) -N(0, /22(L + 12)
Barenblau solutions (MVE) (Serfaty and Väzquez, 2014 : Consider the 3D MVE with the Coulomb interaction kernel 3 with the diffusion coefficient set to zero, i.e. d = 3 and v = O. Let Uniform[2 A]
be the uniform distribution over a set A. Consider the whole domain case where X' RR If po
Uniform[ x S G 1/33 for some tu 2 O, then we have P = Uniform[ x=|| s C (t+ 1/31
Numerical results We present the results of our experiments in Figure I where the first row contains the result for the Lamb-Oseen vortex (2D NSE) and the second row contains the result for the Barenblatt model (3D MVE). The explicit solutions of these models allow to assess the quality of the outputs of the included methods. Specifically, given hypothesis solution pf the ground truth Pe and the interaction kernel K, define the relative & error timestamp as Q(t) def Sa K (pf (=)())) (x()) |dæ, where Ω is some domain where Pt has non-zero density. We are particularly interested in the quality of the convolution term pf since it has physical meanings. In the Biot-Savart kernel case, it the velocity of the fluid, while in the Coulomb case, is the Coulomb field. We set Q to be [-2, 2] for the Lamb-Oseen vortex and to [-0.1, 0.1 3 for the Barenblaut model. For both models, we take v 0.1, to 0.1, and T 1. The neural network that we use is an MLP with 7 hidden layers, each of which has 20 neurons.
From the first column of Figure D we see that the objective loss of all methods has substantially reduced over a training period of 10000 iterations. This excludes the possibility that a baseline has worse performance because the NN is not well-trained, and hence the quality of the solution now solely depends on the efficacy of the method. From the second and third columns, we see that the proposed EINN method significantly outperforms the other two methods in terms of the time- average relative l2 error, i.e. T Jo (!)) and the relative f2 error at the last time stamp Q(T). This shows the advantage of our method.
By employing entropy dissipation of the MVE, we design potential function for hypothesis velocity field such that it controls the KL divergence between the corresponding hypoth- esis solution and the ground truth, for any time stamp within the period of interest. Built on this potential, we proposed the EINN method for MVEs with NN and derived the detailed computa- tion method of the stochastic gradient, using the classic adjoint method. Through empirical studies on examples of the 2D NSE and the 3D MVE with Coulomb interactions, we show the significant advantage of the proposed method, when compared with two SOTA NN based MVE solvers.
Conclusion

---

Acknowledgments and Disclosure of Funding
Zhenfu Wang is supported by the National Key R&D Program of China, Project Number 1021YFA0002000, NSFC grant No.12111009, Young Elite Scientist Sponsorship Program by China Association for Science and Technology (CAST) No. YESS20200028 and the Fundamental Re- search Funds for the Central Universities (the start-up fund), Peking University. Zebang Shen's work is supported by ETH research grant and Swiss National Science Foundation (SNSF) Project Funding No. 200021-207343.
References
N. M. Boffi and E. Vanden Wanden-Eijnden. Probability flow solution of the fokker- Baaker-planck equation. Ma- chine Learning: Science and Technology, 4(3):33012,, 2023.
D. Bresch, P.-E. Jabin, and Z. Wang. On mean-field limits and quantitative estimates with large class of singular kernels: application to the patlak-keller-seee model. Comptes Rendus Mathe- matique, 357(9))708-722) 2019a.
D. Bresch, P.-E. Jabin, and Z. Wang. Modulated free energy and mean field limit, Séminaire Laurent Schwartz- DDP et applications, pages 1-22, 2019b.
S. Cai, Z. Mao, Z. Wang, M. Yin, and G. E. Karniadakis. Physics-informed neural networks (pinns) for fluid mechanies: review. Acta Mechanica Sinica, 37(12):1727-1738, 2021.
J. A. Carrillo, K. Craig, and F. S. Patacchini. blob method for diffusion. Calculus of Variations and Partial Differential Equations, 58:1 1-53, 2019.
J. A. Carrillo, K. Craig, L. Wang, and C. Wei. Primal dual methods for wasserstein gradient flows. Foundations of Computational Mathematics, pages 1-55, 2022.
Cuomo, V. Di Cola, F. Giampaolo, Rozza, M. Raissi, and F. Piccialli. Scientific machine learning through physics-informee informed neural networks: where we are and what's next. Journal of Scientific Computing, 92(3):88, 2022.
T. De Ryck and Mishra. Error analysis for physics-informee informed neural networks (pinns) approximat- ing kolmogorov pdes. Advances in Computational Mathematics, 48(6): 400 2022.
T. De Ryck, S. Lanthaler, and S. Mishra. On the approximation of functions by tanh neural net- works. Neural Nerworks, 143:732 750, 2021. ISSN 0893-6080. doi: http:://do..rrg//00000///
neunet.2021.00 0.002..08..055 URL https: www. wtwpeeeeeeeeeeeeeeeeeeeeeeeeeeee wwwpeeettt.ooccc.eeccecceee
p11/808936082210033
De Ryck, A. Jagtap, and S. Mishra. Error estimates for physics-informed neural networks ap- proximating the aavier-stokes equations. IMA Journal of Numerical Analysis, 2023.
L. Erdōs and H.-T. Yau. A dynamical approach to random matrix theory, volume 28. American Mathematical Soc., 2017.
L. Evans. Partial differential equations, volume 19. American Mathematical Society, 2022.
Golse. On the dynamics of large particle systems in the mean field limit. Macroscopic and large scale phenomena: coarse graining, mean field limits and ergodicity, pages 1-144, 2016.
A. Guillin, P. L.. Bris, and P. Monmarché. Uniform in time propagation of chaos for the 2d vortex model and other singular stochastic systems. arXiv preprint arXiv:2108. 08675, 2021.
G. Gupta, X. Xiao, and Bogdan. Multiwavelet-based operator learning for differential equations. Advances in neural information processing systems, 14:24048-24065 2021.
J. Han, A. Jentzen, and W. E. Solving high-dimensional partial differential equations using deep learning. Proceedings of the National Academy of Sciences, 115(34):855-88500, 2018
P.-E. Jabin and Z., Wang. Quantitative estimates of propagation of chaos for stochastic systems with W- kernels. Inventiones mathematicae, 214(1):523-591, 2018.
10

---

C. Johnson. Numerical solution partial differential equations by the finite element method. Courier Corporation, 2012.
R. Jordan, D. Kinderlehrer, and F. Otto. The variational formulation of the fokker-planck equation. SIAM journal on mathematical analysis, 29(1):1-17, 1998.
G. E. Karniadakis, I. G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, and L. Yang. Physics-informed machine learning. Nature Reviews Physics, 3(6):422 440, 2021.
N. Kovachki, Z. Li, B. Liu, K. Azizzadenesheli, K. Bhattacharya, A. Stuart, and A. Anandkumar. Neural operator: Learning maps between function spaces. arXiv preprint arXiv:2 100.0881,, 2021.
A. Krishnapriyan, A. Gholami, S. Zhe, R. Kirby, and M. W. Mahoney. Characterizing possible failure modes in physics-informet neural networks. Advances in Neural Information Processing Systems, 34:26548-26560, 2021.
L. Li, S. Hurault, and Solomon. Self-consistent velocity matching of probability flows. arXiv preprint arXiv:2301.1372 2023.
Z. Li, Kovachki, K. Arizzadenesheli, B. Liu, Bhattacharya, Stuart, and Anand- kumar. Fourier neural operator for parametric partial differential equations. arXiv preprint arXiv:2010.08895, 2020a.
Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, and A. Anandku- mar. Neural operator: Graph kernel network for partial differential equations. arXiv preprint arXiv:2003.03485 2020b.
D.-G. Long. Convergence of the random vortex method two dimensions. Journal of the American Mathematical Society, 1(4):779-88) 804, 1988.
A. Majda, A. L. Bertozzi, and A. Ogawa. Vorticity and incompressible flow. cambridge texts in applied mathematics. Appl. Mech. Rev., 55(4):B77-B78, 2002.
S. Mishra and R. Molinaro. Estimates on the generalization error of physics-informed neural net- works for approximating pdes. IMA Journal Numerical Analysis, 43(1):1 01 2022, ISSN 0272-4979. doi: 10. 10.0099//manum//rab0933. URL /tttpp////...orr111........ /mm.nmmmmmmmm///////m drab093
F. Moukalled, L. Mangani, M. Darwish, F. Moukalled, L. Mangani, and M. Darwish. The finite volume method. Springer, 2016.
C. Oseen. Uber die wirbelbewegung in einer reibenden flussigkeit. Ark. Mat. Astro. Fys., 7, 1912.
M. Raissi, P. Perdikaris, and G. E. Karniadakis. pyvsics-informed neural networks: A deep learn- ing framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational physics, 378:686-707, 2019.
M. Raissi, A. Yazdani, and G. E. Karniadakis. Hidden fluid mechanics: Learning velocity and pressure fields from flow visualizations. Science, 367(6481) (0066488))10026-441 1026-1030, 2020.
M. Rosenzweig. Mean-field convergence of point vortices the incompressible euler equation with vorticity _ Archive for Rational Mechanics and Analysis, 243(3):1361-1431 2022.
Serfaty. Mean field limit for coulomb-type flows. Duke Math. J., 169(15):2887-2935 2020.
S. Serfaty and J. L. Vázquez. A mean field equation as limit of nonlinear diffusions with fractional laplacian operators. Calculus of Variations and Purtial Differential Equations, 49:1091 120, 2014.
Shen, T. Tang, and I..L.. Wang. Spectral methods: algorithms, analysis and applications, vol-
ume 41. Springer Science & Business Media, 2011.

---

Z. Shen, Z. Wang, S. Kale, A. Ribeiro, A. Karbasi, and H. Hassani. Self seff-consistency of the fokker planck equation. In P.-L. Loh and M. Raginsky, editors, Proceedings of Thirty Fifth Conference on Learning Theory, volume 178 of Proceedings of Machine Learning Research, pages 817-841. PMIR, 02-05 Jul 2022. URL https::/proceed wgg.gg........ wwggsmml...... htt.://g......ees////////e//eeee
html
J. Shin, YeonjongDarbon and G. Em Karniadakis. On the convergence of physics informed neural networks for linear second-order elliptic and parabolic type pdes. Communications in Computa- tional Physics, (8(5):20422207)) 2020.
Y. Shin, Z. Zhang, and G. E. Karniadakis. Error estimates of residual minimization using neural networks for linear pdes. arXiv preprint arXiv:2010.08019 2020.
G. D. Smith, G. D. Smith, and G. D. S. Smith. Numerical solution of partial differential equations: finite difference methods. Oxford university press, 1985.
C. Villani et al. Optimal transport: old and new, volume 338. Springer, 2009.
C. Wang. S. Li, D. He, and L. Wang. Is L physics informed loss always suitable for training physics informed neural network? In Advances in Neural Information Processing Systems, 2022.
E. Weinan, Han, and A. Jentzen. Algorithms for solving high dimensional pdes: from nonlinear monte carlo to machine learning. Nonlinearity, 35(1):278, 2021.
X. Xiao, D. Cao, R. Yang, G. Gupta, G. Liu, C. Yin, R. Balan, and P. Bogdan. Coupled mul- tiwavelet neural operator learning for coupled partial differential equations. arXiv preprint arXiv:2303.02304, 2023.
L. Zhang, J. Han, H. Wang, and R. Car. Deep potential molecular dynamics: a scalable model with the accuracy of quantum mechanics. Physical review letters, 120(14);14300) 2018.
R. Zhang, P. Hu, Q Meng, Y. Wang. R. Zhu, B. Chen, Z.-M. Ma, and T.-Y. Liu. Drvn (deep random vortex network): A new physics-informed machine learning method for simulating and inferring incompressible fluid flows. Physics of Fluids, 34(10): 107112 2022.
12


Reviews:
[{"review_id": "q4hAvbNNPl", "reviewer": "Reviewer_AJKM", "summary": "The paper proposed an entropy-dissipation informed neural network for solving certain PDE systems. Using energy dissipation, the paper designed a special structure that propagates information from previous dynamics. They provided a theoretical guarantee for this framework, including the control of the KL divergence between the target and the estimator and a discussion of the approximation error of the NN. They also did some numerical experiments to verify their claims.", "questions": "1. In Figure 1, first row,  at iteration 0, the value of the objective functions are not the same. \n2. In Figure 1, it seems DRVN and EINN perform similarly. Please comment on this.\n3. Please comment on how EINN compares with more recent neural-net-based PDE solvers, such as Fourier neural networks.\n4. I wonder if more experiments can be done, other than 2D NSE and 3D MVE.", "limitations": "yes", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. The idea of using energy dissipation to design a framework is very interesting. EINN also takes the previous dynamics precisely and provides a uniform-in-time solution. \n\n2. Mathematically speaking, the proofs are solid and well-written.\n\n3. McKean-Vlasov type PDEs is an important class of PDEs for people from both math and physics communities. The presence of the singular kernel does bring challenges, which are handled by EINN.", "weaknesses": "1. The authors specifically take advantage of the special structure of this type of PDEs. On one hand, it performs well for those PDEs; on the other hand, this limits the application of this method to other PDEs. \n2. Numerical experiments are done with only comparison with PINN and DRVN.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for carefully reading our paper and for making helpful suggestions. In the following, we address the comments/questions one by one. The responses will be incorporated in our revision.\n\n> The authors specifically take advantage of the special structure of this type of PDEs. On one hand, it performs well for those PDEs; on the other hand, this limits the application of this method to other PDEs.\n\nWe emphasize that this work is mainly theory-oriented and we focus on the special instance of the McKean-Vlasov equation (Eq.2) so that we can establish rigorous error estimations, e.g. Theorems 2 and 3.\n\nIn terms of solving PDEs numerically, we emphasize that the EINN framework can be applied to a very broad class of evolutionary PDEs that satisfy the preservation of mass: For any PDE that can be written as a continuity equation (Eq.6) so that the operator $\\mathcal{A}$ (Eq.5) is known, one can construct the EINN loss (see Eq.8) accordingly since it only requires the knowledge of $\\mathcal{A}$ and $\\rho_0$. To solve these PDEs numerically, one simply uses the standard stochastic optimization methods, for example ADAM, to minimize the EINN loss. An important subclass of PDEs that we can solve is the PDEs that admit the structure of Wasserstein gradient flow. This includes the highly general diffusion–advection–interaction equation (see for example the equation below Eq.4.14 in [Santambrogio, F. {Euclidean, metric, and Wasserstein} gradient flows: an overview. Bull. Math. Sci. 7, 87–154 (2017).](https://doi.org/10.1007/s13373-017-0101-1)). \nHowever, we did not emphasize these numerical examples since we are not able to establish a general error estimation guarantee for all these PDEs, for example our current framework cannot handle the porous medium equation.\n\nOf course, we acknowledge that our method cannot solve numerically all PDEs. One example is the Poisson equation. However, philosophically, producing highly accurate solutions for PDEs, especially with rigorous error estimation guarantees, is a sophisticated task, which needs to be accomplished by exploiting the underlying structures. Solving all PDEs with a versatile approach like the PINN might not lead to the most accurate solution, as demonstrated in our experiment section. \n\n> In Figure 1, first row, at iteration 0, the value of the objective functions are not the same.\n\nSince the objective functions are different for different methods, it is not surprising that they have a different value initially. In fact, we point out that comparing the function values of different methods may not be meaningful since their objectives are of drastically different meanings.\n\nInstead, we emphasize that the purpose of presenting the function values of different methods is to exclude the possibility that the difference of their solution accuracy (the relative l2 error) is due to insufficient training. Following this logic, to interpret the first column of Figure 1, one should focus on the change of the function value of a single method across different training iterations. Since we observe that all the methods are well trained (the function value sufficiently decreases over a long training period), it is reasonable to infer that the difference of their accuracy is due to the design of the method and the superior solution accuracy of EINN implies that it is a better approach. Please see a similar discussion in lines 352-355 of our paper.\n\n> In Figure 1, it seems DRVN and EINN perform similarly. Please comment on this.\n\nBy \"perform similarly\", maybe the reviewer is referring to the objective value. For this question, please check the discussion above. \n\nHowever, if the reviewer is referring to the solution accuracy, the experiments (the second and third columns) in Figure 1 show that EINN has a much better solution accuracy compared to DRVN. For example, in the 2D-NSE task, the average relative l2 error of EINN is an order of magnitude smaller than that of DRVN and the final-time-stamp error of EINN is around 1/5 of that in DRVN. Consequently, EINN performs much better.\n\n> Please comment on how EINN compares with more recent neural-net-based PDE solvers, such as Fourier neural networks.\n\nTo highlight the major difference between EINN and the approach of Fourier Neural Operator, EINN considers a completely different problem setting as neural operator: **EINN requires no pre-existing data** and the goal is obtain the solution to a PDE by solely exploiting the structure of the equation itself. In contrast, **the neural operator approach is data-driven**, i.e. it relies on the existence of configuration-solution pairs. Here, by configuration-solution pairs, we mean the correspondence between some configurations that determine the PDE, e.g. the initial condition or the viscosity parameter in the fluid dynamics problems, and the pre-existing solution to the PDE given the aforementioned configurations. Consequently, the neural operator approach is more like a regression problem where a neural network is trained to learn the abstract map between the configuration and the solution. In contrast, EINN is more like a numerical PDE solver.\n\nConsequently, EINN and the approach of neural operators are two related but quite distinct research directions. They are related in the sense that EINN can provide the data (configuration-solution pairs) required by the neural operator approach. They are distinct since EINN requires no data a priori, while the neural operator approach is built on the supervised learning paradigm.\n\n> I wonder if more experiments can be done, other than 2D NSE and 3D MVE.\n\nWe have conducted an additional experiment on the Porous Medium Equation. While we could not provide rigorous error analysis, the numerical result seems promising. Please see the results in the attached pdf file."}, {"review_id": "D2UfHlpRbu", "reviewer": "Reviewer_HFPv", "summary": "The authors propose a novel method for solving the McKean-Vlasov Equation (MVE), a Partial Differential Equation (PDE) that models the dynamics of stochastic particle systems with mean-field interactions. The new approach, named Entropy-dissipation Informed Neural Network (EINN), leverages neural networks to estimate the velocity field underlying the MVE and minimize a proposed potential function. This method is particularly designed to tackle singular interaction kernels, where traditional methods struggle. The authors provide empirical comparisons with state-of-the-art (SOTA) neural network-based MVE solvers and demonstrate the superior performance of the EINN framework. Despite the advancements, a noticeable weakness in this work is the omission of a discussion on the neural operator framework for solving PDEs, which represents a significant progression in the field. This lack of consideration limits the comprehensiveness and depth of the study.", "questions": "1. While the authors aim to utilize the Lyapunov function to study the deviation of the hypothesis solution from the true solution, they might have overlooked the possibility of overfitting. There seems to be no clear mention of measures taken to avoid overfitting in the model. How does the model ensure that it is generalizable and not just fitting the data it was trained on?\n2. How can the EINN framework be adapted or expanded to deal with other types of PDEs? \n3. How does the performance of EINN scale with the complexity of the problem, particularly in high dimensional cases? \n4. Given the reliance on NNs, how sensitive is the performance of the EINN approach to the choice of NN architecture and parameters?\n5. The authors should refer to and discuss the neural operator work for solving PDEs as related work. Here are some references: [Li, Zongyi, et al. \"Fourier neural operator for parametric partial differential equations.\"]; [Kovachki, Nikola, et al. \"Neural operator: Learning maps between function spaces.\"]; [Li, Zongyi, et al. \"Neural operator: Graph kernel network for partial differential equations.\"]; [Gupta, Gaurav, Xiongye Xiao, and Paul Bogdan. \"Multiwavelet-based operator learning for differential equations.\"]; [Xiao, Xiongye, et al. \"Coupled Multiwavelet Neural Operator Learning for Coupled Partial Differential Equations.\"].", "limitations": "The authors have specifically designed this approach for solving MVEs, which might limit its general applicability for solving various types of PDEs. ", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. The paper provides a novel approach to solving MVEs. The idea of using entropy dissipation in the underlying system is inventive and beneficial to the community. \n2. The research is solid and thorough. The authors have established theoretical guarantees for their EINN framework and provided clear derivations and proofs. \n3. The empirical comparison of EINN with state-of-the-art methods and the positive results strengthen the claim of the proposed method's effectiveness. The authors ensured a fair comparison by using NNs with the same complexity for all methods tested.", "weaknesses": "1. The lack of discussion on the neural operator framework for solving PDEs. This approach has been a considerable advancement in the field, and the authors' failure to address or refer to it is an oversight. Further exploration or comparison of the neural operator approach could enhance the comprehensiveness and relevance of the paper.\n2. While the authors did a commendable job providing rigorous mathematical proofs, the paper could be difficult for those without a strong background in the relevant mathematical fields. Providing more intuitive explanations or graphical illustrations might have made the methodology more accessible. \n3. The empirical testing seems to be limited in scope. More extensive experimentation, including larger-scale problems and more varied examples, would have strengthened the empirical evidence. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for pointing out interesting related works. In the following, we address the concerns one by one. The following discussion will be incorporated in our revision.\n\n> The lack of discussion on the neural operator framework for solving PDEs. This approach has been a considerable advancement in the field, and the authors' failure to address or refer to it is an oversight. Further exploration or comparison of the neural operator approach could enhance the comprehensiveness and relevance of the paper.\n\nWe thank the reviewer for pointing out this interesting research direction. However, to highlight the major difference between EINN and the approach of the Neural Operator, it's worth noting that they consider **completely different problem settings**: **EINN requires no pre-existing data** and the goal is to obtain the solution to a PDE by solely exploiting the structure of the equation itself. In contrast, **the neural operator approach is data-driven**, i.e. it relies on the existence of configuration-solution pairs. Here, by configuration-solution pairs, we mean the correspondence between some configurations that determine the PDE, e.g. the initial condition or the viscosity parameter in the fluid dynamics problems, and the pre-existing solution to the PDE given the aforementioned configurations. Consequently, the neural operator approach is more like a regression problem where a neural network is trained to learn the abstract map between the configuration and the solution. In contrast, EINN is more like a numerical PDE solver.\n\nConsequently, EINN and the approach of neural operator are two related but quite distinct research directions. They are related in the sense that EINN can provide the data (configuration-solution pairs) required by the neural operator approach. They are distinct since EINN requires no data a priori, while the neural operator approach is built on the supervised learning paradigm.\n\n> While the authors aim to utilize the Lyapunov function to study the deviation of the hypothesis solution from the true solution, they might have overlooked the possibility of overfitting. There seems to be no clear mention of measures taken to avoid overfitting in the model. How does the model ensure that it is generalizable and not just fitting the data it was trained on?\n\nMaybe the reviewer has some misunderstanding of the problem formulation and the proposed approach of our paper. \nAs mentioned in the above response, **EINN is not data-driven, i.e. it requires no data** and the goal of our research to solve PDEs by exploiting only its structure. Hence there is no concern of overfitting. According to our theory, a small EINN loss can directly be translated to a high solution accuracy.\n\n> How can the EINN framework be adapted or expanded to deal with other types of PDEs?\n\nWe emphasize that this work is mainly theory-oriented and we focus on the special instance of the McKean-Vlasov equation (Eq.2) so that we can establish rigorous error estimations, e.g. Theorems 2 and 3.\n\nIn terms of solving PDEs numerically, we emphasize that the EINN framework can be applied to a very broad class of evolutionary PDEs that satisfy the preservation of mass: For any PDE that can be written as a continuity equation (Eq.6) so that the operator $\\mathcal{A}$ (Eq.5) is known, one can construct the EINN loss (see Eq.8) accordingly since it only requires the knowledge of $\\mathcal{A}$ and $\\rho_0$. To solve these PDEs numerically, one simply uses the standard stochastic optimization methods, for example ADAM, to minimize the EINN loss. An important subclass of PDEs that we can solve is the PDEs that admit the structure of Wasserstein gradient flow. This includes the highly general diffusion–advection–interaction equation (see for example the equation below Eq.4.14 in [Santambrogio, F. {Euclidean, metric, and Wasserstein} gradient flows: an overview. Bull. Math. Sci. 7, 87–154 (2017).](https://doi.org/10.1007/s13373-017-0101-1)). \nHowever, we did not emphasize these numerical examples since we are not able to establish a general error estimation guarantee for all these PDEs, for example our current framework cannot handle the porous medium equation.\n\n> How does the performance of EINN scale with the complexity of the problem, particularly in high dimensional cases?\n\nThis is an interesting problem and we intend to test it out in our future work. Currently, this work is mainly theory-oriented and the goal is to establish high quality error estimation for PDEs that involve singular terms. \n\n> Given the reliance on NNs, how sensitive is the performance of the EINN approach to the choice of NN architecture and parameters?\n\nIn our experience, EINN is pretty robust to the NN architecture. Some preliminary results (where we vary the width of the NN) are provided in the PDF file. However, for more complicated tasks, e.g. when the solution has worse regularity, we would expect an NN with stronger representation capability would lead to better performance.\n\n> The authors should refer to and discuss the neural operator work for solving PDEs as related work...\n\nWe thank the reviewer for pointing out these interesting works and will discuss their relation with EINN in our revision. However, we reemphasize that EINN and the neural operator approach consider completely different settings."}, {"review_id": "TrmL2DxOtt", "reviewer": "Reviewer_Jhg5", "summary": "This paper proposes a new method for learning solutions to partial differential equations of the McKean-Vlasov type, and presents a theoretical analysis of the proposed method. In particular, for a certain class of equations, it is shown that the error bound does not increase exponentially in time. This is a significant result that takes advantage of the characteristics of the equations.", "questions": "(1) What are the practical applications of the McKean-Vlasov-type partial differential equations?\n(2) Because computationally expensive methods such as the adjoint method and the Monte Carlo method are required, the computational cost of the proposed method may be larger than the classical approaches such as the finite difference method. So, what are the advantages of the proposed neural-network solver compared to the classical numerical methods (e.g, the finite difference method)?", "limitations": "No potential negative societal impact is expected.", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "The uniform bound presented by the authors is a very strong theoretical result. Although I do not understand the details of the proofs, the proofs seem to be reliable because several references are provided. In my opinion, this paper is an excellent paper as a paper in the research area of numerical analysis.", "weaknesses": "The most important weakness is that it is not clear what practical applications the McKean-Vlasov type of partial differential equations have.  Unlike numerical analysis, which is a research area in mathematics, there are many people who are interested in applications in the machine learning community. So these results may not be of much interest to them. \n\nAnother significant weakness is that the adjoint method and the Monte Carlo method are required in the proposed approach, which are computationally expensive. Considering the computational cost, the proposed method is not considered superior to classical numerical methods, such as the finite difference method.\n\nIn addition, the uniform bound seems to hold only for a class of equations, but this is not clearly written in Introduction. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for pointing out the need to highlight the application of our research and to compare with classical method. In the following we address the concerns one by one. Our response will be incorporated in our revision.\n\n> What are the practical applications of the McKean-Vlasov-type partial differential equations?\n\nMcKean-Vlasov equation describes the evolution of the law of a large interacting particle system. It plays a crucial role in various scientific fields, including physics, control, and social sciences. Just to list a few example applications and how EINN can be utilized accordingly:\n\n- [physics] In the research of plasma physics, MVE with column interaction can be used to model the dynamics of charged electrons in the over-damped limit. In fluid dynamics, the 2D-NSE is a simplification of the 3D Navier-Stokes equation, which is one of the Millennium Prize problems.\n- [control] The mean-field game is formulated as a system of the Hamilton–Jacobi–Bellman equation and the McKean-Vlasov equation. Solving MVE accurately is a crucial part of this research direction. When the interaction between agents is bounded, EINN can directly be used to derive error estimation for solving the MVE part. Please see Theorem 6 in the appendix. \n- [social science] The flocking model was originally designed to model the flocking behavior of birds but now finds application in social sciences, such as in the study of opinion dynamics. When the flocking interaction is bounded, applying the EINN method directly leads to solutions with error estimations.\n\n> What are the advantages of the proposed neural-network solver compared to the classical numerical methods (e.g, the finite difference method)?\n\nIn the following, we list some key advantages of EINN compared to the classical methods. Since there are various existing methods, we specifically focus on the finite difference method as pointed out by the reviewer.\n\n- One of the major advantages of EINN and in general NN based PDE solvers is that they do not suffer from the curse of dimensionality (at least not explicitly). Specifically, the finite different method requires to discritize the space so that the solution to the PDE can be represented using a finite dimensional vector. Consequently, the number of grid points grows exponentially with the ambient dimension. For NN based PDE solver, the solution is represented via an NN, which enjoys a significantly stronger representation capability. Mathematically speaking, to approximate a function with sufficiently good regularity, NN is a universal function approximator with width and depth polynomially grows with the ambient dimension. Please also see Appendix A of our paper where we discussed the relation of the classical methods with NN based approaches.\n\n- Further, to maintain the same level of granularity, the number of grid points also grows quickly with the size of the domain. Consequently, if the problem is defined in an unbounded domain, the implementation of the finite difference method requires domain cut-off to ensure the boundedness. In contrast, EINN can directly handle the unbounded cases.\n\n- Another important advantage of EINN over the finite difference method is that, once trained, EINN allows one to access value and gradient of the solution at any spatial coordinate while the finite difference method only provide information on the grid points.\n\n- Recall that at any given time stamp, the solution to MVE is a probability distribution. With EINN, one can sample from the solution at any given time, and our KL-type error estimation ensures that this sampling is meaningful. In contrast, even when the positivity and the preservation of mass is ensured in the finite difference method, it would not generate sample beyond the grid points.\n\n> In addition, the uniform bound seems to hold only for a class of equations, but this is not clearly written in Introduction.\n\nWe are not sure what the reviewer means by \"this is not clearly written in Introduction\", could you be more specific? If the reviewer is referring to the limitations on the interaction kernel, we have clearly stated in Theorem 1 of the introduction that our theoretical guarantee holds for bounded interaction kernel or the singular interaction kernel in the Coulomb or Biot-Savart cases (please see lines 82-83). \n\n\n\n\n\n"}, {"review_id": "mQ72XbFXK0", "reviewer": "Reviewer_A2vP", "summary": "This work focuses on deriving a uniform-in-time bound to the KL divergence between the NN-derived solution $\\rho_t^f$ and the true solution $\\bar{\\rho}_t$, where $\\rho_t^f$ is obtained via a push forward of samples of $\\bar{\\rho}_0$ with velocity parameterized by the neural network $f$. The neural network $f$ is trained to match the operator $\\mathcal{A}[\\rho]$ via the the Entropy-dissipation Informed Neural Network (EINN) loss given as an expectation over $\\bar{\\rho}_0$ of the integral from time $0$ to $T$ of the difference between $f$ and $\\mathcal{A}[\\rho^f]$. The method is then demonstrated to perform better than SOTA on McKean-Vlasov equations with singular interaction kernels: Coulomb interaction and Biot-Savart interaction.", "questions": "1. Maybe I am misunderstanding, how is using a NN to parameterize $\\mathcal{A}[\\rho^f]$ helping in computation? $\\rho_t$ is still computed via push forward of $\\rho_0$. Would using the true $\\mathcal{A}[\\rho]$ be better? Or is the goal to first compute $f$ as a surrogate to $\\mathcal{A}[\\rho]$ then compute $\\rho^f$ from the PDE with $f$ as a surrogate for $\\mathcal{A}[\\rho]$?\n2. Why is $K_c$ set to $0$ instead of clipping some large value (line 173, equation (19))?\n", "limitations": "1. Addressed in Section 3, approximation error of neural network that if a neural network can represent the true $\\mathcal{A}[\\rho]$, EINN loss is zero. Perhaps reference some papers on neural networks designed to represent $f$ for MVE. ", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. Uniform-in-time bound to the KL divergence for finite duration $T$.\n2. Method was shown to work on singular interaction kernels: Coulomb interaction,\\ and Biot-Savart interaction.\n3. Singularity for Biot-Savart kernel was removed in gradient computation during training. ", "weaknesses": "1. Addressed in the last paragraph of section 4, the similarity in methodology to Shen et al. [2022]. Though the implementation is similar, the theoretical contribution is different. Perhaps include a comparison of their result in Figure 1 with and without generalization to MVE. And perhaps discuss the generalization to MVE.\n2. Also addressed in the second last paragraph of section 4, the similar in loss function to Random Deep Vortex Network (RDVN) of Zhang et al. [2022]. Though the losses are similar, the authors mention a lack of discrepancy control in Zhang et al. [2022]. Perhaps provide more discussion.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for making valuable comments. In the following, we address them one by one. Our response will be incorporated in our revision. \n\n> Addressed in the last paragraph of section 4, the similarity in methodology to Shen et al. [2022] ... And perhaps discuss the generalization to MVE.\n\nAs acknowledged by the reviewer, EINN can handle the McKean-Vlasov equation which is more general than the Fokker-Planck equation considered in Shen et al. [2022]. For the purpose of comparison, let us now focus on the Fokker-Planck equation. In this special instance, the difference between EINN and Shen et al. [2022] is that the former uses functional $\\mathcal{L}^2$ norm of the perturbation $\\delta^f$ while the latter requires the Sobolev norm, making their method much more complicated to implement. We believe the additional terms of the higher order derivatives in the Sobolev norm are artifacts for establishing the theoretical guarantee therein. In fact, with our theoretical analysis, we can directly recover a stronger guarantee (convergence in terms of KL-divergence) than that of Shen et al. [2022] (convergence in terms of Wasserstein-2 distance), by setting $K=0$ (please see Theorem 6 in the appendix). \n\nIn this regard, EINN can be regarded as an improvement (and simplification) of Shen et al. [2022] for the Fokker-Planck equation and a generalization of their work to the MVE case. \n\nIn terms of directly generalizing the Sobolev norm loss of Shen et al. [2022] to the MVE, we are not sure how the higher order gradient of the perturbation $\\delta^f$ can be computed (or even defined) for the (possibly singular) interaction term. To be specific, in the instance of MVE with Coulomb interaction, we have that $K(x) = c \\frac{x}{\\|x\\|^d}$ for $d\\geq 2$ where $c$ is some constant and $d$ is the ambient dimension. The Sobolev norm of $\\delta^f$ involves the gradient of $K(x)$, which contains terms of the order $ \\frac{1}{\\|x\\|^d}$. This term is no longer locally integrable, i.e. for any $\\epsilon > 0$, $\\int_{\\|x\\|\\leq\\epsilon} \\frac{1}{\\|x\\|^d} \\mathrm{d} x = \\infty$. We will discuss this in our revision.\n\n> Also addressed in the second last paragraph of section 4, the similar in loss function to Random Deep Vortex Network (RDVN) of Zhang et al. [2022] ... Perhaps provide more discussion.\n\nWhile both RDVN and EINN share a similar form of minimizing functional $\\mathcal{L}^2$ norm of some perturbation function, we emphasize that perturbation functions therein are of different meanings: The perturbation function of RDVN, represented by $u_t^\\theta - K \\ast \\rho_t^\\theta$, captures the error in estimating the interaction term in the 2D-NSE and the diffusion term ($\\Delta \\rho$ in the MVE in Eq.2) is incorporated in the SDE (see line 313). In contrast, the perturbation function of EINN (defined in Eq.11) is derived from the sensitive analysis of some Lyapunov function and should be understood as the error of estimating the **entire velocity field** of the 2D-NSE. Importantly, the diffusion term is also included in the EINN loss. Consequently, the EINN loss captures the **full information** of the velocity of the 2D-NSE while the RDVN only captures a part of it (the interaction part). We believe this is the reason why the authors of RDVN did not provide rigorous solution recover guarantee in the work Zhang et al. [2022].\n\n> Maybe I am misunderstanding how is using a NN to parameterize $\\mathcal{A}[\\rho^f]$ helping in computation? Would using the true  $\\mathcal{A}[\\rho]$ be better?\n\nThe EINN objective (Eq.8) is derived from the sensitivity analysis, as discussed in section 2 of the paper. One of the practical benefits of using $\\mathcal{A}[\\rho^f]$ in the EINN objective is that $\\mathcal{A}[\\rho^f]\\circ X_t$ can be computed using **only the initial distribution $\\rho_0$ and the hypothesis velocity field $f$** (please see a detailed discussion in section 2.1 in our paper). It is possible that replacing $\\mathcal{A}[\\rho^f]$ with the true velocity field $\\mathcal{A}[\\rho]$ in the EINN loss would lead to a similar theoretical guarantee, but since $\\mathcal{A}[\\rho]$ is **unknown** (in fact this is the quantity we want to learn in the first place), the objective function after this replacement does not yield a practical algorithm.\n\n> Or is the goal to first compute $f$ as a surrogate to $\\mathcal{A}[\\rho]$ then compute $\\rho^f$ from the PDE with $f$ as a surrogate for $\\mathcal{A}[\\rho]$?\n\nYes, the reviewer is correct.\n\n> Why is $K_c$ set to $0$ instead of clipping some large value (line 173, equation (19))?\n\nWe pick $K_c = 0$ for the sake of simplicity in our discussion, a choice that is also common practice in the literature on the mean-field limit. However, we also point out that the analysis below line 173 remains the same if we pick $K_c$ to be any function with the level of singularity less than $\\frac{1}{\\|x\\|^{d-1}}$. Here by \"same\", we mean that the dependence of the cut-off error on the threshold $c$ remains $O(c)$ in Eq.20, when the cut-off threshold is small ($c\\rightarrow 0$). Consequently, one can try larger values for $K_c$ with $\\|x\\|\\leq c$, which might leads to better practical performance.\n\n> Addressed in Section 3, approximation error of neural network that if a neural network can represent the true $\\mathcal{A}[\\rho]$, EINN loss is zero. Perhaps reference some papers on neural networks designed to represent $f$ for MVE.\n\nOne of the neural network structures that meets the requirement in Assumption 2 is [T. De Ryck, S. Lanthaler, and S. Mishra. On the approximation of functions by tanh neural net- works. Neural Networks, 143:732–750, 2021](https://www.sciencedirect.com/science/article/pii/S0893608021003208), as discussed in line 274 of our paper."}]

Author rebuttal:
This pdf file contains the additional experiments on the Porous Medium Equation and a preliminary study of the influence of the width of the NN on EINN's performance.

Meta-review: The reviewers unanimously find the paper interesting and result worth publication. The area chair agrees with the evaluation after reading the discussion and the manuscript.
Decision: Accept (poster)