# Learning Mean-Field Equations from Particle Data Using WSINDy

Daniel A. Messenger, David M. Bortz<EMAIL_ADDRESS><EMAIL_ADDRESS>of Applied Mathematics, University of
Colorado Boulder, 11 Engineering Dr., Boulder, CO 80309, USA.

###### Abstract.

We develop a weak-form sparse identification method for interacting particle
systems (IPS) with the primary goals of reducing computational complexity for
large particle number $N$ and offering robustness to either intrinsic or
extrinsic noise. In particular, we use concepts from mean-field theory of IPS
in combination with the weak-form sparse identification of nonlinear dynamics
algorithm (WSINDy) to provide a fast and reliable system identification scheme
for recovering the governing stochastic differential equations for an IPS when
the number of particles per experiment $N$ is on the order of several thousand
and the number of experiments $M$ is less than 100. This is in contrast to
existing work showing that system identification for $N$ less than 100 and $M$
on the order of several thousand is feasible using strong-form methods. We
prove that under some standard regularity assumptions the scheme converges
with rate ${\mathcal{O}}(N^{-1/2})$ in the ordinary least squares setting and
we demonstrate the convergence rate numerically on several systems in one and
two spatial dimensions. Our examples include a canonical problem from
homogenization theory (as a first step towards learning coarse-grained
models), the dynamics of an attractive-repulsive swarm, and the IPS
description of the parabolic-elliptic Keller-Segel model for chemotaxis.

Keywords: data-driven model selection, interacting particle systems, weak
form, mean-field limit, sparse recovery.

##  1\. Problem Statement

Consider a particle system
$\mathbf{X}_{t}=(X^{(1)}_{t},\dots,X^{(N)}_{t})\in\mathbb{R}^{Nd}$ where on
some fixed time window $t\in[0,T]$, each particle
$X^{(i)}_{t}\in\mathbb{R}^{d}$ evolves according to the overdamped dynamics

$dX_{t}^{(i)}=\left(-\nabla K*\mu^{N}_{t}\left(X^{(i)}_{t}\right)-\nabla
V\left(X^{(i)}_{t}\right)\right)dt+\sigma(X_{t}^{(i)})\,dB_{t}^{(i)}$ (1.1)

with initial data $X_{0}^{(i)}$ each drawn independently from some probability
measure $\mu_{0}\in{\mathcal{P}}_{p}(\mathbb{R}^{d})$, where
${\mathcal{P}}_{p}(\mathbb{R}^{d})$ is the space probability measures on
$\mathbb{R}^{d}$ with finite $p$th moment222We define the $p$th moment of a
probability measure $\mu$ for $p>0$ by $\int_{\mathbb{R}^{d}}|x|^{p}d\mu(x)$..
Here, $K$ is the interaction potential defining the pairwise forces between
particles, $V$ is the local potential containing all exogenous forces,
$\sigma$ is a diffusivity, and $\left(B^{(i)}_{t}\right)_{i=1,\dots,N}$ are
independent Brownian motions each adapted to the same filtered probability
space $(\Omega,{\mathcal{B}},\mathbb{P},({\mathcal{F}}_{t})_{t\geq 0})$. The
empirical measure is defined

$\mu^{N}_{t}:=\frac{1}{N}\sum_{i=1}^{N}\delta_{X^{(i)}_{t}},$

and the convolution $\nabla K*\mu^{N}_{t}$ is defined

$\nabla
K*\mu^{N}_{t}(x)=\nabla\int_{\mathbb{R}^{d}}K(x-y)\,d\mu^{N}_{t}(y)=\frac{1}{N}\sum_{i=1}^{N}\nabla
K\left(x-X^{(i)}_{t}\right)$

where we set $\nabla K(0)=0$ whenever $\nabla K(0)$ is undefined. The recovery
problem we wish to solve is the following.

(P) Let
$\boldsymbol{\mathbb{X}}=(\mathbf{X}_{\mathbf{t}}^{(1)},\dots,\mathbf{X}_{\mathbf{t}}^{(M)})$
be discrete-time data at $L$ timepoints $\mathbf{t}:=(t_{1},\dots,t_{L})$ for
$M$ i.i.d. trials of the process (1.1) with $K=K^{\star}$, $V=V^{\star}$, and
$\sigma=\sigma^{\star}$ and let
$\boldsymbol{\mathbb{Y}}=\boldsymbol{\mathbb{X}}+\varepsilon$ be a corrupted
dataset. For some fixed compact domain ${\mathcal{D}}\subset\mathbb{R}^{d}$
containing $\text{supp}\left(\boldsymbol{\mathbb{Y}}\right)$, and finite-
dimensional hypothesis spaces333The set ${\mathcal{D}}-{\mathcal{D}}$ is
defined ${\mathcal{D}}-{\mathcal{D}}=\\{x-y\ :\
(x,y)\in{\mathcal{D}}\times{\mathcal{D}}\\}$. ${\mathcal{H}}_{K}\subset
L^{2}({\mathcal{D}}-{\mathcal{D}})$, ${\mathcal{H}}_{V}\subset
L^{2}({\mathcal{D}})$, and ${\mathcal{H}}_{\sigma}\subset
L^{2}({\mathcal{D}})$, solve

$\left(\widehat{K},\widehat{V},\widehat{\sigma}\right)=\text{argmin}_{K\in{\mathcal{H}}_{K},V\in{\mathcal{H}}_{V},\sigma\in{\mathcal{H}}_{\sigma}}\left\|{\nabla
K-\nabla
K^{\star}}\right\|_{L^{2}({\mathcal{D}}-{\mathcal{D}})}+\left\|{\nabla
V-\nabla
V^{\star}}\right\|_{L^{2}({\mathcal{D}})}+\left\|{\sigma-\sigma^{\star}}\right\|_{L^{2}({\mathcal{D}})}.$

The problem (P) is clearly intractable because we do not have access to
$K^{\star}$, $V^{\star}$, or $\sigma^{\star}$, and moreover the interactions
between these terms render simultaneous identification of them ill-posed. We
consider two cases: (i) $\varepsilon\neq 0$ and $\sigma^{\star}=0$,
corresponding to purely extrinsic noise, and (ii) $\varepsilon=0$ and
$\sigma^{\star}\neq 0$, corresponding to purely intrinsic noise. The extrinsic
noise case is important for many applications, such as cell tracking, where
uncertainty is present in the position measurements. In this case we examine
$\varepsilon$ representing i.i.d. Gaussian noise with mean zero and
variance444By $\mathbf{I}_{d}$ we mean the identity in $\mathbb{R}^{d}$.
$\epsilon^{2}\mathbf{I}_{d}$ added to each particle position in
$\boldsymbol{\mathbb{X}}$. In the case of purely intrinsic noise,
identification of the diffusivity $\sigma^{\star}$ is required as well as the
deterministic forces on each particle as defined by $K^{\star}$ and
$V^{\star}$. A natural next step is to consider the case with both extrinsic
and intrinsic noise. However, this is a topic for future work and thus beyond
the scope of this article.

##  2\. Background

Interacting particle systems (IPS) such as (1.1) are used to describe physical
and artificial phenomena in a range of fields including astrophysics [51, 25],
molecular dynamics [32], cellular biology [45, 50, 2], and opinion dynamics
[6]. In many cases the number of particles $N$ is large, with cell migration
experiments often tracking $10^{3}$-$10^{6}$ cells and simulations in physics
(molecular dynamics, particle-in-cell, etc.) requiring $N$ in the range
$10^{6}$-$10^{12}$. Inference of such systems from particle data thus requires
efficient means of computing pairwise forces from ${\mathcal{O}}(N^{2})$
interactions at each timestep for multiple candidate interaction potentials
$K$. Frequently, so-called mean-field equations at the continuum level are
sufficient to describe the evolution of the system, however in many cases
(e.g. chemotaxis in biology [29]) only phenomenological mean-field equations
are available. Moreover, it is often unclear how many particles $N$ are needed
for a mean-field description to suffice. Many fields are now developing
machine learning techniques to extract coarse-grained dynamics from high-
fidelity simulations (see [23] for a recent review in molecular dynamics). In
this work we provide a means for inferring governing mean-field equations from
particle data assumed to follow the dynamics (1.1) that is highly efficient
for large $N$, and is effective in learning mean-field equations when $N$ is
in range $10^{3}$-$10^{5}$.

Inference of the drift and diffusion terms for stochastic differential
equations (SDEs) is by now a mature field, with the primary method being
maximum-likelihood estimation, which uses Girsanov’s theorem together with the
Radon-Nykodym derivative to arrive at a log-likelihood function for
regression. See [3, 34] for some early works and [4] for a textbook on this
approach. More recently, sparse regression approaches using the Kramers-Moyal
expansion have been developed [10, 12, 33] and the authors of [43] use sparse
regression to learn population level ODEs from agent-based modeling
simulations. In addition, a neural network-based algorithm was developed in
[15].

Only in the last few years have significant strides been made towards
parameter inference of interacting particle systems such as (1.1) from data.
Apart from some exceptions, such as a Gaussian process regression algorithm
recently developed in [17], applications of maximum likelihood theory are by
far the most frequently studied. An early but often overlooked work by Kasonga
[28] extends the maximum-likelihood approach to inference of IPS, assuming
full availability of the continuous particle trajectories and the diffusivity
$\sigma$. Two decades later, Bishwal [5] further extended this approach to
discrete particle observations in the specific context of linear particle
interactions. In both cases, a sequence of finite-dimensional subspaces is
used to approximate the interaction function, and convergence is shown as the
dimension of the subspace $J$ and number of particles $N$ both approach
infinity. More recently, the maximum likelihood approach has been carried out
in [9, 35] in the case of radial interactions and in [14] in the case of
linear particle interactions and single-trajectory data (i.e. one instance of
the particle system). The authors of [46] recently developed an online maximum
likelihood method for inference of IPS, and in [24] maximum likelihood is
applied to parameter estimation in an IPS for pedestrian flow. It should also
be noted that parameter estimation for IPS is common in biological sciences,
with the most frequently used technique being nonlinear least squares with a
cost function comprised of summary statistics [36, 45].

Problem (P) is made challenging by the coupled effects of $K$, $V$ and
$\sigma$. In each of the previously mentioned algorithms, the assumption is
made that $\sigma$ is known and/or that $K$ takes a specific form (radial or
linear). In addition, the maximum likelihood-based approach approximates the
differential $dX_{t}^{(i)}$ of particle $i$ using a 1st-order finite
difference: $dX_{t}^{(i)}\approx X_{t+\Delta t}^{(i)}-X_{t}^{(i)}$, which is
especially ill-suited to problems involving extrinsic noise in the particle
positions. Our primary goal is to show that the weak-form sparse regression
framework allows for identification of the full model $(K,V,\sigma)$, with
significantly reduced computational complexity, when $N$ is on the order of
several thousands or more. The feasibility of this approach is grounded in the
convergence of IPS to associated mean-field equations. The reduction in
computational complexity follows from the reduction in evaluation of candidate
potentials (as discussed in Section 4.2), as well as the convolutional form of
the weak-form algorithm.

To the best of our knowledge, we present here the first weak-form sparse
regression approach for inference of interacting particle systems. We use a
two-step process: the density of particles is approximated using a density
kernel $G$ and then the WSINDy algorithm (weak-form sparse identification of
nonlinear dynamics) is applied in the PDE setting [40, 39]. WSINDy is a
modified version of the original SINDy algorithm [11, 44] where the weak
formulation of the dynamics is enforced using a family of test functions that
offers reduced computational complexity, high-accuracy recovery in low-noise
regimes, and increased robustness to high-noise scenarios. There are two works
that are most closely related to the current work. In [48], the authors learn
local hydrodynamic equations from active matter particle systems using the
SINDy algorithm in the strong-form PDE setting. In contrast to [48], our
approach learns nonlocal equations using the weak-form, however similarly to
[48] we perform model selection and inference of parameters using sparse
regression at the continuum level. The weak form provides an advantage because
no smoothness is required on the particle density (for requisite smoothness
the authors of [48] use a Gaussian kernel, which is more expensive to compute
than simple particle binning as done here). In [31], the authors apply the
maximum likelihood approach in the continuum setting on the underlying
nonlocal Fokker-Planck equation and learn directly the nonlocal PDE using
strong-form discretizations of the dynamics. While we similarly use the
continuum setting for inference (albiet in weak form), our approach differs
from [31] in that it is designed for the more realistic setting of discrete-
time particle data.

###  2.1. Contributions

The purpose of the present article is to show that the weak form provides an
advantage in speed and accuracy compared with existing inference methods for
particle systems when the number of particles is sufficiently large (on the
order of several thousand or more). The key points of this article include:

  1. (I)

Formulation of a weak-form sparse recovery algorithm for simultaneous
identification of the particle interaction force $K$, local potential $V$, and
diffusivity $\sigma$ from discrete-time particle data

  2. (II)

$L^{1}$ convergence of the resulting full-rank least-squares solution as the
number of particles $N\to\infty$ and timestep $\Delta t\to 0$

  3. (III)

Numerical illustration of (i) theoretical convergence rates in $N$ and (ii)
robustness to either intrinsic randomness (e.g. Brownian motion) or extrinsic
randomness (e.g. additive measurement noise)

###  2.2. Paper Organization

In Section 3 we review results from mean-field theory used to show convergence
of the weak-form method. In Section 4 we introduce the WSINDy algorithm
applied to interacting particles, including hyperparameter selection,
computational complexity, and convergence of the method under suitable
assumptions in the limit of large $N$. Section 5 contains numerical examples
exhibiting the convergence rates of the previous section and examining the
robustness of the algorithm to various sources of corruption, and Section 6
contains a discussion of extension and future directions.

##  3\. Review of mean-field theory

Our weak-form approach utilizes that under fairly general assumptions the
empirical measure $\mu^{N}_{t}$ of the process $\mathbf{X}_{t}$ defined in
(1.1) converges weakly to $\mu_{t}$, the distribution of the associated mean-
field process $X_{t}$ defined in (3.2). Specifically, under suitable
assumptions on $V,K,\sigma$ and $\mu_{0}$, there exists $T>0$ such that for
all $t\in[0,T]$, the mean-field limit555We use the notation $t\to\mu_{t}$ to
denote the evolution of probability measures. Subscripts will not be used to
denote differentiation.

