# Learning Collective Behaviors from Observation

Jinchao Feng111Great Bay University, Dongguan, Guangdong, China, email:
<EMAIL_ADDRESS>and Ming Zhong222Illinois Institute of Technology, Chicago,
IL, USA, email<EMAIL_ADDRESS>
###### Abstract

We present a review of a series of learning methods used to identify the
structure of dynamical systems, aiming to understand emergent behaviors in
complex systems of interacting agents. These methods not only offer
theoretical guarantees of convergence but also demonstrate computational
efficiency in handling high-dimensional observational data. They can manage
observation data from both first- and second-order dynamical systems,
accounting for observation/stochastic noise, complex interaction rules,
missing interaction features, and real-world observations of interacting agent
systems. The essence of developing such a series of learning methods lies in
designing appropriate loss functions using the variational inverse problem
approach, which inherently provides dimension reduction capabilities to our
learning methods.

##  1 Introduction

Accurate discovery of governing differential equations for biological/physical
systems from observations can be traced back to at least Gauss, Lagrange, and
Laplace [1]. Kepler’s laws of basic planetary motion, Newton’s universal law
of gravitation, and Einstein’s general relativity provide excellent examples
for the discovery of physical laws to understand motion from data. The capture
and subsequent description of such physical/biological laws have always been
important research topics. We review a series of learning methods to provide
mathematical insight into observational data; moreover, these methods can also
offer candidate theories to explain the data.

We are interested in cases where the data is assumed to follow certain
dynamical systems (possibly with stochastic noise/forcing). To simplify our
discussion, we use the following autonomous dynamical system:

$\dot{\mathbf{Y}}=\mathbf{f}(\mathbf{Y}),\quad\mathbf{Y}\in\mathbb{R}^{D},$

where the time-dependent variable $\mathbf{Y}$ describes a certain state in
the system (such as position, velocity, temperature, phase, opinion, or a
combination of these), the function
$\mathbf{f}:\mathbb{R}^{D}\rightarrow\mathbb{R}^{D}$ provides the change
needed to update $\mathbf{Y}$. Given the observation of
$\\{\mathbf{Y}(t)\\}_{t\in[0,T]}$, is it possible to identify $\mathbf{f}$ to
explain the data within the given time $[0,T]$? Many system identification
methods have been developed, such as Sparse Identification of Nonlinear
Dynamics (SINDy, [2]), Physics Informed Neural Network (PINN, [3]), Neural ODE
[4], and many more. However, in the case of high-dimensional data, i.e. $D\gg
1$, the identification task becomes computationally prohibitive.

To overcome the curse of dimensionality, we employ a special and effective
dimension reduction technique by constructing our learning methods based on
the unique structure of rhs. This unique structure arises from the fact that
dynamical systems are used to model inspiring natural phenomena, namely
collective behaviors, also known as self-organization. These behaviors can
manifest as clustering, flocking, milling, swarming, and synchronization,
demonstrating certain global patterns from initially random configurations
through local interactions. They have applications in Biology, Chemistry,
Computer Science, Economics, Physics, Social Science, and more (see references
in [5, 6, 7]). These behaviors can be described using dynamical systems. In
fact, first-order systems can be modeled as follows:

$\dot{\mathbf{y}}_{i}=\frac{1}{N}\sum_{i^{\prime}=1,i^{\prime}\neq
i}^{N}\phi(|\mathbf{y}_{i^{\prime}}-\mathbf{y}_{i}|)(\mathbf{y}_{i^{\prime}}-\mathbf{y}_{i}),\quad
i=1,\cdots,N.$

Here, $\mathbf{y}_{i}\in\mathbb{R}^{d}$ describes the state of the $i^{th}$
agent, and $\phi:\mathbb{R}^{+}\rightarrow\mathbb{R}$ gives the interaction
law on how agent $i^{\prime}$ influences the change of state for agent $i$.
This system is a gradient flow of certain system energy with
rotation/permutation invariances and symmetry. Such a first-order system can
be used to model opinion dynamics (formation of consensus), crystal structure,
and other pattern formations in skin pigmentation [8]. By using vector
notation, i.e.

$\mathbf{Y}=\begin{bmatrix}\vdots\\\ \mathbf{y}_{i}\\\
\vdots\end{bmatrix}\in\mathbb{R}^{D=Nd}\quad\text{and}\quad\mathbf{f}_{\phi}(\mathbf{Y})=\begin{bmatrix}\vdots\\\
\frac{1}{N}\sum_{i^{\prime}=1,i^{\prime}\neq
i}^{N}\phi(|\mathbf{y}_{i^{\prime}}-\mathbf{y}_{i}|)(\mathbf{y}_{i^{\prime}}-\mathbf{y}_{i})\\\
\vdots\end{bmatrix}\in\mathbb{R}^{D},$

we end up with a usual autonomous dynamical system in the form of
$\dot{\mathbf{Y}}=\mathbf{f}_{\phi}(\mathbf{Y})$ with
$\mathbf{Y}\in\mathbb{R}^{D}$ and $D=Nd\gg 1$. For example, even for a simple
system with $5$ agents where $\mathbf{x}_{i}\in\mathbb{R}^{2}$, the system
state variable lives in $D=10$ dimensions. How to efficiently deal with such
high dimensionality has become a crucial aspect of developing our learning
algorithms [5, 9, 10, 11, 12, 13, 14, 15]. Our methods are not only proven to
have theoretical guarantees but are also verified through rigorous numerical
testing. Even when learning from real observational data, our method can offer
valuable insights into how the framework of collective behaviors can be used
to explain complex systems.

The remaining sections of the paper are organized as follows: in section 2, we
discuss the general model form used for our learning paradigm. In section 3,
we go through various learning scenarios, from first-order to second-order,
from homogeneous agents to heterogeneous agents, and with Gaussian Process
priors. We compare our methods to other well-established methods in section
3.1.1. Finally, in section 4, we conclude our review with a highlight of
several future directions.

##  2 Model Equation

In order to provide a more complete picture of collective dynamics, we
consider the following state variable $\mathbf{y}_{i}$ for the $i^{th}$ agent
in a system of $N$ interacting agents. Moreover, $\mathbf{y}_{i}$ is expressed
in the form of $\mathbf{y}_{i}=\begin{bmatrix}\mathbf{x}_{i}\\\
\mathbf{v}_{i}\\\ \zeta_{i}\end{bmatrix}$, where
$\mathbf{x}_{i}\in\mathbb{R}^{d}$ represents the position of agent $i$,
$\mathbf{v}_{i}\in\mathbb{R}^{d}$ is its corresponding velocity (hence
$\mathbf{v}_{i}=\dot{\mathbf{x}}_{i}$), and $\zeta_{i}\in\mathbb{R}$ is an
additional state that can be used to describe opinions, phases, excitation
levels (e.g. towards light sources), emotions, etc. The coupled state variable
$\mathbf{y}_{i}\in\mathbb{R}^{2d+1}$ satisfies the following coupled dynamical
system (for $i=1,\cdots,N$):

$\begin{cases}\dot{\mathbf{x}}_{i}&=\mathbf{v}_{i}\\\
m_{i}\dot{\mathbf{v}}_{i}&=F^{\mathbf{v}}(\mathbf{y}_{i})+\sum_{i^{\prime}=1,i^{\prime}\neq
i}^{N}\frac{1}{N_{\mathit{k}_{i}}}\Big{[}\Phi^{E}_{\mathit{k}_{i},\mathit{k}_{i^{\prime}}}(\mathbf{y}_{i},\mathbf{y}_{i^{\prime}})\mathbf{w}^{E}(\mathbf{x}_{i},\mathbf{x}_{i^{\prime}})\\\
&\quad+\Phi^{A}_{\mathit{k}_{i},\mathit{k}_{i^{\prime}}}(\mathbf{y}_{i},\mathbf{y}_{i^{\prime}})\mathbf{w}^{A}(\mathbf{v}_{i},\mathbf{v}_{i^{\prime}})\Big{]}\\\
\dot{\zeta}_{i}&=F^{\zeta}(\mathbf{y}_{i})+\sum_{i^{\prime}=1,i^{\prime}\neq
i}^{N}\frac{1}{N_{\mathit{k}_{i}}}\Phi^{\zeta}_{\mathit{k}_{i},\mathit{k}_{i^{\prime}}}(\mathbf{y}_{i},\mathbf{y}_{i^{\prime}})(\zeta_{i^{\prime}}-\zeta_{i})\end{cases}.$
(1)

Here, we consider the system to be partitioned into $K$ types, i.e.
$C_{k_{1}}\cap C_{k_{2}}=\emptyset$ for $1\leq k_{1}\neq k_{2}=K$ and
$\cup_{k=1}^{K}C_{K}=[N]=\\{1,\cdots,N\\}$, and $N_{k}$ is the agent of agents
in type $k$. Moreover, the types of the agents do not change over time333Such
a restriction can be relaxed, as long as we have the type information for the
agents at all times., and there is a type function
$\mathit{k}_{\cdot}:\mathbb{N}\rightarrow\mathbb{N}$, which returns the type
index for an agent index. Furthermore, the external force
$F^{\mathbf{v}}:\mathbb{R}^{2d+1}\rightarrow\mathbb{R}^{d}$ gives the
environmental effect on the velocity of the agent $i$, similarly
$F^{\zeta}:\mathbb{R}^{2d+1}\rightarrow\mathbb{R}$ gives the environmental
effect on $\xi_{i}$, $\Phi^{E}:\mathbb{R}^{4d+2}\rightarrow\mathbb{R}$ gives
the energy-based interaction,
$\Phi^{A}:\mathbb{R}^{4d+2}\rightarrow\mathbb{R}$ gives the alignment-based
interaction, $\Phi^{\zeta}:\mathbb{R}^{4d+2}\rightarrow\mathbb{R}$ gives the
$\zeta$-based interaction,
$\mathbf{w}^{E},\mathbf{w}^{A}:\mathbb{R}^{2d}\rightarrow\mathbb{R}^{d}$ give
the directions for $\Phi^{E},\Phi^{A}$ respectively.

###### Remark 1

The inclusion of $(\mathbf{x}_{i},\mathbf{v}_{i})$ into the system also makes
it a second-order system (second-order time derivative for $\mathbf{x}_{i}$).
We include the additional state $\zeta_{i}$ to provide a more realistic
modeling of the emergent behaviors. This formulation can also include the
first-order systems as special cases. When
$F^{\mathbf{v}}(\mathbf{y}_{i})=-\nu_{i}\mathbf{v}_{i}$ (the usual friction)
with $\nu_{i}\gg m_{i}$, $\Phi^{A}=0$ and $\Phi^{E}$ depends on
$(\mathbf{x}_{i},\zeta_{i},\mathbf{x}_{i^{\prime}},\zeta_{i^{\prime}})$, the
second-order systems become:

$\mathbf{0}\approx-\mathbf{v}_{i}+\frac{1}{\nu_{i}}\sum_{i^{\prime}=1,i^{\prime}\neq
i}^{N}\frac{1}{N_{\mathit{k}_{i}}}\Phi^{E}_{\mathit{k}_{i},\mathit{k}_{i^{\prime}}}(\mathbf{x}_{i},\zeta_{i},\mathbf{x}_{i^{\prime}},\zeta_{i^{\prime}})\mathbf{w}^{\mathbf{x}}(\mathbf{x}_{i},\mathbf{x}_{i^{\prime}}),$

we end up with a first-order system of $\mathbf{x}_{i}$ coupled with
$\zeta_{i}$ (see the swarmalator model [16]).

Our second-order coupled formulation (1) can include a rather extensive list
of behaviors, such as flocking with external potential [17], anticipated
flocking dynamics [18], fish swarming dynamics [19], swarmalator dynamics
(concurrent swarming and synchronization) [20, 21, 16], line alignment
dynamics [22], where the interaction kernel depends on pairwise state
variables. It is possible to include more state variables, i.e., we can also
consider $(\mathbf{x}_{i},\mathbf{v}_{i},\zeta_{i},\eta_{i})$, where
$\eta_{i}$ describes a different state than $\zeta_{i}$. However, the learning
of such a system only technically involves more terms than what we have
already proposed, and we will leave the details for future projects.

