remarkRemark assumptionAssumption hypothesisHypothesis claimClaim Data-Driven
Model Selection for Second-Order Particle Based DynamicsJinchao Feng, Charles
Kulick, Sui Tang

#  Data-Driven Model Selections of Second-Order Particle Dynamics via
Integrating Gaussian Processes with Low-Dimensional Interacting Structures
††thanks: This work was partially supported by NSF DMS-2111303.

Jinchao Feng School of Sciences, Great Bay University, Dongguan, Guangdong,
China<EMAIL_ADDRESS>Charles Kulick Departments of Mathematics,
University of California, Santa Barbara, Isla Vista, CA ().
<EMAIL_ADDRESS>Sui Tang Departments of Mathematics, University of
California, Santa Barbara, Isla Vista, CA<EMAIL_ADDRESS>
###### Abstract

In this paper, we focus on the data-driven discovery of a general second-order
particle-based model that contains many state-of-the-art models for modeling
the aggregation and collective behavior of interacting agents of similar size
and body type. This model takes the form of a high-dimensional system of
ordinary differential equations parameterized by two interaction kernels that
appraise the alignment of positions and velocities. We propose a Gaussian
Process-based approach to this problem, where the unknown model parameters are
marginalized by using two independent Gaussian Process (GP) priors on latent
interaction kernels constrained to dynamics and observational data. This
results in a nonparametric model for interacting dynamical systems that
accounts for uncertainty quantification. We also develop acceleration
techniques to improve scalability. Moreover, we perform a theoretical analysis
to interpret the methodology and investigate the conditions under which the
kernels can be recovered. We demonstrate the effectiveness of the proposed
approach on various prototype systems, including the selection of the order of
the systems and the types of interactions. In particular, we present
applications to modeling two real-world fish motion datasets that display
flocking and milling patterns up to 248 dimensions. Despite the use of small
data sets, the GP-based approach learns an effective representation of the
nonlinear dynamics in these spaces and outperforms competitor methods.

###### keywords:

Particle-based system, data-driven methods, Gaussian process, kernel ridge
regression, inverse problems, randomized numerical linear algebra

##  1 Introduction

Interacting particle/agent systems are a broad spectrum of complex systems
with multiple components interacting with each other and co-evolving with
time. Individual interactions yield a wide variety of collective behaviors at
different scales and levels of complexity such as clustering, alignment,
swarming, synchronization, or dancing equilibrium. There are numerous real-
world examples of such systems, including the orbits of planets, motion of
self-propelled particles, flocking of birds, schooling of fish, aggregation of
cells, consensus of opinions, and synchronization of oscillators over
networks. Understanding the link between individual interactions and global-
scale collective behaviors is one of the most fundamental problems in various
disciplines.

Modeling interacting agents by differential equations has played a crucial
role in exploring the emergence of collective behaviors from individual
interactions. However, such systems are often high-dimensional and exhibit
many possible dynamical couplings of components that contribute to the
dynamics, making them challenging to study [1, 2, 3, 4]. Despite these
challenges, recent work has made impressive progress in developing a general
physical model derived from Newton’s second law that can capture a wide range
of collective behaviors [5, 6, 7, 8, 9]. This model describes a system of $N$
agents interacting according to a set of ODEs, where each agent’s motion is
influenced by self-propulsion, friction, and interactions with other agents,
represented by energy and alignment-based radial interaction kernels: for
$i=1,\cdots,N$

(1)
$m_{i}\ddot{\bm{x}}_{i}=F_{i}({\bm{x}}_{i},\dot{\bm{x}}_{i},\bm{\alpha}_{i})+\sum_{i^{\prime}=1}^{N}\frac{1}{N}\Big{[}\phi^{E}(||{\bm{x}}_{i^{\prime}}-{\bm{x}}_{i}||)({\bm{x}}_{i^{\prime}}-{\bm{x}}_{i})+{\phi^{A}}(||{{\bm{x}}_{i^{\prime}}-{\bm{x}}_{i}}||)(\dot{\bm{x}}_{i^{\prime}}-\dot{\bm{x}}_{i})\Big{]},$