$\lim_{N\to\infty}\mu^{N}_{t}=\mu_{t}$

holds in the weak topology of measures666Meaning that for all continuous
bounded functions $\phi:\mathbb{R}^{d}\to\mathbb{R}$,
$\int_{\mathbb{R}^{d}}\phi(x)d\mu^{N}_{t}(x)\to\int_{\mathbb{R}^{d}}\phi(x)d\mu_{t}(x)$.,
where $\mu_{t}$ is a weak-measure solution to the mean-field dynamics

$\partial_{t}\mu_{t}=\nabla\cdot\left(\left(\nabla K*\mu_{t}+\nabla
V\right)\mu_{t}\right)+\frac{1}{2}\sum_{i,j=1}^{d}\frac{\partial^{2}}{\partial
x_{i}\partial
x_{j}}\left(\sigma\sigma^{T}\mu_{t}\right),\quad\mu_{0}\in{\mathcal{P}}_{p}(\mathbb{R}^{d}).$
(3.1)

Equation 3.1 describes the evolution of the distribution of the McKean-Vlasov
process

$dX_{t}=\left(-\nabla V\left(X_{t}\right)-\nabla
K*\mu_{t}\left(X_{t}\right)\right)dt+\sigma(X_{t})\,dB_{t}.$ (3.2)

This implies that as $N\to\infty$, an initially correlated particle system
driven by pairwise interaction becomes uncorrelated and only interacts with
its mean field $\mu_{t}$. In particular, the following theorem summarizes
several mean-field results taken from the review article [26] with proofs in
[49, 38]. (Note that for a function $f:\mathbb{R}^{d}\to Y$, where $Y$ is a
metric space with metric $\rho$, we define $\text{Lip}_{S}(f)$ by

$\text{Lip}_{S}(f):=\sup_{\\{x,y\in
S\\}}\frac{\rho(f(x),f(y))}{\left|x-y\right|}$

with $\text{Lip}(f):=\text{Lip}_{\mathbb{R}^{d}}(f)$. Throughout we use
$|\cdot|$ to denote the Euclidean norm.)

######  Theorem 3.1.

Assume that $K$ is globally Lipschitz, $V=0$, and
$\sigma(x)=\sigma=\text{const.}$ In addition assume that
$\mu_{0}\in{\mathcal{P}}_{2}(\mathbb{R}^{d})$. Then for any $T>0$, for all
$t\leq T$ it holds that

  1. (i)

There exists a unique solution $(X_{t},\mu_{t})$ where $X_{t}$ is a strong
solution to (3.2) and $\mu_{t}$ is a weak-measure solution to (3.1).

  2. (ii)

For any $\phi\in C^{1}_{b}(\mathbb{R}^{d})$,

$\mathbb{E}\left|\frac{1}{N}\sum_{i=1}^{N}\phi(X_{i}(t))-\int_{\mathbb{R}^{d}}\phi(x)d\mu_{t}(x)\right|\leq\frac{C}{\sqrt{N}}$
(3.3)

with $C$ depending on $\text{Lip}(\nabla K)$, $\left\|{\phi}\right\|_{C^{1}}$
and $T$.

  3. (iii)

For any $k\in\mathbb{N}$, $a.e.$-$t<T$, the $k$-particle marginal

$\rho_{t}^{(k)}(x_{1},\dots,x_{k}):=\int_{\mathbb{R}^{d(N-k)}}F_{t}^{N}(x_{1},\dots,x_{k},x_{k+1},\dots,x_{N})\,dx_{k+1}\cdots
dx_{N}$

converges weakly to $\mu_{t}^{\otimes k}$ as $N\to\infty$, where
$F^{N}_{t}\in{\mathcal{P}}(\mathbb{R}^{Nd})$ is the distribution of
$\mathbf{X}_{t}$.

Theorem 3.1 immediately extends to the case of $V$ and $\sigma$ both globally
Lipschitz and has been extended to $K$ locally-Lipschitz in [8], $K$ with
Coulomb-type singularity at the origin in [7], and domains with boundaries in
[19, 18]. Analysis of the model (3.1) continues to evolve in various contexts,
including with analysis of equilibria [41, 21, 13] and connections to deep
learning [1].

###  3.1. Weak form

Despite the ${\mathcal{O}}(N^{-1/2})$ convergence of the empirical measure in
Theorem 3.1, it is unclear at what particle number $N$ the mean-field
equations become a suitable framework for inference using particle data, due
to the complex variance structure at any finite $N$. A key piece of the
present work is to show that the weak form of the mean-field equations does
indeed provide a suitable setting when $N$ is at least several thousand.
Moreover, since in many cases (3.1) can only be understood in a weak sense,
the weak form is the natural framework for identification. We say that
$\mu_{t}$ is a weak solution to (3.1) if for any $\psi\in
C^{2}(\mathbb{R}^{d}\times(0,T))$ compactly supported it holds that

$\int_{0}^{T}\int_{\mathbb{R}^{d}}\partial_{t}\psi(x,t)\,d\mu_{t}(x)dt=\int_{0}^{T}\int_{\mathbb{R}^{d}}\left(\nabla\psi(x,t)\cdot\left(\nabla
K*\mu_{t}(x)+\nabla
V(x)\right)-\frac{1}{2}\text{Tr}\left(\nabla^{2}\psi(x,t)\sigma(x)\sigma^{T}(x)\right)\right)d\mu_{t}(x)dt,$
(3.4)

where $\nabla^{2}\psi$ denotes the Hessian of $\psi$ and
$\text{Tr}(\mathbf{A})$ is the trace of the matrix $\mathbf{A}$. Our method
requires discretizing (3.4) for all $\psi\in\Psi$ where
$\Psi=(\psi_{1},\dots,\psi_{n})$ is a suitable test function basis, and
approximating the mean field density $\mu_{t}$ with a discrete density $U_{t}$
constructed from particle data at time $t$. We then find $K,V$, and $\sigma$
within specified finite-dimensional function spaces.

##  4\. Algorithm

We propose the general algorithm 4.1 for discovery of mean-field equations
from particle data. The inputs are a discrete-time sample
$\boldsymbol{\mathbb{Y}}$ containing $M$ experiments each with $N$ particle
positions over $L$ timepoints $\mathbf{t}=(t_{1},\dots,t_{L})$, and the
following hyperparameters are defined by the user: (i) a kernel $G$ used to
map the empirical measure $\mu^{N}_{t}$ to an approximate density $U_{t}$,
(ii) a spatial grid $\mathbf{C}$ over which to evaluate the approximate
density $\mathbf{U}_{t}=U_{t}(\mathbf{C})$, (iii) a library of trial functions
$\mathbb{L}=\\{\mathbb{L}_{K},\mathbb{L}_{V},\mathbb{L}_{\sigma}\\}=\\{(K_{j})_{j=1}^{J_{K}},(V_{j})_{j=1}^{J_{V}},(\sigma_{j})_{j=1}^{J_{\sigma}}\\}$,
(iv) a basis of test functions $\Psi=(\psi_{k})_{k=1}^{n}$, and (v) a
quadrature rule over $(\mathbf{C},\mathbf{t})$ denoted by an inner product
$\left\langle\cdot,\cdot\right\rangle$, and (vi) sparsity factors
$\boldsymbol{\lambda}$ for the modified sequential thresholding algorithm
(MSTLS) reviewed below. We discuss the choices of these hyperparameters in
Section 4.1, the computational complexity of the algorithm in Section 4.2, and
convergence of the algorithm in Section 4.3.

Algorithm 4.1 WSINDy for Particle Systems
$({\widehat{\mathbf{w}}},\,\hat{\lambda})=$
WSINDy$(\boldsymbol{\mathbb{Y}},\mathbf{t}\,;\
G,\,\mathbf{C},\,\mathbb{L},\,\Psi,\,\left\langle\cdot,\cdot\right\rangle,\,\boldsymbol{\lambda})$

1: for $\ell=1:L$ do

2: for $m=1:M$ do

3:
$\mathbf{U}_{\ell}^{(m)}=\int_{\mathbb{R}^{d}}G(\mathbf{C},y)d\mu^{N}_{t_{\ell}}(y)$
where $\mu^{N}_{t_{\ell}}$ is the empirical measure for
$\mathbf{Y}_{t_{\ell}}^{(m)}$

4: end for

5: $\mathbf{U}_{\ell}=\frac{1}{M}\sum_{m=1}^{M}\mathbf{U}_{\ell}^{(m)}$

6: end for

7:

8: for $j=1:J_{K}$ do

9: for $k=1:n$ do

10: $\mathbf{G}_{kj}^{K}=\left\langle\nabla\psi_{k},\mathbf{U}\nabla
K_{j}*\mathbf{U}\right\rangle$

11: end for

12: end for

13:

14: for $j=1:J_{V}$ do

15: for $k=1:n$ do

16: $\mathbf{G}_{kj}^{V}=\left\langle\nabla\psi_{k},\mathbf{U}\nabla
V_{j}\right\rangle$

17: end for

18: end for

19:

20: for $j=1:J_{\sigma}$ do

21: for $k=1:n$ do

22:
$\mathbf{G}_{kj}^{\sigma}=\frac{1}{2}\sum_{p,q=1}^{d}\left\langle\partial_{x_{p}x_{q}}\psi_{k},(\sigma_{j}\sigma_{j}^{T})_{pq}\mathbf{U}\right\rangle$

23: end for

24: end for

25: $\mathbf{G}=[\mathbf{G}^{K}\ \mathbf{G}^{V}\ \mathbf{G}^{\sigma}]$

26:

27: for $k=1:n$ do

28: $\mathbf{b}_{k}=\left\langle\partial_{t}\psi_{k},\mathbf{U}\right\rangle$

29: end for

30:

31:
$(\widehat{\mathbf{w}},\widehat{\lambda})=\text{MSTLS}(\mathbf{G},\mathbf{b};\,\boldsymbol{\lambda})$

###  4.1. Hyperparameter Selection

#### Quadrature

We assume that the set of gridpoints $\mathbf{C}$ in Algorithm 4.1 is chosen
from some compact domain ${\mathcal{D}}\subset\mathbb{R}^{d}$ containing
$\text{supp}\left(\boldsymbol{\mathbb{Y}}\right)$. The choice of $\mathbf{C}$
(and ${\mathcal{D}}$) must be chosen in conjunction with the quadrature
scheme, which includes integration in time using the given timepoints
$\mathbf{t}$ as well as space. For completeness, the inner products in lines
10, 16, 22 and 27 of Algorithm 4.1 are defined in the continuous setting by

$\left\langle
f,g\right\rangle=\int_{0}^{T}\int_{{\mathcal{D}}}f(x,t)g(x,t)dxdt,$

and the convolution in line 10 is defined by

$\nabla K_{j}*U_{t}(x)=\int_{{\mathcal{D}}}\nabla K_{j}(x-y)U_{t}(y)dy.$

In the present work we adopt the scheme used in the application of WSINDy for
local PDEs [39], which includes the trapezoidal rule in space and time with
test functions $\psi$ compactly supported in ${\mathcal{D}}\times(0,T)$. We
take ${\mathcal{D}}$ to be an equally-spaced rectangular grid enclosing
$\text{supp}\left(\boldsymbol{\mathbb{Y}}\right)$ in order to efficiently
evaluate convolution terms. In what follows we denote by
$\left\langle\cdot,\cdot\right\rangle$ the continuous inner product,
$\left\langle\cdot,\cdot\right\rangle_{h}$ the inner product over
${\mathcal{D}}\times[0,T]$ evaluated using the composite trapezoidal rule in
space with meshwidth $h$ and Lebesgue integration in time, and by
$\left\langle\cdot,\cdot\right\rangle_{h,\Delta t}$ the trapezoidal rule in
both space and time, with meshwidth $h$ in space and $\Delta t$ in time. With
some abuse of notation, $f*g$ will denote the convolution of $f$ and $g$,
understood to be discrete or continuous by the context. Note also that we
denote $\mu^{N}$, $\mu$ and $U$ the measures over $\mathbb{R}^{d}\times[0,T]$
defined by $\mu^{N}_{t}\lambda_{[0,T]}$, $\mu_{t}\lambda_{[0,T]}$ and
$U_{t}\lambda_{[0,T]}$, respectively, where $\lambda_{[0,T]}$ is the Lebesgue
measure on $[0,T]$.

#### Density Kernel

Having chosen the domain ${\mathcal{D}}\subset\mathbb{R}^{d}$ containing the
particle data $\boldsymbol{\mathbb{Y}}$, let
$P^{h}=\left\\{B_{k}\right\\}_{k}$ be a partition of ${\mathcal{D}}$
($\cup_{k}B_{k}={\mathcal{D}}$) with $h$ indicating the size of the atoms
$B_{k}$. For the remainder of the paper we take $B_{k}$ to be hypercubes of
equal side length $h$ in order to minimize computation time for integration,
although this is by no means necessary. For particle positions
$\mathbf{X}_{t}$, we define the histogram777The indicator function is defined
$\mathbbm{1}_{A}(x):=\begin{cases}1,&x\in A\\\ 0,&x\notin A\end{cases}$.

$U_{t}=\sum_{k}\frac{1}{|B_{k}|}\mathbbm{1}_{B_{k}}(x)\left(\frac{1}{N}\sum_{i}\mathbbm{1}_{B_{k}}(X^{(i)}_{t})\right)=\int_{\mathcal{D}}G(x,y)d\mu^{N}_{t}(y).$
(4.1)

Here the density kernel is defined

$G(x,y)=\sum_{k}\frac{1}{|B_{k}|}\mathbbm{1}_{B_{k}}(x)\mathbbm{1}_{B_{k}}(y),$