Using the vector notation, we can obtain a simplified system denoted as
$\dot{\mathbf{Y}}=\mathbf{f}_{\Phi^{E},\Phi^{A},\Phi^{\zeta}}(\mathbf{Y})$,
where
$\mathbf{Y}=\begin{bmatrix}\cdots&\mathbf{y}_{i}^{\top}&\cdots\end{bmatrix}^{\top}\in\mathbb{R}^{D=N(2d+1)}$
and

$\mathbf{f}(\mathbf{Y})=\begin{bmatrix}\vdots\\\ \mathbf{v}_{i}\\\
F^{\mathbf{v}}(\mathbf{y}_{i})+\sum_{i^{\prime}=1,i^{\prime}\neq
i}^{N}\frac{1}{N_{\mathit{k}_{i}}}\Big{[}\Phi^{E}_{\mathit{k}_{i},\mathit{k}_{i^{\prime}}}(\mathbf{y}_{i},\mathbf{y}_{i^{\prime}})\mathbf{w}^{E}(\mathbf{x}_{i},\mathbf{x}_{i^{\prime}})\\\
\quad+\Phi^{A}_{\mathit{k}_{i},\mathit{k}_{i^{\prime}}}(\mathbf{y}_{i},\mathbf{y}_{i^{\prime}})(\mathbf{v}_{i^{\prime}}-\mathbf{v}_{i})\Big{]}\\\
F^{\zeta}(\mathbf{y}_{i})+\sum_{i^{\prime}=1,i^{\prime}\neq
i}^{N}\frac{1}{N_{\mathit{k}_{i}}}\Phi^{\zeta}_{\mathit{k}_{i},\mathit{k}_{i^{\prime}}}(\mathbf{y}_{i},\mathbf{y}_{i^{\prime}})(\zeta_{i^{\prime}}-\zeta_{i})\\\
\vdots\end{bmatrix}\in\mathbb{R}^{N(2d+1)}.$

For such a coupled system, we seek to identify
$\mathbf{f}=\mathbf{f}_{\\{\Phi^{E}_{k_{1},k_{2}},\Phi^{A}_{k_{1},k_{2}},\Phi^{\zeta}_{k_{1},k_{2}}\\}_{k_{1},k_{2}}^{K}}$
in terms of the interaction kernels
$\\{\Phi^{E}_{k_{1},k_{2}},\Phi^{A}_{k_{1},k_{2}},\Phi^{\zeta}_{k_{1},k_{2}}\\}_{k_{1},k_{2}}^{K}$
from observations. In a simple system of $5$ agents with $d=2$, we are facing
a learning problem of inferring $25$-dimensional functions of $25$-variables,
which make the problem computationally prohibitive for traditional methods to
infer. We will discuss the details of how to handle such high-dimensional
learning in section 3.

##  3 Learning Methods

We are now prepared to discuss in detail the learning framework developed for
the efficient learning of such high-dimensional dynamical systems. We will
begin our discussion with the simplest case: first-order systems, and then
gradually extend the method to more complex scenarios, such as systems with
multiple types of agents, geometry-constrained dynamics, missing feature maps,
coupled systems (for high-order systems), and learning with a Gaussian Process
prior.

###  3.1 First Order

We use a simple first-order system as the starting point to build our learning
framework. The first-order system considered here is given as follows

$\dot{\mathbf{x}}_{i}=\frac{1}{N}\sum_{i^{\prime}=1,i^{\prime}\neq
i}^{N}\Phi(\mathbf{x}_{i},\mathbf{x}_{i^{\prime}})(\mathbf{x}_{i^{\prime}}-\mathbf{x}_{i}),\quad
i=1,\cdots,N.$ (2)

Here $\mathbf{x}_{i}$ is the state variable of the $i^{th}$ agent, and
$\Phi=\Phi^{E}:\mathbb{R}^{2d}\rightarrow\mathbb{R}$ gives the energy-based
(short-range repulsion and long-range attraction, we omit the superscript $E$
here to streamline the notations) interaction. We only require the interaction
kernel $\Phi$ to be symmetric, i.e.
$\Phi(\mathbf{x}_{i},\mathbf{x}_{i^{\prime}})=\Phi(\mathbf{x}_{i^{\prime}},\mathbf{x}_{i})$.
We shall further assume that $\Phi$ can be written in the following
composition form

$\Phi(\mathbf{x}_{i},\mathbf{x}_{i^{\prime}})=\phi(\boldsymbol{\xi}(\mathbf{x}_{i},\mathbf{x}_{i})),\quad\phi:\mathbb{R}^{d_{\phi}}\rightarrow\mathbb{R}\quad\boldsymbol{\xi}:\mathbb{R}^{2d}\rightarrow\mathbb{R}^{d_{\phi}}.$
(3)

Here $\phi$ is the reduced interaction kernel, and $\boldsymbol{\xi}$ is the
reduced interaction variable. In most of the collective behavior’s models,
$\boldsymbol{\xi}$ is assumed to be the pairwise distance variable, i.e.

$\boldsymbol{\xi}(\mathbf{x}_{i},\mathbf{x}_{i^{\prime}})=|\mathbf{x}_{i^{\prime}}-\mathbf{x}_{i}|,\quad|\mathbf{x}|=\sqrt{\sum_{j=1}^{d}|(\mathbf{x})_{j}|^{2}},\quad\forall\mathbf{x}\in\mathbb{R}^{d}.$

(a) $\mathbf{X}(t)$ vs $\hat{\mathbf{X}}(t)$

(b) $\phi$ vs $\hat{\phi}$

Figure 1: Opinion Dynamics with $N=20$ agents and
$\phi(r)=\chi_{[0,\frac{1}{\sqrt{2}})}(r)+0.1*\chi_{[\frac{1}{\sqrt{2}},1]}(r)$,
learned on $[0,10]$, for other parameters see [5].

Hence we arrive at a simpler and more familiar form of (2)

$\dot{\mathbf{x}}_{i}=\frac{1}{N}\sum_{i^{\prime}=1,i^{\prime}\neq
i}^{N}\phi(|\mathbf{x}_{i^{\prime}}-\mathbf{x}_{i}|)(\mathbf{x}_{i^{\prime}}-\mathbf{x}_{i}),\quad
i=1,\cdots,N.$ (4)

###### Remark 2

(4) can be viewed as a gradient flow of the following form

$\dot{\mathbf{x}}_{i}=\nabla_{\mathbf{x}_{i}}\Big{(}\frac{1}{2N}\sum_{1\leq
i\neq i\leq N}U(|\mathbf{x}_{i^{\prime}}-\mathbf{x}_{i}|)\Big{)},\quad
i=1,\cdots,N.$

Here $U:\mathbb{R}^{+}\rightarrow\mathbb{R}$ is an energy potential that
depends on pairwise distance.

We are now ready to discuss the learning framework based on the $1$
dimensional $\phi$, however, our learning method can be easily extended to the
case for $\Phi$. Given observation data in the form,
$\\{\mathbf{x}^{m}_{i}(t_{l}),\dot{\mathbf{x}}^{m}_{i}(t_{l})\\}_{i,l,m=1}^{N,L,M}$
for $0=t_{1}<t_{2}<\cdots<t_{L}=T$ and $\mathbf{x}_{i}^{m}(0)$s are i.i.d
sample of certain probability distribution $\mu_{0}$, we find the estimator to
$\phi$ from minimizing the following loss functional,

$\mathcal{E}_{L,M,\mathcal{H}}(\varphi)=\frac{1}{LM}\sum_{l,m=1}^{L,M}|\dot{\mathbf{X}}_{l}^{m}-\mathbf{f}_{\varphi}(\mathbf{X}_{l}^{m})|^{2}_{\mathcal{S}},\quad\varphi\in\mathcal{H}.$
(5)

Here
$\mathbf{X}_{l}^{m}=\begin{bmatrix}\cdots&\mathbf{x}_{i}^{\top}&\cdots\end{bmatrix}^{\top}$
and

$\mathbf{f}_{\varphi}(\mathbf{X}_{l}^{m})=\begin{bmatrix}\vdots\\\
\frac{1}{N}\sum_{i^{\prime}=1,i^{\prime}\neq
i}^{N}\phi(|\mathbf{x}_{i^{\prime}}^{m}(t_{l})-\mathbf{x}_{i}|^{m}(t_{l}))(\mathbf{x}_{i^{\prime}}^{m}(t_{l})-\mathbf{x}_{i}^{m}(t_{l}))\\\
\vdots\end{bmatrix}.$

and the norm $|\cdot|_{\mathcal{S}}$ is defined as follows

$|\mathbf{U}-\mathbf{V}|_{\mathcal{S}}^{2}=\frac{1}{N}\sum_{i=1}^{N}|\mathbf{u}_{i}-\mathbf{v}_{i}|^{2}_{\ell^{2}(\mathbb{R}^{d})},\quad\mathbf{U}=\begin{bmatrix}\vdots\\\
\mathbf{u}_{i}\\\ \vdots\end{bmatrix},\quad\mathbf{V}=\begin{bmatrix}\vdots\\\
\mathbf{v}_{i}\\\
\vdots\end{bmatrix},\quad\mathbf{u}_{i},\mathbf{v}_{i}\in\mathbb{R}^{d}.$

Apparently
$\mathbf{X}_{l}^{m},\mathbf{f}_{\varphi}(\mathbf{X}_{l}^{m})\in\mathbb{R}^{D=Nd}$.
We also use a compact and convex hypothesis space $\mathcal{H}$ so that the
existence and uniqueness of the minimization is guaranteed. We denote the the
unique minimizer of $\mathcal{E}_{L,M,\mathcal{H}}$ as
$\hat{\phi}_{L,M,\mathcal{H}}$, i.e.
$\hat{\phi}_{L,M,\mathcal{H}}=\underset{\varphi\in\mathcal{H}}{\operatorname{arg}\operatorname{min}}\;\mathcal{E}_{L,M,\mathcal{H}}(\varphi)$.
In order to show the convergence of $\hat{\phi}_{L,M,\mathcal{H}}$ to $\phi$
when $M\rightarrow\infty$, we will need to consider the expected version of
the loss functional, i.e.

$\mathcal{E}_{L,\mathcal{H}}(\varphi)=\frac{1}{L}\mathbb{E}\Big{[}\sum_{l}^{L}|\dot{\mathbf{X}}_{l}^{m}-\mathbf{f}_{\varphi}(\mathbf{X}_{l}^{m})|_{\mathcal{S}}^{2}\Big{]}.$

We denote the minimizer $\hat{\phi}_{L,\mathcal{H}}$ from minimzing
$\mathcal{E}_{L,\mathcal{H}}$ over the same hypothesis space $\mathcal{H}$,
i.e.
$\hat{\phi}_{L,\mathcal{H}}=\underset{\varphi\in\mathcal{H}}{\operatorname{arg}\operatorname{min}}\;\mathcal{E}_{L,\mathcal{H}}(\varphi)$.
Using law of the large number,
$\hat{\phi}_{L,M,\mathcal{H}}\rightarrow\hat{\phi}_{L,\mathcal{H}}$ as
$M\rightarrow\infty$. The only thing left is to show the distance between
$\hat{\phi}_{L,\mathcal{H}}$ and $\phi$ due to the approximation power of
$\mathcal{H}$. Before we show the lemma needed to prove convergence, we need
to define a weighted $L^{2}(\rho)$ norm where the weight function $\rho$ is
given as follows

$\rho(r)=\frac{1}{{N\choose
2}T}\int_{t=0}^{T}\mathbb{E}_{\mathbf{X}_{0}\sim\mu_{0}}\Big{[}\sum_{i,i^{\prime}=1,i<i^{\prime}}^{N}\delta_{r_{i,i^{\prime}}(t)}(r)\,dt\Big{]}.$
(6)