where $m_{i}\geq 0$ is the mass of the agent $i$;
$\ddot{\bm{x}}_{i}\in\mathbb{R}^{d}$ is the acceleration,
$\dot{\bm{x}}_{i}\in\mathbb{R}^{d}$ is the velocity, and
${\bm{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 with scalar
parameters $\bm{\alpha}_{i}$ describing their strength;
$||{\bm{x}}_{j}-{\bm{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. We summarize the relevant notations in Table
1.

Particular examples of (1) include the first-order systems ($m_{i}\equiv
0,{\phi^{A}}\equiv 0$) that model clustering and aggregation of agents with
application to opinion dynamics [10], the second-order Cucker-Smale model
(${\phi^{E}}\equiv 0$) [11] for the flocking behavior of animals and robots,
the second-order self-propelling particle model (${\phi^{A}}\equiv 0$) that is
shown to reproduce (double) milling, ring, escaping or swarming behaviors of
biological motors [12], and the anticipation dynamics [13]
(${\phi^{A}},{\phi^{E}}\neq 0$) that describes the velocity alignment and
spatial concentration of animal groups. For simplicity of description, we
assume the masses of all agents are the same and equal to $\bm{m}$. We write
the second-order model (1) in a compact form:

(2)
$\bm{m}\bm{Z}(t)=F_{\bm{\alpha}}(\bm{Y}(t))+\mathbf{f}_{\bm{\phi}}(\bm{Y}(t))$

where $\bm{Y}(t):=\begin{bmatrix}\bm{X}(t)\\\
\bm{V}(t)\end{bmatrix}\in\mathbb{R}^{2dN}$ represents the state variable for
the system, $\bm{Z}(t)=\dot{\bm{V}}(t)=\ddot{\bm{X}}(t)$, and
$\mathbf{f}_{\bm{\phi}}(\bm{Y}(t))=\mathbf{f}_{{\phi^{E}},{\phi^{A}}}(\bm{Y}(t))$
represents the sum of energy and alignment-based interactions as in (1).

Table 1: Notations for second-order systems Variable | Definition
---|---
$N$ | number of agents
$m_{i}$ |  mass of agent $i$
${\bm{x}}_{i}(t)\in{}^{d}$ |  position vector of agent $i$ at time $t$
$\dot{{\bm{x}}}_{i}(t)\in{}^{d}$ |  velocity vector of agent $i$ at time $t$
$\ddot{{\bm{x}}}_{i}(t)\in{}^{d}$ |  acceleration vector of agent $i$ at time $t$
$F$ | non-collective force
$\bm{\alpha}$ |  parameters of $F$
${\phi^{E}},{\phi^{A}}$ | energy and alignment-based interaction kernels respectively
$\|\cdot\|$ |  Euclidean norm in $\mathbb{R}^{d}$

###  1.1 Data-driven model selection problem

Recent advancements in data information technology, such as digital imaging,
high-resolution lightweight GPS devices, and particle tracking methods, have
allowed for the gathering of high-resolution trajectory data of individual
particles in various applications. However, a significant issue that remains
scarcely addressed is how to select models that match the observational data.
For example, while there are many theoretical models known to reproduce
flocking patterns, it is challenging to determine which one generates the
pattern observed in the data. Previous theoretical and numerical studies
cannot address this problem, as predetermined governing equations are needed,
and the aim is often to reproduce qualitative rather than quantitative
dynamics.

To address this issue, we consider the data-driven model selection problem,
aiming to select possible models from a general form to match the
observational data. For instance, given the motion data of a school of fish,
we aim to determine whether to use first-order or second-order models and
which types of interactions, such as alignment versus energy-based or both,
contribute to collective patterns. These are challenging questions that
practitioners typically address based on their expertise in the field. In this
paper, we seek to develop data-driven methods to automate this step by
considering a general model that incorporates many classical models as special
cases.

Mathematically, we formulate the problem as follows. Given approximate
observations of multiple trajectory data
$\mathcal{D}_{M,L}:=\\{\bm{Y}^{(m)}(t_{l}),\bm{Z}^{(m)}(t_{l})\\}_{m,l=1}^{M,L}$,
where the observation time instances are denoted by $0=t_{1}<\dots<t_{L}=T$
and $m$ denotes the trial number of experiments starting from different
initial conditions, the goal is to infer the interaction kernels
${\bm{\phi}}={{\phi^{E}},{\phi^{A}}}$ as well as the unknown scalar parameters
$\bm{\alpha}$ and possibly $\bm{m}$ from the trajectory data
$\mathcal{D}_{M,L}$. Subsequently, we use the learned governing equations to
make predictions about future events or simulate new datasets.

###  1.2 Scalable Model Selection by Gaussian processes

The field of data-driven model selection faces two primary practical
challenges. Firstly, there is often limited information available on the
parametric forms of interaction kernels, making it difficult to select a
suitable approximation dictionary. Secondly, datasets may be scarce and noisy.
Gaussian process (GP) based approaches in machine learning offer a solution to
these challenges, as they are known for their ability to learn a rich class of
nonlinear functions without making assumptions about their parametric form and
for quantifying the associated uncertainty. However, the challenge of
scalability to large-scale problems remains a significant hurdle for specific
applications.

This paper proposes a novel approach to address these challenges by leveraging
the inference power of Gaussian processes and developing efficient techniques
to improve scalability. Computationally,

  * •

We propose a novel method by modeling interaction kernels as two independent
Gaussian processes to learn (1) from data with uncertainty quantification. We
investigate whether types of interaction kernels and order information (first
versus second order) of the system can be learned from scarce noisy data. We
conduct intensive numerical experiments on various prototypical systems
exhibiting clustering, milling, and flocking behaviors that demonstrate the
effectiveness.

  * •

We propose effective acceleration techniques based on the recent progress from
randomized numerical linear algebra.

  * •

Our method is applied to modeling two real-world fish motion sets that display
flocking and milling patterns up to 248 dimensions and outperforms competitor
methods that use SINDy and feed-forward neural networks.

Theoretically,

  * •

We derive a Representer theorem that connects the GP-based estimators with the
kernel ridge regression estimators, shedding light on the role of the
hyperparameters in learning. It also provides a basis representation for the
estimators of interaction kernels, which enables efficient trajectory
prediction using learned models over larger time intervals.

  * •

We study the well-posedness of the inverse problem for learning interaction
kernels in a statistical setting.

###  1.3 Relevant works

Integrating machine learning techniques into the data-driven discovery of
dynamical systems (see e.g. [14, 5, 15, 16, 6, 17, 18, 7, 8, 9, 19, 20, 21])
has become a hot topic in scientific machine learning, as it provides powerful
models to represent the complex functional data. In terms of parametric
methods, one can refer to [22] (and references therein) for the most recent
survey on deep learning techniques and [15, 23, 24, 25] for sparse regression
techniques.

Gaussian process regression (GPR) is a non-parametric Bayesian machine
learning technique for supervised learning with a built-in quantification of
uncertainty framework. As such, GPs have been applied to learn ODEs, SDEs, and
PDEs [26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36] and lead to more accurate
and robust models of dynamical systems. Because of the distinctive nature of
dynamical data, it necessitates novel ideas and nontrivial efforts tailored to
particular types of dynamical systems and data regimes. We model the latent
interaction kernels as GPs and imbue them with the structure of our governing
equations (translation and rotational invariance). This makes our work
distinguishable from most works, which model state variables as GPs.

In the context of interaction kernel learning in interacting particle systems,
least square estimators derived from the maximum likelihood method are the
most frequently used, where a challenge lies in the selection of the basis to
represent the interaction kernels. One can refer to [1, 2, 3, 4] for the usage
of a piecewise polynomial basis. The random feature method together with
sparse regression techniques is recently proposed in [37]. One can also refer
to the recent methodology development of interaction kernel/potential learning
in mean-field systems such as [38, 39, 40, 41, 42].

In particular, [3, 4] considered learning theory for heterogeneous systems and
showed that learning multiple interaction kernels simultaneously is
challenging and regularization is necessary.

For scarce noisy data, our GP method leverages the underlying statistical
inference power to select the best basis to represent the observational
dynamics and provides effective regularization. It yields accurate recovery of
the governing equation beyond learning only interaction kernels. It is well-
known that the non-collective force also plays an important role in
determining the collective behaviors. The governing equation recovery makes
our method more practical than previous work that only focuses on interaction
kernels. Further, we analyze the well-posedness of the inverse problems, which
complements the missing analysis in [4]. This work is an extension of our
recent work [43] on a single kernel case where we assumed ${\phi^{A}}\equiv 0$
and the focus was the theoretical framework for error analysis. Here, we
consider a more generalized model involving two types of kernels and consider
the model selection problems. The focus is shifted to the computational
aspects concerning scalability and uncertainty quantification, and real data
applications.

###  1.4 Notation and preliminaries

##### Notation

Let $\rho$ be a Borel positive measure on $D$ dimensional Eucliean space
$\mathbb{R}^{D}$. We use $L^{2}(\mathbb{R}^{D};\rho;\mathbb{R}^{n})$ to denote
the set of $L^{2}(\rho)$-integrable vector-valued functions that map
$\mathbb{R}^{D}$ to $\mathbb{R}^{n}$. For a function $\bm{f}\in
L^{2}(\mathbb{R}^{D};\rho;\mathbb{R}^{n})$, and a vector
$\bm{X}=[{\bm{x}}_{1}^{\top},\cdots,{\bm{x}}_{m}^{\top}]^{T}\in\mathbb{R}^{mD}$
with ${\bm{x}}_{i}\in\mathbb{R}^{D}$, we use the notation $\bm{f}(\bm{X})$ to
represent the image of the vector under the function of $\bm{f}$
componentwisely, namely,
$\bm{f}(\bm{X})=[\bm{f}({\bm{x}}_{1})^{\top},\cdots,\bm{f}({\bm{x}}_{m})^{\top}]^{\top}\in\mathbb{R}^{mn}$.
Let $\mathcal{S}_{1}$ be a measurable subset of $\mathbb{R}^{m}$, the
restriction of the measure $\rho$ on $\mathcal{S}_{1}$, denote by
$\rho\mathbin{\vrule height=5.16663pt,depth=0.0pt,width=0.55974pt\vrule
height=0.55974pt,depth=0.0pt,width=3.87495pt}\mathcal{S}_{1}$, is defined as
$\rho\mathbin{\vrule height=5.16663pt,depth=0.0pt,width=0.55974pt\vrule
height=0.55974pt,depth=0.0pt,width=3.87495pt}\mathcal{S}_{1}(\mathcal{S}_{2})=\rho(\mathcal{S}_{1}\cap\mathcal{S}_{2})$
for any measurable subset $\mathcal{S}_{2}$ of $\mathbb{R}^{D}$. We used
$\mathcal{N}(0,I_{d\times d})$ to denote the standard multivariate Gaussian
distribution in $\mathbb{R}^{d}$.

##### Preliminaries on operator algebras

Let $\mathcal{H}_{1},\mathcal{H}_{2}$ be Hilbert spaces. We use
$\langle\cdot,\cdot\rangle_{\mathcal{H}_{1}}$ to denote the inner product over
$\mathcal{H}_{1}$, and still use $\langle\cdot,\cdot\rangle$ to denote the
inner product on the Euclidean space. We denote by
$\mathcal{B}(\mathcal{H}_{1},\mathcal{H}_{2})$ the set of bounded linear
operators mapping $\mathcal{H}_{1}$ to $\mathcal{H}_{2}$. Let
$A\in\mathcal{B}(\mathcal{H}_{1},\mathcal{H}_{2})$, we use $\mathrm{Im}(A)$ to
denote its range and $\|A\|$ to denote its operator norm. $A$ is a compact
operator if $A$ maps bounded subsets of $\mathcal{H}_{1}$ to relatively
compact subsets of $\mathcal{H}_{2}$ (subsets with compact closure in
$\mathcal{H}_{2}$). We use $A^{*}:\mathcal{H}_{2}\rightarrow\mathcal{H}_{1}$
to denote the adjoint operator of $A$, that is, $\forall f\in\mathcal{H}_{1}$,
$g\in\mathcal{H}_{2}$, $\langle Af,g\rangle_{\mathcal{H}_{2}}=\langle
f,A^{*}g\rangle_{\mathcal{H}_{1}}$.

For $d,N,M,L\in\mathbb{N}^{+}$, let
$\bm{w}=(\bm{w}_{m,l,i})_{m,l,i=1}^{M,L,N},\bm{z}=(\bm{z}_{m,l,i})_{m,l,i=1}^{M,L,N}\in\mathbb{R}^{dNML}$
with $\bm{w}_{m,l,i},\bm{z}_{m,l,i}\in\mathbb{R}^{d}$, we define

(3)
$\langle\bm{w},\bm{z}\rangle=\frac{1}{MLN}\sum_{m,l,i=1}^{M,L,N}\langle\bm{w}_{m,l,i},\bm{z}_{m,l,i}\rangle,$

where $\langle\bm{w}_{m,l,i},\bm{z}_{m,l,i}\rangle$ is the canonical inner
product on $\mathbb{R}^{d}$.

Then for vectors $\bm{y}\in\mathbb{R}^{mdN}$ and functions
$\bm{g}:\mathbb{R}^{mdN}\rightarrow\mathbb{R}^{dn}$ for some
$m,n\in\mathbb{N}^{+}$, and let $\rho$ be a measure at a Borel subset
$\mathcal{Y}$ of $\mathbb{R}^{mdN}$, we have the norm:

(4)
$\lVert\bm{g}(\bm{y})\rVert_{L^{2}(\rho)}^{2}=\frac{1}{n}\int_{\mathcal{Y}}\sum_{i=1}^{n}\lVert\bm{g}_{i}(\bm{y})\rVert^{2}\rho(d\bm{y})$

where $\bm{g}(x)$ is componentwise denoted by
$\bm{g}_{i}(\bm{y}):\mathbb{R}^{mdN}\rightarrow\mathbb{R}^{d}$.

For two Borel positive measures $\rho_{1},\rho_{2}$ defined on
$\mathbb{R}^{D}$, $\rho_{1}$ is said to be absolutely continuous with respect
to $\rho_{2}$, $\rho_{1}\ll\rho_{2}$, if $\rho_{1}(\mathcal{S})=0$ for every
set $\rho_{2}(\mathcal{S})=0$, $\mathcal{S}\subset\mathbb{R}^{D}$. $\rho_{1}$
and $\rho_{2}$ are called equivalent iff $\rho_{1}\ll\rho_{2}$ and
$\rho_{2}\ll\rho_{1}$. The product measure $\rho_{1}\times\rho_{2}$ is defined
to be a measure on $\mathbb{R}^{2D}$ satisfying the property
$(\rho_{1}\times\rho_{2})(\mathcal{S}_{1}\times\mathcal{S}_{2})=\rho_{1}(\mathcal{S}_{1})\rho_{2}(\mathcal{S}_{2})$
for all subsets $\mathcal{S}_{i}\subset\mathbb{R}^{D}$, $i=1,2$.

##### Preliminaries on GPs (Gaussian Processes) Prior

We say $\phi\sim\mathcal{GP}(u,K)$ to denote our prior on $\phi$. In
particular, this means that for any $r\in\mathbb{R}$, the random variable
$\phi(r)$ is Gaussian: $\phi(r)\sim\mathcal{N}(u(r),K(r,r))$, where
$\mathcal{N}$ denotes the normal or multivariable normal distributions,
$u:\mathbb{R}\to\mathbb{R}$ is the mean function, and
$K:\mathbb{R}\times\mathbb{R}\to\mathbb{R}$ is the covariance function.
Similarly, for any $(r,r^{\prime})\in\mathbb{R}^{2}$, the joint distribution
of $\begin{bmatrix}\phi(r)\\\ \phi(r^{\prime})\end{bmatrix}$ is multivariate
Gaussian: $\begin{bmatrix}\phi(r)\\\
\phi(r^{\prime})\end{bmatrix}\sim\mathcal{N}\left(\begin{bmatrix}(u(r)\\\
(u(r^{\prime})\end{bmatrix},\begin{bmatrix}K(r,r)&K(r,r^{\prime})\\\
K(r^{\prime},r)&K(r^{\prime},r^{\prime})\end{bmatrix}\right)$. This extends in
a natural way to any finite set $(r_{1},\dots,r_{N})\in\mathbb{R}^{N}$.

##### Preliminaries on RKHSs

Let $\mathcal{D}$ be a compact subset of $\mathbb{R}^{D}$. We say that
$K:\mathcal{D}\times\mathcal{D}\rightarrow\mathbb{R}$ is a Mercer kernel if it
is continuous, symmetric, and positive semidefinite, i.e., for any finite set
of distinct points $\\{x_{1},\cdots,x_{M}\\}\subset\mathcal{D},$ the matrix
$(K(x_{i},x_{j}))_{i,j=1}^{M}$ is positive semidefinite. For
$x\in\mathbb{R}^{D}$, $K_{x}$ is a function defined on $\mathcal{D}$ such that
$K_{x}(y)=K(x,y)$, $y\in\mathcal{D}$. The Moore–Aronszajn theorem proves that
there is an RKHS $\mathcal{H}_{{K}}$ associated with the kernel $K$, which is
defined to be the closure of the linear span of the set of functions
$\\{K_{x}:x\in\mathcal{D}\\}$ with respect to the inner product
$\langle\cdot,\cdot\rangle_{\mathcal{H}_{{K}}}$ satisfying $\langle
K_{x},K_{y}\rangle_{\mathcal{H}_{{K}}}=K(x,y)$. Let ${K}$ be a Mercer kernel
that is defined on $[0,R]\times[0,R]$ and use $\mathcal{H}_{{K}}$ to denote
the RKHS associated with ${K}$. For two RKHS
$\mathcal{H}_{K_{1}},\mathcal{H}_{K_{2}}$, with
$K_{1},K_{2}:\mathcal{D}\times\mathcal{D}\rightarrow\mathbb{R}$, the product
RKHS $\mathcal{H}_{K_{1}}\times\mathcal{H}_{K_{2}}$ is defined to be the
closure of the linear span of the set of functions
$\\{(K_{1,x_{1}},K_{2,x_{2}}):x_{1},x_{2}\in\mathcal{D}\\}$ with respect to
the inner product
$\langle\cdot,\cdot\rangle_{\mathcal{H}_{K_{1}}\times\mathcal{H}_{K_{2}}}$
satisfying
$\langle(K_{1,x_{1}},K_{2,x_{2}}),(K_{1,y_{1}},K_{2,y_{2}})\rangle_{\mathcal{H}_{K_{1}}\times\mathcal{H}_{K_{2}}}=K_{1}(x_{1},y_{1})+K_{2}(x_{2},y_{2})$.

##  2 Methodology

In this section, we propose a learning approach based on GPs for the model
selection problem.

###  2.1 Two independent Gaussian process priors

We start by modeling the interaction kernel functions ${\phi^{E}}$ and
${\phi^{A}}$ with the priors as two independent Gaussian processes

(5)
${\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})),$

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

###  2.2 Training of hyperparameters via maximum likelihood estimation

In real-world modeling, it is possible that some other parameters such as
$\bm{\alpha}$ and noise level in the data are unknown. In this section, we
detail how to perform the estimation of these physical parameters in the
governing equation via a data-driven hyperparameter tuning process induced by
the Gaussian process. This flexible training procedure distinguishes the
Gaussian process from other kernel-based methods [44, 45, 46] and
regularization-based approaches [47, 48, 49].

We organize the training data into the vector format
$\mathbb{Y}=[\bm{Y}^{(1,1)},\dots,\bm{Y}^{(M,L)}]^{T}\in\mathbb{R}^{dNML}$,
and $\mathbb{Z}=[\bm{Z}^{(1,1)},\dots,\bm{Z}^{(M,L)}]^{T}\in\mathbb{R}^{dNML}$
where

(6)
$\bm{Y}^{(m,l)}=\bm{Y}^{(m)}(t_{l}),\quad\bm{Z}^{(m,l)}=\bm{Z}^{(m)}(t_{l}).$

To model the noise, we assume
$\mathbb{Z}=[\bm{Z}^{(1,1)}_{\sigma^{2}},\dots,\bm{Z}^{(M,L)}_{\sigma^{2}}]^{T}\in\mathbb{R}^{dNML}$
where

(7)
$\bm{m}\bm{Z}^{(m,l)}_{\sigma^{2}}=F_{\bm{\alpha}}(\bm{Y}^{(m,l)})+\mathbf{f}_{{\bm{\phi}}}(\bm{Y}^{(m,l)})+\epsilon^{(m,l)},$

with i.i.d (independent and identically distributed) noise
$\epsilon^{(m,l)}\sim\mathcal{N}(0,\sigma^{2}I_{dN})$ that is also independent
of the Gaussian processes. Later, we will show the role of $\sigma$ in the
prediction step is equivalent to the role of the regularization constant in a
Tikhonov regularization problem.

Therefore, based on the properties of Gaussian processes, with the priors of
${\phi^{E}}$, ${\phi^{A}}$, we have

(8)
$\bm{m}\mathbb{Z}\sim\mathcal{N}(F_{\bm{\alpha}}(\mathbb{Y}),K_{\mathbf{f}_{{\bm{\phi}}}}(\mathbb{Y},\mathbb{Y};\theta)+\sigma^{2}I_{dNML}),$

with the mean vector
$F_{\bm{\alpha}}(\mathbb{Y})=\mathrm{Vec}(\\{F_{\bm{\alpha}}(\bm{Y}^{(m,l)})\\}_{m,l=1}^{M,L})\in\mathbb{R}^{dNML}$,
and
$K_{\mathbf{f}_{{\bm{\phi}}}}(\mathbb{Y},\mathbb{Y};\theta)\in\mathbb{R}^{dNML\times
dNML}$ is the covariance matrix between $\mathbf{f}_{{\bm{\phi}}}(\mathbb{Y})$
and $\mathbf{f}_{{\bm{\phi}}}(\mathbb{Y})$, which can be computed elementwise
based on the covariance functions $K_{\theta^{E}}$, $K_{\theta^{A}}$, see
Appendix Appendix A for detailed formulas.

Thus, for the hyperparameters $\bm{\alpha}$, $\bm{\theta}$, and $\sigma$, we
can train by maximizing the probability of the observational data, which is
equivalent to minimizing the negative log marginal likelihood (NLML) (see
Chapter 4 in [50])

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

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. To solve for the hyperparameters
$(\bm{\alpha},\bm{\theta},\sigma)$, we can apply the conjugate gradient (CG)
optimization (see Chapter 5 in [50]) to minimize the negative log marginal
likelihood. More details are shown in Appendix Appendix A.

Table 2: Notations for covariances

Variable | Definition
---|---
$K_{\bm{\theta}}(\cdot,\cdot)$ | covariance kernel function with parameters $\bm{\theta}$
$K_{{\theta^{E}}}(\cdot,\cdot),K_{{\theta^{A}}}(\cdot,\cdot)$ | covariance kernels for modeling ${\phi^{E}}$, ${\phi^{A}}$
$K_{\mathbf{f}_{\phi}}(\cdot,\cdot)$ |  covariance matrix between $\mathbf{f}_{\phi}(\cdot)$ and $\mathbf{f}_{\phi}(\cdot)$
$K_{\mathbf{f}_{\phi},{\phi^{E}}}(\cdot,\cdot):=K_{{\phi^{E}},\mathbf{f}_{\phi}}(\cdot,\cdot)^{T}$ | covariance matrix between $\mathbf{f}_{\bm{\phi}}(\cdot)$ and ${\phi^{E}}(\cdot)$
$K_{\mathbf{f}_{\phi},{\phi^{A}}}(\cdot,\cdot):=K_{{\phi^{A}},\mathbf{f}_{\phi}}(\cdot,\cdot)^{T}$ | covariance matrix between $\mathbf{f}_{\bm{\phi}}(\cdot)$ and ${\phi^{A}}(\cdot)$

####  2.2.1 Parameter $\bf{m}$\- Model selection of the order for the
dynamical system

When modeling real-world dynamics, sometimes we are not sure whether to use
first-order or second-order systems. From the parameter estimation
perspective, it is equivalent to determining if the mass of particles is equal
to zero. We can train $\bm{m}$ via minimizing (9). If the estimation of
$\bm{m}$ is close to zero, we can consider identifying the dynamics as a
first-order system. One can refer to Section 4.3.

###  2.3 Learning interaction kernels

Table 3: Notations for second-order systems Variable | Definition
---|---
$\bm{X}\in{}^{dN}$ |  vectorization of position vectors $({\bm{x}}_{i})_{i=1}^{N}$
$\bm{V}\in{}^{dN}$ |  vectorization of velocity vectors $({\bm{v}}_{i})_{i=1}^{N}=(\dot{{\bm{x}}}_{i})_{i=1}^{N}$
$\bm{Y}\in{}^{2dN}$ | $\bm{Y}=(\bm{X},\bm{V})^{T}$
$\bm{Z}\in{}^{dN}$ |  vectorization of $(\ddot{{\bm{x}}}_{i})_{i=1}^{N}$
$\bm{r}^{\bm{x}}_{ij},\bm{r}^{{\bm{x}}^{\prime}}_{ij}\in{}^{d}$ |  $\bm{X}(t)_{j}-\bm{X}(t)_{i}$, $\bm{X}(t^{\prime})_{j}-\bm{X}(t^{\prime})_{i}$
$\bm{r}^{\bm{v}}_{ij},\bm{r}^{{\bm{v}}^{\prime}}_{ij}\in{}^{d}$ |  $\bm{V}(t)_{j}-\bm{V}(t)_{i}$, $\bm{V}(t^{\prime})_{j}-\bm{V}(t^{\prime})_{i}$
$r^{\bm{x}}_{ij},r^{{\bm{x}}^{\prime}}_{ij}\in{}^{+}$ | $r^{\bm{x}}_{ik}=\|\bm{r}^{\bm{x}}_{ik}\|,r^{{\bm{x}}^{\prime}}_{ij}=\|\bm{r}^{{\bm{x}}^{\prime}}_{ij}\|$
$\mathbf{f}_{{\phi^{E}}},\mathbf{f}_{{\phi^{A}}}$ |  energy and alignment-based interaction force field
$\mathbf{f}_{{\bm{\phi}}}$ |  interaction force field with ${\bm{\phi}}=({\phi^{E}},{\phi^{A}})$

We plug the estimators of hyperparameters obtained in Section 2.2 into the
system and assume they are known. In this subsection, we show how to learn
interaction kernels. For any $r^{\ast}\in\mathbb{R}$ and the corresponding
values of the kernel functions, $\phi^{\mathrm{type}}(r^{*})$,
$\mathrm{type}=E\text{ or }A$, since we have

(10) $\begin{bmatrix}\bm{m}\mathbb{Z}-F_{\mathbf{\alpha}}(\mathbb{Y})\\\
\phi^{\mathrm{type}}(r^{\ast})\end{bmatrix}\sim\mathcal{N}\left(0,\begin{bmatrix}K_{\mathbf{f}_{{\bm{\phi}}}}(\mathbb{Y},\mathbb{Y})+\sigma^{2}I_{dNML}&K_{\mathbf{f}_{\bm{\phi}},\phi^{\mathrm{type}}}(\mathbb{Y},r^{\ast})\\\
K_{\phi^{\mathrm{type}},\mathbf{f}_{\bm{\phi}}}(r^{\ast},\mathbb{Y})&K_{\theta^{\mathrm{type}}}(r^{\ast},r^{\ast})\end{bmatrix}\right),$

where
$K_{\mathbf{f}_{\bm{\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}_{{\bm{\phi}}}(\mathbb{Y})$
and $\phi^{\mathrm{type}}(r^{*})$. Conditioning on
$\mathbf{f}_{{\bm{\phi}}}(\mathbb{Y})$, we obtain the posterior/predictive
distribution for the kernel function value at $r^{\ast}$,
$\phi^{\mathrm{type}}(r^{\ast})$ (see Lemma D.6 in Appendix for detailed
derivation), i.e.

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

where

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

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.

Input: $(\mathbb{Y},\mathbb{Z})$ (training data), $r^{\ast}$ (test point),
$K_{\bm{\theta}}$ (covariance functions), $\mathbf{f}_{{\bm{\phi}}}$
(interaction function), $F_{\bm{\alpha}}$ (force function)

1:
$(\hat{\bm{\alpha}},\hat{\bm{\theta}},\hat{\sigma}^{2})=\underset{\bm{\alpha},\bm{\theta},\sigma^{2}}{\operatorname{arg}\operatorname{min}}\;-\log
p(\bm{m}\mathbb{Z}|\mathbb{Y},\bm{\alpha},\bm{\theta},\sigma^{2})${solve for
parameters by minimizing NLML (9) using CG}

2:
$L:=\textrm{cholesky}(K_{\mathbf{f}_{{\bm{\phi}}}}(\mathbb{Y},\mathbb{Y})+\hat{\sigma}^{2}I)$

3:
$\gamma:=L^{T}\backslash(L\backslash(\bm{m}\mathbb{Z}-F_{\bm{\hat{\alpha}}}(\mathbb{Y})))$

4:
$K_{E}^{\ast}:=K_{\mathbf{f}_{\bm{\phi}},{\phi^{E}}}(\mathbb{Y},r^{\ast})$$K_{A}^{\ast}:=K_{\mathbf{f}_{\bm{\phi}},{\phi^{A}}}(\mathbb{Y},r^{\ast})$
{compute covariances between $\mathbf{f}_{{\bm{\phi}}}(\mathbb{Y})$ and
${\phi^{E}}(r^{\ast})$, ${\phi^{A}}(r^{\ast})$}

5:
$\bar{\phi}^{E\ast}:=(K_{E}^{\ast})^{T}\gamma$$\bar{\phi}^{A\ast}:=(K_{A}^{\ast})^{T}\gamma$
{predictive mean Eq. 12}

6: ${\bm{v}}_{E}=L\backslash K_{E}^{\ast}$, ${\bm{v}}_{A}=L\backslash
K_{A}^{\ast}$

7:
$var(\phi^{E\ast}):=K_{\hat{\theta}^{E}}(r^{\ast},r^{\ast})-{\bm{v}}_{E}^{T}{\bm{v}}_{E}$$var(\phi^{A\ast}):=K_{\hat{\theta}^{A}}(r^{\ast},r^{\ast})-{\bm{v}}_{A}^{T}{\bm{v}}_{A}$
{predictive variance Eq. 13}

Output: $\bar{\phi}^{E\ast},\bar{\phi}^{A\ast}$ (mean),
$var(\phi^{E\ast}),var(\phi^{A\ast})$ (variance)

Algorithm 1 Learning kernels

###  2.4 Prediction of trajectories and its uncertainty quantification

We use the posterior mean estimators of ${\bm{\phi}}$ in trajectory prediction
by performing numerical simulations of the equations

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

We can also perform uncertainty quantification for the trajectory prediction
via the uncertainty band of $\bm{\hat{\phi}}$. We adopted a Monte Carlo
method, where we used $\bm{\phi}$ sampled from the posterior distribution in
each simulation. Then the predictions of the trajectories are given by the
mean of the trajectories’ samples and the uncertainty band of each trajectory
is given by the standard deviation, with the results of experiments shown in
Section 4. Another possible alternative is to use step-wise uncertainty
quantification based on the numerical integrator scheme such as the one-step
Euler method. In this case, it is easy to compute the variance of the solution
from the posterior distribution of ${\bm{\phi}}$, since the vector field
$\hat{\mathbf{f}}_{\bar{{\bm{\phi}}}}(\bm{Y}(t))$ is a linear combination of
$\hat{{\bm{\phi}}}$ by its definition, which suggests it also follows a
Gaussian distribution and the uncertainty band can be derived from its
covariance matrix.

###  2.5 Acceleration of the Computation

While the full GP methods described above yield extremely accurate predictions
in our empirical examples, a well-known limitation is the computational
complexity; calculating the log determinant of $K_{\mathbf{f}_{{\bm{\phi}}}}$
and inverting the kernel matrices in the maximum likelihood estimation and
prediction steps scales cubically with the matrix dimension, which is
$\mathcal{O}((NdML)^{3})$. Therefore, the naive approach can quickly become
infeasible for large-scale problems. Below, we describe our integrated
approach to the scalable estimation of hyperparameters in maximum likelihood
estimation and scalable kernel prediction.

####  2.5.1 Efficient Hyperparameter Optimization

There are many recent advancements in accelerating the hyperparameter learning
computations in the full GP methods for regression tasks. Our problem,
however, presents many numerical difficulties that dampen runtime gains from
traditional computational methods and must be addressed:

  * •

Lack of sparsity. Many classical acceleration techniques rely on the sparsity
of the kernel matrix $K_{\mathbf{f}_{{\bm{\phi}}}}$. As our kernel depends on
pairwise distance and our modeling is nonlocal, we do not have a sparse kernel
matrix in our formulation. Our method must be able to operate on dense
$K_{\mathbf{f}_{{\bm{\phi}}}}$.

  * •

Extreme ill-conditioning and higher accuracy requirements. The $L^{2}$
condition number of a matrix is the ratio of its maximum and minimum singular
values. When much larger than $1$, the condition number indicates that a
matrix is nearly singular, and thus accuracy-reducing errors in computation
will occur. For many problems, such as those addressed in Section 5, the
kernel matrix $K_{\mathbf{f}_{{\bm{\phi}}}}$ has observed $L^{2}$ condition
number above $10^{15}$. These extremely high condition numbers result in slow
and inaccurate computation when using traditional methods. Our problem is also
an inverse problem while learning our hyperparameters for
$K_{\mathbf{f}_{{\bm{\phi}}}}$ (see Section 3.1). This is very sensitive to
perturbations, especially as our optimization problem for the hyperparameters
is generally not convex. We must carefully balance the tradeoff between
computational time and accuracy.

We empirically observed the approximately low-rank structure of
$K_{\mathbf{f}_{{\bm{\phi}}}}$ in various examples. This motivated us to adapt
two main classes of algorithms in [51] for acceleration (see pseudocode and
additional details in Appendix Appendix B):

  * •

Preconditioned conjugate gradient (PCG) algorithm. The PCG algorithm allows us
to avoid explicit computation of the inverse matrix
$(K_{\mathbf{f}_{{\bm{\phi}}}}(\mathbb{Y},\mathbb{Y};\theta)+\sigma^{2}I)^{-1}$
in both MLE and prediction, as well as compute the coefficients needed in the
stochastic Lanczos quadrature below. Using preconditioners, a classical
numerical technique to lower condition numbers, is a necessity for variance
reduction. In addition, we must maintain a low error tolerance for PCG to
preserve our accuracy throughout learning. Finding effective preconditioners
that are suitable to the unique structure of our kernel matrices is a
challenge. We propose using the Random Gaussian Nystrom preconditioner [52] to
ensure favorable tradeoffs in running time and accuracy. In our practical
implementation, this preconditioner outperformed other low-rank approximation
preconditioners and has low construction and inversion costs, see Section 5.

  * •

Stochastic trace estimation for log determinant acceleration. We utilize the
identity:

$\log\det(K_{\mathbf{f}_{{\bm{\phi}}}}(\mathbb{Y},\mathbb{Y};\theta)+\sigma^{2}I)=\log\det(P)+\log\det(P^{-\frac{1}{2}}(K_{\mathbf{f}_{{\bm{\phi}}}}(\mathbb{Y},\mathbb{Y};\theta)+\sigma^{2}I)P^{-\frac{1}{2}})$

When $P$ is chosen to be a preconditioner, this identity can prove highly
useful. The Random Gaussian Nystrom preconditioner allows us to efficiently
compute $\log\det(P)$, and for the remainder, we use the recently developed
variance reduced Hutchinson’s Estimator [53] combined with stochastic Lanczos
quadrature [51].

##### Analysis of new computational complexity

PCG can reduce explicit inversion complexity from $\mathcal{O}((NdML)^{3})$ to
$\mathcal{O}(t(NdML)^{2})$, where $t$ is the number of iterations. The
stochastic Lanczos quadrature improves log determinant complexity from
$\mathcal{O}((NdML)^{3})$ to
$\mathcal{O}(t\ell(NdML)^{2})+\mathcal{O}(\log\det(P))$, where $t$ is both the
number of eigenvalues and the number of iterations, $\ell$ is the number of
runs of stochastic Lanczos, and $\mathcal{O}(\log\det(P))$ is the complexity
of computing the log determinant of the preconditioner. In practice, we chose
$\ell,t<<NdML$. This lowers the theoretical complexity of these steps to the
quadratic $\mathcal{O}(t\ell(NdML)^{2})$. For the Random Gaussian Nystrom
preconditioner $P$ with rank $r$, we have construction in
$\mathcal{O}(r^{2}(NdML)+r^{3})$ time, inversion in $\mathcal{O}(r^{3})$ time
and log determinant in $\mathcal{O}(r)$ time.

##  3 Theoretical analysis

In this section, we are concerned with two theoretical problems regarding
learning interaction kernels in the prediction step. The first one is to
understand the role of hyperparameters in the prediction step of the Gaussian
process, i.e., $\theta^{E},\theta^{A}$, and the Gaussian noise $\sigma$. The
second one is to study well-posedness as an inverse problem.

As in the prediction step, interaction kernels are the only unknown terms in
the equations. We make the following simplification on the form of equations
to avoid unnecessary technical hurdles:

(15) $\displaystyle\ddot{\bm{X}}(t)$
$\displaystyle=\mathbf{f}_{{\bm{\phi}}}(\bm{Y}(t))=\mathbf{f}_{{\phi^{E}}}(\bm{X}(t))+\mathbf{f}_{{\phi^{A}}}(\bm{Y}(t)),$

where the masses of the agents are assumed to be one and non-collective forces
are assumed to be zero. Our analysis can be extended to general second-order
systems (1) with known mass and non-collective force terms with slight
modifications.

###  3.1 The Representer theorem

In the classical regression setting [50], there is an interesting link between
GP regression and 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}_{M}},r_{\mathbb{V}_{M}},\mathbb{Z}_{\sigma^{2},M}\\}$
instead of the pairs
$\\{r_{\mathbb{X}_{M}},{\phi^{E}}(r_{\mathbb{X}_{M}}),{\phi^{A}}(r_{\mathbb{X}_{M}})\\}$,
where $r_{\mathbb{X}_{M}},r_{\mathbb{V}_{M}}\in\mathbb{R}^{MLN^{2}}$ are the
sets contains all the pairwise distances in $\mathbb{X}_{M}$, and
$\mathbb{V}_{M}$, i.e.

(16)
$r_{\mathbb{X}_{M}}=\\{r_{ij}^{\bm{X}^{(m,l)}}\\}_{i,j,m,l=1}^{N,N,M,L},\quad
r_{\mathbb{V}_{M}}=\\{r_{ij}^{\bm{V}^{(m,l)}}\\}_{i,j,m,l=1}^{N,N,M,L},\quad\mathbb{Z}_{\sigma^{2},M}=\\{\bm{Z}^{(m,l)}_{\sigma^{2}}\\}_{m,l=1}^{M,L},$

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 [54] that helps clarify the role of the hyperparameters.

We assume that $K^{E}$ and $K^{A}$ are two Mercer kernels defined on
$[0,R]\times[0,R]$ for some $R>0$. The true interaction functions
${\phi^{E}}\in\mathcal{H}_{{K}^{E}}$, ${\phi^{A}}\in\mathcal{H}_{{K}^{A}}$,
and

$\kappa^{2}_{E}=\mathrm{sup}_{r\in[0,R]}{{K}^{E}}(r,r)<\infty,$
$\kappa^{2}_{A}=\mathrm{sup}_{r\in[0,R]}{{K}^{A}}(r,r)<\infty.$

######  Theorem 3.1 (Representer theorem).

Let $K^{E}$ and $K^{A}$ be two Mercer kernels that satisfy Assumption (3.1).
Given the training data $\\{\mathbb{Y}_{M},\mathbb{Z}_{\sigma^{2},M}\\}$, 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{\bm{\phi}}=(\bar{\phi}^{E},\bar{\phi}^{A})$ in (12) coincides with the
minimizer 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

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

where $\lambda=\\{\lambda^{E},\lambda^{A}\\}$ and the estimator
$\bar{\bm{\phi}}\in\mathcal{H}_{{K}^{E}}\times\mathcal{H}_{{K}^{A}}$ can also
be represented by

(18) $\bar{\bm{\phi}}=(\sum_{r^{x}\in
r_{\mathbb{X}_{M}}}\hat{c}_{r^{x}}K_{r^{x}}^{E},\sum_{(r^{x},r^{v})\in(r_{\mathbb{X}_{M}}\times
r_{\mathbb{V}_{M}})}\hat{c}_{r^{v}}K_{r^{x}}^{A}),$

with

$\displaystyle\hat{\mathbf{c}}_{r^{x}}=\frac{1}{N}\bm{r}_{\mathbb{X}_{M}}^{T}\cdot(K_{\mathbf{f}_{\bm{\phi}}}(\mathbb{Y}_{M},\mathbb{Y}_{M})+\lambda^{E}NMLI_{dNML})^{-1}\mathbb{Z}_{\sigma^{2},M},$
(19)
$\displaystyle\hat{\mathbf{c}}_{r^{v}}=\frac{1}{N}\bm{r}_{\mathbb{V}_{M}}^{T}\cdot(K_{\mathbf{f}_{\bm{\phi}}}(\mathbb{Y}_{M},\mathbb{Y}_{M})+\lambda^{A}NMLI_{dNML})^{-1}\mathbb{Z}_{\sigma^{2},M},$

where $\hat{\mathbf{c}}_{r^{x}}$, $\hat{\mathbf{c}}_{r^{v}}$ are the
vectorizations of $(\hat{c}_{r^{x}})_{r^{x}\in r_{\mathbb{X}_{M}}}$ and
$(\hat{c}_{r^{v}})_{r^{v}\in r_{\mathbb{V}_{M}}}$ respectively,
$r_{\mathbb{X}_{M}}\times r_{\mathbb{V}_{M}}\in\mathbb{R}^{MLN^{2}\times
MLN^{2}}$ is the set containing all the pairwise distances in $\mathbb{X}_{M}$
and their associated pairwise distances in $\mathbb{V}_{M}$ as defined in
(16), $\bm{r}_{\mathbb{X}_{M}}$ is the block-diagonal matrix defined by
$\mathrm{diag}(\bm{r}_{\bm{X}^{(m,l)}})\in\mathbb{R}^{MLdN\times MLN^{2}}$ and
$\bm{r}_{\bm{X}^{(m,l)}}=\mathrm{diag}(\\{[\bm{r}_{i1}^{\bm{X}^{(m,l)}},\dots,\bm{r}_{iN}^{\bm{X}^{(m,l)}}]\\}_{i=1}^{N})\in\mathbb{R}^{dN\times
N^{2}}$, similarly for $\bm{r}_{\mathbb{V}_{M}}$.

Detailed proof of Theorem 3.1 is shown in Appendix Appendix D. From the
theorem, it is clear how hyperparameters affect the prediction of interaction
kernels: $\theta^{E}$, $\theta^{A}$, and $\sigma$ jointly affect the choice of
Mercer kernels and regularization constant, which becomes quite crucial in
real data applications (see Figure 5). In (3.1), we also see that the
posterior mean estimator $\phi^{E}$ lies in the span of basis functions with
indices determined by the pairwise distances, and their coefficients are
correlated with the basis functions. This is an effect imposed by the
structure of the governing equation encoded in $\mathbf{f}_{{\bm{\phi}}}$.

###  3.2 Well-posedness

We are concerned with the nonparametric learning of interaction kernels. That
is, we do not assume the parametric form of interaction kernels. In this case,
one can not expect to recover the true interaction kernels from finite data as
they live in infinite dimensional spaces. Therefore, it is important to ensure
one can asymptotically identify the true interaction kernels as the number of
observational data snapshots goes to infinity. Otherwise, the empirical
estimators from finite data will have limited value as a scientific and
predictive tool. Mathematically, we study the well-posedness under a
statistical inverse problem setting. We introduce a linear operator
$A:\mathcal{H}_{{K}^{E}}\times\mathcal{H}_{{K}^{A}}\rightarrow
L^{2}(\mathbb{R}^{2dN};\rho_{\bm{Y}};\mathbb{R}^{dN})$ defined by

(20) $\displaystyle A{\bm{\varphi}}=\mathbf{f}_{{\bm{\varphi}}},$

where $\mathbf{f}_{{\bm{\varphi}}}$ is the right hand side of system (15) by
replacing ${\bm{\phi}}$ with ${\bm{\varphi}}$, and $\rho_{\bm{Y}}$ is the
limiting measure on $\mathbb{R}^{2dN}$ that we assume the observational data
are sampled i.i.d from. For example, if we assume that the initial condition
of each trial is sampled i.i.d from a measure, then

(21)
$\displaystyle\rho_{\bm{Y}}(S)=\lim_{M\to\infty}\frac{1}{M}\sum_{m,l=1}^{M,L}\mathbb{1}_{\bm{Y}^{(m,l)}\in
S}$

for any Borel set $S\subset\mathbb{R}^{2dN}$ and the limit does exist in the
weak sense by the law of large numbers. We denote the marginal probability
measures for $\bm{X}$ and $\bm{V}$ by $\rho_{\bm{X}}$, $\rho_{\bm{V}}$
respectively.

Then the well-posedness of (20) is reduced to studying under which conditions
$A$ has a bounded inverse.

####  3.2.1 Well-posedness on an $L^{2}$ space

We first consider the embedding of
$\mathcal{H}_{{K}^{E}}\times\mathcal{H}_{{K}^{A}}$ to a suitable $L^{2}$ space
and consider the well-posedness in a weaker $L^{2}$-norm. Motivated by (18) in
the Representer theorem, we consider the measures $\tilde{\rho}_{r}^{E}$,
$\tilde{\rho}_{r}^{A}$ for ${\phi^{E}}$, ${\phi^{A}}$ based on the structure
of $\mathbf{f}_{{\bm{\varphi}}}$,

(22)
$\displaystyle\tilde{\rho}_{r}^{E}(Q)=\int_{Q}\int_{\mathbb{R}^{dN}}\frac{1}{N(N-1)}\sum_{i\neq
j}\delta_{r_{ij}^{{\bm{x}}^{\ast}}}(r)\cdot(r_{ij}^{{\bm{x}}^{\ast}})^{2}d\rho_{\bm{X}}(\bm{X}^{\ast})dr$
(23)
$\displaystyle\tilde{\rho}_{r}^{A}(Q)=\int_{Q}\int_{\mathbb{R}^{2dN}}\frac{1}{N(N-1)}\sum_{i\neq
j}\delta_{r_{ij}^{{\bm{x}}^{\ast}}}(r)\cdot(r_{ij}^{{\bm{v}}^{\ast}})^{2}d\rho_{\bm{Y}}(\bm{X}^{\ast},\bm{V}^{\ast})dr$

for any set $Q\subset[0,R]$, and $\delta(\cdot)$ is the Dirac $\delta$
distribution. By the continuity,
$\mathcal{H}_{{K}^{E}}\times\mathcal{H}_{{K}^{A}}$ can be naturally embedded
as a subspace of
$L^{2}([0,R]\times[0,R];\tilde{\rho}_{r};\mathbb{R}\times\mathbb{R})$ with
$\tilde{\rho}_{r}=\tilde{\rho}_{r}^{E}\times\tilde{\rho}_{r}^{A}$. One can
follow the proof of Proposition 9 in [43] to show that $A$ is a bounded linear
operator from
$L^{2}([0,R]\times[0,R];\tilde{\rho}_{r};\mathbb{R}\times\mathbb{R})$ to
$L^{2}(\mathbb{R}^{2dN};\rho_{\bm{Y}};\mathbb{R}^{dN})$.

Now we can introduce a sufficient condition to guarantee the existence of a
bounded inverse of $A$ on
$L^{2}([0,R]\times[0,R];\tilde{\rho}_{r};\mathbb{R}\times\mathbb{R})$, called
the coercivity condition:

######  Definition 3.2.

We say that the system (15) satisfies the coercivity condition if
$\forall{\bm{\varphi}}\in\mathcal{H}_{{K}^{E}}\times\mathcal{H}_{{K}^{A}}$,

(24)
$\displaystyle\|A{\bm{\varphi}}\|^{2}_{L^{2}(\rho_{\bm{Y}})}=\|\mathbf{f}_{{\bm{\varphi}}}\|^{2}_{L^{2}(\rho_{\bm{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})}$

for some constants $c_{\mathcal{H}_{{K}^{E}}},c_{\mathcal{H}_{{K}^{A}}}>0$.

Here we show one example to support the coercivity condition.

######  Theorem 3.3.

Consider $\rho_{\bm{Y}}=\begin{bmatrix}\rho_{\bm{X}}\\\
\rho_{\bm{V}}\end{bmatrix}$, where $\rho_{\bm{X}}$ is the product of $N$
independent and identical measures with compact support on $\mathbb{R}^{d}$,
and $\rho_{\bm{V}}$ is defined in the same way and is independent of
$\rho_{\bm{X}}$. Then we have

(25)
$\displaystyle\|\mathbf{f}_{{\bm{\varphi}}}\|^{2}_{L^{2}(\rho_{\bm{Y}})}\geq\frac{N-1}{N^{2}}\|\varphi^{E}\|^{2}_{L^{2}(\tilde{\rho}_{r}^{E})}+\frac{N-1}{N^{2}}\|\varphi^{A}\|^{2}_{L^{2}(\tilde{\rho}_{r}^{A})}$

Detailed proof of Theorem 3.3 is shown in Appendix Appendix C. In [4], the
identifiability of a structured sum of ${\phi^{E}}$ and ${\phi^{A}}$ is
studied. Here we consider a stronger version of identifiability as we want to
individually recover ${\phi^{E}}$ and ${\phi^{A}}$. Note that it is also
possible for distributions on $\mathbb{R}^{dN}$ with non-i.i.d
$\mathbb{R}^{d}$ components that satisfy the coercivity condition. Finally, we
remark that the coercivity condition (24) holds on measure pairs
$(\rho_{1},\rho_{2})$ equivalent to
$(\tilde{\rho}_{r}^{E},\tilde{\rho}_{r}^{A})$. This can provide us with many
nontrivial examples from the special case in Theorem 3.3. We conjecture that
the coercivity condition is generally satisfied and leave further
investigation as future work.

####  3.2.2 Well-posedness on
$\mathcal{H}_{{K}^{E}}\times\mathcal{H}_{{K}^{A}}$ and the convergence
analysis

Now we turn to study the well-posedness on
$\mathcal{H}_{{K}^{E}}\times\mathcal{H}_{{K}^{A}}$ with the stronger RKHS
norm, and we make the following assumption.

We assume that $\tilde{\rho}_{r}$ is non-degenerate on $[0,R]\times[0,R]$.

We remark that the above assumption is mild. For example, we can pick
$\rho_{\bm{Y}}$ to be a uniform measure supported on a large enough cube, then
$\tilde{\rho}_{r}$ satisfies the assumption.

It is straightforward to see that the coercivity condition implies injectivity
of $A$ on $\mathcal{H}_{{K}^{E}}\times\mathcal{H}_{{K}^{A}}$:
${\bm{\varphi}}=0$ everywhere on $[0,R]$ when $A{\bm{\varphi}}=0$ for
${\bm{\varphi}}\in\mathcal{H}_{{K}^{E}}\times\mathcal{H}_{{K}^{A}}$. This is
due to the non-degeneracy of $\tilde{\rho}_{r}$ on $[0,R]\times[0,R]$ and the
continuity of ${\bm{\varphi}}$. Therefore, $A$ is injective. However, showing
$A$ has a bounded inverse on
$\mathcal{H}_{{K}^{E}}\times\mathcal{H}_{{K}^{A}}$ is impossible when it is
infinitely dimensional, as $A$ is a compact operator. Suppose the coercivity
condition (24) holds, then following the theoretical framework developed in
[43], one could prove the well-posedness on a suitable subspace determined by
the source conditions on ${\phi^{E}},{\phi^{A}}$ following inverse problem
literature. In this case, it is possible to prove one could recover both
kernels with a statistically optimal rate under the corresponding RKHS norm.
We obtained the result for the single-kernel case in our recent work [43], and
we leave the work for the double-kernel case for the future investigation.

##  4 Numerical examples

In this section, we investigate the performance of the algorithm proposed in
Section 2 to show the effectiveness of model selection in (1). Specific
instances of (1) have found many applications in modeling the clustering,
swarming, and alignment behaviors of collective agents. The examples include
(1) Cucker-Smale dynamics (CS) with friction force ($m_{i}\equiv 1$,
${\phi^{E}}\equiv 0$, ${\phi^{A}}\neq 0$) in Section 4.2.1, (2) fish milling
dynamics (FM) with friction force ($m_{i}\equiv 1$, ${\phi^{E}}\neq 0$,
${\phi^{A}}\equiv 0$) in Section 4.2.2, (3) anticipation dynamics (AD)
($m_{i}\equiv 1$, ${\phi^{E}},{\phi^{A}}\neq 0$) in Section 4.2.3 and (4)
opinion dynamics (OD) with stubborn agents ($m_{i}\equiv 0$, ${\phi^{E}}\neq
0$, ${\phi^{A}}\equiv 0$) in Section 4.3. In (1)-(3), the mass of agents is
known in advance, i.e, they are second-order systems. We are interested in
learning ${\phi^{E}}$, ${\phi^{A}}$, and other hyperparameters $\bm{\alpha}$
from data, resulting in the selection of types of interactions (energy versus
alignment interactions). In (4), we used the prior knowledge that
${\phi^{A}}\equiv 0$ and investigate if the true zero mass of the opinions and
${\phi^{E}}$ can be learned from data, resulting in the selection of the order
of the system (first versus second order).

The detailed setups of each dynamic are shown in Table 4. We applied the
strategies proposed in Section 2 to learn $\bm{\alpha}$ in $F_{\bm{\alpha}}$,
and the interaction kernels ${\phi^{E}}(r)$, ${\phi^{A}}(r)$. We initialize
the parameters in $\bm{\alpha}$ randomly from the uniform distribution
$\mathcal{U}([0,1])$, and the same for $\sigma$ in the cases with noisy data.
In each experiment, we run 10 independent trials and report the errors of the
estimations for $\bm{\alpha}$, the estimation errors for ${\phi^{E}}$,
${\phi^{A}}$ in the (relative) $L^{\infty}([0,R])$-norm, and compare the
discrepancy between the true trajectories (evolved using $\alpha$,
${\phi^{E}}$, ${\phi^{A}}$) and predicted trajectories (evolved using
$\hat{\alpha}$, $\hat{{\phi^{E}}}$, $\hat{{\phi^{A}}}$) on both the training
time interval $[0,T]$ and on the future time interval $[T,T_{f}]$, over two
different sets of initial conditions (IC) – one taken from the training data,
and one consisting of new samples from the same initial distribution.

##### Real data application

We also apply our method to two real datasets of fish in Section 4.4, where
one shows a flocking behavior and another shows a milling behavior. We fit
them into the Cucker-Smale and fish milling dynamics respectively and perform
comparisons with two other classical approaches: SINDy [55] and feed-forward
neural networks.

##### Numerical Setup.

We simulate the trajectory data $(\mathbb{Y},\mathbb{Z})$ on the time interval
$[0,T]$ with given i.i.d initial conditions generated from the probability
measures specified for each system as shown in Table 4. For the training data
sets, we generate $M$ trajectories and observe each trajectory at $L$
equidistant times $0=t_{1}<t_{2}<\cdots<t_{L}=T$ and add Gaussian noise to
$\mathbb{Z}$ with level $\sigma$. We construct an empirical approximation to
the probability measure $\tilde{\rho}_{r}$, with $2000$ trajectories and let
$[0,R]$ be its support. All ODE systems are evolved using ode$15$s in MATLAB®
with a relative tolerance at $10^{-5}$ and absolute tolerance at $10^{-6}$.
For noise-free training data, we add a jitter constant $\approx 10^{-6}$ as a
way of regularization. We apply the minimize function in the GPML
package***Carl Edward Rasmussen & Hannes Nickisch
(http://gaussianprocess.org/gpml/code) to train the parameters using conjugate
gradient optimization with the partial derivatives shown in Section 2, and set
the maximum number of function evaluations to 400.

In almost all examples, we use the full GP methods, as we use scarce data and
there is no need for acceleration. However, we show the effectiveness of our
acceleration techniques in Fish milling dynamics in section 5 when we have a
larger scale of data.

Table 4: System parameters in the dynamics

System | CS | FM | AD | OD
---|---|---|---|---
$d$ | 2 | 2 | 2 | 1
$N$ | 10 | 10 | 10 | 5
$m_{i}$ | 1 | 1 | 1 | 0
$[0;T;T_{f}]$ | $[0;10;20]$ | $[0;5;10]$ | $[0;10;20]$ | $[0;2;20]$
$\mu_{0}^{\bm{x}}$ | $\mathrm{Unif}([-2,2]^{2})$ | $\mathrm{Unif}([-0.5,0.5]^{2})$ | $\mathrm{Unif}([0,5]^{2})$ | $\mathrm{Unif}([-1,1])$
$\mu_{0}^{\bm{v}}$ | $\mathrm{Unif}([-1,1]^{2})$ | $\mathrm{Unif}([0,0]^{2})$ | $\mathrm{Unif}([0,5]^{2})$ | -
${\phi^{E}}$ | 0 | $\frac{1}{r}\bigg{[}-e^{-2r}+e^{-\frac{r}{4}}\bigg{]},$ | $\frac{0.1}{(1+r)^{2.5}}+\frac{1}{(1+r)^{0.5}}$ | Eq. 29
${\phi^{A}}$ | $\frac{1}{(1+r^{2})^{1/4}}$ | 0 | $\frac{0.1}{(1+r^{2})^{0.5}}$ | 0
$F({\bm{x}}_{i},\dot{\bm{x}}_{i},\bm{\alpha})$ | $\kappa\dot{\bm{x}}_{i}(1-\|\dot{\bm{x}}_{i}\|^{p})$ | $(\gamma-\beta\|\dot{\bm{x}}_{i}\|^{2})\dot{\bm{x}}_{i}$ | 0 | Eq. 30
$\bm{\alpha}$ | $(\kappa,p)=(1,2)$ | $(\gamma,\beta)=(1.5,0.5)$ | - | $(P_{1},\kappa)=(1,10)$

##### Choice of the covariance function.

We choose the Matérn covariance function defined on $[0,R]\times[0,R]$ for the
Gaussian process priors in our numerical experiments, i.e.,

(26)
$K_{\theta}(r,r^{\prime})=s_{\phi}^{2}\frac{2^{1-\nu}}{\Gamma(\nu)}(\frac{\sqrt{2\nu}\|r-r^{\prime}\|}{\omega_{\phi}})^{\nu}B_{\nu}(\frac{\sqrt{2\nu}\|r-r^{\prime}\|}{\omega_{\phi}}),$

where the parameter $\nu>0$ determines the smoothness; $\Gamma(\nu)$ is the
Gamma function; $B_{\nu}$ is the modified Bessel function of the second kind;
and the hyperparameters $\theta=\\{s_{\phi}^{2},\omega_{\phi}\\}$ quantify the
amplitude and scale. In our numerical examples, we choose $\nu=p+1/2$ with
$p=0\text{ or }1$.

The Reproducing Kernel Hilbert Space (RKHS), $\mathcal{H}_{Mat\acute{e}rn}$,
associated with this Matérn kernel is norm-equivalent to the Sobolev space
$W_{2}^{\nu+1/2}([0,R])$ defined by

(27) $W_{2}^{\nu+1/2}([0,R]):=\Big{\\{}f\in
L_{2}([0,R]):\|f\|^{2}_{W_{2}^{\nu+1/2}}:=\sum_{\beta\in\mathbb{N}_{0}^{1}:|\beta|\leq\nu+1/2}\|D^{\beta}f\|_{L_{2}}^{2}<\infty\Big{\\}}.$

That is to say, $\mathcal{H}_{Mat\acute{e}rn}=W_{2}^{s}([0,R])$ as a set of
functions, and there exists constants $c_{1},c_{2}>0$ such that

(28)
$c_{1}\|f\|_{W_{2}^{\mu+\frac{1}{2}}}\leq\|f\|_{\mathcal{H}_{Mat\acute{e}rn}}\leq
c_{2}||f||_{W_{2}^{\mu+\frac{1}{2}}},\quad\forall
f\in\mathcal{H}_{Mat\acute{e}rn}.$

In other words, $\mathcal{H}_{Mat\acute{e}rn}$ consists of functions that are
differentiable up to order $\nu$ and weak differentiable up to order
$s=\nu+\frac{1}{2}$.

(a) CS: ${\phi^{E}}$ (b) CS: ${\phi^{A}}$ (c) CS: trajectory prediction (d)
FM: ${\phi^{E}}$

(e) FM: ${\phi^{A}}$ (f) FM: trajectory prediction (g) AD: ${\phi^{E}}$ (h)
AD: ${\phi^{A}}$

(i) AD: trajectory prediction

Figure 1: Results of learning different dynamics using the Matérn kernel. Top:
Learning CSF ($\\{N,M,L,\sigma\\}=\\{10,6,3,0.1\\}$); Middle: Learning FM
($\\{N,M,L,\sigma\\}=\\{10,3,3,0.1\\}$); Bottom: Learning AD
($\\{N,M,L,\sigma\\}=\\{10,12,3,0.01\\}$). Left, Center: Predictive mean
$\hat{{\phi^{E}}}$ and $\hat{{\phi^{A}}}$ of the true kernels, and two-
standard-deviation band (light blue color) around the means. The grey bars
represent the empirical density of the $\tilde{\rho}_{r}$; Right: the true
(left) versus predicted (right) trajectories using $\hat{\bm{\alpha}}$ and
$\hat{\phi}$ with initial conditions of training data (top) and testing data
(bottom)

###  4.1 Summary of the numerical experiments

  * •

The proposed learning approach performs simultaneous precise model selections
from small amounts of noisy observation data. The numerical results in all
different dynamics show that the algorithm can accurately identify the
existence of energy-based/alignment-based interactions and can learn order
information of dynamics between agents in the systems.

  * •

The GP method selects a kernel basis to represent the underlying sparse
dynamics that generalizes remarkably well in larger time prediction with new
initial conditions. The occasional larger prediction errors that occur in a
larger time interval may be caused by the propagation of estimation errors. We
believe the performance is satisfactory since we only have very limited and
noisy training data. Even in cases where the prediction errors are relatively
large, the estimators can predict remarkably accurate collective behaviors of
the agents, e.g. the consensus in the opinion dynamics, the flocking behavior
in the Cucker-Smale dynamics, and the milling pattern in the fish milling
dynamics.

  * •

In synthetic experiments, the uncertainty quantification band for the
trajectories is rather small ($\mathcal{O}(10^{-3})$), resulting from the
narrow uncertainty bands of ${\bm{\phi}}$. In real data experiments, we found
models using interaction kernels sampled from uncertainty bands all reproduced
the true dynamics very well.

  * •

The real data experiments show that the proposed GP approach combined with the
particle-based models is practically applicable, and outperformed two other
competitors in preserving the physics of the true dynamics.

###  4.2 Model selection for types of interaction kernels

####  4.2.1 Cucker-Smale dynamics with friction force

The Cucker-Smale system [56, 57, 58] is used to model collective behaviors in
a system of agents that follow a prescribed protocol of communication, such as
wedges of bird flocks, lattices in cell organization, or bee hives [59, 60,
61]. We consider the system of $N$ agents in the form (1) with components
defined in Table 4, where ${\phi^{A}}$ is a communication kernel, or influence
function, that makes the agents flock, and $F_{\bm{\alpha}}$ a Rayleigh-type
friction force that pushes all magnitudes of the velocities $\|v_{i}\|$
towards the same value 1 and counteracts the directional alignment forces
governed by ${\phi^{A}}$ to produce a rich variety of collective dynamics
depending on the relative strengths of the involved forces.

In this example, the unknown parameters $\bm{\alpha}=(\kappa,p)$ are nonlinear
with respect to the system. We show the errors of our estimation for
$\bm{\alpha}$ and ${\phi^{E}}$, ${\phi^{A}}$ in Table 5. Note that for this
model, ${\phi^{A}}$ is in the RKHS generated by the Matérn kernel we pick. The
estimated interaction kernel $\hat{\phi}^{A}$ can recover the true
${\phi^{A}}(r)$ almost perfectly in the region within the support of the
empirical $\tilde{\rho}_{r}$ from both noise-free and noisy training data.
Moreover, the true interaction kernel ${\phi^{A}}(r)$ is fully covered in the
uncertainty region we constructed using the posterior variances. Table 5 also
shows that our method can identify the nonexistence of the energy-based
interaction well with small errors (at most $O(10^{-4})$) from zero in
$L^{\infty}([0,R])$-norm. See also in Fig. 1(a),(b). The errors for the
predicted trajectories are shown in Table 6. We can see that in both the
training time interval $[0,10]$ and future time interval $[10,20]$, the
estimators can produce accurate approximations of the true trajectories and
the performance becomes better when we increase the size of training data ($M$
or $L$).

Table 5: Means and standard deviations of the errors of $\hat{\bm{\alpha}}$
(including $\hat{\sigma}$ when noise exists) and $\hat{\phi}$ for different
settings of the CS dynamics.

$\\{N,M,L,\sigma\\}$ | $\|\hat{\bm{\alpha}}-\bm{\alpha}\|_{\infty}$ | $\|\hat{{\phi^{E}}}-0\|_{\infty}$ | $\|\hat{{\phi^{A}}}-{\phi^{A}}\|_{\infty}/\|{\phi^{A}}\|_{\infty}$
---|---|---|---
$\\{10,1,3,0\\}$ | $1.9\cdot 10^{-3}\pm 1.0\cdot 10^{-3}$ | $2.1\cdot 10^{-5}\pm 4.0\cdot 10^{-5}$ | $5.6\cdot 10^{-2}\pm 1.5\cdot 10^{-2}$
$\\{10,3,3,0\\}$ | $\mathbf{1.1\cdot 10^{-3}\pm 7.9\cdot 10^{-4}}$ | $2.6\cdot 10^{-5}\pm 6.5\cdot 10^{-5}$ | $4.5\cdot 10^{-2}\pm 2.0\cdot 10^{-2}$
$\\{10,6,3,0\\}$ | $1.3\cdot 10^{-3}\pm 2.5\cdot 10^{-3}$ | $\mathbf{1.1\cdot 10^{-5}\pm 1.3\cdot 10^{-5}}$ | $\mathbf{3.2\cdot 10^{-2}\pm 1.0\cdot 10^{-2}}$
$\\{10,6,3,0.05\\}$ | $1.1\cdot 10^{-1}\pm 1.1\cdot 10^{-1}$ | $1.2\cdot 10^{-4}\pm 1.6\cdot 10^{-4}$ | $1.6\cdot 10^{-1}\pm 8.6\cdot 10^{-2}$
$\\{10,6,3,0.1\\}$ | $2.3\cdot 10^{-1}\pm 2.3\cdot 10^{-1}$ | $1.4\cdot 10^{-4}\pm 2.9\cdot 10^{-4}$ | $1.8\cdot 10^{-1}\pm 8.0\cdot 10^{-2}$

Table 6: The trajectory prediction errors for different settings.

$\\{N,M,L,\sigma\\}$ | Training IC $[0,10]$ | Training IC $[10,20]$ | new IC $[0,10]$ | new IC $[10,20]$
---|---|---|---|---
$\\{10,1,3,0\\}$ | $4.9\cdot 10^{-4}\pm 4.2\cdot 10^{-4}$ | $6.7\cdot 10^{-4}\pm 1.3\cdot 10^{-3}$ | $1.8\cdot 10^{-3}\pm 4.4\cdot 10^{-3}$ | $1.4\cdot 10^{-2}\pm 4.2\cdot 10^{-2}$
$\\{10,3,3,0\\}$ | $2.5\cdot 10^{-4}\pm 2.0\cdot 10^{-4}$ | $1.5\cdot 10^{-4}\pm 1.3\cdot 10^{-4}$ | $4.9\cdot 10^{-4}\pm 4.9\cdot 10^{-4}$ | $8.7\cdot 10^{-3}\pm 1.7\cdot 10^{-2}$
$\\{10,6,3,0\\}$ | $\mathbf{1.5\cdot 10^{-4}\pm 1.2\cdot 10^{-4}}$ | $\mathbf{9.4\cdot 10^{-5}\pm 9.2\cdot 10^{-5}}$ | $\mathbf{2.7\cdot 10^{-4}\pm 4.1\cdot 10^{-4}}$ | $\mathbf{2.3\cdot 10^{-4}\pm 4.6\cdot 10^{-4}}$
$\\{10,6,3,0.05\\}$ | $2.3\cdot 10^{-2}\pm 1.3\cdot 10^{-2}$ | $1.9\cdot 10^{-2}\pm 1.3\cdot 10^{-2}$ | $2.7\cdot 10^{-2}\pm 1.9\cdot 10^{-2}$ | $2.5\cdot 10^{-2}\pm 2.0\cdot 10^{-2}$
$\\{10,6,3,0.1\\}$ | $4.2\cdot 10^{-2}\pm 2.6\cdot 10^{-2}$ | $3.8\cdot 10^{-2}\pm 2.8\cdot 10^{-2}$ | $4.9\cdot 10^{-2}\pm 3.4\cdot 10^{-2}$ | $4.5\cdot 10^{-2}\pm 3.9\cdot 10^{-2}$

####  4.2.2 Fish-Milling dynamics with friction force

In this subsection, we consider another type of cohesive collective system
that produces milling patterns [62, 63]. A special instance of such systems is
the D’Orsogna model [12, 61, 64], which describes the motion of $N$ self-
propelled particles powered by biological or mechanical motors, that
experience a frictional force, and can produce a rich variety of collective
patterns. We consider the system of $N$ agents of the form (1) with components
defined in Table 4, where the interaction kernel ${\phi^{E}}$ is derived from
the Morse-type potential. Since it is singular at $r=0$, we truncate it at
$r_{0}=0.05$ with a function of the form $ae^{-br}$ to ensure the new function
has a continuous derivative. The force function $F^{{\bm{v}}}$ includes self-
propulsion with strength $\gamma$ and nonlinear drag with strength $\beta$.

The errors of the estimations for $\bm{\alpha}$ after our training procedure
and the learned ${\phi^{E}}$, ${\phi^{A}}$ are shown in Table 7. In this
model, ${\phi^{E}}$ is in the RKHS generated by the chosen Matérn kernel. We
can see that our estimators produced faithful approximations to the true
kernel based on the results we report in Table 7 and Fig. 1(d),(e). They also
show that we can identify the nonexistence of the alignment-based interaction
with very small errors and select the correct model. The discrepancy between
the true trajectories and the predicted trajectories on both the training time
interval $[0,5]$ and future time interval $[5,10]$ are shown in Table 8. Even
if the trajectory prediction errors can go up to $O(10^{-1})$ with the
presence of a relatively large noise for the systems with $N=10$, our
estimators provided faithful predictions to most of the agents in the system
and the milling pattern as shown in Fig. 1(f).

Table 7: Means and standard deviations of the errors of $\hat{\bm{\alpha}}$
(including $\hat{\sigma}$ when noise exists) and $\hat{\phi}$ for different
settings of FM dynamics.

$\\{N,M,L,\sigma\\}$ | $\|\hat{\bm{\alpha}}-\bm{\alpha}\|_{\infty}$ | $\|\hat{{\phi^{E}}}-{\phi^{E}}\|_{\infty}/\|{\phi^{E}}\|_{\infty}$ | $\|\hat{{\phi^{A}}}-0\|_{\infty}$
---|---|---|---
$\\{10,1,3,0\\}$ | $7.9\cdot 10^{-4}\pm 1.0\cdot 10^{-3}$ | $3.6\cdot 10^{-2}\pm 4.3\cdot 10^{-3}$ | $6.6\cdot 10^{-4}\pm 6.9\cdot 10^{-4}$
$\\{10,1,9,0\\}$ | $6.4\cdot 10^{-5}\pm 6.2\cdot 10^{-5}$ | $3.9\cdot 10^{-2}\pm 2.7\cdot 10^{-3}$ | $1.6\cdot 10^{-4}\pm 1.3\cdot 10^{-4}$
$\\{10,3,3,0\\}$ | $\mathbf{4.7\cdot 10^{-5}\pm 5.0\cdot 10^{-5}}$ | $3.8\cdot 10^{-2}\pm 5.4\cdot 10^{-3}$ | $1.2\cdot 10^{-4}\pm 1.7\cdot 10^{-4}$
$\\{10,3,3,0.01\\}$ | $3.4\cdot 10^{-3}\pm 1.9\cdot 10^{-3}$ | $\mathbf{2.9\cdot 10^{-2}\pm 5.7\cdot 10^{-3}}$ | $2.9\cdot 10^{-3}\pm 4.3\cdot 10^{-3}$
$\\{10,3,3,0.05\\}$ | $1.4\cdot 10^{-2}\pm 8.5\cdot 10^{-3}$ | $4.9\cdot 10^{-2}\pm 1.5\cdot 10^{-2}$ | $\mathbf{4.6\cdot 10^{-5}\pm 7.0\cdot 10^{-5}}$
$\\{10,3,3,0.1\\}$ | $3.5\cdot 10^{-2}\pm 7.2\cdot 10^{-2}$ | $7.1\cdot 10^{-2}\pm 2.0\cdot 10^{-2}$ | $2.9\cdot 10^{-2}\pm 9.0\cdot 10^{-2}$

Table 8: The trajectory prediction errors for different settings of FM
dynamics.

$\\{N,M,L,\sigma\\}$ | Training IC $[0,5]$ | Training IC $[5,10]$ | new IC $[0,5]$ | new IC $[5,10]$
---|---|---|---|---
$\\{10,1,3,0\\}$ | $2.1\cdot 10^{-3}\pm 2.0\cdot 10^{-3}$ | $1.0\cdot 10^{-2}\pm 8.7\cdot 10^{-3}$ | $1.9\cdot 10^{-3}\pm 1.9\cdot 10^{-3}$ | $5.4\cdot 10^{-3}\pm 4.4\cdot 10^{-3}$
$\\{10,1,9,0\\}$ | $\mathbf{3.4\cdot 10^{-4}\pm 2.9\cdot 10^{-4}}$ | $\mathbf{1.4\cdot 10^{-3}\pm 1.2\cdot 10^{-3}}$ | $\mathbf{4.7\cdot 10^{-4}\pm 4.2\cdot 10^{-4}}$ | $\mathbf{1.3\cdot 10^{-3}\pm 1.2\cdot 10^{-3}}$
$\\{10,3,3,0\\}$ | $8.1\cdot 10^{-4}\pm 8.0\cdot 10^{-4}$ | $2.2\cdot 10^{-3}\pm 2.0\cdot 10^{-3}$ | $8.8\cdot 10^{-4}\pm 8.8\cdot 10^{-4}$ | $3.5\cdot 10^{-3}\pm 2.8\cdot 10^{-3}$
$\\{10,3,3,0.01\\}$ | $8.3\cdot 10^{-3}\pm 3.8\cdot 10^{-3}$ | $1.8\cdot 10^{-2}\pm 1.2\cdot 10^{-2}$ | $6.6\cdot 10^{-3}\pm 3.2\cdot 10^{-3}$ | $1.4\cdot 10^{-2}\pm 9.3\cdot 10^{-3}$
$\\{10,3,3,0.05\\}$ | $3.4\cdot 10^{-2}\pm 2.1\cdot 10^{-2}$ | $7.1\cdot 10^{-2}\pm 4.7\cdot 10^{-2}$ | $3.7\cdot 10^{-2}\pm 1.9\cdot 10^{-2}$ | $7.0\cdot 10^{-2}\pm 4.7\cdot 10^{-2}$
$\\{10,3,3,0.1\\}$ | $8.0\cdot 10^{-2}\pm 9.8\cdot 10^{-2}$ | $1.5\cdot 10^{-1}\pm 1.9\cdot 10^{-1}$ | $9.5\cdot 10^{-2}\pm 1.3\cdot 10^{-1}$ | $1.5\cdot 10^{-1}\pm 2.3\cdot 10^{-1}$

####  4.2.3 Anticipation Dynamics

In this subsection, we consider a more complicated model where the
interactions depend on both the pairwise distance and the differences in
velocities, i.e. both ${\phi^{E}}$ and ${\phi^{A}}$ are nonzero. The
anticipation dynamics (AD) models in [13] are suitable candidates, and we
consider the system of $N$ agents in the form (1) with components defined in
Table 4.

The errors of the estimations for $\bm{\alpha}$ and the learned ${\phi^{E}}$,
${\phi^{A}}$ are shown in Table 9. In this model, both ${\phi^{E}}$ and
${\phi^{A}}$ are in the RKHS generated by the chosen Matérn kernel. We can see
that our estimators produced faithful approximations to both true kernels
based on the results we report in Table 9 and Fig. 1(g),(h). The comparisons
between the true trajectories and the predicted trajectories on both the
training time interval $[0,10]$ and future time interval $[10,20]$ are shown
in Table 10. The estimators can produce accurate approximations of the true
trajectories with errors at most $O(10^{-1})$, see also Fig. 1(i).

Table 9: Means and standard deviations of the errors of $\hat{\sigma}$ (when
noise exists) and $\hat{\phi}$ for different settings of AD dynamics.

$\\{N,M,L,\sigma\\}$ | $\|\hat{\sigma}-\sigma\|_{\infty}$ | $\|\hat{{\phi^{E}}}-{\phi^{E}}\|_{\infty}/\|{\phi^{E}}\|_{\infty}$ | $\|\hat{{\phi^{A}}}-{\phi^{A}}\|_{\infty}/\|{\phi^{A}}\|_{\infty}$
---|---|---|---
$\\{10,3,3,0\\}$ | - | $9.2\cdot 10^{-2}\pm 7.4\cdot 10^{-3}$ | $4.5\cdot 10^{-2}\pm 1.0\cdot 10^{-2}$
$\\{10,6,3,0\\}$ | - | $7.9\cdot 10^{-2}\pm 6.7\cdot 10^{-3}$ | $4.3\cdot 10^{-2}\pm 5.1\cdot 10^{-3}$
$\\{10,12,3,0\\}$ | - | $\mathbf{7.4\cdot 10^{-2}\pm 6.1\cdot 10^{-3}}$ | $\mathbf{3.6\cdot 10^{-2}\pm 7.0\cdot 10^{-3}}$
$\\{10,12,3,0.005\\}$ | $8.8\cdot 10^{-5}\pm 5.1\cdot 10^{-5}$ | $1.3\cdot 10^{-1}\pm 1.7\cdot 10^{-2}$ | $7.3\cdot 10^{-2}\pm 3.2\cdot 10^{-2}$
$\\{10,12,3,0.01\\}$ | $1.8\cdot 10^{-4}\pm 9.9\cdot 10^{-5}$ | $1.6\cdot 10^{-1}\pm 1.9\cdot 10^{-2}$ | $9.3\cdot 10^{-2}\pm 4.1\cdot 10^{-2}$

Table 10: The trajectory prediction errors for different settings of AD
dynamics.

$\\{N,M,L,\sigma\\}$ | Training IC $[0,10]$ | Training IC $[10,20]$ | new IC $[0,10]$ | new IC $[10,20]$
---|---|---|---|---
$\\{10,3,3,0\\}$ | $\mathbf{4.2\cdot 10^{-4}\pm 3.8\cdot 10^{-4}}$ | $\mathbf{2.3\cdot 10^{-4}\pm 2.1\cdot 10^{-4}}$ | $6.1\cdot 10^{-4}\pm 8.4\cdot 10^{-4}$ | $3.5\cdot 10^{-4}\pm 5.0\cdot 10^{-4}$
$\\{10,6,3,0\\}$ | $6.6\cdot 10^{-4}\pm 7.4\cdot 10^{-4}$ | $3.8\cdot 10^{-4}\pm 4.1\cdot 10^{-4}$ | $7.1\cdot 10^{-4}\pm 9.2\cdot 10^{-4}$ | $3.9\cdot 10^{-4}\pm 5.2\cdot 10^{-4}$
$\\{10,12,3,0\\}$ | $6.2\cdot 10^{-4}\pm 6.8\cdot 10^{-4}$ | $3.3\cdot 10^{-4}\pm 3.7\cdot 10^{-4}$ | $\mathbf{3.7\cdot 10^{-4}\pm 5.2\cdot 10^{-4}}$ | $\mathbf{2.1\cdot 10^{-4}\pm 3.1\cdot 10^{-4}}$
$\\{10,12,3,0.005\\}$ | $1.9\cdot 10^{-3}\pm 2.1\cdot 10^{-3}$ | $1.1\cdot 10^{-3}\pm 1.2\cdot 10^{-3}$ | $1.1\cdot 10^{-3}\pm 1.2\cdot 10^{-3}$ | $6.8\cdot 10^{-4}\pm 7.1\cdot 10^{-4}$
$\\{10,12,3,0.01\\}$ | $3.4\cdot 10^{-3}\pm 4.3\cdot 10^{-3}$ | $1.9\cdot 10^{-3}\pm 2.4\cdot 10^{-3}$ | $1.9\cdot 10^{-3}\pm 2.1\cdot 10^{-3}$ | $1.2\cdot 10^{-3}\pm 1.3\cdot 10^{-3}$

###  4.3 Model selection for the order of systems

An example of opinion dynamics is shown below to test the validity of our
method for identifying the order of dynamic systems. This is a first-order
system of $N$ interacting agents, and each agent $i$ is characterized by a
continuous opinion variable $x_{i}\in\mathbb{R}$. The dynamics of opinion
exchange are governed by the first-order equation mentioned in Section 2.2.1
with

(29) ${\phi^{E}}(r)=\begin{cases}25r&\textrm{if }0\leq r<0.4,\\\ 10&\textrm{if
}0.4\leq r<0.6,\\\ 25-25r&\textrm{if }0.6\leq r<1,\\\ 0&\textrm{if }r\geq
1.\end{cases}$

The interaction kernel ${\phi^{E}}$ encodes the non-repulsive interactions
between agents: all agents aim to align their opinions to their connected
neighbors according to distance-based attractive influences. We consider the
case where there is no non-collective force, i.e.
$F_{i}({\bm{x}}_{i},\bm{\alpha})\equiv 0$. We also consider a more complicated
case where there exist stubborn agents, i.e.

(30)
$F_{i}({\bm{x}}_{i},\bm{\alpha})=\begin{cases}-\kappa({\bm{x}}_{i}-P_{i})&\textrm{if
agent $i$ is stubborn with bias $P_{i}$}\\\ 0&\textrm{otherwise}\end{cases}$

where $F_{i}({\bm{x}}_{i},\bm{\alpha})$ describes the additional influence
induced by the stubbornness: the stubborn agents have strong desires to follow
their bias $P_{i}$, and $\kappa$ controls the rate of convergence towards
their bias. The stubborn agents may cause a major effect on the collective
opinion formation process. If $\kappa=0$, then stubborn agents do not follow
their biases and behave as regular agents.

Table 11 shows the errors of the estimations for $m$, $\bm{\alpha}$, and
${\phi^{E}}(r)$ in 10 independent trails of experiments. It shows our method
can identify the order of dynamics with the estimation of $m\approx 0$ and
learn the interaction kernel $\phi$ simultaneously, see also Fig. 2.

Table 11: Means and standard deviations of the errors of $\hat{m}$ (including
$\hat{\sigma}$ when noise exists) and $\hat{\phi}$ for different settings of
OD dynamics.

Model | $\\{N,M,L,\sigma\\}$ | $\|\hat{m}-0\|_{\infty}$ | $\|\hat{\bm{\alpha}}-\bm{\alpha}\|_{\infty}$ | $\|\hat{\phi}-\phi\|_{\infty}/\|\phi\|_{\infty}$
---|---|---|---|---
OD | $\\{5,6,3,0\\}$ | $8.5\cdot 10^{-4}\pm 9.0\cdot 10^{-4}$ | - | $3.8\cdot 10^{-3}\pm 1.1\cdot 10^{-3}$
OD | $\\{5,6,3,0.1\\}$ | $4.8\cdot 10^{-3}\pm 5.2\cdot 10^{-4}$ | $3.2\cdot 10^{-2}\pm 1.6\cdot 10^{-2}$ | $1.1\cdot 10^{-2}\pm 5.6\cdot 10^{-3}$
ODS | $\\{10,3,3,0\\}$ | $5.5\cdot 10^{-4}\pm 2.8\cdot 10^{-4}$ | $7.2\cdot 10^{-2}\pm 4.1\cdot 10^{-2}$ | $5.2\cdot 10^{-2}\pm 4.4\cdot 10^{-2}$
ODS | $\\{10,3,3,0.1\\}$ | $3.8\cdot 10^{-3}\pm 1.8\cdot 10^{-3}$ | $9.0\cdot 10^{-1}\pm 1.1\cdot 10^{0}$ | $3.3\cdot 10^{-2}\pm 1.9\cdot 10^{-2}$

Figure 2: Model selection of OD ($\\{N,M,L\\}=\\{5,6,3\\}$) and $\sigma=0,0.1$
using the Matérn kernel. Left: Predictive mean $\hat{\phi}$ of the true
kernel, and two-standard-deviation band (light blue color) around the mean.
The grey bars represent the empirical density of the $\tilde{\rho}_{r}$.
Right: the true (left) versus predicted (right) trajectories using
$\hat{\bm{\alpha}}$ and $\hat{\phi}$ with initial conditions of training data
(top) and testing data (bottom) when $\sigma=0.1$.

###  4.4 Real fish data

Finally, we test the performance of our method using two real datasets of
swimming fish by Couzin et al., which are available at ScholarsArchive of
Oregon State University†††Katz, Yael, Kolbjorn Tunstrom, Christos C Ioannou,
Cristian Huepe, and Iain D Couzin, 2021. The URL address is
https://ir.library.oregonstate.edu/concern/datasets/zk51vq07c. The
experimental arena consisted of a white shallow tank of size $2.1\times 1.2$ m
($7\times 4$ ft) surrounded by a floor-to-ceiling white curtain. Water depth
was chosen to be $4.5$-$5$ cm so the schools would be approximately 2D. We
consider two data sets, one is from frame 2201 to frame 2296 which consists of
50 fish and forms a flocking behavior, and another one is from frame 4601 to
frame 4798 which consists of 124 fish and forms a milling behavior. We relabel
them as frame 0 to frame 95 and frame 0 to frame 198 respectively, refer to
more details of the dataset in the supplementary information of [65]. We first
normalize the position data into the region [0,1], and then we smooth the data
by using a moving window average with a window size of 10 frames and apply the
finite difference method to calculate the velocities and accelerations.

##### Flocking behavior example

For the first data set, as shown in Fig. 3(e), the fish will eventually follow
approximately the same direction as time evolves. In this case, the magnitude
of the velocity data of fish is relatively small. The velocities can be
considered the same, as long as their normalized direction vectors are very
close. So the fish exhibit approximate flocking behavior (i.e.
$\|\bm{v}_{i}-\bm{v}_{c}\|\approx 0$ for all $i$ and some common velocity
$\bm{v}_{c}$). Therefore, we use the Cucker-Smale system shown in Section
4.2.1 to model the flocking behavior, i.e. considering the governing equation
(1) with corresponding interaction kernel ${\phi^{A}}$ and force
$F({\bm{x}}_{i},\dot{\bm{x}}_{i},\bm{\alpha})$, $\bm{\alpha}=(\kappa,p)$ shown
in Table 4 for CS dynamics.

The training data consists of frame 0 and frame 28. In the training procedure,
we first use a subset of data with two selected agents, the initialization of
hyperparameters for $\theta^{E},\theta^{A},\sigma$, and $\bm{\alpha}$ are
$(1,1),(1,1),0.001,(1,1)$, and we set the length of runs in the minimizer
solver to be 100. The results shown in Fig. 3(a),(b) suggest there only exist
alignment-based interactions since the estimated energy-based interaction
$\hat{\phi}^{E}\equiv 0$. Therefore, we use all data to learn the system with
only ${\phi^{A}}$. After we obtain the estimators, we run the learned
dynamical system on the time interval [0,20] with frame 0 as the initial
condition. We find that the simulated position data at $t=19$ matches the
position data at frame 95 very well. We then compare the original position
data set with the simulated ones at $t=0:0.2:19$.

##### Milling behavior example

For the second data set, as shown in Fig. 3(g), the fish will eventually
follow approximately a milling pattern. Therefore, we use the Fish-Milling
system shown in Section 4.2.2 to model the milling behavior, i.e. considering
the governing equation (1) with interaction kernel ${\phi^{E}}$ and force
$F({\bm{x}}_{i},\dot{\bm{x}}_{i},\bm{\alpha})$, $\bm{\alpha}=(\gamma,\beta)$
shown in Table 4 for FM dynamics.

The training data consists of frame 0 and frame 28. In the training procedure,
we first use a subset of data with two selected agents, the initialization of
hyperparameters for $\theta^{E},\theta^{A},\sigma$, and $\bm{\alpha}$ are
$(1,1),(1,1),0.001,(1,1)$, and we set the length of runs in the minimizer
solver to be 100. With the estimated parameters, we obtained the estimators
for ${\phi^{E}}$ and ${\phi^{A}}$ using all data, and run the learned
dynamical system at the time interval [0,38] with the frame 0 as the initial
condition.

(a) Real data CS: ${\phi^{E}}$

(b) Real data CS: ${\phi^{A}}$

(c) Real data FM: ${\phi^{E}}$

(d) Real data FM: ${\phi^{E}}$

(e) Real data CS trajectory prediction (f) Real data CS comparisons

(g) Real data FM trajectory prediction (h) Real data FM comparisons

Figure 3: (a)(b): Fitting into a Cucker-Smale system (dim=100), estimated
${\phi^{E}}$ and ${\phi^{A}}$ with all data; (c)(d): Fitting into a Fish-
Milling system (dim=248), estimated ${\phi^{E}}$ and ${\phi^{A}}$ with data of
2 agents; (e)(g): true dynamics (left) v.s. the predicted dynamics using our
proposed approach, SINDy model, and FNN model, with frame 0 as the initial
condition; (f)(h): baseline comparisons using the group polarisation parameter
$M(t)$.

##### Measure of Performance

To evaluate the performance at the group level, we consider the group
polarisation $M(t)$ [66], which is a vector order parameter that encapsulates
both the direction and degree of the fish alignment, which is defined by

(31) $M(t)=\frac{1}{N}\sum_{i=1}^{N}\hat{{\bm{v}}}_{i}(t)$

where $\hat{{\bm{v}}}_{i}(t)={\bm{v}}_{i}(t)/\|{\bm{v}}_{i}(t)\|$ is the
direction of motion of the i-th fish (at time t). When $|M|$ is close to 1,
the fish are moving in a coherent direction, whereas when $|M|$ is close to
zero, there is no prevailing direction and individual motion is effectively
isotropic.

##### Baseline Comparisons

We perform comparisons with approaches that learn the right-hand side function
of (1) directly from trajectory data: the first one is SINDy [55], which aims
at finding a sparse representation for each row of governing equations in a
(typically large) dictionary; the second one is regression using feed-forward
neural networks, for which we use the MATLAB® 2021a Deep Learning Toolbox™.

For the SINDy model, we apply a reasonably large dictionary consisting of
monomials up to order 2, sines, and cosines of frequencies
$\\{k\\}_{k=1}^{10}$. For the neural network model, we consider a three-layer
FNN (Feed-Forward Neural Network) with $[50,50,25]$ hidden units for the
flocking behavior example, and a two-layer FNN with $[40,20]$ hidden units for
the milling behavior example.

The predictive trajectories for the flocking behavior example using different
models are shown in Fig. 3(e). We compare the performances in terms of group
polarisation $M(t)$ in Fig. 3(f), and the order-1 Wasserstein distances
between the empirical distributions of $M(t)$ in true data and the predicted
dynamics are $W(p_{M}^{true},p_{M}^{gp})=0.0112$,
$W(p_{M}^{true},p_{M}^{SINDy})=0.5497$, and
$W(p_{M}^{true},p_{M}^{NN})=0.5869$, see Fig. 4 (Left). The results for the
milling behavior example are shown in Fig. 3(g). We compare the performances
in Fig. 3(h), and the order-1 Wasserstein distances between the empirical
distributions of $M(t)$ in the true data and the predicted dynamics are
$W(p_{M}^{true},p_{M}^{gp})=0.0087$, $W(p_{M}^{true},p_{M}^{SINDy})=0.0530$,
and $W(p_{M}^{true},p_{M}^{NN})=0.1055$, see Fig. 4 (Right). In both examples,
we can see that although both predictions using the SINDy and FNN models look
similar to the true trajectories, based on the group polarisation parameter
$M(t)$ and comparing the changes of $M(t)$ in $t$ or the empirical
distributions of $M(t)$, only our model using GP captures the group behaviors.

Figure 4: Baseline comparisons using the group polarisation parameter $M(t)$.
Left: the empirical distribution of $M(t)$ in the flocking behavior example,
the order-1 Wasserstein distances between true data and the predicted dynamics
are $W(p_{M}^{true},p_{M}^{gp})=0.0112$,
$W(p_{M}^{true},p_{M}^{SINDy})=0.5497$, and
$W(p_{M}^{true},p_{M}^{NN})=0.5869$. Right: the empirical distribution of
$M(t)$ in the milling behavior example, the order-1 Wasserstein distances
between true data and the predicted dynamics are
$W(p_{M}^{true},p_{M}^{gp})=0.0087$, $W(p_{M}^{true},p_{M}^{SINDy})=0.0530$,
and $W(p_{M}^{true},p_{M}^{NN})=0.1055$.

We also compare our result with two other GP models in the milling behavior
example, where the parameters $\bm{\alpha}$ are not estimated properly: (1) we
use the initial values of hyperparameters, i.e. let
$\theta^{E},\theta^{A},\sigma$ and $\bm{\alpha}$ equal
$(1,1),(1,1),0.001,(1,1)$, and do not train those hyperparameters; (2) we
apply the noise-free model, i.e. do not consider noise and let $\sigma\equiv
0$. The results of these two GP models are shown in Fig. 5.

Figure 5: Fitting into a Fish-Milling system (dim=248). Top: estimated
${\phi^{E}}$, ${\phi^{A}}$, and predictive trajectory with initial parameters,
i.e. no training (NT) for hyperparameters $\theta^{E},\theta^{A},\sigma$, and
$\bm{\alpha}$; Bottom: estimated ${\phi^{E}}$, ${\phi^{A}}$, and predictive
trajectory using noise-free (NF) model, i.e. $\hat{\sigma}=0$.

##  5 Acceleration Result Comparison

We now present our acceleration (see Section 2.5) results for a 20-dimensional
Fish Milling (FM) system with increasing observational data. When we have a
larger amount of observational data, we will focus on learning
$\sigma,\gamma,$ and $\beta$, and use $\theta^{E}=\theta^{A}=1$ for a default
prior with the Matérn kernel. We will show the impact on kernel predictions is
minimal.

We use $\nu=\frac{3}{2}$ for all examples below. All results shown are
averaged over $10$ complete runs with standard deviation included where we
used the same training data but with initialized hyperparameters uniformly at
random from an interval centered at the ground truth with radius 0.5 in each
trial. We use the randomized Gaussian Nystrom preconditioner [52] for all
tests with rank the floor of $\frac{30}{\log(12)}\cdot\log(\frac{NdML}{10})$.
While we would ideally use the effective rank of our kernel matrix, this is
expensive to compute in practice and we resort to empirical approximation for
our trials.

Figure 6: Log plot of absolute error of learned hyperparameters $\gamma$
(red), $\beta$ (yellow), and $\sigma$ (blue) for the FM system
($\\{N,M,L\\}=\\{20,M,6\\}$) with varying $M$. True values are
$\sigma=0.01,\gamma=1.5,\beta=0.5$. Shown also is a runtime comparison with
Full GP.

Figure 7: Accelerated kernel prediction error for the experiments above. Left
is the $L^{\infty}$ error, relative for $\phi^{E}$ and absolute for $\phi^{A}$
as $\phi^{A}=0$ is the ground truth. Right is the $L_{\rho}^{2}$ error,
relative for $\phi^{E}$ and absolute for $\phi^{A}$.

Figure 8: Relative trajectory prediction error on testing data, FM system ($\\{N,M,L\\}=\\{20,M,6\\}$) with varying $M$. This plot uses test error on the full interval of $[0,10]$. Prediction error quickly goes to zero in testing. Shown also is a runtime comparison with full GP. Table 12: System parameters in Fish Milling above $d$ | $N$ | $[0;T;T_{f}]$ | $\alpha=(\gamma,\beta)$ | $\mu_{0}^{\bm{x}}$ | $\mu_{0}^{\bm{v}}$
---|---|---|---|---|---
2 | 20 | $[0,5,10]$ | $(1.5,0.5)$ | $\mathrm{Unif}([-1,1]^{2})$ | $\mathrm{Unif}([0,0]^{2})$

Figure 8 shows that our hyperparameter learning method is able to accurately
recover the hyperparameters $\sigma,\gamma,\beta$ with greatly improved
runtime compared to the full GP method. Once these hyperparameters are
learned, our acceleration can also be utilized for the prediction of the
kernel, which also has a low observed error, see Figure 8. Most errors of the
kernel prediction occur away from the support of observed data in the FM
system and do not affect the trajectory prediction of our system. This is
quantified in the very low relative $L^{2}$ error of the predicted
trajectories of the FM system using our predicted kernel as shown in Figure 8.

These results provide clear evidence of successful acceleration options while
maintaining highly acceptable accuracy. While the running time of accelerated
MLE can still be expensive for prohibitively large data, the accelerated
method scales much better than the fully explicit method and opens up exciting
possibilities in modeling large datasets. We note that prediction also scales
quite well and relies only upon PCG and preconditioner choice, allowing the
usage of efficient cross-validation techniques for hyperparameter choice in
certain classes of problems.

Our central findings are the following:

  * •

We have discovered that accurate hyperparameter recovery can be achieved using
a small set of observational data, and using more training data does not
necessarily improve the accuracy. This is due to the lack of consistency in
the training of MLE, which is a well-known result in Gaussian process
regression. We recommend that one should split a small subset for
hyperparameter tuning and then use the full dataset for kernel learning. We
have seen empirical success with this method.

  * •

In 10 trials with small $M$, we often observed one or two trials with
relatively large recovery errors in hyperparameters. We removed these outliers
from our data before plotting above. We attribute this to the instability of
the Lanczos algorithm or the non-convexity of the optimization problem, as in
these cases, we observed that the minimization of MLE stopped very early.
Nonetheless, we would like to point out that even in these cases, we obtained
very satisfying performance in kernel learning and trajectory prediction.

  * •

There are additional opportunities for acceleration in kernel learning that
depend on the specific problem and infrastructure available. For instance, in
the case of $\nu=\frac{1}{2}$, we may exploit sparsity in a decomposition of
the kernel matrix
$K_{\mathbf{f}_{{\bm{\phi}}}}(\mathbb{Y},\mathbb{Y};\theta)$, as developed in
[67], while maintaining desired exactness. We leave the extension of this
method to all half-integer $\nu$ values for future work. Furthermore, there
are avenues for accelerating GP learning using modern hardware. With access to
GPUs, one can parallelize the explicit construction of kernel matrices and the
Lanczos algorithm calculations. These steps are embarrassingly parallel and
allow for demanding much greater accuracy.

##  6 Final remarks and future work

In this work, we present an approach based on Gaussian processes to perform
the model selection of particle/agent-based models from scarce and noisy data.
We propose efficient acceleration techniques to improve the scalability. The
methodology is extendable to cover heterogeneous systems with multiple types
of agents and external potentials. It is also possible to extend the learning
approach to the mean-field limits of the particle models. Another line of
future work is to apply the quantitative framework developed in this paper to
design a data acquisition plan (active learning). The goal is to optimize the
kernel learning using the least amount of trajectory data by looking at their
marginal pairwise distance distributions. We leave it as future work.

## Acknowledgments

Charles Kulick was partially supported by NSF DMS-2111303. S.T. was partially
supported by Hellman Family Faculty Fellowship, and the NSF DMS-2111303. S.T.
would like to thank Hengrui Luo and Didong Li for their helpful discussions.

##  7 Appendix

##  Appendix A Learning approach for model selection

Our learning approach is a generalization of the methodology proposed in [43];
to be self-contained, we state the detailed formulation here.

######  Lemma A.1.

Let ${\bm{\phi}}=({\phi^{E}},{\phi^{A}})$ be two Gaussian processes with mean
zero and covariance function
$K_{\theta^{E}},K_{\theta^{A}}:[0,R]\times[0,R]\to\mathbb{R}$ respectively,
i.e.,
$\phi^{\mathrm{type}}\sim\mathcal{GP}(0,K_{\theta^{\mathrm{type}}}(r,r^{\prime}))$,
type = $E$ or $A$, and
$\bm{m}\bm{Z}(t)=F_{\bm{\alpha}}(\bm{Y}(t))+\mathbf{f}_{\bm{\phi}}(\bm{Y}(t))$
as defined in (2). Then for any $t,t^{\prime}\in[0,T]$, we have that,

(32) $\begin{bmatrix}\bm{m}\bm{Z}(t)\\\
\bm{m}\bm{Z}(t^{\prime})\end{bmatrix}\sim\mathcal{N}\left(\begin{bmatrix}F_{\bm{\alpha}}(\bm{Y}(t))\\\
F_{\bm{\alpha}}(\bm{Y}(t^{\prime}))\end{bmatrix},K_{\mathbf{f}_{{\bm{\phi}}}}(\bm{Y}(t),\bm{Y}(t^{\prime}))\right),$

where $K_{\mathbf{f}_{\phi}}(\bm{X}(t),\bm{X}(t^{\prime})))$ is the covariance
matrix
$\mathrm{Cov}(\mathbf{f}_{\phi}(\bm{X}(t)),\mathbf{f}_{\phi}(\bm{X}(t^{\prime})))$
with $(i,j)$th block

$\displaystyle\mathrm{Cov}([\mathbf{f}_{{\bm{\phi}}}(\bm{Y})]_{i},[\mathbf{f}_{{\bm{\phi}}}(\bm{Y}^{\prime})]_{j})=\frac{1}{N^{2}}\sum_{k\neq
i,k^{\prime}\neq
j}\big{(}K_{\theta^{E}}(r^{\bm{x}}_{ik},r^{{\bm{x}}^{\prime}}_{jk^{\prime}})\bm{r}^{\bm{x}}_{ik}{\bm{r}^{{\bm{x}}^{\prime}}_{jk^{\prime}}}^{T}+K_{\theta^{A}}(r^{\bm{x}}_{ik},r^{{\bm{x}}^{\prime}}_{jk^{\prime}})\bm{r}^{\bm{v}}_{ik}{\bm{r}^{{\bm{v}}^{\prime}}_{jk^{\prime}}}^{T}\big{)},$

see Table 3 for the definitions.

######  Proof A.2.

For $\phi\sim\mathcal{GP}(0,K_{\theta}(r,r^{\prime}))$, and any
$r,r^{\prime}\in[0,R]$, we have that,

(34) $\displaystyle\mathbb{E}[\phi(r)]$ $\displaystyle=$ $\displaystyle 0,$
(35) $\displaystyle\mathrm{Cov}[\phi(r),\phi(r^{\prime})]$ $\displaystyle=$
$\displaystyle K_{\theta}(r,r^{\prime}).$

Therefore, for any collection of states $\\{r_{i}\\}_{i=1}^{n}\subset[0,R]$,
and $\\{a_{i}\\}_{i=1}^{n},\\{b_{i}\\}_{i=1}^{n}\subset\mathbb{R}$, the linear
operator on function values
$\mathcal{L}(\\{\phi(r_{i})\\}_{i=1}^{n}):=(a_{i}\phi(r_{i})+b_{i})_{i=1}^{n}$
satisfies

(36)
$\mathcal{L}(\\{\phi(r_{i})\\}_{i=1}^{n})\sim\mathcal{N}(\mathrm{vec}(\\{b_{i}\\}_{i=1}^{n}),\Sigma_{\mathcal{L}(\phi)}),$

where $\mathcal{N}$ denotes the Gaussian distribution,
$\mathrm{vec}(\\{b_{i}\\}_{i=1}^{n})\in\mathbb{R}^{n}$ is the vectorization of
$\\{b_{i}\\}_{i=1}^{n}$, and the covariance matrix
$\Sigma_{\mathcal{L}(\phi)}=\\{a_{i}a_{j}K_{\theta}(r_{i},r_{j})\\}_{i,j=1}^{n}\in\mathbb{R}^{n\times
n}$.

Therefore, since ${\phi^{E}}$, ${\phi^{A}}$ are independent, and
$\mathbf{f}_{{\bm{\phi}}}(\bm{Y}(t))$ is linear in ${\bm{\phi}}$, for any $t$,
$t^{\prime}$, we have that

(37) $\begin{bmatrix}\mathbf{f}_{{\bm{\phi}}}(\bm{Y}(t))\\\
\mathbf{f}_{{\bm{\phi}}}(\bm{Y}(t^{\prime}))\end{bmatrix}\sim\mathcal{N}(\bm{0},K_{\mathbf{f}_{{\bm{\phi}}}}(\bm{Y}(t),\bm{Y}(t^{\prime}))),$

where $K_{\mathbf{f}_{{\bm{\phi}}}}(\bm{Y}(t),\bm{Y}(t^{\prime})))$ is the
covariance matrix

(38)
$\displaystyle\mathrm{Cov}(\mathbf{f}_{{\bm{\phi}}}(\bm{Y}(t)),\mathbf{f}_{{\bm{\phi}}}(\bm{Y}(t^{\prime})))=\big{(}\mathrm{Cov}([\mathbf{f}_{{\bm{\phi}}}(\bm{Y}(t))]_{i},[\mathbf{f}_{{\bm{\phi}}}(\bm{Y}(t^{\prime})]_{j}))\big{)}_{i,j=1}^{N,N},$

with $(i,j)$th block

$\displaystyle\mathrm{Cov}([\mathbf{f}_{{\bm{\phi}}}(\bm{Y})]_{i},[\mathbf{f}_{{\bm{\phi}}}(\bm{Y}^{\prime})]_{j})=\frac{1}{N^{2}}\sum_{k\neq
i,k^{\prime}\neq
j}\big{(}K_{\theta^{E}}(r^{\bm{x}}_{ik},r^{{\bm{x}}^{\prime}}_{jk^{\prime}})\bm{r}^{\bm{x}}_{ik}{\bm{r}^{{\bm{x}}^{\prime}}_{jk^{\prime}}}^{T}+K_{\theta^{A}}(r^{\bm{x}}_{ik},r^{{\bm{x}}^{\prime}}_{jk^{\prime}})\bm{r}^{\bm{v}}_{ik}{\bm{r}^{{\bm{v}}^{\prime}}_{jk^{\prime}}}^{T}\big{)},$

Thus, by (2), the observation $\bm{Z}$ in the model follows the Gaussian
distribution

(39) $\begin{bmatrix}\bm{m}\bm{Z}(t)\\\
\bm{m}\bm{Z}(t^{\prime})\end{bmatrix}\sim\mathcal{N}(\begin{bmatrix}F_{\bm{\alpha}}(\bm{Y}(t))\\\
F_{\bm{\alpha}}(\bm{Y}(t^{\prime}))\end{bmatrix},K_{\mathbf{f}_{{\bm{\phi}}}}(\bm{Y}(t),\bm{Y}(t^{\prime}))).$

Then suppose that the training data consists of

(40)
$\displaystyle\\{\mathbb{Y}_{M},\mathbb{Z}_{\sigma^{2},M}\\}=\\{\mathbb{X}_{M},\mathbb{V}_{M},\mathbb{Z}_{\sigma^{2},M}\\}$

with

$\displaystyle\mathbb{X}_{M}$
$\displaystyle=\mathrm{Vec}\big{(}\\{\bm{X}^{(m,l)}\\}_{m,l=1}^{M,L}\big{)}\in\mathbb{R}^{dNML},$
$\displaystyle\mathbb{V}_{M}$
$\displaystyle=\mathrm{Vec}\big{(}\\{\bm{V}^{(m,l)}\\}_{m,l=1}^{M,L}\big{)}=\mathrm{Vec}\big{(}\\{\dot{\bm{X}}^{(m,l)}\\}_{m,l=1}^{M,L}\big{)}\in\mathbb{R}^{dNML},$
$\displaystyle\mathbb{Z}_{\sigma^{2},M}$
$\displaystyle=\mathrm{Vec}\big{(}\\{\bm{Z}^{(m,l)}_{\sigma^{2}}\\}_{m,l=1}^{M,L}\big{)}=\mathrm{Vec}\big{(}\\{\ddot{\bm{X}}^{(m,l)}+\sigma^{2}\bm{\epsilon}^{(m,l)}\\}_{m,l=1}^{M,L}\big{)}\in\mathbb{R}^{dNML}$

where we observe the dynamics at $0=t_{1}<t_{2}<\cdots<t_{L}=T$; $m$ indexes
trajectories corresponding to different initial conditions at $t_{1}=0$;
$\bm{X}^{(m,1)}\stackrel{{\scriptstyle i.i.d}}{{\sim}}\mu_{0}^{\bm{x}}$,
$\bm{V}^{(m,1)}\stackrel{{\scriptstyle i.i.d}}{{\sim}}\mu_{0}^{\bm{v}}$,
$(\mu_{0}^{\bm{x}},\mu_{0}^{\bm{v}})$ are two independent probability measure
on $\mathbb{R}^{dN}$; the noise term
$\bm{\epsilon}^{(m,l)}\stackrel{{\scriptstyle
i.i.d}}{{\sim}}\mathcal{N}(\bm{0},I_{dN})$; we assume that
$\mu_{0}=(\mu_{0}^{\bm{x}},\mu_{0}^{\bm{v}})$ is independent of the
distribution of noise.

Applying Lemma A.1, we now derive the negative log marginal likelihood for
training parameters $\bm{\alpha}$, $\bm{\theta}$, and $\sigma$, with given
observational data as specified above.

######  Proposition A.3.

Denote $\bm{Y}^{(m,l)}=\bm{Y}^{(m)}(t_{l})$ and
$\bm{Z}^{(m,l)}_{\sigma^{2}}=\bm{Z}^{(m)}(t_{l})+\epsilon^{(m,l)}$ with i.i.d
noise $\epsilon^{(m,l)}\sim\mathcal{N}(0,\sigma^{2}I_{dN\times dN})$. Suppose
we are given the training data set
$(\mathbb{Y}_{M},\mathbb{Z}_{\sigma^{2},M}):=$
$\\{(\bm{Y}^{(m,l)},\bm{Z}^{(m,l)}_{\sigma^{2}})\\}_{m,l=1}^{M,L}$ for
$M,L\in\mathbb{N}$, such that

(41)
$\bm{Z}^{(m,l)}_{\sigma^{2}}=F_{\bm{\alpha}}(\bm{Y}^{(m,l)})+\mathbf{f}_{{\bm{\phi}}}(\bm{Y}^{(m,l)})+\epsilon^{(m,l)},$

with $F_{\bm{\alpha}}$, $\mathbf{f}_{{\bm{\phi}}}$ defined in Table 3. Then
the negative log marginal likelihood of $\mathbb{Z}_{\sigma^{2},M}$ given
$\mathbb{Y}_{M}$ and parameters $\bm{\alpha}$, $\theta$, $\sigma$ satisfies

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

where $I_{dNML}$ is the identity matrix of consistent size.

######  Proof A.4.

Using Lemma A.1, since $\epsilon^{(m,l)}$ is i.i.d Gaussian noise and is
independent of the initial distributions, we have that

(44)
$\bm{m}\mathbb{Z}_{\sigma^{2},M}\sim\mathcal{N}(F_{\bm{\alpha}}(\mathbb{Y}_{M}),K_{\mathbf{f}_{{\bm{\phi}}}}(\mathbb{Y}_{M},\mathbb{Y}_{M};\bm{\theta})+\sigma^{2}I_{dNML}),$

where the mean vector
$F_{\bm{\alpha}}(\mathbb{Y}_{M})=\mathrm{Vec}((F_{\bm{\alpha}}(\bm{Y}^{(m,l)}))_{m,l=1}^{M,L})\in\mathbb{R}^{dNML}$,
and the covariance matrix
$K_{\phi}(\mathbb{Y}_{M},\mathbb{Y}_{M};\theta)=\big{(}\mathrm{Cov}(\mathbf{f}_{{\bm{\phi}}}(\bm{Y}^{(i,j)}),\mathbf{f}_{{\bm{\phi}}}(\bm{Y}^{(i^{\prime},j^{\prime})}))\big{)}_{i,i^{\prime},j,j^{\prime}=1}^{M,M,L,L}\in\mathbb{R}^{dNML\times
dNML}$ can be computed by using (38). According to the properties of the
Gaussian distribution, given $\mathbb{Y}_{M}$ and parameters $\bm{\alpha}$,
$\theta$, $\sigma$, we have the negative log marginal likelihood function as
shown in (43).

As mentioned in the main text, we can apply the gradient-based method [68], to
minimize the negative log marginal likelihood and solve for the
hyperparameters $(\bm{\alpha},\bm{\theta},\sigma)$.

######  Proposition A.5.

Let
$\bm{\gamma}=(K_{\mathbf{f}_{{\bm{\phi}}}}(\mathbb{Y}_{M},\mathbb{Y}_{M};\bm{\theta})+\sigma^{2}I)^{-1}(\bm{m}\mathbb{Z}_{\sigma^{2},M}-F_{\bm{\alpha}}(\mathbb{Y}_{M}))$.
The partial derivatives of the marginal likelihood w.r.t. the parameters
$\bm{\alpha}$, $\bm{\theta}$, and $\sigma$ can be computed as follows:

(45) $\displaystyle\frac{\partial}{\partial\bm{\alpha}_{i}}\log
p(\bm{m}\mathbb{Z}_{\sigma^{2},M}|\mathbb{Y}_{M},\bm{\alpha},\bm{\theta},\sigma^{2})$
$\displaystyle=\bm{\gamma}^{T}\frac{\partial
F_{\bm{\alpha}}(\mathbb{Y}_{M})}{\partial\bm{\alpha}_{i}}.$ (46)
$\displaystyle\frac{\partial}{\partial\bm{\theta}_{j}}\log
p(\bm{m}\mathbb{Z}_{\sigma^{2},M}|\mathbb{Y}_{M},\bm{\alpha},\bm{\theta},\sigma^{2})$
$\displaystyle=\frac{1}{2}\mathrm{Tr}\left((\bm{\gamma}\bm{\gamma}^{T}-(K_{\mathbf{f}_{{\bm{\phi}}}}(\mathbb{Y}_{M},\mathbb{Y}_{M};\bm{\theta})+\sigma^{2}I)^{-1})\frac{\partial
K_{\mathbf{f}_{{\bm{\phi}}}}(\mathbb{Y}_{M},\mathbb{Y}_{M};\bm{\theta})}{\partial\bm{\theta}_{j}}\right).$
(47) $\displaystyle\frac{\partial}{\partial\sigma}\log
p(\bm{m}\mathbb{Z}_{\sigma^{2},M}|\mathbb{Y}_{M},\bm{\alpha},\bm{\theta},\sigma^{2})$
$\displaystyle=\mathrm{Tr}\left((\bm{\gamma}\bm{\gamma}^{T}-(K_{\mathbf{f}_{{\bm{\phi}}}}(\mathbb{Y}_{M},\mathbb{Y}_{M};\bm{\theta})+\sigma^{2}I)^{-1})\right)\sigma.$

With the updated prior $\phi$ from $\bm{\theta}$, and the parameters
$\bm{\alpha}$, $\sigma$, we show the detailed derivation of our estimators for
the prediction $\phi(r^{*})$ at $r^{\ast}\in[0,R]$.

######  Theorem A.6.

Suppose we are given the training data set
$(\mathbb{Y}_{M},\mathbb{Z}_{\sigma^{2},M}):=\\{(\bm{Y}^{(m,l)},\bm{Z}^{(m,l)}_{\sigma^{2}})\\}_{m,l=1}^{M,L}$
defined in Proposition A.3, and the hyperparameters
$(\bm{\alpha},\bm{\theta},\sigma)$ are known. Then for any $r^{\ast}\in[0,R]$,
type = $E$ or $A$, $\phi^{\mathrm{type}}(r^{\ast})$ satisfies

(48)
$p(\phi^{\mathrm{type}}(r^{\ast})|\mathbb{Y}_{M},\mathbb{Z}_{\sigma^{2},M})\sim\mathcal{N}(\bar{\phi}^{\mathrm{type}},var(\bar{\phi}^{\mathrm{type}})),$

where

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

and
$K_{\mathbf{f}_{\bm{\phi}},\phi^{\mathrm{type}}}(\mathbb{Y}_{M},r^{*})=K_{\phi^{\mathrm{type}},\mathbf{f}_{\bm{\phi}}}(r^{*},\mathbb{Y}_{M})^{T}$
denotes the covariance matrix between
$\mathbf{f}_{{\bm{\phi}}}(\mathbb{Y}_{M})$ and $\phi^{\mathrm{type}}(r^{*})$.

######  Proof A.7.

Since $\mathbf{f}_{{\bm{\phi}}}(\mathbb{Y}_{M})$ is defined componentwisely as
in (2), for any $r^{\ast}\in[0,R]$, we have that

(51) $\begin{bmatrix}\mathbf{f}_{{\bm{\phi}}}(\mathbb{Y}_{M})\\\
\phi^{\mathrm{type}}(r^{\ast})\end{bmatrix}\sim\mathcal{N}\left(0,\begin{bmatrix}K_{\mathbf{f}_{{\bm{\phi}}}}(\mathbb{Y}_{M},\mathbb{Y}_{M})&K_{\mathbf{f}_{\bm{\phi}},\phi^{\mathrm{type}}}(\mathbb{Y}_{M},r^{\ast})\\\
K_{\phi^{\mathrm{type}},\mathbf{f}_{\bm{\phi}}}(r^{\ast},\mathbb{Y}_{M})&K_{\theta^{\mathrm{type}}}(r^{\ast},r^{\ast})\end{bmatrix}\right),$

where $K_{\mathbf{f}_{\bm{\phi}}}(\mathbb{Y}_{M},\mathbb{Y}_{M})$ is the
covariance matrix between $\mathbf{f}_{{\bm{\phi}}}(\mathbb{Y}_{M})$ and
$\mathbf{f}_{{\bm{\phi}}}(\mathbb{Y}_{M})$ as we defined in Proposition A.3,
and
$K_{\mathbf{f}_{\bm{\phi}},\phi^{\mathrm{type}}}(\mathbb{Y}_{M},r^{*})=K_{\phi^{\mathrm{type}},\mathbf{f}_{\bm{\phi}}}(r^{*},\mathbb{Y}_{M})^{T}$
is the covariance matrix between $\mathbf{f}_{{\bm{\phi}}}(\mathbb{Y}_{M})$
and $\phi^{\mathrm{type}}(r^{*})$, i.e.,
$K_{\mathbf{f}_{\bm{\phi}},\phi^{\mathrm{type}}}(\mathbb{Y}_{M},r^{*})=(\mathrm{Cov}(\mathbf{f}_{{\bm{\phi}}}(\bm{Y}^{(m,l)}),\phi^{\mathrm{type}}(r^{\ast})))_{m,l=1}^{M,L}$
and the i-th component of
$\mathrm{Cov}(\mathbf{f}_{{\bm{\phi}}}(\bm{Y}^{(m,l)}),\phi^{\mathrm{type}}(r^{\ast}))$
is computed by

(52)
$\displaystyle\mathrm{Cov}([\mathbf{f}_{{\bm{\phi}}}(\bm{Y}^{(m,l)})]_{i},{\phi^{E}}(r^{\ast}))$
$\displaystyle=\frac{1}{N}\sum_{k\neq
i}K_{\theta^{E}}(r_{ik}^{\bm{X}^{(m,l)}},r^{\ast})\bm{r}_{ij}^{\bm{X}^{(m,l)}},$
(53)
$\displaystyle\mathrm{Cov}([\mathbf{f}_{{\bm{\phi}}}(\bm{Y}^{(m,l)})]_{i},{\phi^{A}}(r^{\ast}))$
$\displaystyle=\frac{1}{N}\sum_{k\neq
i}K_{\theta^{A}}(r_{ik}^{\bm{X}^{(m,l)}},r^{\ast})\bm{r}_{ij}^{\bm{V}^{(m,l)}}.$

Note that
$\bm{m}\bm{Z}^{(m,l)}_{\sigma^{2}}=F_{\bm{\alpha}}(\bm{Y}^{(m,l)})+\mathbf{f}_{\phi}(\bm{X}^{(m,l)})+\epsilon^{(m,l)}$
with i.i.d noise $\epsilon^{(m,l)}\sim\mathcal{N}(0,\sigma^{2}I_{dN})$ for all
$(m,l)$, so we have

(54)
$\begin{bmatrix}\bm{m}\mathbb{Z}_{\sigma^{2},M}-F_{\bm{\alpha}}(\mathbb{Y}_{M})\\\
\phi^{\mathrm{type}}(r^{\ast})\end{bmatrix}\sim\mathcal{N}\left(0,\begin{bmatrix}K_{\mathbf{f}_{\phi}}(\mathbb{Y}_{M},\mathbb{Y}_{M})+\sigma^{2}I_{dNML}&K_{\mathbf{f}_{\phi},\phi^{\mathrm{type}}}(\mathbb{Y}_{M},r^{\ast})\\\
K_{\phi^{\mathrm{type}},\mathbf{f}_{\bm{\phi}}}(r^{\ast},\mathbb{Y}_{M})&K_{\theta^{\mathrm{type}}}(r^{\ast},r^{\ast})\end{bmatrix}\right),$

Therefore, based on the properties of the joint Gaussian distribution (see
Lemma D.6), conditioning on $(\mathbb{Y}_{M},\mathbb{Z}_{\sigma^{2},M})$, we
have that

(55)
$p(\phi^{\mathrm{type}}(r^{\ast})|\mathbb{Y}_{M},\mathbb{Z}_{\sigma^{2},M},r^{\ast})\sim\mathcal{N}(\bar{\phi}^{\mathrm{type}},var(\bar{\phi}^{\mathrm{type}})),$

where $\bar{\phi}^{\mathrm{type}}$ and $var(\bar{\phi}^{\mathrm{type}})$ are
defined as in (49) and (50).

##  Appendix B Psuedocode for Acceleration

In this section, we discuss in detail the acceleration of the computations
used in our GP framework. We first review the bottleneck in the computation:
our likelihood function evaluation is very slow, as it involves inverting the
kernel matrix
$(K_{\mathbf{f}_{{\bm{\phi}}}}(\mathbb{Y},\mathbb{Y};\theta)+\sigma^{2}I)^{-1}$
and computing the log determinant
$\log\det(K_{\mathbf{f}_{{\bm{\phi}}}}(\mathbb{Y},\mathbb{Y};\theta)+\sigma^{2}I)$.
We also require evaluation of the gradient for exact optimization, which
further requires evaluation of the trace
$\text{Tr}((K_{\mathbf{f}_{{\bm{\phi}}}}(\mathbb{Y},\mathbb{Y};\theta)+\sigma^{2}I)^{-1}\frac{\partial
K_{\mathbf{f}_{{\bm{\phi}}}}(\mathbb{Y},\mathbb{Y};\theta)}{\partial\theta_{i}})$
for each parameter $\theta_{i}$ as shown in Proposition A.5.

Our primary goal is to avoid explicit inversion of the kernel matrix entirely
by utilizing the Preconditioned Conjugate Gradient (PCG) algorithm, see
Algorithm 2. PCG is an iterative method that can solve systems $Ax=b$ for $x$
without explicitly inverting $A$ through clever choices of update at each
step. This algorithm is central for scalability when solving large-scale
linear systems with positive definite matrices in the numerical linear algebra
literature. Note that the standard CG method is unlikely to work, as our
kernel matrix is likely to be very ill-conditioned. An efficient
preconditioner will be necessary to avoid extremely slow convergence. As
mentioned in the main paper, we recommend the Randomized Gaussian Nystrom
preconditioner for improving performance.

Input:
$\color[rgb]{0,0,0}K_{\mathbf{f}_{{\bm{\phi}}}}(\mathbb{Y},\mathbb{Y};\theta)+\sigma^{2}I\color[rgb]{0,0,0}$
(matrix-vector multiplication of kernel),
$\color[rgb]{0,0,0}P\color[rgb]{0,0,0}$ (preconditioner), $\bm{b}$ (target
vector), $\bm{x_{0}}$ (initial guess), $\mathrm{errorTol}$ (error tolerance),
$t$ (iterations)

1:
$r_{0}:=\bm{b}-\color[rgb]{0,0,0}(K_{\mathbf{f}_{{\bm{\phi}}}}(\mathbb{Y},\mathbb{Y};\theta)+\sigma^{2}I)\color[rgb]{0,0,0}\bm{x_{0}}$

2: $z_{0},d_{0}:=\color[rgb]{0,0,0}P\color[rgb]{0,0,0}r_{0}$

3: while $\left\|r_{n}\right\|>\mathrm{errorTol}$ and $n<t$

4:
$v_{n}:=\color[rgb]{0,0,0}(K_{\mathbf{f}_{{\bm{\phi}}}}(\mathbb{Y},\mathbb{Y};\theta)+\sigma^{2}I)\color[rgb]{0,0,0}d_{n-1}$

5: $\alpha_{n}:=\frac{r_{n-1}^{T}z_{n-1}}{d_{n-1}^{T}v_{n}}$

6: $x_{n}:=x_{n-1}+\alpha_{n}d_{n-1}$

7: $r_{n}:=r_{n-1}-\alpha_{n}v_{n}$

8: $z_{n}:=\color[rgb]{0,0,0}P\color[rgb]{0,0,0}r_{n}$

9: $\beta_{n}:=\frac{z_{n}^{T}r_{n}}{z_{n-1}^{T}r_{n-1}}$

10: $d_{n}:=z_{n}+\beta_{n}d_{n-1}$

Output: $x_{n}$ (solution to
$\color[rgb]{0,0,0}(K_{\mathbf{f}_{{\bm{\phi}}}}(\mathbb{Y},\mathbb{Y};\theta)+\sigma^{2}I)\color[rgb]{0,0,0}\bm{x}=\bm{b}$),

$\\{(\alpha_{i},\beta_{i})\text{ for all }i\leq n\\}$ (exclusively for
constructing Lanczos weights)

Algorithm 2 Preconditioned CG for solving
$(K_{\mathbf{f}_{{\bm{\phi}}}}(\mathbb{Y},\mathbb{Y};\theta)+\sigma^{2}I)\bm{x}=\bm{b}$

Now we can solve the problem of slow likelihood function evaluations. Instead
of inversion, our proper preconditioner will allow us to apply PCG and reduce
the computational complexity from cubic for inversion to quadratic, see
Algorithm 2. Note that PCG is only limited by the runtime of matrix-vector
multiplication, and in the presence of sparsity or other structural features
that allow for linear time matrix-vector multiplication, the complexity of PCG
will also reduce to linear time. This can be accomplished in the
$\nu=\frac{1}{2}$ case using [69].

Then we consider the log determinant evaluation. Using stochastic Lanczos
quadrature, we can instead compute an estimator for
$\text{Tr(log(}\color[rgb]{0,0,0}(K_{\mathbf{f}_{{\bm{\phi}}}}(\mathbb{Y},\mathbb{Y};\theta)+\sigma^{2}I)\color[rgb]{0,0,0}))$,
see Algorithm 3. This algorithm requires quadrature weights, but these can be
efficiently recovered by running the PCG algorithm and arranging
$\alpha,\beta$ in a tridiagonal matrix, as seen in [70]. Then we apply
stochastic trace estimation, as developed in [51]. These methods also extend
to gradient calculations.

Input:
$\color[rgb]{0,0,0}(K_{\mathbf{f}_{{\bm{\phi}}}}(\mathbb{Y},\mathbb{Y};\theta)+\sigma^{2}I)\color[rgb]{0,0,0}$
(matrix-vector multiplication for kernel matrix),
$\color[rgb]{0,0,0}P\color[rgb]{0,0,0}$ (preconditioner), $n$ (number of test
vectors), $m$ (number of Lanczos coefficients)

1: for $i$ from $1$ to $n$

2: $v_{n}\sim\text{Rademacher}$ {draw from Rademacher distribution}

3:
$T:=PCG(\color[rgb]{0,0,0}(K_{\mathbf{f}_{{\bm{\phi}}}}(\mathbb{Y},\mathbb{Y};\theta)+\sigma^{2}I)\color[rgb]{0,0,0},\color[rgb]{0,0,0}P\color[rgb]{0,0,0},v_{n},\ell=m)$
{get Lanczos coefficients}

4: $[W,\lambda]:=eig(T)$

5: for $j$ from $1$ to $m$

6: $\gamma_{i}:=\gamma_{i}+W_{1,j}^{2}\log(\lambda_{j})$

7:
$tr_{est}:=\log\det(\color[rgb]{0,0,0}P\color[rgb]{0,0,0})+\frac{NdML}{n}\sum_{i=1}^{n}\gamma_{i}$

Output: $tr_{est}$ (estimated trace of
$\log(\color[rgb]{0,0,0}(K_{\mathbf{f}_{{\bm{\phi}}}}(\mathbb{Y},\mathbb{Y};\theta)+\sigma^{2}I)\color[rgb]{0,0,0})$)

Algorithm 3 Stochastic Trace Estimation with Lanczos Quadrature

When choosing a preconditioner, we must have a method for fast and accurate
computation of matrix-vector multiplication by
$\color[rgb]{0,0,0}P^{-1}\color[rgb]{0,0,0}$ and evaluation of $\log\det(P)$,
as these are necessary operations in the above algorithms.

One widely applicable class of preconditioners for positive semi-definite
matrices is the low-rank Nystrom approximation. The central idea is to create
a low-rank approximation $P$ of a matrix $A$ of interest, with the expectation
that $P^{-1}A$ will have a condition number close to $1$. One common
implementation is to subsample $r$ columns of the matrix and use these to
construct an approximation for the missing entries with rank at most $r$.

The randomized Gaussian Nystrom preconditioner builds on this idea. Written in
a general form, we have
$P=A\Omega(\Omega^{T}A\Omega)^{\dagger}\Omega^{T}A^{T}$ for a chosen matrix
$\Omega\in\mathbb{R}^{NdML\times r}$. Column subsampling is a special case
where columns of the matrix $\Omega$ have a single non-zero entry of the unit
$1$. However, $\Omega$ can also be populated with randomized Gaussian entries.
This idea, developed in [52], has resulted in better empirical performance and
enjoys theoretical support. For an implementation see Algorithm 4.

Input:
$\color[rgb]{0,0,0}K_{\mathbf{f}_{{\bm{\phi}}}}(\mathbb{Y},\mathbb{Y};\theta)\color[rgb]{0,0,0}$
(matrix-vector multiplication for kernel matrix), $\sigma^{2}$ (noise
hyperparameter), $r$ (rank of preconditioner)

1: $\Omega\sim$ Standard Gaussian $\in\mathbb{R}^{NdML\times r}$

2: $R=\text{QR}(\Omega)$ using economy QR

3:
$Y=\color[rgb]{0,0,0}K_{\mathbf{f}_{{\bm{\phi}}}}(\mathbb{Y},\mathbb{Y};\theta)\color[rgb]{0,0,0}R$

4: $\nu=\text{eps}(||Y||_{F})$

5: $Y_{\nu}=Y+\nu R$

6: $C=\text{chol}(R^{T}Y_{\nu})$

7: $B=Y_{\nu}/C$

8: $U,\Sigma=\text{svd}(B)$ using economy svd

9: $\Lambda=\text{max}(0,\Sigma^{2}-\nu I)$

10: $P^{-1}=(\Lambda(-1)+\sigma^{2})U(\Lambda+\sigma^{2}I)^{-1}U^{T}+I-UU^{T}$

11: $\log\det(P)=\text{sum(sum(}\log(\frac{1}{\Lambda+\sigma^{2}})))$

Output: $P^{-1}$, $\log\det(P)$ (needed preconditioner quantities)

Algorithm 4 Randomized Gaussian Nystrom Preconditioner

##  Appendix C Proof of the Coercivity condition in Section 3.2

######  Theorem C.1.

Consider $\rho_{\bm{Y}}=\begin{bmatrix}\rho_{\bm{X}}\\\
\rho_{\bm{V}}\end{bmatrix}$, where $\rho_{\bm{Y}}$ is the product of $N$
independent and identical measures with compact support on $\mathbb{R}^{d}$
and $\rho_{\bm{V}}$ is defined in the same way and is independent of
$\rho_{\bm{X}}$. Then we have

(56)
$\displaystyle\|\mathbf{f}_{{\bm{\varphi}}}\|^{2}_{L^{2}(\rho_{\bm{Y}})}\geq\frac{N-1}{N^{2}}\|\varphi^{E}\|^{2}_{L^{2}(\tilde{\rho}_{r}^{E})}+\frac{N-1}{N^{2}}\|\varphi^{A}\|^{2}_{L^{2}(\tilde{\rho}_{r}^{A})}$

######  Proof C.2.

Following the definition of measure $\rho_{\bm{Y}}$ and the norm in (4), we
have

$\displaystyle\|\mathbf{f}_{{\bm{\varphi}}}\|^{2}_{L^{2}(\rho_{\bm{Y}})}$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}\left\|\sum_{i^{\prime}=1}^{N}\frac{1}{N}\Big{[}\varphi^{E}(|{\bm{x}}_{i^{\prime}}-{\bm{x}}_{i}|)({\bm{x}}_{i^{\prime}}-{\bm{x}}_{i})+\varphi^{A}(|{{\bm{x}}_{i^{\prime}}-{\bm{x}}_{i}}|)({\bm{v}}_{i^{\prime}}-{\bm{v}}_{i})\Big{]}\right\|^{2}_{L^{2}(\rho_{\bm{Y}})}$
$\displaystyle=\frac{1}{N^{3}}\sum_{i=1}^{N}\left(\left(\sum_{j=k=1}^{N}+\sum_{j\neq
k=1}^{N}\right)C_{i,j,k}^{E}+C_{i,j,k}^{A}+D_{i,j,k}\right)$ (57)
$\displaystyle=\frac{N-1}{N^{2}}(\|\varphi^{E}\|_{L^{2}(\tilde{\rho}_{r}^{E})}^{2}+\|\varphi^{A}\|_{L^{2}(\tilde{\rho}_{r}^{A})}^{2})+\mathcal{R}$

where

$\displaystyle C_{i,j,k}^{E}$
$\displaystyle=\langle\varphi^{E}(\|{\bm{x}}_{j}-{\bm{x}}_{i}\|)({\bm{x}}_{j}-{\bm{x}}_{i}),\varphi^{E}(\|{\bm{x}}_{k}-{\bm{x}}_{i}\|)({\bm{x}}_{k}-{\bm{x}}_{i})\big{\rangle}_{L^{2}(\rho_{\bm{Y}})},$
$\displaystyle C_{i,j,k}^{A}$
$\displaystyle=\langle\varphi^{A}(\|{\bm{x}}_{j}-{\bm{x}}_{i}\|)({\bm{v}}_{j}-{\bm{v}}_{i}),\varphi^{A}(\|{\bm{x}}_{k}-{\bm{x}}_{i}\|)({\bm{v}}_{k}-{\bm{v}}_{i})\big{\rangle}_{L^{2}(\rho_{\bm{Y}})},$
$\displaystyle D_{i,j,k}$
$\displaystyle=\langle\varphi^{E}(\|{\bm{x}}_{j}-{\bm{x}}_{i})\|({\bm{x}}_{j}-{\bm{x}}_{i}),\varphi^{A}(\|{\bm{x}}_{k}-{\bm{x}}_{i}\|)({\bm{v}}_{k}-{\bm{v}}_{i})\big{\rangle}_{L^{2}(\rho_{\bm{Y}})}$
$\displaystyle\quad\quad+\langle\varphi^{A}(\|{\bm{x}}_{j}-{\bm{x}}_{i}\|)({\bm{v}}_{j}-{\bm{v}}_{i}),\varphi^{E}(\|{\bm{x}}_{k}-{\bm{x}}_{i}\|)({\bm{x}}_{k}-{\bm{x}}_{i})\big{\rangle}_{L^{2}(\rho_{\bm{Y}})}=0,$
$\displaystyle\mathcal{R}$
$\displaystyle=\frac{1}{N^{3}}\sum_{i=1}^{N}\sum_{j\neq k,j\neq i,k\neq
i}(C_{ijk}^{A}+C_{ijk}^{E}).$

By the property of $\rho_{\bm{Y}}$, when $i,j,k$ are distinct, we have

$\displaystyle C_{ijk}^{E}$
$\displaystyle=\mathbb{E}\big{[}\varphi^{E}(\|X_{1}-X_{2}\|)\varphi^{E}(\|X_{1}-X_{3}\|)\left\langle
X_{2}-X_{1},X_{3}-X_{1}\right\rangle\big{]}$ $\displaystyle C_{ijk}^{A}$
$\displaystyle=\mathbb{E}\big{[}\varphi^{A}(\|X_{1}-X_{2}\|)\varphi^{A}(\|X_{1}-X_{3}\|)\big{]}\mathbb{E}\big{[}\left\langle
V_{2}-V_{1},V_{3}-V_{1}\right\rangle\big{]},$

for all $(i,j,k)$, where $X_{i}$s and $V_{i}$s are identical copies of the
position and velocity variables ${\bm{x}}_{i},{\bm{v}}_{i}$s. From the Lemma
C.3 below,

$\displaystyle C_{ijk}^{E}\geq 0,C_{ijk}^{A}\geq 0$

and we used the fact

$\mathbb{E}\big{[}\left\langle
V_{2}-V_{1},V_{3}-V_{1}\right\rangle\big{]}=\mathbb{E}(\|V_{1}\|^{2})-\|\mathbb{E}(V_{1})\|^{2}\geq
0$

Therefore,

$\displaystyle\|\mathbf{f}_{{\bm{\varphi}}}\|^{2}_{L^{2}(\rho_{\bm{Y}})}\geq\frac{N-1}{N^{2}}(\|\varphi^{E}\|_{L^{2}(\tilde{\rho}_{r}^{E})}^{2}+\|\varphi^{A}\|_{L^{2}(\tilde{\rho}_{r}^{A})}^{2})$

The proof of Theorem C.1 uses the following lemma.

######  Lemma C.3.

If $X,Y,Z$ are i.i.d random vectors, then for any measurable function $g$ on
$\mathbb{R}^{d}$, we have that

$\mathbb{E}[g(X-Y)g(X-Z)\langle{X-Y,X-Z}\rangle]\geq 0,$
$\mathbb{E}[g(X-Y)g(X-Z)]\geq 0,$

provided the expectation exists.

######  Proof C.4.

Without loss of generality, suppose the probability density function of $X$ is
$p(x)$. (The discrete distribution case follows from the same argument). Let
$(U,V)=(X-Y,X-Z)$. By the independence of $X,Y,Z$, the pdf of $(U,V)$ is

$p(u,v)=\int p(x)p(x-u)p(x-v)dx.$

Since

$\sum_{i=1}^{N}\sum_{j=1}^{N}c_{i}\bar{c}_{j}p(u_{i},u_{j})=\int
p(x)|\sum_{i=1}^{N}c_{i}p(x-u_{i})|^{2}dx\geq 0,$

which means $p(u,v)$ is positive definite (p.d.) As $\langle{u,v}\rangle$ is
p.d and $g(u)g(v)$ is p.d. [71], we get $g(u)g(v)\langle{u,v}\rangle p(u,v)$
is p.d.. Note that

(58)
$\displaystyle\mathbb{E}[g(X-Y)g(X-Z)\langle{X-Y,X-Z}\rangle]=\int_{\mathbb{R}^{2d}}g(u)g(v)\langle{u,v}\rangle
p(u,v)dudv,$

if the function $g(u)g(v)\langle{u,v}\rangle p(u,v)$ is measurable and
integrable. Then the inequality holds by p.d. property. Similarly, one can
prove the second inequality.

##  Appendix D Proof of Representer Theorem

We prove the Representer Theorem (Theorem 3.1 in main text Section 3.1) by
using an operator-theoretic approach.

######  Proposition D.1.

Given the empirical noisy trajectory data
$(\mathbb{Y}_{M},\mathbb{Z}_{\sigma^{2},M})=\\{\mathbb{X}_{M},\mathbb{V}_{M},\mathbb{Z}_{\sigma^{2},M}\\}$.
We define the sampling operator
$A_{M}:\mathcal{H}_{{K}^{E}}\times\mathcal{H}_{{K}^{A}}\rightarrow\mathbb{R}^{dNML}$
by

(59) $\displaystyle
A_{M}{\bm{\varphi}}=\mathbf{f}_{{\bm{\varphi}}}(\mathbb{X}_{M})$
$\displaystyle:=\mathrm{Vec}(\\{\mathbf{f}_{{\bm{\varphi}}}(\bm{Y}^{(m,l)})\\}_{m,l=1}^{M,L})=\mathrm{Vec}(\\{\mathbf{f}_{\varphi^{E}}(\bm{Y}^{(m,l)})+\mathbf{f}_{\varphi^{A}}(\bm{Y}^{(m,l)})\\}_{m,l=1}^{M,L}),$

where $\mathbb{R}^{dNML}$ is equipped with the inner product defined in (3).

  * 1.

The adjoint operator $A_{M}^{*}$ is a finite rank operator. For any noise
vector $\mathbb{W}$ in $\mathbb{R}^{dNML}$, let
$\mathbb{W}_{m,l,i}\in\mathbb{R}^{d}$ denote the $i$-th component of $(m,l)$th
block of $\mathbb{W}$ as the same way in $\mathbb{Y}_{M}$, then we have

$\displaystyle A^{*}_{M}\mathbb{W}=$
$\displaystyle\bigg{(}\frac{1}{LM}\sum_{l,m=1}^{L,M}\sum_{i=1,i^{\prime}\neq
i}^{N}\frac{1}{N^{2}}K_{r_{ii^{\prime}}^{\bm{X}^{(m,l)}}}^{E}\langle\bm{r}_{ii^{\prime}}^{\bm{X}^{(m,l)}},\mathbb{W}_{m,l,i}\rangle,$
(60)
$\displaystyle\qquad\qquad\frac{1}{LM}\sum_{l,m=1}^{L,M}\sum_{i=1,i^{\prime}\neq
i}^{N}\frac{1}{N^{2}}K_{r_{ii^{\prime}}^{\bm{X}^{(m,l)}}}^{A}\langle\bm{r}_{ii^{\prime}}^{\bm{V}^{(m,l)}},\mathbb{W}_{m,l,i}\rangle\bigg{)}.$

For any function
${\bm{\varphi}}\in\mathcal{H}_{{K}^{E}}\times\mathcal{H}_{{K}^{A}}$, we have
that

(61) $\displaystyle B_{M}{\bm{\varphi}}:=A^{*}_{M}A_{M}{\bm{\varphi}}=$
$\displaystyle\bigg{(}\frac{1}{LM}\sum_{l,m=1}^{L,M}\sum_{i=1,i^{\prime},i^{\prime\prime}\neq
i}^{N}\frac{1}{N^{3}}K_{r_{ii^{\prime}}^{\bm{X}^{(m,l)}}}^{E}(\langle\varphi^{E},K_{r_{ii^{\prime\prime}}^{\bm{X}^{(m,l)}}}^{E}\rangle_{\mathcal{H}_{{K}^{E}}}\langle\bm{r}_{ii^{\prime}}^{\bm{X}^{(m,l)}},\bm{r}_{ii^{\prime\prime}}^{\bm{X}^{(m,l)}}\rangle$
$\displaystyle\hskip
72.26999pt+\langle\varphi^{A},K_{r_{ii^{\prime\prime}}^{\bm{X}^{(m,l)}}}^{A}\rangle_{\mathcal{H}_{{K}^{A}}}\langle\bm{r}_{ii^{\prime}}^{\bm{X}^{(m,l)}},\bm{r}_{ii^{\prime\prime}}^{\bm{V}^{(m,l)}}\rangle),$
(62)
$\displaystyle\frac{1}{LM}\sum_{l,m=1}^{L,M}\sum_{i=1,i^{\prime},i^{\prime\prime}\neq
i}^{N}\frac{1}{N^{3}}K_{r_{ii^{\prime}}^{\bm{X}^{(m,l)}}}^{A}(\langle\varphi^{A},K_{r_{ii^{\prime\prime}}^{\bm{X}^{(m,l)}}}^{A}\rangle_{\mathcal{H}_{{K}^{A}}}\langle\bm{r}_{ii^{\prime}}^{\bm{V}^{(m,l)}},\bm{r}_{ii^{\prime\prime}}^{\bm{V}^{(m,l)}}\rangle$
(63) $\displaystyle\hskip
72.26999pt+\langle\varphi^{E},K_{r_{ii^{\prime\prime}}^{\bm{X}^{(m,l)}}}^{E}\rangle_{\mathcal{H}_{{K}^{E}}}\langle\bm{r}_{ii^{\prime}}^{\bm{V}^{(m,l)}},\bm{r}_{ii^{\prime\prime}}^{\bm{X}^{(m,l)}}\rangle)\bigg{)}.$

  * 2.

If $\bm{\lambda}=(\lambda^{E},\lambda^{A})>0$, a unique minimizer
$\phi_{\mathcal{H}_{{K}^{E}}\times\mathcal{H}_{{K}^{A}}}^{\bm{\lambda},M}$
that solves

$\underset{{\bm{\varphi}}\in\mathcal{H}_{{K}^{E}}\times\mathcal{H}_{{K}^{A}}}{\operatorname{arg}\operatorname{min}}\;\mathcal{E}^{\bm{\lambda},M}({\bm{\varphi}}):=\|A_{M}{\bm{\varphi}}-\mathbb{Z}_{\sigma^{2},M}\|^{2}+\|\sqrt{\bm{\lambda}}\cdot{\bm{\varphi}}\|_{\mathcal{H}_{{K}^{E}}\times\mathcal{H}_{{K}^{A}}}^{2}$

exists and is given by

(64)
$\displaystyle\phi_{\mathcal{H}_{{K}^{E}}\times\mathcal{H}_{{K}^{A}}}^{\bm{\lambda},M}=(B_{M}+\bm{\lambda})^{-1}A_{M}^{*}\mathbb{Z}_{\sigma^{2},M}.$

where we interpret the map $\bm{\lambda}({\bm{\phi}})$ by
$\bm{\lambda}\cdot{\bm{\phi}}=(\lambda^{E}{\phi^{E}},\lambda^{A}{{\phi^{A}}})$.

######  Proof D.2.

The part 1 of Proposition D.1 can be derived by using the identity $\langle
A_{M}{\bm{\varphi}},\bm{w}\rangle=\langle{\bm{\varphi}},A_{M}^{*}\bm{w}\rangle_{\mathcal{H}_{{K}^{E}}\times\mathcal{H}_{{K}^{A}}}$.
Part 2 of Proposition D.1 is straightforward by solving the normal equation.

Now we derive a basis representation formula for the empirical minimizer of
(64)

######  Theorem D.3.

If $\bm{\lambda}>0$, then the minimizer of the regularized empirical risk
functional $\mathcal{E}^{\bm{\lambda},M}(\cdot)$ has the form

(65)
$\phi_{\mathcal{H}_{{K}^{E}}\times\mathcal{H}_{{K}^{A}}}^{\lambda,M}=(\sum_{r^{x}\in
r_{\mathbb{X}_{M}}}\hat{c}_{r^{x}}K_{r^{x}}^{E},\sum_{(r^{x},r^{v})\in(r_{\mathbb{X}_{M}}\times
r_{\mathbb{V}_{M}})}\hat{c}_{r^{v}}K_{r^{x}}^{A}),$

where $r_{\mathbb{X}_{M}}\in\mathbb{R}^{MLN^{2}}$ is the set contains all the
pair distances in $\mathbb{X}_{M}$, i.e.

(66)
$r_{\mathbb{X}_{M}}=\begin{bmatrix}r_{11}^{(1,1)},\dots,r_{1N}^{(1,1)},\dots,r_{N1}^{(1,1)},\dots,r_{NN}^{(1,1)},\dots,r_{11}^{(M,L)},\dots,r_{1N}^{(M,L)},\dots,r_{N1}^{(M,L)},\dots,r_{NN}^{(M,L)}\end{bmatrix}^{T},$

and $r_{\mathbb{X}_{M}}\times r_{\mathbb{V}_{M}}\in\mathbb{R}^{MLN^{2}\times
MLN^{2}}$ is the set contains all the pair distances in $\mathbb{X}_{M}$ and
their associated pair distances in $\mathbb{V}_{M}$.
Moreover, we have

$\displaystyle\hat{c}_{r^{x}}=\frac{1}{N}\bm{r}_{\mathbb{X}_{M}}^{T}\cdot(K_{\mathbf{f}_{\bm{\phi}}}(\mathbb{Y}_{M},\mathbb{Y}_{M})+\lambda^{E}NMLI)^{-1}\mathbb{Z}_{\sigma^{2},M},$
(67)
$\displaystyle\hat{c}_{r^{v}}=\frac{1}{N}\bm{r}_{\mathbb{V}_{M}}^{T}\cdot(K_{\mathbf{f}_{\bm{\phi}}}(\mathbb{Y}_{M},\mathbb{Y}_{M})+\lambda^{A}NMLI)^{-1}\mathbb{Z}_{\sigma^{2},M},$

where the block-diagonal matrix
$\bm{r}_{\mathbb{X}_{M}}=\mathrm{diag}(\bm{r}_{\bm{X}^{(m,l)}})\in\mathbb{R}^{MLdN\times
MLN^{2}}$ and $\bm{r}_{\bm{X}^{(m,l)}}\in\mathbb{R}^{dN\times N^{2}}$ defined
by

(68)
$\bm{r}_{\bm{X}^{(m,l)}}=\begin{bmatrix}\bm{r}_{11}^{(m,l)},\dots,\bm{r}_{1N}^{(m,l)}&\bm{0}&\cdots&\bm{0}\\\
\bm{0}&\bm{r}_{21}^{(m,l)},\dots,\bm{r}_{2N}^{(m,l)}&\cdots&\bm{0}\\\
\vdots&\vdots&\ddots&\vdots\\\
\bm{0}&\bm{0}&\cdots&\bm{r}_{N1}^{(m,l)},\dots,\bm{r}_{NN}^{(m,l)}\end{bmatrix}\
,$

and same for $\bm{r}_{\mathbb{V}_{M}}$.

######  Proof D.4.

Let $\mathcal{H}_{K^{E},M}$ be the subspace of $\mathcal{H}_{{K}^{E}}$ spanned
by the set of functions $\\{K_{r}^{E}:r\in r_{\mathbb{X}_{M}}\\}$, and
similarly for $\mathcal{H}_{K^{A},M}$. By Proposition Proposition D.1, we know
that
$B_{M}(\mathcal{H}_{K,M}^{E}\times\mathcal{H}_{K,M}^{A})\subset\mathcal{H}_{K,M}^{E}\times\mathcal{H}_{K,M}^{A}$.
Since $B_{M}$ is self-adjoint and compact, by the spectral theory of self-
adjoint compact operator (see [72]),
$\mathcal{H}_{K,M}^{E}\times\mathcal{H}_{K,M}^{A}$ is also an invariant
subspace for the operator $(B_{M}+\lambda I)^{-1}$. Then by (64), there exists
vectors $\hat{c}_{r^{x}}$, $\hat{c}_{r^{v}}$ such that

(69) $\phi_{\mathcal{H}_{K^{E}\times K^{A}}}^{\lambda,M}=(\sum_{r^{x}\in
r_{\mathbb{X}_{M}}}\hat{c}_{r^{x}}K_{r^{x}}^{E},\sum_{(r^{x},r^{v})\in(r_{\mathbb{X}_{M}}\times
r_{\mathbb{X}_{M}})}\hat{c}_{r^{v}}K_{r^{x}}^{A}).$

Then, multiplying $(B_{M}+\bm{\lambda})$ on both sides of (64) and plugging in
(69), we can obtain

$\begin{cases}\big{(}\bm{r}_{\mathbb{X}_{M}}^{T}\bm{r}_{\mathbb{X}_{M}}K^{E}(r_{\mathbb{X}_{M}},{r_{\mathbb{X}_{M}}})+\lambda^{E}N^{3}MLI\big{)}\hat{c}_{r^{x}}+\bm{r}_{\mathbb{X}_{M}}^{T}\bm{r}_{\mathbb{V}_{M}}K^{A}(r_{\mathbb{X}_{M}},{r_{\mathbb{X}_{M}}})\hat{c}_{r^{v}}&=N\bm{r}_{\mathbb{X}_{M}}^{T}\mathbb{Z}_{\sigma^{2},M}\\\
\big{(}\bm{r}_{\mathbb{V}_{M}}^{T}\bm{r}_{\mathbb{V}_{M}}K^{A}(r_{\mathbb{X}_{M}},{r_{\mathbb{X}_{M}}})+\lambda^{A}N^{3}MLI\big{)}\hat{c}_{r^{v}}+\bm{r}_{\mathbb{V}_{M}}^{T}\bm{r}_{\mathbb{X}_{M}}K^{E}(r_{\mathbb{X}_{M}},{r_{\mathbb{X}_{M}}})\hat{c}_{r^{x}}&=N\bm{r}_{\mathbb{V}_{M}}^{T}\mathbb{Z}_{\sigma^{2},M}\end{cases}$

using the matrix representation of $(B_{M}+\bm{\lambda})$ with respect to the
spanning sets $\\{K_{r}^{E}:r\in r_{\mathbb{X}_{M}}\\}$ and $\\{K_{r}^{A}:r\in
r_{\mathbb{X}_{M}}\\}$.

Recall that we have
$K^{E}(r_{\mathbb{X}_{M}},{r_{\mathbb{X}_{M}}})=(K^{E}(r_{ij},r_{i^{\prime}j^{\prime}}))_{r_{ij},r_{i^{\prime}j^{\prime}}\in
r_{\mathbb{X}_{M}}}$, $K^{A}(r_{\mathbb{X}_{M}},{r_{\mathbb{X}_{M}}})=$
$(K^{A}(r_{ij},r_{i^{\prime}j^{\prime}}))_{r_{ij},r_{i^{\prime}j^{\prime}}\in
r_{\mathbb{X}_{M}}}$ and
$K_{\mathbf{f}_{\bm{\phi}}}(\mathbb{Y}_{M},\mathbb{Y}_{M})=\mathrm{Cov}(\mathbf{f}_{{\bm{\phi}}}(\mathbb{Y}_{M}),\mathbf{f}_{{\bm{\phi}}}(\mathbb{Y}_{M}))$,
so using the identity

(70)
$\displaystyle\bm{r}_{\mathbb{X}_{M}}K^{E}(r_{\mathbb{X}_{M}},{r_{\mathbb{X}_{M}}})\bm{r}_{\mathbb{X}_{M}}^{T}+\bm{r}_{\mathbb{V}_{M}}K^{A}(r_{\mathbb{X}_{M}},{r_{\mathbb{X}_{M}}})\bm{r}_{\mathbb{V}_{M}}^{T}=N^{2}K_{\mathbf{f}_{\bm{\phi}}}(\mathbb{Y}_{M},\mathbb{Y}_{M})$

and the fact that the matrices
$\big{(}\bm{r}_{\mathbb{X}_{M}}^{T}\bm{r}_{\mathbb{X}_{M}}K^{E}(r_{\mathbb{X}_{M}},{r_{\mathbb{X}_{M}}})+\lambda^{E}N^{3}MLI\big{)}$,
$\big{(}\bm{r}_{\mathbb{V}_{M}}^{T}\bm{r}_{\mathbb{V}_{M}}K^{A}(r_{\mathbb{X}_{M}},{r_{\mathbb{X}_{M}}})+\lambda^{A}N^{3}MLI\big{)}$
are invertible, one can verify that

(71)
$\displaystyle\begin{cases}\hat{c}_{r^{x}}&=\frac{1}{N}\bm{r}_{\mathbb{X}_{M}}^{T}\cdot(K_{\mathbf{f}_{\bm{\phi}}}(\mathbb{Y}_{M},\mathbb{Y}_{M})+\lambda^{E}NMLI)^{-1}\mathbb{Z}_{\sigma^{2},M},\\\
\hat{c}_{r^{v}}&=\frac{1}{N}\bm{r}_{\mathbb{V}_{M}}^{T}\cdot(K_{\mathbf{f}_{\bm{\phi}}}(\mathbb{Y}_{M},\mathbb{Y}_{M})+\lambda^{A}NMLI)^{-1}\mathbb{Z}_{\sigma^{2},M},\end{cases}$

is the solution.

Now we are ready to finish the proof of the Representer theorem.

######  Proof D.5.

Let $\tilde{K}^{E}=\frac{\sigma^{2}K^{E}}{MNL\lambda^{E}}$,
$\tilde{K}^{A}=\frac{\sigma^{2}K^{A}}{MNL\lambda^{A}}$.

Since ${\phi^{E}}\sim\mathcal{GP}(0,\tilde{K}^{E})$,
${\phi^{A}}\sim\mathcal{GP}(0,\tilde{K}^{A})$, the posterior mean in
$\eqref{eq:estimated phi}$ will then become

$\displaystyle\bar{\phi}_{M}^{E}(r^{\ast})$
$\displaystyle=\tilde{K}_{{\phi^{E}},\mathbf{f}_{{\bm{\phi}}}}(r^{\ast},\mathbb{X}_{M})(\tilde{K}_{\mathbf{f}_{\bm{\phi}}}(\mathbb{Y}_{M},\mathbb{Y}_{M})+\sigma^{2}I)^{-1}\mathbb{Z}_{\sigma^{2},M}$
$\displaystyle=\frac{1}{N}\tilde{K}_{r_{\mathbb{X}_{M}}^{T}}^{E}(r^{\ast})\bm{r}_{\mathbb{X}_{M}}^{T}(\tilde{K}_{\mathbf{f}_{\bm{\phi}}}(\mathbb{Y}_{M},\mathbb{Y}_{M})+\sigma^{2}I)^{-1}\mathbb{Z}_{\sigma^{2},M}$
$\displaystyle=\frac{1}{N}K_{r_{\mathbb{X}_{M}}^{T}}^{E}(r^{\ast})\bm{r}_{\mathbb{X}_{M}}^{T}(K_{\mathbf{f}_{\bm{\phi}}}(\mathbb{Y}_{M},\mathbb{Y}_{M})+NML\lambda^{E}I)^{-1}\mathbb{Z}_{\sigma^{2},M}$
$\displaystyle=K_{{\phi^{E}},\mathbf{f}_{{\bm{\phi}}}}(r^{\ast},\mathbb{X}_{M})(K_{\mathbf{f}_{\bm{\phi}}}(\mathbb{Y}_{M},\mathbb{Y}_{M})+NML\lambda^{E}I)^{-1}\mathbb{Z}_{\sigma^{2},M}$
$\displaystyle=\sum_{r\in r_{\mathbb{X}_{M}}}\hat{c}_{r}K_{r}^{E},$

where $\hat{c}_{r}$ is defined in (D.3) and we used the identity
$K_{{\phi^{E}},\mathbf{f}_{{\bm{\phi}}}}(r^{\ast},\mathbb{X}_{M})=\frac{1}{N}K_{r_{\mathbb{X}_{M}}^{T}}^{E}(r^{\ast})\bm{r}_{\mathbb{X}_{M}}^{T}$
(also for $\tilde{K}$) in the proof. Similarly, we can get the posterior mean
for $\bar{\phi}_{M}^{A}(r^{\ast})$ .

######  Lemma D.6.

Let ${\bm{x}}$ and $\bm{y}$ be jointly Gaussian random vectors

(72) $\begin{bmatrix}{\bm{x}}\\\
\bm{y}\end{bmatrix}\sim\mathcal{N}(\begin{bmatrix}\mu_{{\bm{x}}}\\\
\mu_{\bm{y}}\end{bmatrix},\begin{bmatrix}A&C\\\ C^{T}&B\end{bmatrix}),$

then the marginal distribution of ${\bm{x}}$ and the conditional distribution
of ${\bm{x}}$ given $\bm{y}$ are

(73) ${\bm{x}}\sim\mathcal{N}(\mu_{{\bm{x}}},A),\quad\textrm{and
}{\bm{x}}|\bm{y}\sim\mathcal{N}(\mu_{{\bm{x}}}+CB^{-1}(\bm{y}-\mu_{\bm{y}}),A-CB^{-1}C^{T}).$

######  Proof D.7.

See, e.g. [50], Appendix A.

## References

  * [1] 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.