and in this setting the corresponding spatial grid
$\mathbf{C}=(\mathbf{c}_{k})_{k}$ is the set of center-points of the boxes
$B_{k}$, from which we define the discrete histogram data
$\mathbf{U}_{t}=U_{t}(\mathbf{C})$. The discrete histogram $\mathbf{U}_{t}$
then serves as an approximation to the mean-field distribution $\mu_{t}$.

Pointwise estimation of densities from samples of particles usually requires
large numbers of particles to achieve reasonably low variance, and in general
the variance grows inversely proportional to the bin width $h$. One benefit of
the weak form is that integrating against a histogram $U$ does not suffer from
the same increase in variance with small $h$. In particular,

######  Lemma 4.1.

Let $\mathbf{Y}=(Y^{(1)},\dots,Y^{(N)})$ be i.i.d. samples from
$\mu\in{\mathcal{P}}(\mathbb{R}^{d})$ with associated empirical measure
$\mu^{N}$ and let $U$ be the histogram computed with kernel $G$ using (4.1)
with $n$ bins of equal sidelength $h$. Then for any $\psi$ in $C^{1}$
compactly supported in ${\mathcal{D}}$, we have the root-mean-squared error

$\left(\mathbb{E}\left[\left(\left\langle\psi,U\right\rangle_{{\mathcal{T}}}-\left\langle\psi,\mu\right\rangle\right)^{2}\right]\right)^{1/2}\leq\left\|{\psi}\right\|_{C^{1}}\left(\sqrt{2}^{d-3}h+N^{-1/2}\right).$

###### Proof.

Using the compact support of $\psi$, we have

$\left\langle\psi,U\right\rangle_{{\mathcal{T}}}=\left\langle\psi,\int_{\mathbb{R}^{d}}G(\cdot,y)d\mu^{N}\right\rangle_{h}=\left\langle\psi^{\mathbf{C}},\mu^{N}\right\rangle=\frac{1}{N}\sum_{i=1}^{N}\psi^{\mathbf{C}}(Y^{(i)})$

where

$\psi^{\mathbf{C}}(x)=\sum_{k=1}^{K}\psi(c_{k})\mathbbm{1}_{B_{k}}(x)$ (4.2)

is the midpoint approximation of $\psi$. We then have the squared bias

$\displaystyle\text{bias}(\left\langle\psi,U\right\rangle_{{\mathcal{T}}})^{2}$
$\displaystyle=\left(\mathbb{E}\left[\left\langle\psi,U\right\rangle_{{\mathcal{T}}}\right]-\int_{\mathcal{D}}\psi(x)d\mu(x)\right)^{2}$
$\displaystyle=\left(\int_{\mathcal{D}}\psi^{\mathbf{C}}(x)d\mu(x)-\int_{\mathcal{D}}\psi(x)d\mu(x)\right)^{2}$
$\displaystyle\leq\left\|{\nabla\psi}\right\|^{2}_{\infty}2^{d-3}h^{2}$

and the variance, using the fact that $Y^{(i)}$ are independent,

$\text{Var}\left(\left\langle\psi,U\right\rangle_{{\mathcal{T}}}\right)=\frac{1}{N}\mathbb{E}_{X\sim\mu}\Big{[}\psi^{\mathbf{C}}(X)\left(\psi^{\mathbf{C}}(X)-\mathbb{E}_{Z\sim\mu}\left[\psi^{\mathbf{C}}(Z)\right]\right)\Big{]}\leq\frac{1}{N}\left\|{\psi}\right\|_{\infty}^{2}.$

The result follows since

$\mathbb{E}\left[\left(\left\langle\psi,U\right\rangle_{{\mathcal{T}}}-\left\langle\psi,\mu\right\rangle\right)^{2}\right]=\text{bias}(\left\langle\psi,U\right\rangle_{{\mathcal{T}}})^{2}+\text{Var}\left(\left\langle\psi,U\right\rangle_{{\mathcal{T}}}\right).$

∎

The previous lemma in particular shows that small bin-width $h$ does not
negatively impact $\left\langle\psi,U\right\rangle_{h}$ as an estimator of
$\left\langle\psi,\mu\right\rangle$, which is in contrast to $U(x)$ as a
pointwise estimator of $\mu(x)$. For example, if we assume that $\mathbf{Y}$
is sampled from a $C^{1}$ density $\mu$, it is well known that the mean-square
optimal bin width is $h={\mathcal{O}}(N^{-1/3})$ [22]. Summarizing this
result, elementary computation reveals the pointwise bias for $x\in B_{k}$,

$\text{bias}(U(x))=\mathbb{E}\left[U(x)\right]-\mu(x)=\frac{\mu(B_{k})}{|B_{k}|}-\mu(x):=\mu(\xi)-\mu(x)$

for some $\xi\in B_{k}$. Letting $L_{k}=\max_{x\in B_{k}}|\nabla\mu(x)|$, we
have

$\text{bias}(U(x))^{2}\leq L_{k}^{2}2^{d-1}h^{2}.$

For the variance we get

$\text{Var}\left(U(x)\right)=\frac{1}{N}\frac{\mu(B_{k})(1-\mu(B_{k}))}{|B_{k}|^{2}}=\frac{\mu(\xi)}{N}\left(1-\mu(B_{k})\right)\frac{1}{\sqrt{2}^{d-1}h},$

and hence a bound for the mean-squared error

$\mathbb{E}\left[\left(U(x)-\mu(x)\right)^{2}\right]\leq
L_{k}^{2}2^{d-1}h^{2}+\frac{\mu(\xi)}{N\sqrt{2}^{d-1}}h^{-1}.$

Minimizing the bound over $h$ we find an approximately optimal box width

$h^{*}=\left(\frac{\rho(\xi)}{2^{\frac{3d-1}{2}}L_{k}^{2}}\right)^{1/3}N^{-1/3}={\mathcal{O}}(N^{-1/3}),$

which provides an overall pointwise root-mean-squared error of
${\mathcal{O}}(N^{-1/3})$. Hence, not only does the weak-form remove the
inverse $h$ dependence in the variance, but fewer particles are needed to
accurately approximate integrals of the density $\mu$.

#### Trial Function Library

The general Algorithm 4.1 does not impose a radial structure for the
interaction potential $K$, nor does it assume any prior knowledge that the
particle system is in fact interacting. In the examples below we use monomial
and/or trigonometric libraries for $K$, $V$, and $\sigma$ to show that sparse
regression is effective in selecting the correct interaction terms from a
library also containing local drift terms. (Details of the libraries used in
examples can be found in Appendix A).

#### Test Function Basis

For the test functions $(\psi_{k})_{k\in[n]}$ we use the same approach as the
PDE setting [39], namely we fix a reference test function $\psi$ and set

$\psi_{k}(x,t)=\psi(\mathbf{x}_{k}-x,t_{k}-t)$

where ${\mathcal{Q}}:=\\{(\mathbf{x},t_{k})\\}_{k\in[n]}$ is a fixed set of
query points. This, together with a separable representation

$\psi(x,t)=\phi_{1}(x_{1})\cdots\phi_{d}(x_{d})\phi_{d+1}(t),$

enables construction of the linear system $\mathbf{G}$, $\mathbf{b}$ using the
FFT. We choose $\phi_{j}$, $1\leq j\leq d+1$, of the form

$\phi_{m,p}(v;\,\Delta):=\max\left\\{1-\left(\frac{v}{m\Delta}\right)^{2},0\right\\}^{p}$
(4.3)

where $m$ is the integer support parameter such that $\phi_{m,p}$ is supported
on $2m+1$ points of spacing $\Delta\in\\{h,\Delta t\\}$ and $p\geq 1$ is the
degree of $\phi_{m,p}$. Using the convergence analysis below, we need
$\phi_{j}$ to be three times continuously differentiable, yet with Lipschitz
constant as low as possible, hence for $\phi_{d+1}$ (along the time axis) we
pick $p=p_{t}=3$, while for all spatial dimensions we set $p=p_{x}=5$. We
choose the support parameters $m_{t}$ and $m_{x}$ using the changepoint
algorithm in [39, Appendix A], enforcing that each $\phi_{j}$ is supported on
a minimum of 6 gridpoints and a maximum of half of the width of the domain
along the respective coordinate. For $Q$, we sample points uniformly within
${\mathcal{C}}$ of frequency $s_{x}$ or $s_{t}$ depending on if the coordinate
in spatial or temporal. The resulting values
$m_{x},m_{t},p_{x},p_{t},s_{x},s_{t}$ for each example can be found in
Appendix A.

#### Sparsity Regression

As in [39], we enforce sparsity using a modified sequential thresholding
least-squares algorithm. With888The Moore-Penrose inverse
$\mathbf{A}^{\dagger}$ is defined for a rank-$r$ matrix $\mathbf{A}$ using the
reduced SVD $\mathbf{A}=U_{r}\Sigma_{r}V^{*}_{r}$ as
$\mathbf{A}^{\dagger}:=V_{r}\Sigma^{-1}_{r}U_{r}^{*}$. The subscript $r$
denotes restriction to the first $r$ columns.
$\mathbf{w}^{0}=\mathbf{G}^{\dagger}\mathbf{b}$, we define

$\text{MSTLS}(\mathbf{G},\mathbf{b};\lambda\,)\qquad\begin{dcases}\hskip
12.23447pt{\mathcal{I}}^{\ell}=\\{1\leq i\leq SJ\ :\
L^{\lambda}_{i}\leq|\mathbf{w}^{\ell}_{i}|\leq U^{\lambda}_{i}\\}\\\
\mathbf{w}^{\ell+1}=\text{argmin}_{\text{supp}\left(\mathbf{w}\right)\subset{\mathcal{I}}^{\ell}}\left\|{\mathbf{G}\mathbf{w}-\mathbf{b}}\right\|_{2}^{2}.\end{dcases}$
(4.4)

where the bounds are defined

$\begin{dcases}L_{i}^{\lambda}=\lambda\max\left\\{1,\
\frac{\left\|{\mathbf{b}}\right\|}{\left\|{\mathbf{G}_{i}}\right\|}\right\\}\\\
U_{i}^{\lambda}=\frac{1}{\lambda}\min\left\\{1,\
\frac{\left\|{\mathbf{b}}\right\|}{\left\|{\mathbf{G}_{i}}\right\|}\right\\}\end{dcases},\qquad
1\leq i\leq SJ.$ (4.5)

We then select the sparsity threshold $\widehat{\lambda}$ as the smallest
minimizer of the cost function

${\mathcal{L}}(\lambda)=\frac{\left\|{\mathbf{G}(\mathbf{w}^{\lambda}-\mathbf{w}^{0})}\right\|_{2}}{\left\|{\mathbf{G}\mathbf{w}^{0}}\right\|_{2}}+\frac{\left\|{\mathbf{w}^{\lambda}}\right\|_{0}}{SJ}$
(4.6)

over all $\lambda$ in a specified finite set $\boldsymbol{\lambda}$, where
$\mathbf{w}^{\lambda}:=\text{MSTLS}(\mathbf{G},\mathbf{b};\lambda\,)$. We set
the final model coefficients to
${\widehat{\mathbf{w}}}:=\mathbf{w}^{\widehat{\lambda}}$. The bounds (4.5)
enforce a quasi-dominant balance rule, such that
$\left\|{\mathbf{w}_{i}\mathbf{G}_{i}}\right\|_{2}$ is within
$-\log_{10}(\lambda)$ orders of magnitude from
$\left\|{\mathbf{b}}\right\|_{2}$ and $|\mathbf{w}_{i}|$ is within
$-\log_{10}(\lambda)$ orders of magnitude from $1$ (the coefficient of time
derivative $\partial_{t}\mu_{t}$). Minimizers of the cost function
${\mathcal{L}}$ then equally weight the accuracy and sparsity of
$\mathbf{w}^{\widehat{\lambda}}$. By choosing $\widehat{\lambda}$ to be the
smallest minimizer of ${\mathcal{L}}$ over $\boldsymbol{\lambda}$, we identify
the thresholds $\lambda\in\boldsymbol{\lambda}$ such that
$\lambda<\widehat{\lambda}$ as those resulting in an overfit model.

###  4.2. Computational Complexity

To compute convolutions against each $\nabla K_{j}$, we evaluate
$\partial_{x_{i}}K_{j}$, $1\leq i\leq d$, at the grid $\mathbf{C}-\mathbf{C}$
defined by

$\mathbf{C}-\mathbf{C}:=\\{x\in\mathbb{R}^{d}\ :\
x=(i_{1}h,\dots,i_{d}h),\quad-m_{\ell}\leq i_{\ell}\leq m_{\ell}\\}$

where $h$ is the spacing of $\mathbf{C}$ and $m_{\ell}$, $1\leq\ell\leq d$, is
the number of points in $\mathbf{C}$ along the $\ell$th spatial dimension.
Then $\mathbf{C}-\mathbf{C}$ discretizes the set

${\mathcal{D}}-{\mathcal{D}}:=\\{x-y\in\mathbb{R}^{d}\ :\
(x,y)\in{\mathcal{D}}\times{\mathcal{D}}\\}$

which contains all observed interparticle distances. (In words, to form
$\mathbf{C}-\mathbf{C}$ we shift $\mathbf{C}$ to lie in the positive orthant
$\\{x\in\mathbb{R}^{d}\ :\ x_{\ell}\geq 0,\quad 1\leq\ell\leq d\\}$, and then
reflect $\mathbf{C}$ through each coordinate plane $x_{\ell}=0$,
$1\leq\ell\leq d$.) In this way $\partial_{x_{i}}K_{j}$ is evaluated at
$2^{d}|\mathbf{C}|$ points, where $|\mathbf{C}|$ is the number of points in
the grid $\mathbf{C}$. Define
$\partial_{x_{i}}\mathbf{K}_{j}:=\partial_{x_{i}}K_{j}(\mathbf{C}-\mathbf{C})$.
Since $\mathbf{C}$ is equally spaced, we use the $d$-dimensional FFT to
compute the convolutions

$\partial_{x_{i}}\mathbf{K}_{j}*\mathbf{U}_{t}\approx\partial_{x_{i}}K_{j}*U_{t}(\mathbf{C}),\quad
t\in\mathbf{t},$