Here $\mathbf{r}_{i,i^{\prime}}=\mathbf{x}_{i^{\prime}}-\mathbf{x}_{i}$ and
$r_{i,i^{\prime}}=|\mathbf{r}_{i,i^{\prime}}|$. $\rho$ is used to illustrate
the distribution of pairwise distance data through the interaction of agents
within the system; such data is also used to learn $\phi$, however, we do not
have direct access to the data points of $\phi(r)$, rather we have the
indirect observation of $\phi$. With this definition of $\rho$, we can define
our weighted $L^{2}(\rho)$ norm as follows

$|\phi_{1}(\cdot)\cdot-\phi_{2}(\cdot)\cdot|_{L^{2}(\rho)}^{2}=\int_{r\in\text{supp}(\rho)}|\phi_{1}(r)-\phi_{2}(r)|r\,d\rho(r),$
(7)

for any $\phi_{1},\phi_{2}\in L^{2}(\text{supp}(\rho),\rho)$. Now we are ready
to show the lemma needed to prove the convergence theorem.

###### Definition 1

The dynamical system given by (4) with IC sampled from $\mu_{0}$ on
$\mathbb{R}^{D=Nd}$, satisfies the coercivity condition on a set $\mathcal{H}$
if there exists a constant $C_{L,N,\mathcal{H}}$ such that for all
$\varphi\in\mathcal{H}$ with $\varphi(\cdot)\cdot\in L^{2}(\rho)$,

$C_{L,N,\mathcal{H}}|\varphi(\cdot)\cdot|^{2}_{L^{2}(\rho)}\leq\frac{1}{LN}\sum_{i,l=1}^{N,L}\mathbb{E}|\frac{1}{N}\sum_{i^{\prime}=1,i^{\prime}\neq
i}^{N}\varphi(r_{i,i^{\prime}}(t_{l})\mathbf{r}_{i,i^{\prime}}(t_{l})|^{2}.$

The coercivity condition given by definition 1 basically states that
minimizing $\mathcal{E}_{L,\mathcal{H}}$ would also minimize the distance
between $\hat{\phi}_{L,\mathcal{H}}$ and $\phi$ under weighted $L^{2}(\rho)$
norm, i.e.

$C_{L,N,\mathcal{H}}|\phi(\cdot)\cdot-\varphi(\cdot)\cdot|^{2}_{L^{2}(\rho)}\leq\mathcal{E}_{L,\mathcal{H}}(\varphi)$

To show this, simply use the fact that $\phi$ is the interaction kernel which
gives the observation data, i.e.
$\dot{\mathbf{x}}_{i}=\frac{1}{N}\sum_{i^{\prime}=1,i^{\prime}\neq
i}^{N}\phi(r_{i,i^{\prime}})\mathbf{r}_{i,i^{\prime}}$.

###### Theorem 1

Assume that $\phi\in\mathcal{K}_{R,S}$ (an admissible set
$\mathcal{K}_{R,S}=\\{\phi\in
C^{1}(\mathbb{R}_{+}):\text{supp}(\phi)\subset[0,R],\sup_{r\in[0,R]}|\phi(r)|+|\phi^{\prime}(r)|\leq
S\\}$ for some $R,S>0$). Let $\\{\mathcal{H}_{n}\\}_{n}$ be a sequence of
subspaces of $L^{\infty}([0,R])$, with $\text{dim}(\mathcal{H}_{n})\leq
c_{0}n$ and
$\inf_{\varphi\in\mathcal{H}_{n}}|\varphi-\phi|_{L^{\infty}([0,R])}\leq
c_{1}n^{-s}$, for some constants $c_{0},c_{1},s>0$. Assume that the coercivity
condition holds on
$\mathcal{H}\coloneqq\overline{\cup_{n=1}^{\infty}\mathcal{H}_{n}}$. Such a
sequence exists, for example, if $\phi$ is $s$-Hölder regular, and can be
chosen so that $\mathcal{H}$ is compact in $L^{2}(\rho)$. Choose
$n_{*}=(\frac{M}{\log(M)})^{\frac{1}{2s+1}}$. Then, there exists a constant
$C=C(c_{0},c_{1},R,S)$ such that

$\mathbb{E}\Big{[}|\hat{\phi}_{L,M,\mathcal{H}_{n_{*}}}(\cdot)\cdot-\phi(\cdot)\cdot|_{L^{2}(\rho)}\Big{]}\leq\frac{C}{C_{L,N,\mathcal{H}}}(\frac{M}{\log(M)})^{\frac{s}{2s+1}}.$

Our learning rate is as optimal as if it was learned from the regression
setting, however different from regression, we do not have direct access of
the data points $\phi(r_{i})$. For more detailed discussion, see [5].

When we choose a basis for $\mathcal{H}_{n}$, i.e.
$\mathcal{H}_{n}=\text{span}\\{\psi_{1},\cdots,\psi_{n}\\}$ and let
$\varphi=\sum_{\eta=1}^{n}\alpha_{\eta}\psi_{\eta}$, then the aforementioned
minimization problem can be re-written as a linear system as
$A\vec{\alpha}=\vec{b}$, where
$\vec{\alpha}=\begin{bmatrix}\alpha_{1}&\cdots&\alpha_{n}\end{bmatrix}^{\top}$.
Moreover we have theoretical guarantee for the well-conditioning of the
system, hence it is computational tractable to solve $A\vec{\alpha}=\vec{b}$.
Further more, for big data set, one can use parallelization to assemble $A$
and $\vec{b}$, which makes the overall process (near-)linear. See [5] for the
detailed description of the algorithm and
https://github.com/mingjzhong/LearningDynamics for the software pakcage.

####  3.1.1 Other System Identification Methods

There are many different methods of learning dynamical systems, i.e. using the
parametric structure of the right-hand side with Bayesian inference, using the
sparse structure of the right-hand side (SINDy [2]), using neural networks to
estimate the right-hand side (PINN [3] and NeuralODE [4]). We will discuss
SINDy and the Neural Network approach and make a direct comparison of these
two methods to ours.

SINDy stands for Sparse Identification of Nonlinear Dynamics is a data-driven
identification method for dynamical systems. It assumes that the dynamical
system $\dot{\mathbf{Y}}=\mathbf{f}(\mathbf{Y})$ can be approximated using the
following form

$\mathbf{f}(\mathbf{Y})\approx
a_{1}\psi_{1}(\mathbf{Y})+\cdots+a_{n}\psi_{n}(\mathbf{Y}),\quad\mathbf{Y}=\begin{bmatrix}(\mathbf{Y})_{1}\\\
\vdots\\\ (\mathbf{Y})_{D}\end{bmatrix}\in\mathbb{R}^{D}.$

$\psi_{i}$ is applied component wise to $\mathbf{Y}$, i.e.

$\psi_{i}(\mathbf{Y})=\begin{bmatrix}\psi_{i}(y_{1})\\\ \vdots\\\
\psi_{i}(y_{D})\end{bmatrix},\quad\mathbf{Y}=\begin{bmatrix}y_{1}\\\ \vdots\\\
y_{D}\end{bmatrix}.$

Furthermore, $\\{\psi_{1},\cdots,\psi_{n}\\}$ is a set of basis functions from
a predetermined dictionary, e.g. they can be a set of polynomials,
sine/cosine, and/or negative power polynomials. Given the observation data
$\\{\mathbf{Y}(t_{1}),\cdots,\mathbf{Y}(t_{L})\\}$, we assemble the matrix

$\mathbf{Y}=\begin{bmatrix}\mathbf{Y}^{\top}(t_{1})\\\ \vdots\\\
\mathbf{Y}^{\top}(t_{L})\end{bmatrix}=\begin{bmatrix}y_{1}(t_{1})&y_{2}(t_{1})&\cdots&y_{D}(t_{1})\\\
y_{1}(t_{2})&y_{2}(t_{2})&\cdots&y_{D}(t_{2})\\\
\vdots&\vdots&\ddots&\vdots\\\
y_{1}(t_{L})&y_{2}(t_{L})&\cdots&y_{D}(t_{L})\end{bmatrix}\in\mathbb{R}^{L\times
D}.$

We construct a dictionary (or library) $\Theta(\mathbf{Y})$ of nonlinear
candidate functions of $\mathbf{Y}$, i.e.

$\Theta(\mathbf{Y})=\begin{bmatrix}|&|&&|\\\
\psi_{1}(\mathbf{Y})&\psi_{2}(\mathbf{Y})&\cdots&\psi_{n}(\mathbf{Y})\\\
|&|&&|\end{bmatrix}\in\mathbb{R}^{L\times nD}.$

Again $\psi_{i}(\mathbf{Y})$ is applied component wise, i.e.

$\psi_{i}(\mathbf{Y})=\begin{bmatrix}\psi_{i}(y_{1}(t_{1}))&\psi_{i}(y_{2}(t_{1}))&\cdots&\psi_{i}(y_{D}(t_{1}))\\\
\psi_{i}(y_{1}(t_{2}))&\psi_{i}(y_{2}(t_{2}))&\cdots&\psi_{i}(y_{D}(t_{2}))\\\
\vdots&\vdots&\ddots&\vdots\\\
\psi_{i}(y_{1}(t_{L}))&\psi_{i}(y_{2}(t_{L}))&\cdots&\psi_{i}(y_{D}(t_{L}))\end{bmatrix}.$

The set of basis functions can be
$\\{1,y,y^{2},y^{3},\cdots,\sin(y),\cos(y),\cdots\\}$, depending on prior
knowledge or computational capacity. We write the dynamical system in the new
form

$\dot{\mathbf{Y}}=\Theta(\mathbf{Y})\Xi,\quad\Xi=\begin{bmatrix}\boldsymbol{\xi}_{1}&\cdots&\boldsymbol{\xi}_{D}\end{bmatrix},\quad\boldsymbol{\xi}_{i}\in\mathbb{R}^{nD\times
D}.$

SINDy assumes that the expansion of $\Theta(\mathbf{Y})$ is the same for all
time, it simplifies the structure of $\Xi$ down to

$\Xi=\begin{bmatrix}\boldsymbol{\xi}&\boldsymbol{\xi}&\cdots&\boldsymbol{\xi}\end{bmatrix}=\boldsymbol{\xi}\cdot\mathbf{1},\quad\boldsymbol{\xi}\in\mathbb{R}^{nD\times
1}\quad\text{and}\quad\mathbf{1}=\begin{bmatrix}1&1&\cdots&1\end{bmatrix}\in\mathbb{R}^{1\times
D}.$

Furthermore, SINDy assumes that $\mathbf{f}(\mathbf{Y}(t))$ admits a spare
representation in $\Theta(\mathbf{Y})$. Moreover, it finds a parsimonious
model by performing least squares regression with sparsity-promoting
regularization, i.e.

$\hat{\boldsymbol{\xi}}=\underset{\boldsymbol{\xi}\in\mathbb{R}^{nD\times
1}}{\operatorname{arg}\operatorname{min}}\;\\{|\dot{\mathbf{Y}}-\Theta(\mathbf{Y})(\boldsymbol{\xi}\cdot\mathbf{1})|_{\text{Frob}}^{2}+\lambda|\boldsymbol{\xi}|_{1}\\}.$

The SINDy approach has been shown to be effective for observational data from
various kinds of dynamics, especially effective for the Lorenz system (cite).
The assumption that the right hand side function can be approximated well by a
set of predetermined dictionary opens up a question on how to optimally choose
such a dictionary. In the case of collective dynamics, a direct application is
ineffective. Although the form of collective dynamics can be written as
polynomials in components of $\mathbf{y}_{i}$, the coefficients are changing
all the time. In order to use it for collective dynamics, one has to re-work
the approximation scheme, see [23] for an example.

We can also use neural network structures to solve and infer information about
dynamical systems. Given the observation data
$\\{\mathbf{Y}(t_{l}),\dot{\mathbf{Y}}(t_{l})\\}_{l=1}^{L}$, we find the
estimator to $\mathbf{f}$ as in $\dot{\mathbf{Y}}=\mathbf{f}(\mathbf{Y})$ from
minimizing the following loss functional

$\text{Loss}(\tilde{\mathbf{f}})=\frac{1}{L}\sum_{l=1}^{L}|\dot{\mathbf{Y}}(t_{l})-\tilde{\mathbf{f}}(\mathbf{Y}(t_{l}))|_{\ell^{2}(\mathbb{R}^{D})}^{2},\quad\tilde{\mathbf{f}}\in\mathcal{H}.$

Here $\tilde{\mathbf{f}}:\mathbb{R}^{D}\rightarrow\mathbb{R}^{D}$ and
$\mathcal{H}$ is a set of neural networks of the same depth, same neurons on
each hidden layer, and same activation function. The neural network solution
is

$\hat{\mathbf{f}}=\underset{\tilde{\mathbf{f}}\in\mathcal{H}}{\operatorname{arg}\operatorname{min}}\;\text{Loss}(\tilde{\mathbf{f}}).$

Neural network approximation is shown to be effective for high-dimensional
functions estimation. However, the right hand side $\mathbf{f}(\mathbf{Y})$
might have a special structure that an usual neural network might fail to
capture. In this case, NeuralODE [4], which uses the ODE solver and Recurrent
Neural Network structure to train, can be used to capture special properties
of the ODE system. But in the case of collective dynamics, when $\mathbf{f}$
on $\mathbf{Y}$ is learned directly, we restrict ourselves to a fixed system
(fixed number of agents), recall that $\mathbf{Y}$ is a concatenation of all
system state variables. In other words, if we switch to a system where the
number of agents is different from training data, we have to re-train. PINN
(Physics Informed Neural Networks [3]), on the other hand, may be adopted for
such special need. As it can learn both the right hand side and the ODE
solution together.

###  3.2 Heterogeneous Agents

It is natural to consider heterogeneous agents, i.e. agents of different
types, for there are many dynamics that have multiple types of agents
involved, e.g. leader-follower, predator-prey, pedestrian-vehicle, etc. We
consider the scenario that the system of $N$ agents is partitioned into $K$
types, i.e. $C_{k_{1}}\cap C_{k_{2}}=\emptyset$ for $1\leq k_{1}\neq k_{2}=K$
and $\cup_{k=1}^{K}C_{K}=[N]=\\{1,\cdots,N\\}$. Moreover the type of the
agents does not change in time444Such restriction can be relaxed, as long as
we know the type information for the agents at all time., and there is a type
function $\mathit{k}_{\cdot}:\mathbb{N}\rightarrow\mathbb{N}$, which returns
the type index for an agent index.

(a) $\mathbf{X}(t)$ vs $\hat{\mathbf{X}}(t)$

(b) $\phi_{k_{1},k_{2}}$ vs $\hat{\phi}_{k_{1},k_{2}}$

Figure 2: Predator-Preys Dynamics with $1:\text{prey}$, $2:\text{predator}$,
$N_{1}=19$ preys, $N_{2}=1$ predator, learned on $[0,5]$, for other parameters
see [5].

Then we can consider a simple extension of (4) as follows

$\dot{\mathbf{x}}_{i}=\sum_{i^{\prime}=1,i^{\prime}\neq
i}^{N}\frac{1}{N_{\mathit{k}_{i}}}\phi_{\mathit{k}_{i},\mathit{k}_{i^{\prime}}}(|\mathbf{x}_{i^{\prime}}-\mathbf{x}_{i}|)(\mathbf{x}_{i^{\prime}}-\mathbf{x}_{i}),\quad
i=1,\cdots,N.$ (8)

Now the interaction kernel $\phi$ does not only depend on the pairwise
distance variable $r$, but is also different for agents in different types.
Given the same observation data, i.e.
$\\{\mathbf{x}^{m}_{i}(t_{l}),\dot{\mathbf{x}}^{m}_{i}(t_{l})\\}_{i,l,m=1}^{N,L,M}$,
we will learn a total of $K^{2}$ interaction kernels, i.e.
$\\{\phi_{k_{1},k_{2}}\\}_{k_{1},k_{2}=1}^{K}$, all together by minimizing a
slightly updated loss functional,

$\mathcal{E}^{\text{Mul}}_{L,M,\boldsymbol{\mathcal{H}}}(\boldsymbol{\varphi})=\frac{1}{LM}\sum_{l,m=1}^{L,M}|\dot{\mathbf{X}}_{l}^{m}-\mathbf{f}_{\boldsymbol{\varphi}}(\mathbf{X}_{l}^{m})|_{\mathcal{S}^{\prime}}^{2},$
(9)

Here
$\boldsymbol{\varphi}=\\{\varphi_{k_{1},k_{2}}\in\mathcal{H}_{k_{1},k_{2}}\\}_{k_{1},k_{2}=1}^{K}$,
the direct sum space
$\boldsymbol{\mathcal{H}}=\oplus_{k_{1},k_{2}=1}^{K}\mathcal{H}_{k_{1},k_{2}}$,
$\mathbf{X}_{l}^{m}$ is the usual vector notation and

$\mathbf{f}_{\boldsymbol{\varphi}}(\mathbf{X}_{l}^{m})=\begin{bmatrix}\vdots\\\
\sum_{i^{\prime}=1,i^{\prime}\neq
i}^{N}\frac{1}{N_{\mathit{k}_{i}}}\phi_{\mathit{k}_{i},\mathit{k}_{i^{\prime}}}(|\mathbf{x}_{i^{\prime}}^{m}(t_{l})-\mathbf{x}_{i}^{m}(t_{l})|)(\mathbf{x}_{i^{\prime}}^{m}(t_{l})-\mathbf{x}_{i}^{m}(t_{l}))\\\
\vdots\end{bmatrix}.$

The norm, $|\cdot|_{\mathcal{S}^{\prime}}$, is changed slightly as

$|\mathbf{U}-\mathbf{V}|_{\mathcal{S}^{\prime}}^{2}=\sum_{i=1}^{N}\frac{1}{N_{\mathit{k}_{i}}}|\mathbf{u}_{i}-\mathbf{v}_{i}|^{2}_{\ell^{2}(\mathbb{R}^{d})},\quad\mathbf{U},\mathbf{V}\in\mathbb{R}^{D=Nd}.$

For proof of convergence, it follows similarly from the proof presented in
[5].

###  3.3 Riemannian Geometry Constraints

The state variable, i.e. $\mathbf{x}_{i}\in\mathbb{R}^{d}$, might be moving on
a low dimensional manifold as $d$ gets large. We consider the case when
$\mathbf{x}_{i}\in\mathcal{M}\subset\mathbb{R}^{d}$, where $\mathcal{M}$ is a
Riemannian manifold naturally embedded in $\mathbb{R}^{d}$. We further assume
that we know enough information about the manifold, i.e. the pair
$(\mathcal{M},g)$, with $g$ being the Riemannian metric, is given to us. We
also consider that the first order model is updated as

$\mathbf{v}_{i}=\sum_{i^{\prime}=1,i^{\prime}\neq
i}^{N}\frac{1}{N_{\mathit{k}_{i}}}\phi_{\mathit{k}_{i},\mathit{k}_{i^{\prime}}}(|\mathbf{x}_{i^{\prime}}-\mathbf{x}_{i}|)\mathbf{w}^{\mathbf{x}}(\mathbf{x}_{i},\mathbf{x}_{i^{\prime}}),\quad
i=1,\cdots,N.$ (10)

Here
$\mathbf{v}_{i}=\dot{\mathbf{x}}_{i}\in\mathcal{T}_{\mathbf{x}_{i}}\mathcal{M}$
(where $\mathcal{T}_{\mathbf{x}}\mathcal{M}$ is the tangent space of
$\mathcal{M}$ at $\mathbf{x}$), and
$\mathbf{w}^{\mathbf{x}}(\mathbf{x}_{i},\mathbf{x}_{i^{\prime}})$ gives the
unit tangent direction on the geodesic from $\mathbf{x}_{i}$ to
$\mathbf{x}_{i^{\prime}}$ if $\mathbf{x}_{i^{\prime}}$ is not in the cut
locust set of $\mathbf{x}_{i}$, otherwise
$\mathbf{w}^{\mathbf{x}}(\mathbf{x}_{i},\mathbf{x}_{i^{\prime}})=\mathbf{0}$.
We also assume that each $\phi_{k_{1},k_{2}}$ is defined on $[0,R]$, where $R$
is sufficiently small so that length-minimizing geodesics exists uniquely.
Then such first order model is a well defined gradient flow model of a sum of
pairwise potential energy. In order to respect the geometry, we update the
loss functional so that it has the information about the Riemannian metric
information, i.e.

$\mathcal{E}^{\text{Rie}}_{L,M,\boldsymbol{\mathcal{H}}}(\boldsymbol{\varphi})=\frac{1}{LM}\sum_{l,m=1}^{L,M}|\dot{\mathbf{X}}_{l}^{m}-\mathbf{f}_{\boldsymbol{\varphi}}(\mathbf{X}_{l}^{m})|_{g}^{2},$

where the new norm, $|\cdot|_{g}$, is given as

$|\dot{\mathbf{X}}_{l}^{m}-\mathbf{f}_{\boldsymbol{\varphi}}(\mathbf{X}_{l}^{m})|_{g}^{2}=\sum_{i=1}^{N}\frac{1}{N_{\mathit{k}_{i}}}|\dot{\mathbf{x}}_{i,l}^{m}-\sum_{i^{\prime}=1,i^{\prime}\neq
i}^{N}\frac{1}{N_{\mathit{k}_{i^{\prime}}}}\varphi_{\mathit{k}_{i},\mathit{k}_{i^{\prime}}}(r_{i,i^{\prime},l}^{m})\mathbf{w}^{m}_{i,i^{\prime},l}|_{\mathcal{T}_{\mathbf{x}_{i,l}^{m}}\mathcal{M}}^{2}.$

Here $\mathbf{x}_{i,l}^{m}=\mathbf{x}_{i}^{m}(t_{l})$,
$r_{i,i^{\prime},l}^{m}=|\mathbf{x}_{i^{\prime}}^{m}(t_{l})-\mathbf{x}_{i}^{m}(t_{l})|$,
and
$\mathbf{w}^{m}_{i,i^{\prime},l}=\mathbf{w}^{\mathbf{x}}(\mathbf{x}_{i,l}^{m},\mathbf{x}_{i^{\prime},l}^{m})$.
With this new norm, we are able to preserve the convergence rate as if the
data is in Euclidean space, see [10] for details.

###  3.4 Feature Map Learning

Although there are many interacting-agent systems that demonstrate various
complex emergent behaviors that are only involved with the interaction kernel
functions that depend on the pairwise distance as we described above, however,
it may depend on some other unknown but a very small number of variables (e.g.
pairwise distances/angles of velocities, pairwise differences of phases,
etc…), to solve for this case, we consider the interacting-agent systems of
$N$ agents governed by the equations in the form 2, i.e.,

$\dot{\mathbf{x}}_{i}(t)=\frac{1}{N}\sum_{i^{\prime}=1}^{N}\Phi(\mathbf{x}_{i}(t),\mathbf{x}_{i^{\prime}}(t))(\mathbf{x}_{i^{\prime}}(t)-\mathbf{x}_{i}(t))\,,$
(11)

for $i=1,\ldots,N$, $t\in[0,T]$; $\mathbf{x}_{i}\in\mathbb{R}^{d}$ is a state
vector, and $\Phi:\mathbb{R}^{d}\times\mathbb{R}^{d}\rightarrow\mathbb{R}$ is
the interaction kernel, governing how the state of agent $i^{\prime}$
influences the state of agent $i$. Note that while the state space of the
system is $dN$-dimensional, the interaction kernel is a function of $2d$
dimensions. And we consider the case that the interaction kernel $\Phi$ is a
function that depends on a smaller number of natural variables
$\boldsymbol{\xi}\in\mathbb{R}^{d^{\prime}}$, with $d^{\prime}\ll 2d$, which
are functions of pairs
$(\mathbf{x}_{i},\mathbf{x}_{i^{\prime}})\in\mathbb{R}^{d}\times\mathbb{R}^{d}$.

(a) $\mathbf{X}(t)$ vs $\hat{\mathbf{X}}(t)$

(b) $\phi$ vs $\hat{\phi}$ (with or without $B$)

Figure 3: Power Law Dynamics with $N=3$ agents (new algorithm), for other
parameters see [11].

In other words, the interaction kernel $\Phi$ can be factorized (as
composition of functions) as a reduced interaction kernel
$\phi:\mathbb{R}^{d^{\prime}}\rightarrow\mathbb{R}$, and reduced variables
$\boldsymbol{\xi}:\mathbb{R}^{d}\times\mathbb{R}^{d}\rightarrow\mathbb{R}^{d^{\prime}}$,
i.e.

$\Phi(\mathbf{x}_{i},\mathbf{x}_{i^{\prime}})=\phi(\boldsymbol{\xi}(\mathbf{x}_{i},\mathbf{x}_{i^{\prime}})).$
(12)

When
$\boldsymbol{\xi}(\mathbf{x}_{i},\mathbf{x}_{i^{\prime}}):=||\mathbf{x}_{i}-\mathbf{x}_{i^{\prime}}||$,
$\Phi(\mathbf{x}_{i},\mathbf{x}_{i^{\prime}})=\phi(\boldsymbol{\xi}(\mathbf{x}_{i},\mathbf{x}_{i^{\prime}}))=\phi(|\mathbf{x}_{i}-\mathbf{x}_{i^{\prime}}|)$,
with $\phi:\mathbb{R}_{+}\rightarrow\mathbb{R}$, it becomes the simple model
which is a special case of what we discussed above.

Given observations of trajectory data, i.e.
$\\{\mathbf{x}_{i}^{(m)}(t_{l}),\dot{\mathbf{x}}_{i}^{(m)}(t_{l})\\}_{i,l,m=1}^{N,L,M}$,
with $0=t_{1}<\cdots<t_{L}=T$, initial conditions
$\mathbf{x}_{i}^{(m)}(t_{1})\sim\mu_{0}$ (some unknown probability
distribution on the state space $\mathbb{R}^{dN}$) and $m=1,\ldots,M$, we are
interested in estimating the interaction kernel $\Phi$, the only unknown
quantity in (2). Note that this is again not a regression problem, but an
inverse problem, since values of the function $\Phi$ are not observed, but
only averages of such values, as per the right-hand side of (2), are.

We proceed by further assuming that $\boldsymbol{\xi}$ can be factorized as
$\boldsymbol{\xi}(\mathbf{x}_{i},\mathbf{x}_{i^{\prime}}):=B\mathbf{z}(\mathbf{x}_{i},\mathbf{x}_{i^{\prime}})$,
where $\mathbf{z}:\mathbb{R}^{d}\times\mathbb{R}^{d}\rightarrow\mathbb{R}^{D}$
is a known feature map from pairs of states to high-dimensional feature
vectors, and the feature reduction map
$B:\mathbb{R}^{D}\rightarrow\mathbb{R}^{d^{\prime}}$ is an unknown linear map
to be estimated. Since we will estimate the reduced interaction kernel $\phi$
in a non-parametric fashion, what really matters is the range of $B$, and so
we may assume $B$ to be orthogonal. Therefore, the factorization (12) may be
rewritten in the form
$\Phi(\mathbf{x}_{i},\mathbf{x}_{i^{\prime}})=\phi(B\mathbf{z}(\mathbf{x}_{i},\mathbf{x}_{i^{\prime}}))$.
To estimate $\phi$ and $B$ we proceed in two steps.

Step $1$. Estimating the feature reduction map $B$, and the variables
${\boldsymbol{\xi}}$. While the choice of the feature map $\mathbf{z}$ is
typically application-dependent, and can incorporate symmetries, or physical
constraints on the system, a rather canonical choice is the map to polynomials
in the states, and here we restrict ourselves to second-order polynomials, and
therefore assume that:

$\displaystyle\mathbf{z}(\mathbf{x}_{i},\mathbf{x}_{i^{\prime}})$ (13)
$\displaystyle\coloneqq\big{[}\mathbf{x}_{i},\mathbf{x}_{i^{\prime}},((\mathbf{x}_{i})_{j}(\mathbf{x}_{i})_{j^{\prime}})_{\begin{subarray}{c}j\leq
j^{\prime}\end{subarray}},(({\mathbf{x}_{i^{\prime}})_{j}}(\mathbf{x}_{i^{\prime}})_{j^{\prime}})_{\begin{subarray}{c}j\leq
j^{\prime}\end{subarray}},((\mathbf{x}_{i})_{j}(\mathbf{x}_{i^{\prime}})_{j^{\prime}})_{\begin{subarray}{c}j,j^{\prime}=1,\dots,d\end{subarray}}\big{]}^{T}$

Here $\mathbf{z}\in\mathbb{R}^{D=2d^{2}+3d}$. Higher-order polynomials can be
added for complicated interactions and how to optimally choose such polynomial
basis will be for future study.

The techniques of [24], similar to the previously existing multi-index
regression works, are only applicable to the regression setting, and therefore
not directly applicable. However, at this point we note that if we allow the
learning procedure to conduct “experiments” and have access to a simulator of
(2) for different values of $N$, and in particular for $N=2$, it can conduct
observations of trajectories of two agents, and from those equations that
reveal the values taken by the interaction kernels can be extracted, therefore
taking us back to a regression setting, and apply the Multiplicatively
Perturbed Least Squares approach of [24] to obtain an estimate $\hat{B}$ of
the feature reduction map $B$. While the further study on how to optimally add
the higher order terms in (13) is ongoing, the techniques we use to estimate
$\phi$ follow those in [24], which are not cursed by the dimension $D$, making
it possible to increase the dimension of (the range of) the feature map
$\mathbf{y}$ with relatively small additional sampling requirements.

Step $2$. Estimating the reduced interaction kernel $\phi$. Once an estimate
$\hat{B}$ for the feature reduction map $B$ has been constructed, we proceed
to estimate $\phi$. Always in the case $N=2$, we proceed by projecting the
pairs of states, using the estimated feature reduction map $\hat{B}$, to
$\mathbb{R}^{d^{\prime}}$, and use a non-parametric regression technique on
that subspace, aimed at minimizing, over a suitable set of functions
$\psi\in\mathcal{H}$, the error functional

$\mathcal{E}_{L,M,\mathcal{H}}(\varphi)\coloneqq\frac{1}{LM}\sum_{l,m=1}^{L,M}|\dot{\mathbf{X}}_{l}^{m}-\mathbf{f}_{\varphi}(\mathbf{X}_{l}^{m})|_{\mathcal{S}}^{2},$

where

$\mathbf{f}_{\varphi}(\mathbf{X}_{l}^{m})=\begin{bmatrix}\vdots\\\
\frac{1}{N}\sum_{i^{\prime}=1,i^{\prime}\neq
i}^{N}\varphi(\hat{B}\mathbf{z}_{i,i^{\prime},l}^{m})(\mathbf{x}_{i^{\prime}}^{m}(t_{l})-\mathbf{x}_{i}^{m}(t_{l}))\\\
\vdots\end{bmatrix},$

where
$\mathbf{z}_{i,i^{\prime},l}^{m}=\mathbf{z}(\mathbf{x}_{i}^{m}(t_{l}),\mathbf{x}_{i^{\prime}}^{m}(t_{l}))$.
We will choose $\mathcal{H}$ to be a convex and compact (in the $L^{\infty}$
norm) subset of a subspace of functions of the estimated variables
$\hat{B}\mathbf{y}$, e.g. spanned by splines with knots on a grid, or
piecewise polynomials. The dimension of $\mathcal{H}$ will be chosen as a
suitably increasing function of the number of training trajectories $M$,
following the ideas of [5] and [12] (in the case of multiple reduced
variables), see [11] for detailed discussion.

###  3.5 Coupled Systems

We consider two types of coupled systems in order to introduce more involved
sets of states. First, we take the state variable $\mathbf{y}_{i}$ as
$\mathbf{y}_{i}=\begin{bmatrix}\mathbf{x}_{i}\\\ \zeta_{i}\end{bmatrix}$ with
$\mathbf{x}_{i}\in\mathbb{R}^{d}$ and $\zeta_{i}\in\mathbb{R}$ (we can also
take $\zeta_{i}\in\mathbb{R}^{d^{\prime}}$ for any $d^{\prime}\geq 1$), hence
$\mathbf{y}_{i}\in\mathbb{R}^{d+1}$. $\mathbf{x}_{i}$ usually describes the
position and $\zeta_{i}$ for phase/excitation/opinion (e.g. the swarmalator
model in cite three papers). The interaction comes not only from each
$\mathbf{y}_{i}$ on the agents, but also from the intra-agent interaction of
$\mathbf{x}_{i}$ and $\zeta_{i}$, which gives rise to a set more complicated
patterns. The change of $\mathbf{y}_{i}$ is governed by the following ODE
systems

$\begin{cases}\dot{\mathbf{x}}_{i}&=F^{\mathbf{x}}(\mathbf{y}_{i})+\sum_{i^{\prime}=1,i^{\prime}\neq
i}^{N}\frac{1}{N_{\mathit{k}_{i}}}\phi^{E}_{\mathit{k}_{i},\mathit{k}_{i^{\prime}}}(\boldsymbol{\xi}^{E}(\mathbf{y}_{i},\mathbf{y}_{i^{\prime}}))(\mathbf{x}_{i^{\prime}}-\mathbf{x}_{i})\\\
\dot{\zeta}_{i}&=F^{\zeta}(\mathbf{y}_{i})+\sum_{i^{\prime}=1,i^{\prime}\neq
i}^{N}\frac{1}{N_{\mathit{k}_{i}}}\phi^{\zeta}_{\mathit{k}_{i},\mathit{k}_{i^{\prime}}}(\boldsymbol{\xi}^{\zeta}(\mathbf{y}_{i},\mathbf{y}_{i^{\prime}}))(\zeta_{i^{\prime}}-\zeta_{i})\end{cases},\quad
i=1,\cdots,N.$ (14)

Here $F^{\mathbf{x}}:\mathbb{R}^{d+1}\rightarrow\mathbb{R}^{d}$ is a force
which governs how the interaction of the $i^{th}$ agent and its surrounding
environment affects the change of $\mathbf{x}_{i}$, similarly for
$F^{\zeta}:\mathbb{R}^{d+1}\rightarrow\mathbb{R}$ on $\zeta_{i}$,
$\boldsymbol{\xi}^{E}:\mathbb{R}^{2d+2}\rightarrow\mathbb{R}^{d^{E}}$ is an
energy based reduced variable ($1\leq d^{E}\ll 2d+2$),
$\phi^{E}:\mathbb{R}^{d^{E}}\rightarrow\mathbb{R}$ is the energy based reduced
interaction kernel,
$\boldsymbol{\xi}^{\zeta}:\mathbb{R}^{2d+2}\rightarrow\mathbb{R}^{d^{\zeta}}$
is an $\zeta$-based reduced variable and $\phi^{\zeta}$ is its corresponding
reduced interaction kernel. Given the observation
$\\{\mathbf{y}_{i}^{m}(t_{l}),\dot{\mathbf{y}}_{i}^{m}(t_{l})\\}_{i,l,m=1}^{N,L,M}$,
we find the set of reduced interaction kernels555Here we assume that the
reduced variables $\boldsymbol{\xi}^{E}$ and $\boldsymbol{\xi}^{\zeta}$ are
known to us. by minizming the following loss functionals, first for
$\boldsymbol{\varphi}^{E}=\\{\varphi^{E}_{k_{1},k_{2}}\\}_{k_{1},k_{2}}^{K}$,

$\mathcal{E}^{\mathbf{x}}_{L,M,\boldsymbol{\mathcal{H}}^{E}}(\boldsymbol{\varphi}^{E})=\frac{1}{LM}\sum_{l,m=1}^{L,M}|\dot{\mathbf{Y}}_{l}^{m}-\mathbf{f}_{\boldsymbol{\varphi}^{E}}(\mathbf{Y}_{l}^{m})|_{\mathcal{S}^{\prime}}^{2}$

for
$\boldsymbol{\varphi}^{E}=\\{\varphi^{E}_{k_{1},k_{2}}\in\mathcal{H}^{E}_{k_{1},k_{2}}\\}_{k_{1},k_{2}}^{K}$
and
$\boldsymbol{\mathcal{H}}^{E}=\oplus_{k_{1},k_{2}}^{K}\mathcal{H}^{E}_{k_{1},k_{2}}$,

$\mathbf{Y}_{l}^{m}=\begin{bmatrix}\vdots\\\ \mathbf{x}_{i}^{m}(t_{l})\\\
\xi_{i}^{m}(t_{l})\\\ \vdots\end{bmatrix}$

and

$\mathbf{f}_{\boldsymbol{\varphi}^{E}}(\mathbf{X}_{l}^{m})=\begin{bmatrix}\vdots\\\
F^{\mathbf{x}}(\mathbf{y}_{i,l}^{m})+\sum_{i^{\prime}=1,i^{\prime}\neq
i}^{N}\frac{1}{N_{\mathit{k}_{i}}}\phi^{E}_{\mathit{k}_{i},\mathit{k}_{i^{\prime}}}(\boldsymbol{\xi}^{E}(\mathbf{y}_{i,l}^{m},\mathbf{y}_{i^{\prime},l}^{m}))(\mathbf{x}_{i^{\prime},l}^{m}-\mathbf{x}_{i,l}^{m})\\\
\vdots\end{bmatrix}$

Here $\mathbf{x}_{i,l}^{m}=\mathbf{x}_{i}^{m}(t_{l})$ and
$\mathbf{y}_{i,l}^{m}=\mathbf{y}_{i}^{m}(t_{l})$. Similarly for learning
$\boldsymbol{\varphi}^{\zeta}=\\{\varphi^{\zeta}_{k_{1},k_{2}}\\}_{k_{1},k_{2}}^{K}$,

$\mathcal{E}^{\zeta}_{L,M,\boldsymbol{\mathcal{H}}^{\zeta}}(\boldsymbol{\varphi}^{\zeta})=\frac{1}{LM}\sum_{l,m=1}^{L,M}|\dot{\mathbf{Y}}_{l}^{m}-\mathbf{f}_{\boldsymbol{\varphi}^{\zeta}}(\mathbf{Y}_{l}^{m})|_{\mathcal{S}^{\prime}}^{2}$

for
$\boldsymbol{\varphi}^{\zeta}=\\{\varphi^{\zeta}_{k_{1},k_{2}}\in\mathcal{H}^{\zeta}_{k_{1},k_{2}}\\}_{k_{1},k_{2}}^{K}$
and
$\boldsymbol{\mathcal{H}}^{E}=\oplus_{k_{1},k_{2}}^{K}\mathcal{H}_{k_{1},k_{2}}^{\zeta}$
and

$\mathbf{f}_{\boldsymbol{\varphi}^{\zeta}}(\mathbf{Y}_{l}^{m})=\begin{bmatrix}\vdots\\\
F^{\zeta}(\mathbf{y}_{i,l}^{m})+\sum_{i^{\prime}=1,i^{\prime}\neq
i}^{N}\frac{1}{N_{\mathit{k}_{i}}}\phi^{\zeta}_{\mathit{k}_{i},\mathit{k}_{i^{\prime}}}(\boldsymbol{\xi}^{\zeta}(\mathbf{y}_{i,l}^{m},\mathbf{y}_{i^{\prime},l}^{m}))(\zeta_{i^{\prime},l}^{m}-\zeta_{i,l}^{m})\\\
\vdots\end{bmatrix}$

Here $\zeta_{i,l}^{m}=\zeta_{i}^{m}(t_{l})$.

Next, we take the state variable $\mathbf{y}_{i}$ as
$\mathbf{y}_{i}=\begin{bmatrix}\mathbf{x}_{i}\\\ \mathbf{v}_{i}\\\
\zeta_{i}\end{bmatrix}$, here $\mathbf{x}_{i},\mathbf{v}_{i}\in\mathbb{R}^{d}$
and $\zeta_{i}\in\mathbb{R}$ and we specifically require that
$\mathbf{v}_{i}=\dot{\mathbf{x}}_{i}$. The change of $\mathbf{y}_{i}$ is
governed by the following second-order ODE system

$\begin{cases}\dot{\mathbf{x}}_{i}&=\mathbf{v}_{i},\\\
m_{i}\dot{\mathbf{v}}_{i}&=F^{\mathbf{v}}(\mathbf{y}_{i})+\sum_{i^{\prime}=1,i^{\prime}\neq
i}^{N}\frac{1}{N_{\mathit{k}_{i}}}\Big{[}\phi^{E}_{\mathit{k}_{i},\mathit{k}_{i^{\prime}}}(\boldsymbol{\xi}^{E}(\mathbf{y}_{i},\mathbf{y}_{i^{\prime}}))(\mathbf{x}_{i^{\prime}}-\mathbf{x}_{i})\\\
&\quad+\phi^{A}_{\mathit{k}_{i},\mathit{k}_{i^{\prime}}}(\boldsymbol{\xi}^{A}(\mathbf{y}_{i},\mathbf{y}_{i^{\prime}}))(\mathbf{x}_{i^{\prime}}-\mathbf{x}_{i})\Big{]}\\\
\dot{\zeta}_{i}&=F^{\zeta}(\mathbf{y}_{i})+\sum_{i^{\prime}=1,i^{\prime}\neq
i}^{N}\frac{1}{N_{\mathit{k}_{i}}}\phi^{\zeta}_{\mathit{k}_{i},\mathit{k}_{i^{\prime}}}(\boldsymbol{\xi}^{\zeta}(\mathbf{y}_{i},\mathbf{y}_{i^{\prime}}))(\zeta_{i^{\prime}}-\zeta_{i})\end{cases},\quad
i=1,\cdots,N.$ (15)

Similarly $F^{\mathbf{v}}:\mathbb{R}^{2d+1}\rightarrow\mathbb{R}^{d}$ is a
force which governs how the interaction of the $i^{th}$ agent and its
surrounding environment affects the change of $\mathbf{v}_{i}$, similarly for
$F^{\zeta}:\mathbb{R}^{2d+1}\rightarrow\mathbb{R}$ on $\zeta_{i}$,
$\boldsymbol{\xi}^{E}:\mathbb{R}^{4d+2}\rightarrow\mathbb{R}^{d^{E}}$ is an
energy based reduced variable ($1\leq d^{E}\ll 4d+2$),
$\phi^{E}:\mathbb{R}^{d^{E}}\rightarrow\mathbb{R}$ is the energy based reduced
interaction kernel,
$\boldsymbol{\xi}^{A}:\mathbb{R}^{4d+2}\rightarrow\mathbb{R}^{d^{A}}$ is an
alignment based reduced variable ($1\leq d^{A}\ll 4d+2$),
$\phi^{A}:\mathbb{R}^{d^{A}}\rightarrow\mathbb{R}$ is the alignment based
reduced interaction
kernel$\boldsymbol{\xi}^{\zeta}:\mathbb{R}^{4d+2}\rightarrow\mathbb{R}^{d^{\zeta}}$
is an $\zeta$-based reduced variable and $\phi^{\zeta}$ is its corresponding
reduced interaction kernel. Given

In order to learn
$\boldsymbol{\varphi}^{E}=\\{\varphi^{E}_{k_{1},k_{2}}\\}_{k_{1},k_{2}}^{K}$
and
$\boldsymbol{\varphi}^{A}=\\{\varphi^{A}_{k_{1},k_{2}}\\}_{k_{1},k_{2}}^{K}$,
we use a slightly updated loss functional

$\mathcal{E}^{E+A}_{L,M,\boldsymbol{\mathcal{H}}^{E},\boldsymbol{\mathcal{H}}^{A}}(\boldsymbol{\varphi}^{E},\boldsymbol{\varphi}^{A})=\frac{1}{LM}\sum_{l,m=1}^{L,M}|\dot{\mathbf{Y}}_{l}^{m}-\mathbf{f}_{\boldsymbol{\varphi}^{E},\boldsymbol{\varphi}^{A}}(\mathbf{Y}_{l}^{m})|_{\mathcal{S}^{\prime}}^{2}$

where

$\mathbf{f}_{\boldsymbol{\varphi}^{E},\boldsymbol{\varphi}^{A}}(\mathbf{Y}_{l}^{m})=\begin{bmatrix}\vdots\\\
F^{\mathbf{v}}(\mathbf{y}_{i,l}^{m})+\sum_{i^{\prime}=1,i^{\prime}\neq
i}^{N}\frac{1}{N_{\mathit{k}_{i}}}\Big{[}\phi^{E}_{\mathit{k}_{i},\mathit{k}_{i^{\prime}}}(\boldsymbol{\xi}^{E}(\mathbf{y}_{i,l}^{m},\mathbf{y}_{i^{\prime},l}^{m}))(\mathbf{x}_{i^{\prime},l}^{m}-\mathbf{x}_{i,l}^{m})\\\
+\phi^{A}_{\mathit{k}_{i},\mathit{k}_{i^{\prime}}}(\boldsymbol{\xi}^{A}(\mathbf{y}_{i,l}^{m},\mathbf{y}_{i^{\prime},l}^{m}))(\mathbf{v}_{i^{\prime},l}^{m}-\mathbf{v}_{i,l}^{m})\Big{]}\\\
\vdots\end{bmatrix}$

where $\mathbf{v}_{i,l}^{m}=\mathbf{v}_{i}^{m}(t_{l})$.

###  3.6 Learning with Gaussian Priors

Gaussian process regression (GPR) serves as a non-parametric Bayesian machine
learning technique designed for supervised learning, equipped with an inherent
framework for quantifying uncertainty. Consequently, GPR has found application
in the study of ordinary differential equations (ODEs), stochastic
differential equations (SDEs), and partial differential equations (PDEs) [25,
26, 27, 28, 29, 30, 31, 32, 33, 34, 35], resulting in more accurate and robust
models for dynamical systems. Given the unique characteristics of dynamical
data, it necessitates novel concepts and substantial efforts tailored to
specific types of dynamical systems and data regimes. In our work, we model
latent interaction kernels as Gaussian processes, embedding them with the
underlying structure of our governing equations, including translation and
rotational invariance. This distinguishes our approach from most other works,
which model state variables as Gaussian processes.

Despite the challenges involved, recent mathematical advancements have led to
the development of a general physical model based on Newton’s second law, such
as the methods we discussed in section 3.1.1. This model has been demonstrated
to accurately capture a wide range of collective behaviors. Specifically, the
model describes a system of N agents that interact according to a system of
ODEs, where for each agent $i=1,\cdots,N$:

$\displaystyle\begin{aligned}
&m_{i}\ddot{\mathbf{x}}_{i}=F_{i}(\mathbf{x}_{i},\dot{\mathbf{x}}_{i},\boldsymbol{\alpha}_{i})+\sum_{i^{\prime}=1}^{N}\frac{1}{N}\Big{[}\phi^{E}(|\mathbf{x}_{i^{\prime}}-\mathbf{x}_{i}|)(\mathbf{x}_{i^{\prime}}-\mathbf{x}_{i})\\\
&\qquad+{\phi^{A}}(|{\mathbf{x}_{i^{\prime}}-\mathbf{x}_{i}}|)(\dot{\mathbf{x}}_{i^{\prime}}-\dot{\mathbf{x}}_{i})\Big{]},\end{aligned}$
(16)

where $m_{i}\geq 0$ is the mass of the agent $i$;
$\ddot{\mathbf{x}}_{i}\in\mathbb{R}^{d}$ is the acceleration;
$\dot{\mathbf{x}}_{i}\in\mathbb{R}^{d}$ is the velocity;
$\mathbf{x}_{i}\in\mathbb{R}^{d}$ is the position of agent $i$; the first term
$F_{i}$ is a parametric function of position and velocities, modeling self-
propulsion and frictions of agent $i$ with the environment, and the scalar
parameters $\boldsymbol{\alpha}_{i}$ describing their strength;
$|\mathbf{x}_{j}-\mathbf{x}_{i}|$ is the Euclidean distance; and the 1D
functions $\phi^{E},\phi^{A}:\mathbb{R}^{+}\rightarrow\mathbb{R}$ are called
the energy and alignment-based radial interaction kernels respectively. The
$\phi^{E}$ term describes the alignment of positions based on the difference
of positions; the $\phi^{A}$ term describes the alignment of velocities based
on the difference of velocities. Our primary objective is to infer the
interaction kernels ${\boldsymbol{\phi}}=\\{{\phi^{E}},{\phi^{A}}\\}$ as well
as the unknown scalar parameters $\boldsymbol{\alpha}$ and potentially
$\boldsymbol{m}$ from the observed trajectory data. Subsequently, we utilize
the learned governing equations to make predictions regarding future events or
simulate new datasets.

To learn the model given by (16), we initiate the process by modeling the
interaction kernel functions ${\phi^{E}}$ and ${\phi^{A}}$ with the priors as
two independent Gaussian processes

${\phi^{E}}\sim\mathcal{GP}(0,K_{\theta^{E}}(r,r^{\prime})),\qquad{\phi^{A}}\sim\mathcal{GP}(0,K_{\theta^{A}}(r,r^{\prime})),$
(17)

where $K_{\theta^{E}}$, $K_{\theta^{A}}$ are covariance functions with
hyperparameters $\boldsymbol{\theta}=({\theta^{E}},{\theta^{A}})$.
$\boldsymbol{\theta}$ can either be chosen by the modeler or tuned via a data-
driven procedure discussed later.

Then given the noisy observational data
$\mathbb{Y}=[\mathbf{Y}^{(1,1)},\dots,\mathbf{Y}^{(M,L)}]^{T}\in\mathbb{R}^{dNML}$,
and
$\mathbb{Z}=[\mathbf{Z}^{(1,1)}_{\sigma^{2}},\dots,\mathbf{Z}^{(M,L)}_{\sigma^{2}}]^{T}\in\mathbb{R}^{dNML}$
where

$\mathbf{Y}^{(m,l)}=\mathbf{Y}^{(m)}(t_{l}):=\begin{bmatrix}\mathbf{X}(t)\\\
\dot{\mathbf{X}}(t)\end{bmatrix}\in\mathbb{R}^{2dN},$ (18)

and

$\boldsymbol{m}\mathbf{Z}^{(m,l)}_{\sigma^{2}}=\boldsymbol{m}\ddot{\mathbf{X}}^{(m,l)}_{\sigma^{2}}:=F_{\boldsymbol{\alpha}}(\mathbf{Y}^{(m,l)})+\mathbf{f}_{{\boldsymbol{\phi}}}(\mathbf{Y}^{(m,l)})+\epsilon^{(m,l)},$
(19)

with
$\mathbf{f}_{\boldsymbol{\phi}}(\mathbf{Y}(t))=\mathbf{f}_{{\phi^{E}},{\phi^{A}}}(\mathbf{Y}(t))$
represents the sum of energy and alignment-based interactions as in (16), and
$\epsilon^{(m,l)}\sim\mathcal{N}(0,\sigma^{2}I_{dN})$ are i.i.d noise, based
on the properties of Gaussian processes. With the priors of ${\phi^{E}}$,
${\phi^{A}}$, we can train the hyperparameters $\boldsymbol{m}$,
$\boldsymbol{\alpha}$, $\boldsymbol{\theta}$ and $\sigma$ by maximizing the
probability of the observational data, which is equivalent to minimizing the
negative log marginal likelihood (NLML) (see Chapter 4 in [36])

$\displaystyle-\log
p(\boldsymbol{m}\mathbb{Z}|\mathbb{Y},\boldsymbol{\alpha},\boldsymbol{\theta},\sigma^{2})$
$\displaystyle=\frac{1}{2}(\boldsymbol{m}\mathbb{Z}-F_{\boldsymbol{\alpha}}(\mathbb{Y}))^{T}(K_{\mathbf{f}_{{\boldsymbol{\phi}}}}(\mathbb{Y},\mathbb{Y};\boldsymbol{\theta})+\sigma^{2}I_{dNML})^{-1}(\boldsymbol{m}\mathbb{Z}-F_{\boldsymbol{\alpha}}(\mathbb{Y}))$
$\displaystyle\qquad+\frac{1}{2}\log|K_{\mathbf{f}_{{\boldsymbol{\phi}}}}(\mathbb{Y},\mathbb{Y};\boldsymbol{\theta})+\sigma^{2}I_{dNML}|+\frac{dNML}{2}\log
2\pi.$ (20)

and the posterior/predictive distribution for the interaction kernels
$\phi^{\mathrm{type}}(r^{*})$, $\mathrm{type}=E\text{ or }A$, can be obtained
by

$p(\phi^{\mathrm{type}}(r^{\ast})|\mathbb{Y},\mathbb{Z},r^{\ast})\sim\mathcal{N}(\bar{\phi}^{\mathrm{type}},var(\bar{\phi}^{\mathrm{type}})),$
(21)

where

$\bar{\phi}^{\mathrm{type}}=K_{\phi^{\mathrm{type}},\mathbf{f}_{\boldsymbol{\phi}}}(r^{\ast},\mathbb{Y})(K_{\mathbf{f}_{{\boldsymbol{\phi}}}}(\mathbb{Y},\mathbb{Y})+\sigma^{2}I_{dNML})^{-1}(\boldsymbol{m}\mathbb{Z}-F_{\boldsymbol{\alpha}}(\mathbb{Y})),$
(22)
$var(\bar{\phi}^{\mathrm{type}})=K_{\theta^{\mathrm{type}}}(r^{\ast},r^{\ast})-K_{\phi^{\mathrm{type}},\mathbf{f}_{\boldsymbol{\phi}}}(r^{\ast},\mathbb{Y})(K_{\mathbf{f}_{{\boldsymbol{\phi}}}}(\mathbb{Y},\mathbb{Y})+\sigma^{2}I_{dNML})^{-1}K_{\mathbf{f}_{\boldsymbol{\phi}},\phi^{\mathrm{type}}}(\mathbb{Y},r^{\ast}).$
(23)

and
$F_{\boldsymbol{\alpha}}(\mathbb{Y})=\mathrm{Vec}(\\{F_{\boldsymbol{\alpha}}(\mathbf{Y}^{(m,l)})\\}_{m=1,l=1}^{M,L})\in\mathbb{R}^{dNML}$,
$K_{\mathbf{f}_{{\boldsymbol{\phi}}}}(\mathbb{Y},\mathbb{Y};\theta)\in\mathbb{R}^{dNML\times
dNML}$ is the covariance matrix between
$\mathbf{f}_{{\boldsymbol{\phi}}}(\mathbb{Y})$ and
$\mathbf{f}_{{\boldsymbol{\phi}}}(\mathbb{Y})$, which can be compute
elementwisely based on the covariance functions $K_{\theta^{E}}$,
$K_{\theta^{A}}$,
$K_{\mathbf{f}_{\boldsymbol{\phi}},\phi^{\mathrm{type}}}(\mathbb{Y},r^{*})=K_{\phi^{\mathrm{type}},\mathbf{f}_{\phi}}(r^{*},\mathbb{Y})^{T}$
denotes the covariance matrix between
$\mathbf{f}_{{\boldsymbol{\phi}}}(\mathbb{Y})$ and
$\phi^{\mathrm{type}}(r^{*})$.

Note here the marginal likelihood does not simply favor the models that fit
the training data best, but induces an automatic trade-off between data-fit
and model complexity. We can use the posterior mean estimators of
${\boldsymbol{\phi}}$ in trajectory prediction by performing numerical
simulations of the equations

$\boldsymbol{m}\hat{\mathbf{Z}}(t)=F_{\boldsymbol{\hat{\alpha}}}(\mathbf{Y}(t))+\hat{\mathbf{f}}_{\bar{{\boldsymbol{\phi}}}}(\mathbf{Y}(t)).$
(24)

and the posterior variance $var(\bar{\phi}^{\mathrm{type}})$ can be used as a
good indicator for the uncertainty of the estimation
$\bar{\phi}^{\mathrm{type}}$ based on our Bayesian approach, see Figure 4 and
also other examples in [14, 15].

(a) $\phi^{E}$ vs $\hat{\phi}^{E}$

(b) $\phi^{A}$ vs $\hat{\phi}^{A}$

(c) $\mathbf{X}(t)$ vs $\hat{\mathbf{X}}(t)$

Figure 4: Flocking with external potential (FwEP) model [37] with
$\\{N,L,M,\sigma\\}=\\{20,6,3,0.01\\}$, and $\phi^{E}(r)=1$,
$\phi^{A}(r)=\frac{1}{(1+r^{2})^{1/2}}$, for other parameters see [12]. The
light blue regions are two-standard-deviation bands around the means which
indicate the uncertainty of the estimators.

In classical regression setting [36], there is an interesting link between GP
regression with the kernel ridge regression (KRR), where the posterior mean
can be viewed as a KRR estimator to solve a regularized least square empirical
risk functional. In our setting, we have noisy functional observations of the
interaction kernels, i.e., the
$\\{r_{\mathbb{X}},r_{\mathbb{V}},\mathbb{Z}\\}$ instead of the pairs
$\\{r_{\mathbb{X}},{\phi^{E}}(r_{\mathbb{X}}),{\phi^{A}}(r_{\mathbb{X}})\\}$,
where $r_{\mathbb{X}},r_{\mathbb{V}},\in\mathbb{R}^{MLN^{2}}$ are the sets
contains all the pairwise distances in $\mathbb{X}$, and $\mathbb{V}$, so we
face an inverse problem here, instead of a classical regression problem.
Thanks to the linearity of the inverse problem, we can still derive a
Representer theorem [38] that helps understand the role of the
hyperparameters:

######  Theorem 2 (Connection with KRR)

Given the training data $\\{\mathbb{Y},\mathbb{Z}\\}$, if the priors
${\phi^{E}}\sim\mathcal{GP}(0,\tilde{K}^{E})$,
${\phi^{A}}\sim\mathcal{GP}(0,\tilde{K}^{A})$ with
$\tilde{K}^{E}=\frac{\sigma^{2}K^{E}}{MNL\lambda^{E}}$,
$\tilde{K}^{A}=\frac{\sigma^{2}K^{A}}{MNL\lambda^{A}}$ for some
$\lambda^{E},\lambda^{A}>0$, then the posterior mean
$\bar{\boldsymbol{\phi}}=(\bar{\phi}^{E},\bar{\phi}^{A})$ in (22) coincides
with the minimizer,
${\boldsymbol{\phi}}_{\mathcal{H}_{{K}^{E}}\times\mathcal{H}_{{K}^{A}}}^{\lambda,M}$,
of the regularized empirical risk functional $\mathcal{E}^{\lambda,M}(\cdot)$
on $\mathcal{H}_{{K}^{E}}\times\mathcal{H}_{{K}^{A}}$ where
$\mathcal{E}^{\lambda,M}(\cdot)$ is defined by

$\displaystyle\mathcal{E}^{\lambda,M}({\boldsymbol{\varphi}}):$
$\displaystyle=\frac{1}{LM}\sum_{l=1,m=1}^{L,M}\|\mathbf{f}_{{\boldsymbol{\varphi}}}(\mathbf{Y}^{(m,l)})-\mathbf{Z}_{\sigma^{2}}^{(m,l)}\|^{2}+\lambda^{E}\|\varphi^{E}\|_{\mathcal{H}_{{K}^{E}}}^{2}+\lambda^{A}\|\varphi^{A}\|_{\mathcal{H}_{{K}^{A}}}^{2}.$
(25)

From the theorem, it is clear to see how hyperparameters affect the prediction
of interaction kernels: $\theta^{E}$, $\theta^{A}$, and $\sigma$ jointly
affect the choice of Mercer kernels and regularization constant. And we
establish a novel operator-theoretical framework that generalized the analysis
of kernel regression methods [36] and linear inverse problems to interacting
particle systems. Under Hölder type source conditions on $\phi$, we prove the
reconstruction error converges at an upper rate in $M$:

###### Theorem 3

If the system satisfies the coercivity condition on
$\mathcal{H}_{{K}^{E}}\times\mathcal{H}_{{K}^{A}}$, i.e.
$\forall{\boldsymbol{\varphi}}\in\mathcal{H}_{{K}^{E}}\times\mathcal{H}_{{K}^{A}}$,

$\displaystyle\|A{\boldsymbol{\varphi}}\|^{2}_{L^{2}(\rho_{\mathbf{Y}})}=\|\mathbf{f}_{{\boldsymbol{\varphi}}}\|^{2}_{L^{2}(\rho_{\mathbf{Y}})}\geq
c_{\mathcal{H}_{{K}^{E}}}\|\varphi^{E}\|^{2}_{L^{2}(\tilde{\rho}_{r}^{E})}+c_{\mathcal{H}_{{K}^{A}}}\|\varphi^{A}\|^{2}_{L^{2}(\tilde{\rho}_{r}^{A})}$
(26)

for some constants $c_{\mathcal{H}_{{K}^{E}}},c_{\mathcal{H}_{{K}^{A}}}>0$,
then if we choose $\lambda\asymp M^{-\frac{1}{2\gamma+2}}$ for some
$\gamma\in(0,\frac{1}{2}]$, for any $\delta\in(0,1)$, it holds with
probability at least $1-\delta$ that

$\|{\boldsymbol{\phi}}_{\mathcal{H}_{{K}^{E}}\times\mathcal{H}_{{K}^{A}}}^{\lambda,M}-{\boldsymbol{\phi}}\|_{\mathcal{H}_{{K}^{E}}\times\mathcal{H}_{{K}^{A}}}\lesssim
C\log(\frac{8}{\delta})M^{-\frac{\gamma}{2\gamma+1}},$

where $C$ is a constant depends on ${\boldsymbol{\phi}}$,
$\mathcal{H}_{{K}^{E}}\times\mathcal{H}_{{K}^{A}}$, $\sqrt{dNL}$ and $\sigma$.

see [14, 15] for detailed discuss.

##  4 Conclusion

We review a series of methods for learning collective behaviors, ranging from
first-order to second-order systems, from homogeneous types of agents to
heterogeneous types of agents, from dynamics constrained on Riemannian
manifolds, from dynamics where feature variables are unknown, and learning
with a Gaussian prior. Our methods are proven to be convergent and efficient
for handling high-dimensional observational data. There are many potential
future directions to explore with this series of methods. For example, how to
learn type information from data, how to learn geometry and dynamics from
data, how to handle features in second-order systems and beyond, how to deal
with stochastic forces in second-order systems, and how to learn the
interaction kernel when there is only one-time snapshot of the trajectory.
These ideas present challenging and yet intriguing perspectives for future
projects.

## Acknowledgement

JF and MZ contributed equally to research and writing of the manuscript. JF
and MZ would like to thank Mauro Maggioni (Johns Hopkins University) for the
invaluable discussion. MZ is partially supported by NSF-AoF-$2225507$ and
startup fund from Illinois Tech.

## References

  * [1] S. M. Stigler.  The History of Statistics: The Measurement of Uncertainty Before $1900$.  Harvard University Press, Cambridge, MA, 1st edition, 1986.
  * [2] 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.
  * [3] M. Raissi, P. Perdikaris, and G.E. Karniadakis.  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  Journal of Computational Physics, 378:686–707, 2019.
  * [4] Ricky T. Q. Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud.  Neural ordinary differential equations.  In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 31. Curran Associates, Inc., 2018.
  * [5] Fei Lu, Ming Zhong, Sui Tang, and Mauro Maggioni.  Nonparametric inference of interaction laws in systems of agents from trajectory data.  Proceedings of the National Academy of Sciences, 116(29):14424–14433, 2019.
  * [6] Eitan Tadmor.  On the mathematics of swarming: emergent behavior in alignment dynamics.  Notices of the AMS, 68(4):493 – 503, 2021.
  * [7] Eitan Tadmor.  Long time and large crowd dynamics of fully discrete cucker-smale alignment models.  Pure and Applied Functional Analysis, 8(2):603 – 626, 2023.
  * [8] Scott Camazine, Jean-Louis Deneubourg, Nigel R. Franks, James Sneyd, Guy Theraula, and Eric Bonabeau.  Self-Organization in Biological Systems.  Princeton University Press, NJ, 1st edition, 2003.
  * [9] Ming Zhong, Jason Miller, and Mauro Maggioni.  Data-driven discovery of emergent behaviors in collective dynamics.  Physica D: Nonlinear Phenomena, 411:132542, 2020.
  * [10] Mauro Maggioni, Jason J Miller, Hongda Qiu, and Ming Zhong.  Learning interaction kernels for agent systems on riemannian manifolds.  In Marina Meila and Tong Zhang, editors, Proceedings of the 38th International Conference on Machine Learning, volume 139 of Proceedings of Machine Learning Research, pages 7290–7300. PMLR, 18–24 Jul 2021.
  * [11] Jinchao Feng, Mauro Maggioni, Patrick Martin, and Ming Zhong.  Learning interaction variables and kernels from observations of agent-based systems.  IFAC-PapersOnLine, 55(30):162–167, 2022.  25th International Symposium on Mathematical Theory of Networks and Systems MTNS 2022.
  * [12] Jason Miller, Sui Tang, Ming Zhong, and Mauro Maggioni.  Learning theory for inferring interaction kernels in second-order interacting agent systems.  Sampling Theory, Signal Processing, and Data Analysis, 21(1):21, 2023.
  * [13] Ming Zhong, Jason Miller, and Mauro Maggioni.  Machine learning for discovering effective interaction kernels between celestial bodies from ephemerides, 2021.
  * [14] Jinchao Feng, Charles Kulick, Yunxiang Ren, and Sui Tang.  Learning particle models of swarming from data with Gaussian processes.  arXiv preprint arXiv:2106.02735, 2022.
  * [15] Jinchao Feng, Charles Kulick, and Sui Tang.  Data-driven model selections of second-order particle dynamics via integrating gaussian processes with low-dimensional interacting structure.  preprint, 2023.
  * [16] Baoli Hao, Ming Zhong, and Kevin O’Keeffe.  Mixed sign interactions in the 1d swarmalator model, 2023.
  * [17] R. Shu and E. Tadmor.  Flocking hydrodynamics with external potentials.  Archive for Rational Mechanics and Analysis, 238:347 – 381, 2020\.
  * [18] R. Shu and E. Tadmor.  Anticipation breeds alignment.  Archive for Rational Mechanics and Analysis, 240:203 – 241, 2021\.
  * [19] Audrey Filella, Francois Nadal, Clement Sire, Eva Kanso, and Christophe Eloy.  Model of collective fish behavior with hydrodynamic interactions.  Phys. Rev. Lett., 120:198101, 2018.
  * [20] K. P. O’Keeffe, H. Hong, and S. H. Strogatz.  Oscillators that sync and swarm.  Nat Commun, 8:1504, 2017.
  * [21] Trenton Gerew and Ming Zhong.  Concurrent emergence of clustering, flocking and synchronization in systems of interacting agents, 2023.
  * [22] James M Greene, Eitan Tadmor, and Ming Zhong.  The emergence of lines of hierarchy in collective motion of biological systems.  Physical Biology, 20(5):055001, jun 2023.
  * [23] Daniel A. Messenger and David M. Bortz.  Weak sindy for partial differential equations.  Journal of Computational Physics, 443:110525, 2021.
  * [24] M.P. Martin.  Multiplicatively Perturbed Least Squares for Dimension Reduction.  PhD thesis, Johns Hopkins University, 2021.
  * [25] Markus Heinonen, Cagatay Yildiz, Henrik Mannerström, Jukka Intosalmi, and Harri Lähdesmäki.  Learning unknown ode models with Gaussian processes.  In International Conference on Machine Learning, pages 1959–1968. PMLR, 2018.
  * [26] Cedric Archambeau, Dan Cornford, Manfred Opper, and John Shawe-Taylor.  Gaussian process approximations of stochastic differential equations.  In Gaussian Processes in Practice, pages 1–16. PMLR, 2007.
  * [27] Cagatay Yildiz, Markus Heinonen, Jukka Intosalmi, Henrik Mannerstrom, and Harri Lahdesmaki.  Learning stochastic differential equations with Gaussian processes without gradient matching.  In 2018 IEEE 28th International Workshop on Machine Learning for Signal Processing (MLSP), pages 1–6. IEEE, 2018.
  * [28] Zheng Zhao, Filip Tronarp, Roland Hostettler, and Simo Särkkä.  State-space Gaussian process for drift estimation in stochastic differential equations.  In ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 5295–5299. IEEE, 2020.
  * [29] Maziar Raissi, Paris Perdikaris, and George Em Karniadakis.  Machine learning of linear differential equations using Gaussian processes.  Journal of Computational Physics, 348:683–693, 2017.
  * [30] Jiuhai Chen, Lulu Kang, and Guang Lin.  Gaussian process assisted active learning of physical laws.  Technometrics, pages 1–14, 2020.
  * [31] Hongqiao Wang and Xiang Zhou.  Explicit estimation of derivatives from data and differential equations by Gaussian process regression.  International Journal for Uncertainty Quantification, 11(4), 2021\.
  * [32] Yifan Chen, Bamdad Hosseini, Houman Owhadi, and Andrew M Stuart.  Solving and learning nonlinear pdes with gaussian processes.  Journal of Computational Physics, 447:110668, 2021.
  * [33] Seungjoon Lee, Mahdi Kooshkbaghi, Konstantinos Spiliotis, Constantinos I Siettos, and Ioannis G Kevrekidis.  Coarse-scale PDEs from fine-scale observations via machine learning.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 30(1):013141, 2020.
  * [34] J-L Akian, Luc Bonnet, Houman Owhadi, and Éric Savin.  Learning “best” kernels from data in Gaussian process regression. With application to aerodynamics.  Journal of Computational Physics, 470:111595, 2022.
  * [35] Matthieu Darcy, Boumediene Hamzi, Jouni Susiluoto, Amy Braverman, and Houman Owhadi.  Learning dynamical systems from data: a simple cross-validation perspective, part ii: nonparametric kernel flows.  preprint, 2021.
  * [36] Christopher KI Williams and Carl Edward Rasmussen.  Gaussian processes for machine learning, volume 2.  MIT press Cambridge, MA, 2006.
  * [37] Ruiwen Shu and Eitan Tadmor.  Flocking hydrodynamics with external potentials.  Archive for Rational Mechanics and Analysis, 238:347–381, 2020\.
  * [38] Houman Owhadi and Clint Scovel.  Operator-adapted wavelets, fast solvers, and numerical homogenization: from a game theoretic approach to numerical approximation and algorithm design, volume 35.  Cambridge University Press, 2019.