where only entries corresponding to interactions within $\mathbf{C}$ need to
be retained. For $d=1$ this amounts to
${\mathcal{O}}(|\mathbf{C}|\log|\mathbf{C}|)$ flops per timestep. For $d=2$
and higher dimensions, the $d$-dimensional FFT is considerably slower unless
one of the arrays is separable. Trial interaction potentials $K_{j}$ can be
chosen to be a sum of separable functions,

$K(x)=\sum_{q=1}^{Q}k_{1,q}(x_{1})\cdots k_{d,q}(x_{d}),$

in which case only a series of one-dimensional FFTs are needed, and again the
cost is ${\mathcal{O}}(|\mathbf{C}|\log|\mathbf{C}|)$ per timestep. When $K$
is not separable, we propose using a low-rank approximation

$\partial_{x_{i}}\mathbf{K}_{j}=\sum_{q=1}^{Q}\sigma_{q}\mathbf{k}_{1,q}\otimes\cdots\otimes\mathbf{k}_{d,q}$

to exploit the efficiency of FFT in one dimension. For $d=2$, this is
accomplished using the truncated SVD, while for higher dimensions there does
not exist a unique best rank-$Q$ tensor approximation, although several
efficient algorithms are available to compute a sufficiently accurate
decomposition [37, 47, 27] (and the field of fast tensor decompositions is
advancing rapidly). In the examples below we consider only $d=1$ and $d=2$,
and leave extension to higher dimensions to future work.

Using low-rank approximations, the mean-field approach provides a significant
reduction in computational complexity compared to direct evaluations of
particle trajectories when $N$ is sufficiently large. A particle-level
computation of the nonlocal force in weak-form requires evaluating terms of
the form

$\sum_{\ell=1}^{L}\left(\frac{1}{N^{2}}\sum_{i=1}^{N}\sum_{j=1}^{N}\partial_{x}\psi(X^{(i)}_{t_{\ell}},t_{\ell})\partial_{x}K(X^{(i)}_{t_{\ell}}-X^{(j)}_{t_{\ell}})\right)\Delta
t.$

For a single candidate interaction potential $K$, a collection of $J$ test
functions $\psi$, and $M$ experiments, this amounts to $JMNL+MLN^{2}$ function
evaluations in $\mathbb{R}^{d}$ and ${\mathcal{O}}(MJLN^{2})$ flops. If we use
the proposed method, exploiting the convolutional structure of integration
against a separable reference test function $\psi$ and a rank $Q$
approximation of $\partial_{x}\mathbf{K}$, we instead evaluate

$\partial_{x}\psi*\left(U(\partial_{x}K*U)\right)$

using $2^{d}|\mathbf{C}|$ evaluations of $\partial_{x}K$, reused at each of
the $L$ timepoints, and ${\mathcal{O}}(L|\mathbf{C}|Q\log(|\mathbf{C}|))$
flops999Neglecting the cost of computing the histogram $\mathbf{U}$ and
evaluating $\psi(\mathbf{C})$, amounting to an additional
${\mathcal{O}}(NML+|\mathbf{C}|)$ flops, as these terms are reused in each
column of $\mathbf{G}$ and $\mathbf{b}$. Figure 1 provides a visualization of
the reduction in function evaluations for $L=100$ timepoints and $M=10$
experiments over a range of $N$ and $|\mathbf{C}|^{1/d}$ (points along each
spatial dimension when $|\mathbf{C}|$ is a hypercube) in $d=2$ and $d=3$
spatial dimensions. Table 5 in Appendix A lists walltimes for the examples
below, showing that with $N=64,000$ particles the full algorithm implemented
in MATLAB runs in under 10 seconds with all computations in serial on a laptop
with an AMD Ryzen 7 pro 4750u processor and 38.4 GB of RAM. The dependence on
$N$ is only through the ${\mathcal{O}}(N)$ computation of the histrogram,
hence this approach may find applications in physical coarse-graining (e.g. of
molecular dynamics or plasma simulations).

|
---|---
Figure 1. Factor by which the mean-field evaluation of interaction forces
using histograms reduces total function evaluations as a function of particle
number $N$ and average gridpoints per coordinate $|\mathbf{C}|^{1/d}$ for data
with $M=10$ experiments each with $L=100$ timepoints. For example, with $d=2$
spatial dimensions (left) and $N>2000$ particles, the number of function
evaluations is reduced by at least a factor of $10^{4}$.

###  4.3. Convergence

We now show that the estimators $\widehat{K}$, $\widehat{V}$, and
$\widehat{\sigma}$ of the weak-form method converge with a rate
${\mathcal{O}}(h+N^{-1/2}+\Delta t^{\eta})$ when ordinary least squares is
used (i.e. $\boldsymbol{\lambda}=0$) and only $M=1$ experiment is available.
Here $\eta>0$ is the Hölder exponent of the sample paths of the process
$\mathbf{X}_{t}$. We assume that ${\mathcal{D}}$, $\mathbf{C}$, $G$, $P^{h}$
and the resulting histogram $\mathbf{U}=(\mathbf{U}_{t})_{t\leq T}$ are as in
Section 4.1. We make the following assumptions on the true model and resulting
linear system throughout this section.

######  Assumptions 1.

Let $p\geq 1$ be fixed.

  1. (I)

For each $N\geq 2$, $\mathbf{X}_{t}=(X^{(1)}_{t},\dots,X^{(N)}_{t})$ is a
strong solution to (1.1) for $t\in[0,T]$, and for some $\eta>0$ the sample
paths $t\to X^{(i)}_{t}(\omega)$ are almost-surely $\eta$-Hölder continuous,
i.e.

$|X^{(i)}_{t}(\omega)-X^{(i)}_{s}(\omega)|\leq
C_{\eta}|t-s|^{\eta},\quad\forall 0\leq s\leq t\leq T,\quad\forall
i\in[N],\quad\text{for a.e. }\omega\in\Omega.$

  2. (II)

The initial particle distribution $\mu_{0}$ satisfies the moment bound

$\int_{\mathbb{R}^{d}}|x|^{p}d\mu_{0}(x):=M_{p}<\infty.$

  3. (III)

$\nabla K^{\star}$ and $\nabla V^{\star}$ satisfy for some $C_{p}>0$ the
growth bound:

$|\nabla V^{\star}(x)-\nabla V^{\star}(y)|+|\nabla K^{\star}(x)-\nabla
K^{\star}(y)|\leq C_{p}|x-y|(1+\max\\{|x|,|y|\\}^{p-1}),\quad
x,y\in\mathbb{R}^{d}$

  4. (IV)

For the same constant $C_{p}>0$, it holds that101010For
$\mathbf{A}\in\mathbb{R}^{d\times d}$ the Frobenius norm is defined
$\left\|{\mathbf{A}}\right\|_{F}=\sqrt{\text{Tr}(\mathbf{A}^{T}\mathbf{A})}$

$\left\|{\sigma^{\star}(x)-\sigma^{\star}(y)}\right\|_{F}\leq
C_{p}|x-y|^{1/2}(1+\max\\{|x|,|y|\\}^{p/2-1/2}),\quad x,y\in\mathbb{R}^{d}$

  5. (V)

The test functions $(\psi_{k})_{k\in[n]}\subset
C^{3}(\mathbb{R}^{d}\times(0,T))$ are compactly supported and together with
the library $\mathbb{L}$ are such that $\mathbf{G}$ has full column rank.
Moreover $\left\|{\mathbf{G}^{\dagger}}\right\|_{1}\leq C_{\mathbf{G}}$ almost
surely, where $\left\|{\mathbf{G}^{\dagger}}\right\|_{1}$ is the induced
matrix 1-norm of $\mathbf{G}^{\dagger}$.

  6. (VI)

The true functions $K^{\star}$, $V^{\star}$ and $\sigma^{\star}$ are in the
span of $\mathbb{L}$.

######  Remark 4.1.

Some consequences of Assumption 1 are the following: $\eta$-Hölder continuous
sample paths implies that for each $t\in[0,T]$,

$\int_{\mathbb{R}^{d}}|x|^{p}d\mu^{N}_{t}=\frac{1}{N}\sum_{i=1}^{N}|X^{(i)}_{t}|^{p}\leq\frac{2^{p}}{N}\sum_{i=1}^{N}|X^{(i)}_{0}|^{p}+C_{\eta}2^{p}t^{p\eta}.$

Together with the $p$th moment bound on $\mu_{0}$, this implies

$\mathbb{E}\left[\sup_{t\leq
T}\int_{\mathbb{R}^{d}}|x|^{p}d\mu^{N}_{t}\right]\leq
2^{p}(M_{p}+C_{\eta}T^{p\eta}),$ (4.7)

independent of $N$. The growth bounds on $\nabla K^{\star}$, $\nabla
V^{\star}$ and $\sigma^{\star}$ imply that for some $C>0$,

$|\nabla K^{\star}(x)|+|\nabla
V^{\star}(x)|+\left\|{\sigma^{\star}(x)(\sigma^{\star}(x))^{T}}\right\|_{F}\leq
C(1+|x|^{p}),$ (4.8)

where $\left\|{\cdot}\right\|_{F}$ is the Frobenius norm.

We will now define some notation and prove some lemmas. Define the weak-form
operator

$\mathscr{L}(\rho,\psi,\left\langle\cdot,\cdot\right\rangle):=\left\langle\partial_{t}\psi-\nabla\psi\cdot\nabla
K^{\star}*\rho-\nabla\psi\cdot\nabla
V^{\star}+\frac{1}{2}\text{Tr}\left(\nabla^{2}\psi\sigma^{\star}(\sigma^{\star})^{T}\right),\rho\right\rangle,$
(4.9)

where $\rho=(\rho_{t})_{t\leq T}$ is a curve in
${\mathcal{P}}_{p}(\mathbb{R}^{d})$, $\psi$ is a $C^{2}$ function compactly
supported over $\mathbb{R}^{d}\times(0,T)$, and
$\left\langle\cdot,\cdot\right\rangle$ is an inner product over
$\mathbb{R}^{d}\times(0,T)$. If $\rho=(\mu_{t})_{t\leq T}$ is a weak solution
to (3.1) and $\left\langle\cdot,\cdot\right\rangle$ is the
$L^{2}(\mathbb{R}^{d})$ inner product then
$\mathscr{L}(\rho,\psi,\left\langle\cdot,\cdot\right\rangle)=0$. If instead
$\rho=(\mu^{N}_{t})_{t\leq T}$, then by Itô’s formula
$\mathscr{L}(\rho,\psi,\left\langle\cdot,\cdot\right\rangle)$ takes the form
of an Itô integral, and we have the following:

######  Lemma 4.2.

There exists a constant $C>0$ independent of $N$, such that

$\mathbb{E}\left[\left|\mathscr{L}(\mu^{N},\psi,\left\langle\cdot,\cdot\right\rangle)\right|\right]\leq\frac{C}{\sqrt{N}}.$

###### Proof.

Applying Itô’s formula to the process
$\frac{1}{N}\sum_{i=1}^{N}\psi(X^{(i)},t)$, we get that

$\mathscr{L}(\mu^{N},\psi,\left\langle\cdot,\cdot\right\rangle)=\frac{1}{N}\sum_{i=1}^{N}\int_{0}^{T}\nabla\psi(X^{(i)}_{t},t)^{T}\sigma^{\star}(X^{(i)}_{t})dB^{(i)}_{t}.$

Note that each integral on the right-hand side is a local martingale, since
(4.8) ensures boundedness of $\nabla\psi(x,t)^{T}\sigma^{\star}(x)$ over any
compact set in $\mathbb{R}^{d}$, hence has mean zero. By independence of the
Brownian motions $B^{(i)}_{t}$, exchangeability of $X^{(i)}_{t}$, the moment
bound (4.7), and the growth bounds on $\sigma$, the Itô isometry gives us

$\displaystyle\mathbb{E}\left[\mathscr{L}(\mu^{N},\psi,\left\langle\cdot,\cdot\right\rangle)^{2}\right]$
$\displaystyle=\frac{1}{N}\int_{0}^{T}\mathbb{E}_{X\sim\rho^{(1)}}\left[\left|\nabla\psi(X,t)^{T}\sigma^{\star}(X)\right|^{2}\right]dt$
$\displaystyle=\frac{1}{N}\int_{0}^{T}\mathbb{E}\left[\int_{\mathbb{R}^{d}}\left|\nabla\psi(x,t)^{T}\sigma^{\star}(x)\right|^{2}d\mu^{N}_{t}(x)\right]dt$
$\displaystyle\leq\frac{C^{\prime}}{N}\left\|{\nabla\psi}\right\|^{2}_{2,\infty}\int_{0}^{T}\mathbb{E}\left[1+\int_{\mathbb{R}^{d}}|x|^{p}d\mu^{N}_{t}(x)\right]dt$
$\displaystyle\leq CN^{-1}$

where $C$ depends on $M_{p}$, $C_{p}$, $T$ and $\psi$. Above, $\rho^{(1)}$ is
the $X^{(1)}_{t}$-marginal of the process $\mathbf{X}_{t}\in\mathbb{R}^{dN}$,
and $\left\|{f(x)}\right\|_{p,q}$ for vector-valued functions $f$ denotes the
$L^{q}$ norm over $x$ of the $\ell^{p}$ norm of $f(x)$. The result follows
from Jensen’s inequality. ∎

With the following lemma, we can relate the histogram $U$ to the empirical
measure $\mu^{N}$ through $\mathscr{L}$ using the inner product
$\left\langle\cdot,\cdot\right\rangle_{h}$ defined by trapezoidal-rule
integration in space and continuous integration in time.

######  Lemma 4.3.

For $C$ independent of $N$ and $h$, it holds that

$\mathbb{E}\left[|\mathscr{L}(U,\psi,\left\langle\cdot,\cdot\right\rangle_{h})-\mathscr{L}(\mu^{N},\psi,\left\langle\cdot,\cdot\right\rangle)|\right]\leq
Ch.$

###### Proof.

Using the notation $f^{\mathbf{C}}$ from Lemma 4.1 to denote piecewise
constant approximation of a function $f$ over the domain ${\mathcal{D}}$ using
the grid $\mathbf{C}$, we have

$\displaystyle\mathscr{L}(U,\psi,\left\langle\cdot,\cdot\right\rangle_{h})-\mathscr{L}(\mu^{N},\psi,\left\langle\cdot,\cdot\right\rangle)$
$\displaystyle=-\underbrace{\Big{(}\left\langle(\nabla\psi\cdot((\nabla
K^{\star})^{\mathbf{C}}*\mu^{N}))^{\mathbf{C}},\mu^{N}\right\rangle-\left\langle\nabla\psi\cdot\nabla
K^{\star}*\mu^{N},\mu^{N}\right\rangle\Big{)}}_{E_{\text{interact}}}$
$\displaystyle\qquad+\left\langle\partial_{t}\psi^{\mathbf{C}}-\partial_{t}\psi,\mu^{N}\right\rangle-\left\langle((\nabla\psi\cdot\nabla
V^{\star})^{\mathbf{C}}-\nabla\psi\cdot\nabla V^{\star},\mu^{N}\right\rangle$
$\displaystyle\qquad+\frac{1}{2}\left\langle\text{Tr}\left(\nabla^{2}\psi\sigma^{\star}(\sigma^{\star})^{T}\right)^{\mathbf{C}}-\text{Tr}\left(\nabla^{2}\psi\sigma^{\star}(\sigma^{\star})^{T}\right),\mu^{N}\right\rangle$
$\displaystyle=E_{\text{interact}}+E_{\text{linear}}.$

The right-hand side includes an interaction error $E_{\text{interact}}$
followed by a sum $E_{\text{linear}}$ of terms that are linear in the
difference between a locally Lipschitz function and its piecewise constant
approximation. Hence, we can bound $E_{\text{linear}}$ using smoothness of
$\psi$, the moment assumptions on $\mu^{N}_{t}$ and the growth assumptions on
$V$ and $\sigma$. Specifically, for $x\in B_{k}$ with center $\mathbf{c}_{k}$,
the growth assumptions imply

$\displaystyle|\nabla\psi(x)\cdot\nabla
V^{\star}(x)-\nabla\psi(\mathbf{c}_{k})\cdot\nabla V^{\star}(\mathbf{c}_{k})|$
$\displaystyle\leq
Ch\Big{(}(\left\|{\nabla\psi}\right\|_{2,\infty}+\text{Lip}(\nabla\psi))(1+|x|^{p})\Big{)}$
$\displaystyle|\text{Tr}\left(\nabla^{2}\psi(x)\sigma^{\star}(x)(\sigma^{\star}(x))^{T}\right)-\text{Tr}\left(\nabla^{2}\psi(\mathbf{c}_{k})\sigma^{\star}(\mathbf{c}_{k})(\sigma^{\star}(\mathbf{c}_{k}))^{T}\right)|$
$\displaystyle\leq
C^{\prime}h\Big{(}(\left\|{\nabla^{2}\psi}\right\|_{F,\infty}+\text{Lip}(\nabla^{2}\psi))(1+|x|^{p})\Big{)}$

for $C,C^{\prime}$ depending on $p,d$ and $C_{p}$, hence

$|E_{\text{linear}}|\leq C^{\prime\prime}\sup_{|\alpha|\leq
3}\text{Lip}(\partial^{\alpha}\psi)\left(T+\int_{0}^{T}\int_{\mathbb{R}^{d}}|x|^{p}d\mu^{N}_{t}dt\right)h.$
(4.10)

Similarly, for the interaction error we use that for $x\in B_{k}$ and $y\in
B_{j}$ with centers $\mathbf{c}_{k}$ and $\mathbf{c}_{j}$, we have

$\displaystyle\left|\nabla\psi(\mathbf{c}_{k})\cdot\nabla
K^{\star}(\mathbf{c}_{k}-\mathbf{c}_{j})-\nabla\psi(x)\cdot\nabla
K^{\star}(x-y)\right|$
$\displaystyle\leq|\nabla\psi(\mathbf{c}_{k})|\left|\nabla
K^{\star}(\mathbf{c}_{k}-\mathbf{c}_{j})-\nabla K^{\star}(x-y)\right|$
$\displaystyle\qquad+\left|\nabla\psi(\mathbf{c}_{k})-\nabla\psi(x)\right|\left|\nabla
K^{\star}(x-y)\right|$ $\displaystyle\leq
C^{\prime\prime\prime}h\left(\left\|{\nabla\psi}\right\|_{2,\infty}+\text{Lip}(\nabla\psi)\right)(1+|x-y|^{p})$

with $C^{\prime\prime\prime}$ also depending on $p$, $d$, and $C_{p}$. From
this we have

$|E_{\text{interact}}|\leq
C^{\prime\prime\prime\prime}\left(T+\int_{0}^{T}\int_{\mathbb{R}^{d}}\int_{\mathbb{R}^{d}}|x-y|^{p}d\mu^{N}_{t}(y)d\mu^{N}_{t}(x)dt\right)h.$
(4.11)

The result follows from taking expectation and using the moment bound (4.7),
where the final constant $C$ depends on $p,d,C_{p},M_{p},T,\eta$ and $\psi$.

∎

To incorporate discrete effects, we consider the difference between
$\mathscr{L}(U,\psi,\left\langle\cdot,\cdot\right\rangle_{h})$ and
$\mathscr{L}(U,\psi,\left\langle\cdot,\cdot\right\rangle_{h,\Delta t})$, where
recall that $\left\langle\cdot,\cdot\right\rangle_{h,\Delta t}$ denotes
trapezoidal rule integration in space with meshwidth $h$ and in time with
sampling rate $\Delta t$.

######  Lemma 4.4.

For $C$ independent of $N$, $h$ and $\Delta t$, it holds that

$\mathbb{E}\left[|\mathscr{L}(U,\psi,\left\langle\cdot,\cdot\right\rangle_{h})-\mathscr{L}(U,\psi,\left\langle\cdot,\cdot\right\rangle_{h,\Delta
t})|\right]\leq C(h+\Delta t^{\eta}).$

###### Proof.

Again rewriting the spatial trapezoidal-rule integration in the form
$\int_{\mathbb{R}^{d}}\varphi^{\mathbf{C}}(x)d\mu^{N}_{t}$, we see that

$\mathscr{L}(U,\psi,\left\langle\cdot,\cdot\right\rangle_{h})-\mathscr{L}(U,\psi,\left\langle\cdot,\cdot\right\rangle_{h,\Delta
t})$ (4.12)

reduces to four terms of the form

$A(\varphi):=\frac{1}{N}\sum_{i=1}^{N}\left(\int_{0}^{T}\varphi^{\mathbf{C}}(X^{(i)}_{t})dt-\frac{\Delta
t}{2}\sum_{\ell=1}^{L}\left(\varphi^{\mathbf{C}}(X^{(i)}_{t_{\ell+1}})+\varphi^{\mathbf{C}}(X^{(i)}_{t_{\ell}})\right)\right),$

for $\varphi\in\left\\{\partial_{t}\psi,\nabla\psi\cdot\nabla
V^{\star},\text{Tr}(\nabla^{2}\psi\sigma^{\star}(\sigma^{\star})^{T}),\nabla\psi\cdot\nabla
K^{\star}*\mu^{N}_{t}\right\\}$. Similarly to the bounds derived for
$|\varphi(x)-\varphi^{\mathbf{C}}(x)|$ in Lemma 4.3, the growth bounds on
$V^{\star},K^{\star}$ and $\sigma^{\star}$ imply in general that

$|\varphi(x)-\varphi(y)|\leq C|x-y|\left(1+\max\\{|x|,|y|\\}^{p}\right).$

Rewriting the summands in $A(\varphi)$,

$\int_{0}^{T}\varphi^{\mathbf{C}}(X^{(i)}_{t})dt-\frac{\Delta
t}{2}\sum_{\ell=1}^{L}\left(\varphi^{\mathbf{C}}(X^{(i)}_{t_{\ell+1}})+\varphi^{\mathbf{C}}(X^{(i)}_{t_{\ell}})\right)$
$=\sum_{\ell=1}^{L}\underbrace{\int_{t_{\ell}}^{t_{\ell+1}}\left(\frac{t-t_{\ell}}{\Delta
t}\right)(\varphi^{\mathbf{C}}(X^{(i)}_{t})-\varphi^{\mathbf{C}}(X^{(i)}_{t_{\ell+1}}))dt}_{I_{1}}+\underbrace{\int_{t_{\ell}}^{t_{\ell+1}}\left(\frac{t_{\ell+1}-t}{\Delta
t}\right)(\varphi^{\mathbf{C}}(X^{(i)}_{t})-\varphi^{\mathbf{C}}(X^{(i)}_{t_{\ell}}))dt}_{I_{2}},$

and using

$|\varphi^{\mathbf{C}}(x)-\varphi^{\mathbf{C}}(y)|\leq|\varphi(x)-\varphi(\mathbf{c}_{k})|+|\varphi(x)-\varphi(y)|+|\varphi(y)-\varphi(\mathbf{c}_{\ell})|\leq
C(2h+|x-y|)(1+\max\\{|x|,|y|\\}^{p})$

where $x\in B_{k}$ and $y\in B_{\ell}$, we see that for $I_{1}$,

$\left|\int_{t_{\ell}}^{t_{\ell+1}}\left(\frac{t-t_{\ell}}{\Delta
t}\right)(\varphi^{\mathbf{C}}(X^{(i)}_{t})-\varphi^{\mathbf{C}}(X^{(i)}_{t_{\ell+1}}))dt\right|$
$\leq\int_{t_{\ell}}^{t_{\ell+1}}\left(\frac{t-t_{\ell}}{\Delta
t}\right)C(2h+|X_{t}^{(i)}-X^{(i)}_{t_{\ell+1}}|)(1+\max\\{|X^{(i)}_{t}|,|X^{(i)}_{t_{\ell+1}}|\\}^{p})dt$
$\leq\int_{t_{\ell}}^{t_{\ell+1}}\left(\frac{t-t_{\ell}}{\Delta
t}\right)C^{\prime}(2h+|t_{\ell+1}-t|^{\eta}|)(1+\max\\{|X^{(i)}_{t}|,|X^{(i)}_{t_{\ell+1}}|\\}^{p})dt.$

Taking expectation on both sides and using the moment bound (4.7), we get

$\mathbb{E}\left[\left|\int_{t_{\ell}}^{t_{\ell+1}}\left(\frac{t-t_{\ell}}{\Delta
t}\right)(\varphi^{\mathbf{C}}(X^{(i)}_{t})-\varphi^{\mathbf{C}}(X^{(i)}_{t_{\ell+1}}))dt\right|\right]\leq
C\left(\Delta th+\Delta t^{1+\eta}\right).$

We get the same bound for $I_{2}$. Summing over $\ell$, and taking the average
in $i$, we then get

$\mathbb{E}\left[|A(\varphi)|\right]\leq C(h+\Delta t^{\eta}),$

which implies the desired bound on the difference (4.12). ∎

The previous estimates directly lead to the following bound on the model
coefficients ${\widehat{\mathbf{w}}}$:

######  Theorem 4.1.

Let ${\widehat{\mathbf{w}}}=\mathbf{G}^{\dagger}\mathbf{b}$ be the learned
model coefficients and $\mathbf{w}^{\star}$ the true model coefficients. For
$C$ independent of $N,h$ and $\Delta t$ it holds that

$\mathbb{E}\left[\left\|{{\widehat{\mathbf{w}}}-\mathbf{w}^{\star}}\right\|_{1}\right]\leq
C\left(h+N^{-1/2}+\Delta t^{\eta}\right).$

###### Proof.

Using that $K^{\star}$, $V^{\star}$ and $\sigma^{\star}$ are in the span of
$\mathbb{L}$, we have that

$\mathbf{b}_{k}=\left\langle\partial_{t}\psi_{k},\mathbf{U}\right\rangle_{h,\Delta
t}=\mathscr{L}(U,\psi_{k},\left\langle\cdot,\cdot\right\rangle_{h,\Delta
t})+\mathbf{G}^{T}_{k}\mathbf{w}^{\star}:=\mathbf{L}_{k}+\mathbf{G}^{T}_{k}\mathbf{w}^{\star},$

where $\mathbf{G}^{T}_{k}$ is the $k$th row of $\mathbf{G}$. From the previous
lemmas, we have

$\mathbb{E}\left[|\mathbf{L}_{k}|\right]\leq\mathbb{E}\left[|\mathscr{L}(U,\psi_{k},\left\langle\cdot,\cdot\right\rangle_{h,\Delta
t})-\mathscr{L}(U,\psi_{k},\left\langle\cdot,\cdot\right\rangle_{h})|\right]$
$+\mathbb{E}\left[|\mathscr{L}(U,\psi_{k},\left\langle\cdot,\cdot\right\rangle_{h})-\mathscr{L}(\mu^{N},\psi_{k},\left\langle\cdot,\cdot\right\rangle)|\right]+\mathbb{E}\left[|\mathscr{L}(\mu^{N},\psi_{k},\left\langle\cdot,\cdot\right\rangle)|\right]$
$\leq C^{\prime}\left(h+N^{-1/2}+\Delta t^{\eta}\right).$

Using that $\mathbf{G}$ is full rank, it holds that
${\widehat{\mathbf{w}}}=\mathbf{G}^{\dagger}\mathbf{b}=\mathbf{G}^{\dagger}\mathbf{L}+\mathbf{w}^{\star}$,
hence the result follows from the uniform bound on
$\left\|{\mathbf{G}^{\dagger}}\right\|_{1}$:

$\mathbb{E}\left[\left\|{{\widehat{\mathbf{w}}}-\mathbf{w}^{\star}}\right\|_{1}\right]\leq\left\|{\mathbf{G}^{\dagger}}\right\|_{2}\mathbb{E}\left[\left\|{\mathbf{L}}\right\|_{1}\right]\leq
C^{\prime}n\left\|{\mathbf{G}^{\dagger}}\right\|_{1}\left(h+N^{-1/2}+\Delta
t^{\eta}\right).$

∎

Under the assumption that $K^{\star},V^{\star}$ and $\sigma^{\star}$ are
contained in the span of $\mathbb{L}$, an immediate corollary is

$\mathbb{E}\left[\left\|{K^{\star}-\widehat{K}}\right\|_{L^{2}({\mathcal{D}}-{\mathcal{D}})}+\left\|{V^{\star}-\widehat{V}}\right\|_{L^{2}({\mathcal{D}})}+\left\|{\left\|{\sigma^{\star}(\sigma^{\star})^{T}-\widehat{\sigma}(\widehat{\sigma})^{T}}\right\|_{F}}\right\|_{L^{2}({\mathcal{D}})}\right]\leq
C\left(h+N^{-1/2}+\Delta t^{\eta}\right).$

Finally, setting $h=N^{-\alpha}$ for $\alpha>0$ will ensure convergence as
$N\to\infty$ and $\Delta t\to 0$. We now make several remarks about the
practical (algorithmic) implementation with respect to this theoretical
convergence.

######  Remark 4.2.

blah

  * •

An important case of Theorem 4.1 is $\sigma^{\star}=0$, in which case
$\mu^{N}_{t}$ itself is a weak-measure solution to the mean-field equation
(3.1) and $\left\|{{\widehat{\mathbf{w}}}-\mathbf{w}^{\star}}\right\|_{2}\leq
C(h+\Delta t^{\eta})$, with $\eta\geq 2$. Although the examples below only
explore $\sigma^{\star}=0$ with nonzero extrinsic noise $\varepsilon$ (Figures
5 and 9), we note that when $\varepsilon=0$ and $\sigma^{\star}=0$ (not shown)
Algorithm 4.1 recovers systems to high accuracy similarly to WSINDy applied to
local dynamical systems [40, 39].

  * •

Algorithm 4.1 in general implements sparse regression, yet Theorem 4.1 deals
with ordinary least squares. In practice, sparse regression provides an
improvement over least squares since $\mathbf{G}$ typically has a high
condition number (e.g. due to coupled effects of $K$,$V$, and $\sigma$). Since
least squares is a common subroutine of many sparse regression algorithms
(inluding the MSTLS algorithm used here), the result is still relevant to
sparse regression. Lastly, the full-rank assumption on $\mathbf{G}$ implies
that as $N\to\infty$ sequential thresholding reduces to least squares.

  * •

Theorem 4.1 assumes data from a single experiment ($M=1$), while the examples
below show that $M>1$ experiments improves results. For any fixed $M>1$, the
$N\to\infty$ limit results in convergence, however, the $N$-fixed and
$M\to\infty$ limit does not result in convergence, as this does not lead to
the mean-field equations111111Note that the opposite convergence holds for the
algorithm introduced in [35]: $N$-fixed, $M\to\infty$ results in recovery of
$K$.. The examples below show that using $M>1$ has a practical advantage.

  * •

Many interesting examples have non-Lipschitz $\nabla K$, in particular a lack
of smoothness at $x=0$. If $\nabla K$ has a jump discontinuity at $x=0$ and
$\mu^{N}_{t}$ does not concentrate to a singular measure as $N\to\infty$, then
the bound (4.11) may be modified to include another ${\mathcal{O}}(h)$ term
coming from short-range interactions (i.e. within an ${\mathcal{O}}(h)$
distance). The examples below are chosen in part to show that
${\mathcal{O}}(N^{-1/2})$ convergence holds for $\nabla K$ with jumps at the
origin.

##  5\. Examples

We now demonstrate the successful identification of several particle systems
in one and two spatial dimensions as well as the ${\mathcal{O}}(N^{-1/2})$
convergence predicted in Theorem 4.1. In each case we use Algorithm 4.1 to
discover a mean-field equation of the form (3.1) from discrete-time particle
data. For each dataset we simulate the associated interacting particle system
$\mathbf{X}_{t}$ given by (1.1) using the Euler-Maruyama scheme (initial
conditions and timestep are given in each example). We assess the ability of
WSINDy to select the correct model using the true positivity ratio121212For
example, identification of the true model
($\text{supp}\left({\widehat{\mathbf{w}}}\right)=\text{supp}\left(\mathbf{w}^{\star}\right)$)
results in a TPR(${\widehat{\mathbf{w}}})=1$, while identification of only
half of the correct nonzero terms and no additional falsely identified terms
results in TPR(${\widehat{\mathbf{w}}})=0.5$.

$\text{TPR}({\widehat{\mathbf{w}}})=\frac{\text{TP}}{\text{TP}+\text{FN}+\text{FP}}$
(5.1)

where TP is the number of correctly identified nonzero coefficients, FN is the
number of coefficients falsely identified as zero, and FP is the number of
coefficients falsely identified as nonzero [30]. To demonstrate the
${\mathcal{O}}(N^{-1/2})$ convergence, for correctly identified models (i.e.
TPR$({\widehat{\mathbf{w}}})=1$) we compute the relative $\ell_{2}$ error in
the recovered interaction force $\nabla\widehat{K}$, local force
$\nabla\widehat{V}$, and diffusivity $\widehat{\sigma}$ over
$\mathbf{C}-\mathbf{C}$ and $\mathbf{C}$, respectively. Results are averaged
over 100 trials.

For the computational grid $\mathbf{C}$ we first compute the sample standard
deviation $s$ of $\boldsymbol{\mathbb{Y}}$ and we choose ${\mathcal{D}}$ to be
the rectangular grid extending $3s$ from the mean of $\boldsymbol{\mathbb{Y}}$
in each direction. We then set $\mathbf{C}$ to have 128 points in $x$ and $y$
for $d=2$ dimensions, and 256 points in $x$ for $d=1$, noting that these
numbers are fairly arbitrary, and used to show that the grid need not be too
large. We set the sparsity factors so that $\log_{10}(\boldsymbol{\lambda})$
contains 100 equally spaced points from $-4$ to 0. More information on the
specifications of each example can be found in Appendix A.

###  5.1. Two-Dimensional Local Model

The first system we examine is a constant advection / variable diffusivity
model with mean-field equation131313Since the model is local, (5.2) is the
Fokker-Planck equation for the distribution of each particle, rather than only
in the limit of infinite particles.

$\partial_{t}\mu_{t}=-\partial_{x}\mu_{t}-\partial_{y}\mu_{t}+\Delta\left[\left(1+0.95\cos(\omega
x)\cos(\omega y)\right)\mu_{t}\right].$ (5.2)

The purpose of this example is three-fold. First, we are interested in the
ability of Algorithm 4.1 to correctly identify a local model from a library
containing both local and nonlocal terms. Next, we evaluate whether the
${\mathcal{O}}(N^{-1/2})$ convergence is realized. Lastly, we investigate
whether for large $\omega$ the weak-form identifies the associated homogenized
equation (see e.g. [52])

$\partial_{t}\mu_{t}=-\partial_{x}\mu_{t}-\partial_{y}\mu_{t}+\overline{\omega}\Delta\mu_{t},$
(5.3)

where $\overline{\omega}$ is given by the harmonic mean of diffusivity:

$\overline{\omega}=\left(\int_{\mathcal{D}}\frac{dxdy}{1+0.95\cos(x)\cos(y)}\right)^{-1}.$

For $\omega\in\\{1,20\\}$ we evolve the particles from an initial Gaussian
distribution with mean zero and covariance $\mathbf{I}_{2}$ and record
particle positions for $100$ timesteps with $\Delta t=0.02$ (subsampled from a
simulation with timestep $10^{-4}$). We use a rectangular domain
${\mathcal{D}}$ of approximate sidelength $10$ and compute histograms with 128
bins in $x$ and $y$ for a spatial resolution of $\Delta x\approx 0.078$ (see
Figure 2 for solution snapshots), over which $\overline{\omega}\approx 0.62$.
For $\omega=1$ we compare recovered equations with the full model (5.2), while
for $\omega=20$ we compare with (5.3), for comparison computing
$\overline{\omega}$ over each domain ${\mathcal{D}}$ using MATLAB’s integral2.
Figure 3 shows that as the particle number increases, we do in fact recover
the desired equations, with TPR$({\widehat{\mathbf{w}}})$ approaching one as
$N$ increases. For $\omega=1$ we observe ${\mathcal{O}}(N^{-1/2})$ convergence
of the local potential $\widehat{V}$ and the diffusivity $\widehat{\sigma}$.
For $\omega=20$, we observe approximate ${\mathcal{O}}(N^{-1/2})$ convergence
of $\widehat{V}$, and $\widehat{\sigma}$ converging to within $2\%$ of
$\sqrt{2\overline{\omega}}$, the homogenized diffusivity (higher accuracy can
hardly be expected for $\omega=20$ since (5.3) is itself an approximation in
the limit of infinite $\omega$).

|
---|---
|
Figure 2. Snapshots at time $t=2\Delta t=0.06$ (left) and $t=100\Delta t=2$ (right) of histograms computed with 128 bins in $x$ and $y$ from 16,384 particles evolving under (5.2) with $\omega=1$ (top) and $\omega=20$ (bottom). |  |
---|---|---
Figure 3. Convergence of $\widehat{\sigma}$ (left) and $\nabla\widehat{V}$
(middle) for (5.2) with $\omega\in\\{1,20\\}$, as well as
TPR$({\widehat{\mathbf{w}}})$ (right). For $\omega=1$, comparison is made with
respect to the exact model (5.2), while for $\omega=20$ results are compared
to (5.3).

###  5.2. One-Dimensional Nonlocal Model

|  |
---|---|---
|  |
Figure 4. Histograms computed with 256 bins width $h=0.0234$ from 8000
particles in 1D evolving under $K^{\star}=K_{\text{QANR}}(x)$ (5.4). Top left
to top right: $\sigma^{\star}(x)=0$, $\sigma^{\star}(x)=\sqrt{2(0.1)}$,
$\sigma^{\star}(x)=\sqrt{2(0.1)}|x-2|$. Bottom: deterministic particles with
i.i.d. Gaussian noise added to particle positions with resulting noise ratios
(left to right) $\epsilon=0.0316,0.1,0.316$.

We simulate the evolution of particle systems under the quadratic attraction /
Newtonian repulsion potential

$K_{\text{QANR}}(x)=\frac{1}{2}x^{2}-|x|$ (5.4)

with no external potential $(V=0)$. The $-|x|$ portion of $K_{\text{QANR}}$ ,
leading to a discontinuity in $\nabla K$, is the one-dimensional free-space
Green’s function for $-\Delta$. For $d\geq 1$, when replaced by the
corresponding Green’s function in $d$ dimensions, the distribution of
particles evolves under $K_{\text{QANR}}$ into the characteristic of the unit
ball in $\mathbb{R}^{d}$, which has implications for design and control of
autonomous systems [20]. We compare three diffusivity profiles, $\sigma(x)=0$
corresponding to zero intrinsic noise, $\sigma(x)=\sqrt{2(0.1)}$ leading to
constant-diffusivity intrinsic noise, and $\sigma(x)=\sqrt{2(0.1)}|x-2|$
leading to variable-diffusivity intrinsic noise. With zero intrinsic noise
($\sigma(x)=0$), we examine the effect of extrinsic noise on recovery, and
assume uncertainty in the particle positions due to measurement noise at each
timestep, $\boldsymbol{\mathbb{Y}}=\boldsymbol{\mathbb{X}}+\varepsilon$, for
$\varepsilon\sim{\mathcal{N}}(0,\epsilon^{2}\left\|{\mathbf{X}_{\mathbf{t}}}\right\|_{\text{RMS}}^{2})$
i.i.d. and $\epsilon\in\\{0.01,0.0316,0.1,0.316\\}$. In this way $\epsilon$ is
the noise ratio, such that
$\left\|{\varepsilon}\right\|_{F}/\left\|{\boldsymbol{\mathbb{X}}}\right\|_{F}\approx\epsilon$
(computed with $\varepsilon$ and $\boldsymbol{\mathbb{X}}$ stretched into
column vectors).

Measurement data consists of 100 timesteps at resolution $\Delta t=0.01$,
coarsened from simulations with timestep $0.001$. Initial particle positions
are drawn from a mixture of three Gaussians each with standard deviation
$0.005$. Histograms are constructed with 256 bins of width $h=0.0234$. Typical
histograms for each noise level are shown in Figure 4 computed one experiment
with $N=8000$ particles.

For the case of extrinsic noise (Figure 5), we use only one experiment ($M=1$)
and examine the number of particles $N$ and the noise ratio $\epsilon$. We
find that recovery is accurate and reliable for $\epsilon\leq 0.1$, yielding
correct identification of $K_{\text{QANR}}$ with less than $1\%$ relative
error in at least $98/100$ trials. Increasing $N$ from 500 to 8000 leads to
minor improvements in accuracy for $\epsilon\leq 0.1$, but otherwise has
little effect, implying that for low to moderate noise levels the mean field
equations are readily identifiable even from smaller particle systems. For
$\epsilon=10^{-1/2}\approx 0.3162$ (see Figure 4 (bottom right) for an example
histogram), we observe a decrease in TPR(${\widehat{\mathbf{w}}}$) (Figure 5
middle panel) resulting from the generic identification of a linear diffusion
term $\nu\partial_{xx}u$ with $\nu\approx 0.05$. Using that
$\sqrt{2\nu}\approx\sqrt{2(0.05)}=\epsilon$, we can identify this as the best-
fit intrinsic noise model. Furthermore, increases in $N$ lead to reliable
identification of the drift term, as measured by
TPR(${\widehat{\mathbf{w}}}_{drift}$) (rightmost panel Figure 5) which is the
restriction of TPR to drift terms $\mathbb{L}_{K}$ and $\mathbb{L}_{V}$.

For constant diffusivity $\sigma(x)=\sqrt{2(0.1)}$ (Figure 6), the full model
is recovered with less than $3\%$ errors in $\hat{K}$ and $\hat{\sigma}$ in at
least 98/100 trials when the total particle count $NM$ is at least $8000$, and
yields errors less than $1\%$ for $NM\geq$ 16,000. The error trends for
$\hat{K}$ and $\hat{\sigma}$ in this case both strongly agree with the
predicted ${\mathcal{O}}(N^{-1/2})$ rate. For non-constant diffusivity
$\sigma(x)=\sqrt{2(0.1)}|x-2|$ (Figure 7), we also observe robust recovery
(TPR$({\widehat{\mathbf{w}}})\geq 0.95$) for $NM\geq 8000$ with error trends
close to ${\mathcal{O}}(N^{-1/2})$, although the accuracy in $\hat{K}$ and
$\hat{\sigma}$ is diminished due to the strong order $\Delta t^{1/2}$
convergence of Euler-Maruyama applied to diffusivities $\sigma$ that are
unbounded in $x$ [42].

|  |
---|---|---
Figure 5. Recovery of (3.1) in one spatial dimension for $K^{\star}=K_{\text{QANR}}$ and $\sigma^{\star}=0$ under different levels of observational noise $\epsilon$. Left: relative error in learned interaction kernel $\widehat{K}$. Middle: true positivity ratio for full model (3.1). Right: true positivity ratio for drift term. |  |
---|---|---
Figure 6. Recovery of (3.1) in one spatial dimension for $K^{\star}=K_{\text{QANR}}$ and $\sigma^{\star}=\sqrt{2(0.1)}$ |  |
---|---|---
Figure 7. Recovery of (3.1) in one spatial dimension for
$K^{\star}=K_{\text{QANR}}$ and $\sigma^{\star}=\sqrt{2(0.1)|x-2|}$

###  5.3. Two-Dimensional Nonlocal Model

|  |
---|---|---
|  |
Figure 8. Histograms created from 4000 particles evolving under logarithmic
attraction (equation (5.5)) with varying noise levels at times (left to right)
$t=4$, $t=8$, and $t=12$. Top: $\epsilon=0.316$, $\sigma=0$ (extrinsic only).
Bottom: $\epsilon=0$, $\sigma=(4\pi)^{-1/2}\approx 0.28$ (intrinsic only).

We now discuss an example of singular interaction in two spatial dimensions
using the logarithmic potential

$K(x)=\frac{1}{2\pi}\log|x|$ (5.5)

with constant diffusivity $\sigma(x)=\sigma\in\\{0,\frac{1}{\sqrt{4\pi}}\\}$.
This example corresponds to the parabolic-elliptic Keller-Segel model of
chemotaxis, where $\sigma_{c}:=\frac{1}{\sqrt{4\pi}}$ is the critical
diffusivity such that $\sigma>\sigma_{c}$ leads diffusion-dominated spreading
of particles throughout the domain (vanishing particle density at every point
in $\mathbb{R}^{2}$) and $\sigma<\sigma_{c}$ leads to aggregation-dominated
concentration of the particle density to the dirac-delta located at the center
of mass of the initial particle density [16, 13]. For $\sigma=0$ we examine
the affect of additive i.i.d. measurement noise
$\varepsilon\sim{\mathcal{N}}(0,\epsilon^{2}\left\|{\mathbf{X}_{\mathbf{t}}}\right\|_{\text{RMS}}^{2})$
for $\epsilon\in\\{0.01,0.0316,0.1,0.316,1\\}$.

We simulate the particle system with a cutoff potential

$K_{\delta}(x)=\begin{dcases}\frac{1}{2\pi}\left(\log(\delta)-1+\frac{|x|}{\delta}\right),&|x|<\delta\\\
\frac{1}{2\pi}\log|x|,&|x|\geq\delta\end{dcases}$ (5.6)

with $\delta=0.01$, so that $K_{\delta}$ is Lipschitz and $\nabla K_{\delta}$
has a jump discontinuity at the origin. Initial particle positions are
uniformly distributed on a disk of radius 2 and the particle position data
consists of $81$ timepoints recorded at a resolution $\Delta t=0.1$, coarsened
from $0.0025$. Histograms are created with $128\times 128$ bins in $x$ and $y$
of sidelength $h=0.0469$ (see Figure 8 for histogram snapshots over time). We
examine $M=2^{0},\dots,2^{6}$ experiments with $N=2000$ or $N=4000$ particles.

|
---|---
Figure 9. Recovery of (3.1) in two spatial dimensions with $K^{\star}$ given
by (5.5) from deterministic particles ($\sigma^{\star}=0$) with extrinsic
noise $\epsilon$.

In Figure 9 we observe a similar trend in the $\sigma=0$ case as in the 1D
nonlocal example, namely that recovery for $\epsilon\leq 0.1$ is robust with
low errors in $\widehat{K}$ (on the order of $0.0032$), only in this case the
full model is robustly recovered up to $\epsilon=0.316$. At $\epsilon=1$, with
$N=4000$ the method frequently identifies a diffusion term $\nu\Delta u$ with
$\nu\approx 0.5=\epsilon^{2}/2$, and for $N=2000$ the method occasionally
identifies the backwards diffusion equation
$\partial_{t}\mu_{t}=-\alpha\Delta\mu_{t}$, $\alpha>0$. This is easily
prevented by enforcing positivity, which we leave as an extension for future
work.

|  |
---|---|---
Figure 10. Recovery of (3.1) in two spatial dimensions with $K^{\star}$ given
by (5.5) and $\sigma^{\star}=\frac{1}{\sqrt{4\pi}}$.

With diffusivity $\sigma=\frac{1}{\sqrt{4\pi}}$, we obtain
TPR$({\widehat{\mathbf{w}}})$ approximately greater than 0.95 for $NM\geq
16,000$ (Figure 10, right), with an error trend in $\widehat{K}$ following an
${\mathcal{O}}(N^{-1/2})$ rate, and a trend in $\widehat{\sigma}$ of roughly
${\mathcal{O}}(N^{-2/3})$. Since convergence in $M$ for any fixed $N$ is not
covered by the theorem above, this shows that combining multiple experiments
may yield similar accuracy trends for moderately-sized particle systems.

##  6\. Discussion

We have developed a weak-form method for sparse identification of governing
equations for interacting particle systems using the formalism of mean-field
equations. In particular, we have investigating two lines of inquiry, (1) is
the mean-field setting applicable for inference from medium-size batches of
particles? And (2) can a low-cost, low-regularity density approximation such
as a histogram be used to enforce weak-form agreement with the mean-field PDE?
We have demonstrated on several examples that the answer is yes to both
questions, despite the fact that the mean-field equations are only valid in
the limit of infinitely many particles ($N\to\infty$). This framework is
suitable for systems of several thousand particles in one and two spatial
dimensions, and we have proved convergence in $N$ for the associated least-
squares problem using simple histograms as approximate particle densities. In
addition, the sparse regression approach allows one to identify the full
system, including interaction potential $K$, local potential $V$, and
diffusivity $\sigma$.

It was initially unclear whether the mean-field setting could be utilized for
finite particle batches, hence this can be seen as a proof of concept, with
the potential for many improvements and extensions. On the subject of density
estimation, histograms lead to piecewise-constant approximations and resulting
${\mathcal{O}}(h)$ errors, hence choosing a density kernel $G$ to achieve
high-accuracy quadrature without sacrificing the ${\mathcal{O}}(N)$ runtime of
histogram computation seems prudent. The computational grid $\mathbf{C}$ is
also a free parameter, and may be optimized in tandem with a quadrature rule.
The equally-spaced approach combined with the trapezoidal rule, as applied
here, has several advantages, but may need adjustment for higher dimensions.
Another obvious improvement would be to enforce convex constraints in the
regression problem, such as lower bounds on diffusivity, or $K$ with long-
range attraction depending on the distribution
$\rho_{rr}\in{\mathcal{P}}([0,\infty))$ of pairwise distances (see [35] for
further use $\rho_{rr}$). For extensions, the example system (5.2) and
resulting homogenization motivates further study of effective equations for
systems with complex microstructure. In other fields this is described as
coarse-graining. A related line of study is inference of 2nd-order particle
systems, as explored in [48], which often lead to an infinite hierachy of
mean-field equations. Our weak-form approach may provide a principled method
for truncated and closing such hierarchies using particle data.

##  7\. Acknowledgements

This research was supported in part by the NSF Mathematical Biology MODULUS
grant 2054085, in part by the NSF/NIH Joint DMS/NIGMS Mathematical Biology
Initiative grant R01GM126559, and in part by the NSF Computing and
Communications Foundations grant 1815983. This work also utilized resources
from the University of Colorado Boulder Research Computing Group, which is
supported by the National Science Foundation (awards ACI-1532235 and
ACI-1532236), the University of Colorado Boulder, and Colorado State
University. The authors would also like to thank Prof. Vanja Dukić (University
of Colorado at Boulder, Department of Applied Mathematics) for insightful
discussions and helpful suggestions of references.

## References

  * [1] Dyego Araújo, Roberto I Oliveira, and Daniel Yukimura.  A mean-field limit for certain deep neural networks.  arXiv preprint arXiv:1906.00193, 2019.
  * [2] Dapeng Bi, Xingbo Yang, M Cristina Marchetti, and M Lisa Manning.  Motility-driven glass and jamming transitions in biological tissues.  Physical Review X, 6(2):021011, 2016.
  * [3] Bo Martin Bibby and Michael Sørensen.  Martingale estimation functions for discretely observed diffusion processes.  Bernoulli, pages 17–39, 1995.
  * [4] Jaya PN Bishwal.  Parameter estimation in stochastic differential equations.  Springer, 2007.
  * [5] Jaya Prakash Narayan Bishwal et al.  Estimation in interacting diffusions: Continuous and discrete sampling.  Applied Mathematics, 2(9):1154–1158, 2011.
  * [6] Vincent D Blondel, Julien M Hendrickx, and John N Tsitsiklis.  Continuous-time average-preserving opinion dynamics with opinion-dependent communications.  SIAM Journal on Control and Optimization, 48(8):5214–5240, 2010\.
  * [7] Niklas Boers and Peter Pickl.  On mean field limits for dynamical systems.  Journal of Statistical Physics, 164(1):1–16, 2016.
  * [8] François Bolley, José A Canizo, and José A Carrillo.  Stochastic mean-field limit: non-lipschitz forces and swarming.  Mathematical Models and Methods in Applied Sciences, 21(11):2179–2210, 2011.
  * [9] Mattia Bongini, Massimo Fornasier, Markus Hansen, and Mauro Maggioni.  Inferring interaction rules from observations of evolutive systems i: The variational approach.  Mathematical Models and Methods in Applied Sciences, 27(05):909–951, 2017.
  * [10] Lorenzo Boninsegna, Feliks Nüske, and Cecilia Clementi.  Sparse learning of stochastic dynamical equations.  The Journal of chemical physics, 148(24):241723, 2018.
  * [11] Steven L Brunton, Joshua L Proctor, and J Nathan Kutz.  Discovering governing equations from data by sparse identification of nonlinear dynamical systems.  Proceedings of the national academy of sciences, 113(15):3932–3937, 2016.
  * [12] Jared L Callaham, J-C Loiseau, Georgios Rigas, and Steven L Brunton.  Nonlinear stochastic modelling with Langevin regression.  Proceedings of the Royal Society A, 477(2250):20210092, 2021.
  * [13] JA Carrillo, MG Delgadino, and FS Patacchini.  Existence of ground states for aggregation-diffusion equations.  Analysis and applications, 17(03):393–423, 2019.
  * [14] Xiaohui Chen.  Maximum likelihood estimation of potential energy in interacting particle systems from single-trajectory data.  Electronic Communications in Probability, 26:1–13, 2021.
  * [15] Xiaoli Chen, Liu Yang, Jinqiao Duan, and George Em Karniadakis.  Solving inverse stochastic problems from discrete particle observations using the Fokker–Planck equation and physics-informed neural networks.  SIAM Journal on Scientific Computing, 43(3):B811–B830, 2021.
  * [16] Jean Dolbeault and Benoît Perthame.  Optimal critical mass in the two dimensional Keller–Segel model in r2.  Comptes Rendus Mathematique, 339(9):611–616, 2004.
  * [17] Jinchao Feng, Yunxiang Ren, and Sui Tang.  Data-driven discovery of interacting particle systems using gaussian processes.  arXiv preprint arXiv:2106.02735, 2021.
  * [18] Razvan C Fetecau, Hui Huang, Daniel Messenger, and Weiran Sun.  Zero-diffusion limit for aggregation equations over bounded domains.  arXiv preprint arXiv:1809.01763, 2018.
  * [19] Razvan C Fetecau, Hui Huang, and Weiran Sun.  Propagation of chaos for the Keller–Segel equation over bounded domains.  Journal of Differential Equations, 266(4):2142–2174, 2019.
  * [20] Razvan C Fetecau, Yanghong Huang, and Theodore Kolokolnikov.  Swarm dynamics and equilibria for a nonlocal aggregation model.  Nonlinearity, 24(10):2681, 2011.
  * [21] Razvan C Fetecau and Mitchell Kovacic.  Swarm equilibria in domains with boundaries.  SIAM Journal on Applied Dynamical Systems, 16(3):1260–1308, 2017\.
  * [22] David Freedman and Persi Diaconis.  On the histogram as a density estimator: L2 theory.  Zeitschrift für Wahrscheinlichkeitstheorie und verwandte Gebiete, 57(4):453–476, 1981.
  * [23] Paraskevi Gkeka, Gabriel Stoltz, Amir Barati Farimani, Zineb Belkacemi, Michele Ceriotti, John D Chodera, Aaron R Dinner, Andrew L Ferguson, Jean-Bernard Maillet, Hervé Minoux, et al.  Machine learning force fields and coarse-grained variables in molecular dynamics: application to materials and biological systems.  Journal of Chemical Theory and Computation, 16(8):4757–4775, 2020\.
  * [24] Susana N Gomes, Andrew M Stuart, and Marie-Therese Wolfram.  Parameter estimation for macroscopic pedestrian dynamics models from microscopic data.  SIAM Journal on Applied Mathematics, 79(4):1475–1500, 2019.
  * [25] Jiawei Guo.  The progress of three astrophysics simulation methods: Monte-carlo, pic and mhd.  In Journal of Physics: Conference Series, volume 2012, page 012136\. IOP Publishing, 2021.
  * [26] Pierre-Emmanuel Jabin and Zhenfu Wang.  Mean field limit for stochastic particle systems.  In Active Particles, Volume 1, pages 379–402. Springer, 2017.
  * [27] Jun-Gi Jang and U Kang.  D-tucker: Fast and memory-efficient tucker decomposition for dense tensors.  In 2020 IEEE 36th International Conference on Data Engineering (ICDE), pages 1850–1853. IEEE, 2020.
  * [28] Raphael A Kasonga.  Maximum likelihood theory for large interacting systems.  SIAM Journal on Applied Mathematics, 50(3):865–875, 1990.
  * [29] Evelyn F Keller and Lee A Segel.  Model for chemotaxis.  Journal of theoretical biology, 30(2):225–234, 1971.
  * [30] John H. Lagergren, John T. Nardini, G. Michael Lavigne, Erica M. Rutter, and Kevin B. Flores.  Learning partial differential equations for biological transport models from noisy spatio-temporal data.  Proc. R. Soc. A., 476(2234):20190800, February 2020.
  * [31] Quanjun Lang and Fei Lu.  Learning interaction kernels in mean-field equations of 1st-order systems of interacting particles.  arXiv preprint arXiv:2010.15694, 2020.
  * [32] Tony Lelievre and Gabriel Stoltz.  Partial differential equations and stochastic methods in molecular dynamics.  Acta Numerica, 25:681–880, 2016.
  * [33] Yang Li and Jinqiao Duan.  Extracting governing laws from sample path data of non-gaussian stochastic dynamical systems.  arXiv preprint arXiv:2107.10127, 2021.
  * [34] Andrew W Lo.  Maximum likelihood estimation of generalized itô processes with discretely sampled data.  Econometric Theory, 4(2):231–247, 1988.
  * [35] Fei Lu, Mauro Maggioni, and Sui Tang.  Learning interaction kernels in heterogeneous systems of agents from multiple trajectories.  J. Mach. Learn. Res., 22:32–1, 2021.
  * [36] Ryan Lukeman, Yue-Xian Li, and Leah Edelstein-Keshet.  Inferring individual rules from collective behavior.  Proceedings of the National Academy of Sciences, 107(28):12576–12580, 2010.
  * [37] Osman Asif Malik and Stephen Becker.  Low-rank tucker decomposition of large tensors using tensorsketch.  Advances in neural information processing systems, 31:10096–10106, 2018.
  * [38] Sylvie Méléard.  Asymptotic behaviour of some interacting particle systems; mckean-vlasov and boltzmann models.  In Probabilistic models for nonlinear partial differential equations, pages 42–95. Springer, 1996.
  * [39] Daniel A Messenger and David M Bortz.  Weak SINDy for partial differential equations.  Journal of Computational Physics, page 110525, 2021.
  * [40] Daniel A Messenger and David M Bortz.  Weak SINDy: Galerkin-based data-driven model selection.  Multiscale Modeling & Simulation, 19(3):1474–1497, 2021.
  * [41] Daniel A Messenger and Razvan C Fetecau.  Equilibria of an aggregation model with linear diffusion in domains with boundaries.  Mathematical Models and Methods in Applied Sciences, 30(04):805–845, 2020.
  * [42] Grigorii Noikhovich Milstein.  Numerical integration of stochastic differential equations, volume 313.  Springer Science & Business Media, 1994.
  * [43] John T Nardini, Ruth E Baker, Matthew J Simpson, and Kevin B Flores.  Learning differential equation models from stochastic agent-based model simulations.  Journal of the Royal Society Interface, 18(176):20200987, 2021.
  * [44] Samuel H Rudy, Steven L Brunton, Joshua L Proctor, and J Nathan Kutz.  Data-driven discovery of partial differential equations.  Science Advances, 3(4):e1602614, 2017.
  * [45] Néstor Sepúlveda, Laurence Petitjean, Olivier Cochet, Erwan Grasland-Mongrain, Pascal Silberzan, and Vincent Hakim.  Collective cell motion in an epithelial sheet can be quantitatively described by a stochastic interacting particle model.  PLoS computational biology, 9(3):e1002944, 2013.
  * [46] Louis Sharrock, Nikolas Kantas, Panos Parpas, and Grigorios A Pavliotis.  Parameter estimation for the mckean-vlasov stochastic differential equation.  arXiv preprint arXiv:2106.13751, 2021.
  * [47] Yiming Sun, Yang Guo, Charlene Luo, Joel Tropp, and Madeleine Udell.  Low-rank tucker approximation of a tensor from streaming data.  SIAM Journal on Mathematics of Data Science, 2(4):1123–1150, 2020\.
  * [48] Rohit Supekar, Boya Song, Alasdair Hastewell, Alexander Mietke, and Jörn Dunkel.  Learning hydrodynamic equations for active matter from particle simulations and experiments.  arXiv preprint arXiv:2101.06568, 2021.
  * [49] Alain-Sol Sznitman.  Topics in propagation of chaos.  In Ecole d’été de probabilités de Saint-Flour XIX—1989, pages 165–251. Springer, 1991.
  * [50] Paul Van Liedekerke, MM Palm, N Jagiella, and Dirk Drasdo.  Simulating tissue mechanics with agent-based models: concepts, perspectives and some novel results.  Computational particle mechanics, 2(4):401–444, 2015.
  * [51] Michael S Warren and John K Salmon.  Astrophysical n-body simulations using hierarchical tree data structures.  Proceedings of Supercomputing, 1992.
  * [52] E Weinan.  Principles of multiscale modeling.  Cambridge University Press, 2011.

##  Appendix A Notation & Specifications for Examples

Variable |  Definition |  Domain
---|---|---
$K$ |  pairwise interaction potential |  $L^{1}_{loc}(\mathbb{R}^{d},\mathbb{R})$
$V$ |  local potential |  $C(\mathbb{R}^{d},\mathbb{R})$
$\sigma$ |  diffusivity |  $C(\mathbb{R}^{d},\mathbb{R}^{d\times d})$
$N$ |  number of particles per experiment |  $\\{2,3,\dots\\}$
$d$ |  dimension of latent space |  $\mathbb{N}$
$T$ |  final time |  $(0,\infty)$
$(\Omega,{\mathcal{B}},\mathbb{P},({\mathcal{F}}_{t})_{t\geq 0})$ |  filtererd probability space |
$(B^{(i)}_{t})_{i=1}^{N}$ |  independent $\mathbb{R}^{d}$ Brownian motions on $(\Omega,{\mathcal{B}},\mathbb{P},({\mathcal{F}}_{t})_{t\geq 0})$ |
$X^{(i)}_{t}$ |  $i$th particle in the particle system (1.1) at time $t$ |  $\mathbb{R}^{d}$
$\mathbf{X}_{t}$ |  $N$-particle system (1.1) at time $t$ |  $\mathbb{R}^{Nd}$
$\mu^{N}_{t}$ |  empirical measure |  ${\mathcal{P}}(\mathbb{R}^{d})$
$F^{N}_{t}$ |  distribution of the process $\mathbf{X}_{t}$ in $\mathbb{R}^{Nd}$ |
$X_{t}$ |  mean-field process (3.2) at time $t$ |  $\mathbb{R}^{Nd}$
$\mu_{t}$ |  distribution of $X_{t}$ |  ${\mathcal{P}}(\mathbb{R}^{d})$
$\mathbf{t}$ |  $L$ discrete timepoints |  $[0,T]$
$\boldsymbol{\mathbb{X}}_{\mathbf{t}}$ |  Collection of $M$ independent samples of $\mathbf{X}_{t}$ at $\mathbf{t}$ |  $\mathbb{R}^{MLNd}$
$\boldsymbol{\mathbb{Y}}_{\mathbf{t}}$ |  Sample of $\mathbf{X}_{\mathbf{t}}$ corrupted with i.i.d. additive noise |  $\mathbb{R}^{MLNd}$
$U_{t}$ |  approximate density from particle positions |  ${\mathcal{P}}(\mathbb{R}^{d})$
$G$ |  density kernel mapping $\mu^{N}_{t}$ to $U_{t}$ |  $L^{1}(\mathbb{R}^{d}\times\mathbb{R}^{d},\mathbb{R})$
${\mathcal{D}}$ |  spatial support of $U_{t}$, $t\in[0,T]$ |  compact subset of $\mathbb{R}^{d}$
$\mathbf{C}$ |  discretization of ${\mathcal{D}}$ |
$\mathbf{U}_{t}$ |  discrete approximate density $U_{t}(\mathbf{C})$ |
$\left\langle\cdot,\cdot\right\rangle_{h}$ |  semi-discrete inner product, trapezoidal rule over $\mathbf{C}$ |
$\left\langle\cdot,\cdot\right\rangle_{h,\Delta t}$ |  fully-discrete inner product, trapezoidal rule over $\mathbf{C}\times\mathbf{t}$ |
$\mathbb{L}_{K}$ |  library of candidate interaction forces |
$\mathbb{L}_{V}$ |  library of candidate local forces |
$\mathbb{L}_{\sigma}$ |  library of candidate diffusivities |
$\mathbb{L}$ |  $(\mathbb{L}_{K},\mathbb{L}_{V},\mathbb{L}_{\sigma})$ |
$\Psi$ |  set of $n$ test functions $(\psi_{k})_{k=1}^{n}$ |
$\phi_{m,p}(v;\,\Delta)$ |  test functions used in this work (equation (4.3)) |
$\boldsymbol{\lambda}$ |  set of sparsity thresholds |
${\mathcal{L}}$ |  loss function for sparsity thresholds (equation (4.6)) |
Table 1. Notations used throughout. Mean-field Term | Trial Function Library
---|---
$\nabla\cdot(U\nabla K*U)$ |  $\nabla\cdot(U\nabla|x|^{m}*U)$, $m\in\\{1,2,3,4,5,6,7\\}$
$\nabla\cdot(U\nabla V)$ |  $\partial_{x_{i}}\left(U\cos(mx_{1})\cos(nx_{2})\right)$, $(m,n)\in\\{0,1,2,3,4,5\\}$, $i\in\\{1,2\\}$
$\frac{1}{2}\sum_{i=1}^{d}\frac{\partial^{2}(U\sigma\sigma^{T})_{ij}}{\partial x_{i}\partial x_{j}}$ |  $\Delta(U\cos(mx_{1})\cos(nx_{2}))$, $(m,n)\in\\{0,1,2,3,4,5\\}$
Table 2. Trial function library for local 2D example (Section 5.1). Mean-field Term | Trial Function Library
---|---
$\nabla\cdot(U\nabla K*U)$ |  $\partial_{x}\cdot(U\partial_{x}|x|^{m}*U)$, $m\in\\{1,2,3,4,5,6,7\\}$
$\nabla\cdot(U\nabla V)$ |  $\partial_{x}\left(Ux^{m}\right)$, $m\in\\{0,2,3,4,5,6,7,8\\}$
$\frac{1}{2}\sum_{i=1}^{d}\frac{\partial^{2}(U\sigma\sigma^{T})_{ij}}{\partial x_{i}\partial x_{j}}$ |  $\partial_{xx}(Ux^{m})$, $m\in\\{0,1,2,3,4,5,6,7,8\\}$
Table 3. Trial function library for nonlocal 1D example (Section 5.2). Mean-field Term | Trial Function Library
---|---
$\nabla\cdot(U\nabla K*U)$ | $\begin{dcases}\nabla\cdot(U\nabla|x|^{m}*U),&m\in\\{2,3,4,5,6\\}\\\ \nabla\cdot(U\nabla\left[|x|^{1/2}\right]_{\delta}*U)\\\ \nabla\cdot(U\nabla\left[|x|(\log|x|-1)\right]_{\delta}*U)\\\ \nabla\cdot(U\nabla\left[\log|x|\right]_{\delta}*U)\\\ \end{dcases}$
$\nabla\cdot(U\nabla V)$ |  $\nabla\left(U\nabla(x_{1}^{m}x_{2}^{n})\right)$, $(m,n)\in\mathbb{N}\times\mathbb{N},1\leq m+n\leq 6$
$\frac{1}{2}\sum_{i=1}^{d}\frac{\partial^{2}(U\sigma\sigma^{T})_{ij}}{\partial x_{i}\partial x_{j}}$ |  $\frac{\partial^{2}}{\partial x_{i}\partial x_{j}}(U\cos(mx_{1})\cos(nx_{2}))$, $(m,n)\in\\{0,1,2\\}$, $(i,j)\in\\{1,2\\}$
Table 4. Trial function library for nonlocal 2D example (Section 5.3). Interaction potentials $[\,K\,]_{\delta}$ indicate cutoff potentials of the form (5.6) with $\delta=0.01$ such that the resulting potential is Lipschitz. Example | $m_{x}$ | $m_{t}$ | $p_{x}$ | $p_{t}$ | $s_{x}$ | $s_{t}$ | size$(\mathbf{G})$ | $\kappa(\mathbf{G})$ | size$(\mathbf{U})$ | $(h,\Delta t)$ | Walltime
---|---|---|---|---|---|---|---|---|---|---|---
Local 2D | 31 | 16 | 5 | 3 | 10 | 5 | 686$\times 85$ | $3.8\times 10^{7}$ | $128\times 128\times 101$ | $(0.078,0.02)$ | 9.7s
Nonlocal 1D | 29 | 8 | 5 | 3 | 5 | 1 | 3368$\times 24$ | $8.7\times 10^{8}$ | $256\times 101$ | $(0.023,0.01)$ | 2.6s
Nonlocal 2D | 25 | 8 | 5 | 3 | 8 | 1 | 6500$\times 53$ | $4.8\times 10^{7}$ | $128\times 128\times 81$ | $(0.047,0.1)$ | 8.5s
Table 5. Discretization parameters for examples. Note the condition number
$\kappa(\mathbf{G})$ and walltime are listed for representative samples with
64,000 total particles.