text:
\section{Introduction} \label{sec:intro}
Physical, biological, and social systems across all scales of complexity and size can often be described as dynamical systems written in terms of interacting agents (e.g. particles, cells, humans, planets, ...). Rich theories have been developed  to explain the collective behavior of these interacting agents across many fields including astronomy, particle physics, economics, social science, and biology. Examples include predator-prey systems, molecular dynamics, coupled harmonic oscillators, flocking birds or milling fish, human social interactions, and celestial mechanics, to name a few. In order to encompass many of these examples, we will describe a very general second-order, heterogeneous (the agents can be of different types), interacting (the acceleration of an agent is a function of properties of the other agents) agent system that includes external forces, masses of the agents, multivariable interaction kernels, and an additional environment variable that is a dynamical property of the agent (for example, a firefly having its luminescence varying in time). We propose a learning approach that combines machine learning and dynamical systems in order to provide highly accurate dynamical models of the observation data from these systems.

The model and learning framework presented in sections \ref{sec:ModelDesc}-\ref{sec:LT} includes a very large number of relevant systems and allows for their modeling.
Clustering of opinions \cite{Krause2000,CKFL2005, BHT2009, MT2014} is a simple first-order case that exhibits clustering. Flocking of birds \cite{CS2007,CM2008,CD2011} provides a simple example of a second-order system that exhibits an emergent shared velocity of all agents. Milling of fish \cite{Chuang2007,Abaid2010,Albi2014,Chuang2016} is another second-order model and presents both a $2$ and $3$-dimensional milling pattern over long time and introduces a non-collective force from the environment.
A model of oscillators (fireflies) that sync and swarm together, and have their dynamics governed by their positions and a phase variable $\xi$, was studied by \cite{Strogatz2000,OKeeffe2017,OKeeffe2018,Okeeffe2019}.
There are also models that include both energy and alignment interaction kernels, a particular case of this is the anticipation dynamics model from \cite{shu2019anticipation}, which we study in this work.

One can also consider a collection of celestial bodies interacting via the gravitational potential, which was initially studied in \cite{Zhong20} and is further studied in the upcoming \cite{GravLetter20}. All of these models fit into our framework and we have presented detailed studies of these, and others in this work as well as \cite{Zhong20, Tang2019, lu2019nonparametric}.
These dynamics exhibit a wide range of emergent behaviors, and
as shown in \cite{Vicsek_model,TKIHLC2013,CM2008,GC2004,
CHDOB2007,BDT2017,MT2014}, the behaviors can be studied when the governing equations are known. However, if the equations are not known and the data consists of only trajectories, we still wish to develop a model that can both make accurate predictions of the trajectories and discover a dynamical form that accurately reflects their emergent properties. To achieve this, we present a theoretically optimal learning algorithm that is accurate, captures emergent behavior for large time, and, by exploiting the structure of the collective dynamical system, avoids the curse of dimensionality.

Applying machine learning to the sciences has experienced tremendous growth in recent years, a small selection of general applications related to the ideas in this work include: learning PDEs (\cite{Bar-Sinai,Schaeffer6634, Kevri_2}), modeling dynamical systems (\cite{Kevrik_1, Kevri_3, Kevri_4}), governing equations (\cite{Champion, Kevri_5}), biology (\cite{Chiel993}), fluid mechanics (\cite{Raissi1026, Han219}), many-body problems in quantum systems (\cite{Carleo602}), mean-field games (\cite{Ruthotto}), meteorology (\cite{Ham2019}), and dynamical systems (\cite{Costa1501, Brunton3932, YairE, Bongard9943}). These, and the references therein, give a flavor of the diverse range of applications.
A vast literature exists in the context of learning dynamical systems. In the case of a general nonlinear dynamical system, symbolic regression has been developed to learn the underlying form of the equations from data, see \cite{bongard2007automated, schmidt2009distilling}. Sparse regression techniques which use an extremely large collection of functions, often containing most major mathematical functions, are fit to the data with a sparsity condition that only allows a few terms to appear in the final model. Detailed study and development of these approaches can be found for SINDy in (\cite{BPK2016, RBPK2017, BKP2017}), a LASSO-type penalty (\cite{han2015robust, kang2019ident}), and sparse Bayesian regression (\cite{zhang2018robust}). Other approaches consider multiscale methods, statistical mechanics, or force-based models, see \cite{BCCCCGLOPPVZ2008,BCGMSVW2012}.
Deep learning has also been applied to learn dynamical systems, for ODEs see \cite{raissi2018multistep, rudy2019deep} and for PDEs see \cite{raissi2018deep, raissi2018hidden, long2017pde}, as well as the references therein.

The majority of the earliest work in inferring interaction kernels in systems of the type \eqref{eq:intro:second_order}, \eqref{eq:2ndOrder} occurred in the Physics literature, going back to the works of Newton. From the viewpoint of purely data-driven analysis of the equations, requiring limited or no physical reasoning, foundational work was \cite{LLEK2010,KTIHC2011}. In these works, the interaction kernels are assumed to be in the span of a known family of functions and parameters are estimated.
 In statistics, the problem of parameter estimation in dynamical systems from observations is classical, e.g. \cite{varah1982spline,brunel2008parameter,
liang2008parameter,cao2011robust,ramsay2007parameter}. The question of identifiability of the parameter emerges, see e.g. \cite{dattner2015optimal,miao2011identifiability}. Our work is closely related to this viewpoint but our parameter is now infinite-dimensional, with identifiability discusses in section \ref{sec:coercivity}.

Our learning approach is based on exploiting the structure of collective dynamical systems and nonparametric estimation techniques (see \cite{cucker2002mathematical, Tsybakov:2008, Gyorfi06, devore2006approximation, binev2005universal}).
We focus here on second-order models and the form of the equations, generalizing the first order models (see discussion in Appendix \ref{sec:1stOrder}), derived from Newton's second law: for $i = 1, \ldots, N$
\begin{align}\label{eq:intro:second_order}
\hspace{-2cm} m_i\ddot\mbf{x}_i(t) = \force^{\dot\bx}(\mbf{x}_i, \dot\mbf{x}_i) + \frac{1}{N} \sum_{\substack{i' = 1,\\ i' \neq i}}^N &\phi^E(\norm{\mbf{x}_{i'}(t) - \mbf{x}_i(t)})(\mbf{x}_{i'}(t) - \mbf{x}_i(t)) \nonumber \\ &+ \phi^A(\norm{\mbf{x}_{i'}(t) - \mbf{x}_i(t)})(\dot{\mbf{x}}_{i'}(t) - \dot{\mbf{x}}_i(t)) .
\end{align}
Here, $m_i$ is the mass of the $i^{th}$ agent, $\mbf{x}_i$ is its position, $\force^{\dot\bx}$ is a non-collective force, and $\phi^E: \real^+ \rightarrow \real$, $\phi^A: \real^+ \rightarrow \real$ are known as the interaction kernels.  A significant amount of research on modeling collective dynamics is concerned with inducing desired collective behaviors (flocking, clustering, milling, etc.) from relatively simple, local and non-local interaction kernels, often from known or specific parametric families of relatively simple functions. Here we consider a non-parametric, inverse-problem-based approach to infer the interaction kernels from observations of trajectory data, especially within short-time periods.  In \cite{BFHM17}, a convergence study of learning unknown interaction kernels from observation of first-order models of homogeneous agents was done for increasing $N$, the number of agents.
The estimation problem with $N$ fixed, but the number of trajectories $M$ varying, for first-order and second-order models of heterogeneous agents was numerically studied in \cite{lu2019nonparametric} and learning theory on these first-order models was developed in \cite{Tang2019, li2019identifiability}.  Further extension of the model and algorithm to more complicated second-order systems, with particular emphasis on emergent collective behaviors, was discussed in \cite{Zhong20}.  A big data application to real celestial motion ephemerides is developed and discussed in \cite{GravLetter20}.  In this work, we provide a rigorous learning theory covering the models presented in \cite{Zhong20}, as well as the second-order models introduced in \cite{lu2019nonparametric}.  We consider generalizations of the models in \cite{Zhong20}, to include models with higher-dimensional interaction kernels, that do not depend only on pairwise distances. Compared to the theories studied in \cite{Tang2019, li2019identifiability}, our theory focuses on second-order models with interaction kernels of the form $\phi^{E}(r)r + \phi^{A}(r)\dot{r}$ (with $r$ and $\dot r$ representing norms of differences of positions and, respectively, velocities of pairs of agents); additionally, we discuss the identifiability and separability of $\phi^{E}$ and $\phi^{A}$ from the sum.

The overall objective of the algorithm can be stated as: given trajectory data generated from an interacting agent system, we wish to learn the underlying interaction kernels, from which we will understand its long-term and emergent behavior, and ultimately build a highly accurate approximate model that faithfully captures the dynamics. We make minimal assumptions on the form of the interaction kernels and the various forces involved, and in some cases the assumptions are made for purely theoretical reasons and the algorithm can still perform well when they do not hold for a given system.
We offer a learning approach to address these collective systems by first discovering the governing equations from the observational data, and then using the estimated equations for large-time predictions.

The approach in this paper builds on many of these ideas and uses observation data coming from collective dynamical systems of the form \eqref{eq:2ndOrder}, to learn the underlying interaction functions. This variational approach was initially developed in \cite{BFHM17, lu2019nonparametric} and further studied and extended in \cite{Tang2019,Zhong20}.
Our analysis of this system blends differential equations, inverse problems and nonparametric regression, and (statistical) learning theory. A central insight is that we exploit the form of \eqref{eq:2ndOrder} to move the inference task to just the unknown functions ($\phi^E,\phi^A,\phi^{\xi}$) allowing us to avoid the curse of dimensionality incurred if we were to directly perform regression against the high-dimensional phase space and trajectory data of the system, which provides independence of observations where each observation is a different trajectory.

To use the trajectory data to derive estimators, we consider appropriate hypothesis spaces in which to build our estimators, measures adapted to the dynamics, norms, and other performance metrics, and ultimately an inverse problem built from these tools.
Once we have obtained this estimated interaction kernel, we want to study its properties as a function of the amount of trajectory data we receive, which is the $M$ trajectories sampled from different initial conditions from the same underlying system. Here we study properties of the error functional, establish the uniqueness of its minimizers, and use the probability measures to define a dynamics-adapted norm to measure the error of our estimators over the hypothesis spaces. In comparing the estimators to the true interaction kernels, we first establish concentration estimates over the hypothesis space.

Our first main result is the strong consistency of our learned estimators asymptotically. For the relevant definitions see \ref{sec:learn_framework} and for the full theorem see section \ref{sec:consistency}, which for the model \eqref{eq:intro:second_order} yields,
\begin{align}
\lim_{M\rightarrow \infty}\|\widehat{\bm{\phi}}^{EA} - \bm{\phi}^{EA}\|_{\LtwoB(\bm{\rho}_T^{EA,L})} =0 \text{ with probability one.}
\end{align}
where $\bm{\rho}_T^{EA,L}$ is a dynamics-adapted measure on pairwise distances -- and we use a weighted $\LtwoB$ space detailed in section \ref{sec:LT}, particularly \eqref{meas:vectorized}.
Perhaps most importantly, we give a rate of convergence in terms of the $M$ trajectories. We achieve the minimax rate of convergence for any number of variables $|\bm{\mathcal{V}}|$ in the interaction kernels. See section \ref{sec:rateofconv} for the full theorem, (see section \ref{sec:learn_framework} for relevant definitions) which is given by:
\begin{align}
 \mathbb{E}_{\bm{\mu^{\bY}}}\Big[\| \widehat{\bm{\phi}}^{EA}
 -\bm{\phi}^{EA} \|^2_{\LtwoB(\bm{\rho}_T^{EA,L})} \Big] \leq C \left(\frac{\log M}{M}\right)^{\frac{2s}{2s+|\bm{\mathcal{V}}|}}\,.
 \end{align}
In the case of model \eqref{eq:intro:second_order}, $|\bm{\mathcal{V}}|=1$, as in the results for first-order systems \cite{BFHM17,lu2019nonparametric,Tang2019}.

This means that our estimators converge at the same rate in $M$ as the best possible estimator (up to a logarithmic factor) one could construct when the initial conditions are randomly sampled from some underlying initial condition distribution denoted $\bm{\mu^{\bY}}$ throughout this work, see (section \ref{sec:main_results}).

To solve the inverse problem, we give a detailed discussion of an essential link between these three aspects, the notion of coercivity of the system - detailed in section \ref{sec:coercivity}. Coercivity plays a key role in the approximation properties, the ability to approximate the interaction kernel and the learning theory. We also present numerical examples, see the detailed numerical study in \cite{Zhong20}, which help to explain why the particular norms we define are the right choice, as well as show excellent performance on complex dynamical systems in section \ref{sec:numerics}.

Our paper is structured as follows. The first part of the paper describes the model, learning framework, inference problem, and the basic tools needed for the learning theory. These ideas are all explained in detail in sections \ref{sec:ModelDesc}-\ref{sec:LT}. If one wishes to quickly jump to the theoretical sections, and then refer back to the definitions as needed, we have provided table \ref{tab:2ndOrder_def},\ref{tab:2ndOrderlearning_def} which explains the model equations and outlines the definitions and concepts needed for the learning theory and general theoretical results, respectively.
The theoretical part of the paper (sections \ref{sec:coercivity}-\ref{sec:trajectory}) discusses fundamental questions of identifiability and solvability of the inverse problem, consistency, and rate of convergence of the estimators, and the ability to control trajectory error of the evolved trajectories using our estimators. Some key highlights of our theoretical contributions are described in \ref{sec:theory_overview}, with full details in the corresponding sections.
Lastly, we consider applications in section \ref{sec:numerics}, as well as have many additional proofs and details in appendices \ref{s:trajectoryerrorproof}-\ref{sec:algorithm_numeric}.

\subsection{Comparison with existing work}
Our learning approach discovers the governing structure of a particular subset of dynamical systems of the form,
\[
\dot\bY(t) = \mbf{F}_{\bm{\phi}^{EA}, \bm{\phi}^{\xi}}(\bY(t)), \quad t \in [0, T].
\]
from observations $\{\bY(t_l), \dot\bY(t_l)\}_{l = 1}^L$, by implicitly inferring the right hand side, $\mbf{F}_{\bm{\phi}^{EA}, \bm{\phi}^{\xi}}$. The main difficulties in establishing an effective theory of learning $\mbf{F}_{\bm{\phi}^{EA}, \bm{\phi}^{\xi}}$ are the \textit{curse of dimensionality} caused by the size of $\bY$, which is often $N(2d + 1)$, where $N$ is the number of agents, $d$ the dimension of physical space; and the \textit{dependence} of the observation data, for example $\bY(t_{l + 1})$ depends on $\bY(t_{l})$.

There are many techniques which can be used to tackle the high-dimension of the data set: sparsity assumptions, dimension reduction, reduced-order modeling, and machine learning techniques trained using gradient-based optimization. The dependent nature of the data prevents traditional regression-based approaches, see the discussion in \cite{Tang2019}, but many of the approaches above successfully address this. Our work, however, exploits the interacting-agent structure of collective dynamical systems, which is driven by a collection of two-body interactions where each interaction depends only on pairwise data between the states of agents, as in \eqref{eq:intro:second_order}. With this structure in mind, we are able to reduce the ambient dimension of the data $N(2d + 1)$ to the dimension of the variables in the interaction kernels, which is independent of $N$. We also naturally incorporate the dependence in the data in an appropriate manner by considering trajectories generated from different initial conditions.

Our theoretical results focus on the joint learning of ${\bm{\phi}}^{E}, {\bm{\phi}}^{A}$ that takes into account their natural weighted direct sum structure that is described in the following sections, which is different from the learning theory on single ${\bm{\phi}}^{E}$'s considered in \cite{lu2019nonparametric, Tang2019}.  The current theoretical framework is not able to conclusively show that ${\bm{\phi}}^{E}$ and ${\bm{\phi}}^{A}$ can be learned separately; however we demonstrate in various numerical experiments that by learning ${\bm{\phi}}^{E}$ and ${\bm{\phi}}^{A}$ jointly, we still achieve strong performance. Finally, we note that the first-order theory developed in \cite{Tang2019} is a special case of our second-order theory, see details in appendix \ref{sec:1stOrder}.

\section{Model description} \label{sec:ModelDesc}
In order to motivate the choice of second-order models considered in this paper, we begin our discussion with a simple second-order model derived from the classical mechanics point of view.  Let us consider a closed system of $N$ homogeneous agents (or particles) equipped with a certain type of Lagrangian energy, namely $L(t)$ for the whole system, given as follows,
\[
L(t) = \frac{1}{2}m_i\norm{\dot\mbf{x}_i(t)}^2 - \frac{1}{N}\sum_{1 \le i < i' \le N}U(\norm{\mbf{x}_{i'}(t) - \mbf{x}_i(t)}), \quad i = 1, \ldots, N.
\]
Here $U$ is a potential energy depending on pairwise distance.  From the Lagrange equation, $\frac{d}{dt}\Big(\partial_{\dot\mbf{x}_i}L\Big) = \partial_{\mbf{x}_i}L$, we obtain a simple second-order collective dynamics model
\begin{equation}\label{eq:2nd_simple}
m_i\ddot\mbf{x}_i(t) = \frac{1}{N}\sum_{i' = 1, i' \neq i}^N\phi^{E}(\norm{\mbf{x}_{i'}(t) - \mbf{x}_i(t)})(\mbf{x}_{i'}(t) - \mbf{x}_i(t)), \quad i = 1, \ldots, N.
\end{equation}
Here, $\phi^{E}(r) = \frac{U'(r)}{r}$ represents an energy-based interaction between agents.  For example, taking $U(r) = \frac{NGm_{i'}m_i}{r}$, it becomes the celebrated model for Newton's universal gravity.
In order to incorporate more complicated behaviors into the model equation, we consider alignment-based interactions (to align velocities, so that short-range repulsion, mid-range alignment, and long arrange attraction are all present), auxiliary variables describing internal states of agents (emotion, excitation, phases, etc.), and non-collective forces (interaction with the environment).  We also consider a system of $N$ heterogeneous agents, such that the agents belong to $K$ disjoint types $\{C_{k}\}_{k = 1}^{K}$, with $N_{k}$ being the number of agents of type $k$. They interact according to the system of ODEs
\begin{equation}
\begin{dcases}
m_i\ddot\mbf{x}_i(t) &= \force^{\dot\bx}(\mbf{x}_i(t), \dot\mbf{x}_i(t), \xi_i(t)) + \sum_{i'=1}^N \frac{1}{N_{\mathpzc{k}_{i'}}} \Big[\phi^{E}_{\mathpzc{k}_i \mathpzc{k}_{i'}}(r_{ii'}(t), \topwd^E_{i i'}(t))(\mbf{x}_{i'}(t) - \mbf{x}_i(t)) \\
&\qquad\qquad\qquad\qquad\qquad\qquad\qquad + \phi^{A}_{\mathpzc{k}_i \mathpzc{k}_{i'}}(r_{ii'}(t), \topwd^A_{i i'}(t))(\dot\mbf{x}_{i'}(t) - \dot\mbf{x}_i(t))\Big]  \\
\dot\xi_i(t) &= \force^{\xi}(\mbf{x}_i(t), \dot\mbf{x}_i(t), \xi_i(t)) + \sum_{i'=1}^N \frac{1}{N_{\mathpzc{k}_{i'}}} \phi^{\xi}_{\mathpzc{k}_i \mathpzc{k}_{i'}}(r_{ii'}(t), \topwd^{\xi}_{i i'}(t))(\xi_{i'}(t) - \xi_i(t))
\end{dcases}
\label{eq:2ndOrder}
\end{equation}
for $i = 1, \ldots, N$, where $\mathpzc{k}_i, \mathpzc{k}_{i'} \in \{1,\cdots, K\}$ are the indices of the agent type of the agents $i$ and $i'$ respectively; the interaction kernels, $\phi^{E}_{k k'}, \phi^{A}_{k\idxcl'}, \phi^{\xi}_{k\idxcl'}$, are in general different for interacting agents of different types, and they not only depend on pairwise distance $r_{ii'}(t)$, but also on  other features, given by $\topwd^E_{i i'}(t), \topwd^A_{i i'}(t), \topwd^{\xi}_{i i'}(t)$.  For example, the interactions between birds can also depend on the field of vision, not just the distance between pairs of birds. Note that we will often suppress the explicit dependence on time $t$ when it is clear from context.
The unknowns, for which we will construct estimators, in these equations, are the functions $\phi^{E}_{\mathpzc{k}_i \mathpzc{k}_{i'}}$, $\phi^{A}_{\mathpzc{k}_i \mathpzc{k}_{i'}}$ and $\phi^{\xi}_{\mathpzc{k}_i \mathpzc{k}_{i'}}$; everything else is assumed given.

Table \ref{tab:2ndOrder_def} gives a detailed explanation for the definition of the variables used in \eqref{eq:2ndOrder}.
We note that in what follows the notation, ${\{E,A,\xi\}}$, attached to any expression means that there is one of those maps, functions, etc. for each element in the set ${\{E,A,\xi\}}$. It is a convenient way to avoid excessive repetition of similar definitions.
\begin{table}[H]
\centering
\small{
\small{\begin{tabular}{ c | l }
\hline
Variable                    & Definition \\
\hline\hline
$i,i'$                     & index of agent, from $1,\ldots, N$ \\
\hline
$m_i$                       & mass of agent $i$ \\
\hline
$\mbf{x}_i(t),\dot{\mbf{x}}_i(t),\ddot{\mbf{x}}_i(t)\in \real^d$          & position/velocity/acceleration vector of agent $i$ at time $t$ \\
\hline
$\xi_i, \dot{\xi}_i$					& auxiliary variable, and its derivative  \\
\hline
$\|\cdot\|$ & Euclidean norm in $\mathbb{R}^d$ \\
\hline
$K$                    & number of types of agents \\
\hline
$k, k'$            & index of type of agents, from $1,\ldots, K$  \\
\hline
$N_{k}$                & number of agents in type $k$ \\
\hline
$\mathpzc{k}_i$                   & type index of agent $i$ \\
\hline
$C_{k}$              &  set of indices for type $k$ agent, subset of $\{1, \ldots, N\}$   \\
\hline
$\phi_{kk'} $ (or $\phi_{kk''}$) & interaction kernel: influence of agent $k'$ (or $k''$) on agent $k$ \\
\hline
$\force^{\dot\bx}, \force^{\xi}$         & non-collective forces on $\ddot\mbf{x}_i$ and $\dot\xi_i$ respectively \\
\hline
$\phi^{E}, \phi^{A}, \phi^{\xi}$ & energy, alignment, and environment-based interaction kernels respectively \\
\hline
\iffalse
$M$							& number of trajectories \\
\hline
$L$							& number of observations in time \\
\hline
$n$							& number of basis functions used \\
\fi
\hline $\topwdT $      & Features, $\topwdT(\mbf{x}, \dot{\mbf{x}}, \xi, \mbf{x} ', \dot{\mbf{x}}', \xi '): \real^{4d + 2} \rightarrow \mathbb{R}^{p}$ \\

\hline
$\mbf{s}_{(k, k')}^{{\{E,A,\xi\}}}$            & Feature map, $\pi_{k\idxcl'}^{{\{E,A,\xi\}}}\circ \topwdT(\mbf{x}, \dot\mbf{x}, \xi, \mbf{x}', \dot\mbf{x}', \xi'): \real^{4d + 2} \rightarrow \real^{p_{(k,k')}^{{\{E,A,\xi\}}}}$ \\

\hline
$\mbf{s}_{ii'}^{{\{E,A,\xi\}}}$            & Feature evaluation,
$\mbf{s}_{(\mathpzc{k}_i,\mathpzc{k}_{i'})}^{{\{E,A,\xi\}}}(\mbf{x}_i, \dot\mbf{x}_i, \xi_i, \mbf{x}_{i'}, \dot\mbf{x}_{i'}, \xi_{i'}) \in \real^{p_{\mathpzc{k}_{i}\mathpzc{k}_{i'}}^{{\{E,A,\xi\}}}}$ \\

\hline
\end{tabular}}
}
\caption{Notation for the variables in \eqref{eq:2ndOrder}.}
\label{tab:2ndOrder_def}
\end{table}
The specific instances of the feature map $\topwdT$ together with corresponding projections $\pi_{kk'}^{{\{E,A,\xi\}}}$ include a variety of systems that have found a wide range of applications in physics, biology, ecology, and social science; see the examples in the chart below. We assume that the function $\mathcal{V}$ is Lipschitz and known and so are all the $\pi_{kk'}^{{\{E,A,\xi\}}}$'s. The Lipschitz assumption is necessary for us to control the trajectory error and ensure the well-posedness of the system. The function $\mathcal{V}$ is a uniform way to collect all of the different variables (functions of the inputs) used across any of the $(k,k')$ pairs over all of the $E,A,\xi$ functions in the system. This uniformity is helpful when discussing the rate of convergence, among other places. Examples of where this generality matters emerge naturally, say when one has a different number of variables across interaction kernels for different pairs $(k,k')$, or when the energy and alignment kernels depend on $r$ and then additional but distinct other variables.
From this uniform set of variables, we then project (which of course implies that the feature maps are all Lipschitz) to arrive at the relevant function $\mbf{s}_{(k, k')}^{{\{E,A,\xi\}}}$ for each pair and each of the elements of the wildcard. Lastly, we can then evaluate this map at the specific pair of agents $(i,i')$, that leads to the feature evaluation, $\mbf{s}_{ii'}^{{\{E,A,\xi\}}}$ which is the expression used in the model equation \eqref{eq:2ndOrder}.

The models encompassed by the form \eqref{eq:2ndOrder} are quite diverse. For a concrete example, please see section \ref{sec:example_AD}.
We summarize many examples in table \ref{tab:2ndOrder_Examples} with a shaded cell indicating that the model has that characteristic, and an empty cell indicates that the model does not have this characteristic. A numeric value indicates this is the number of unique variables, $\bm{\mathcal{V}}, \bm{\mathcal{V}}_{\xi}$ used within the EA or $\xi$ portions of the system. The number of these unique variables specifies the dimension in the minimax convergence rate, see section \ref{sec:rateofconv}.

\begin{table}[]
\begin{tabular}{|l||l|l|l|l|l|l|l|l|l|l|l|l|}
\hline
      & \multicolumn{12}{l|}{Properties} \\ \hline
Model &$\phi^E$&$\phi^A$&$m_i$&$\force^{\dot\bx}$&$\phi^{\xi}$&$\force^{\xi}$ &$K>1$&$\topwd^E$&$\topwd^A$&$\topwd^{\xi}$&$|\bm{\mathcal{V}}|$&$|\bm{\mathcal{V}}_{\xi}|$\\ \hline
Anticipation Dynamics  &    \cellcolor[HTML]{343434}    &   \cellcolor[HTML]{343434}     &  \cellcolor[HTML]{343434}   &         &            &           &     &    \cellcolor[HTML]{343434}     &     \cellcolor[HTML]{343434}    &          &          2       &  \\ \hline
Celestial Mechanics     &   \cellcolor[HTML]{343434}      &        &  \cellcolor[HTML]{343434}   &         &            &           & \cellcolor[HTML]{343434}     &         &         &          &         1       &  \\ \hline
Cucker-Smale      &        &   \cellcolor[HTML]{343434}      &  \cellcolor[HTML]{343434}   &         &            &           &     &         &         &          &       1          &  \\ \hline
Fish Milling 2D      &    \cellcolor[HTML]{343434}   &        & \cellcolor[HTML]{343434}   &   \cellcolor[HTML]{343434}       &            &           &     &         &         &          &        1         &  \\ \hline
Fish Milling 3D &    \cellcolor[HTML]{343434}     &        &  \cellcolor[HTML]{343434}   &    \cellcolor[HTML]{343434}      &            &           &     &         &         &          &          1       &  \\ \hline
Flocking w. Ext. Poten.  &    \cellcolor[HTML]{343434}     & \cellcolor[HTML]{343434}       & \cellcolor[HTML]{343434}   &    \cellcolor[HTML]{343434}      &            &           &     &         &         &          &          1       &  \\ \hline
Phototaxis     &        &    \cellcolor[HTML]{343434}     & \cellcolor[HTML]{343434}   &  \cellcolor[HTML]{343434}        &           \cellcolor[HTML]{343434}  &       \cellcolor[HTML]{343434}     &      &         &         &          &       1          &  1 \\ \hline
Predator-Swarm ($2^{nd}$ Order)    &    \cellcolor[HTML]{343434}     &        &  \cellcolor[HTML]{343434}   &   \cellcolor[HTML]{343434}       &            &           &    \cellcolor[HTML]{343434}  &         &         &          &      1           &  \\ \hhline{|=||=|=|=|=|=|=|=|=|=|=|=|=|}
Lennard-Jones     &   \cellcolor[HTML]{343434}      &        &     &         &            &    \cellcolor[HTML]{343434}        &     &         &         &          &          1       &  \\ \hline
 Opinion Dynamics  &  \cellcolor[HTML]{343434}      &        &    &         &            &           &     &         &         &          &         1        &  \\ \hline
 Predator-Swarm ($1^{st}$ Order)    &    \cellcolor[HTML]{343434}     &        &     &          &            &           &    \cellcolor[HTML]{343434}  &         &         &          &      1           &  \\ \hline
Synchronized Oscillator    &        &    \cellcolor[HTML]{343434}     &     &   \cellcolor[HTML]{343434}       &    \cellcolor[HTML]{343434}         &   \cellcolor[HTML]{343434}         &     &     \cellcolor[HTML]{343434}     &         &      \cellcolor[HTML]{343434}     &                2 & 2 \\ \hline
\end{tabular}
\caption{Summary of the models studied in this work and in \cite{lu2019nonparametric,Tang2019,Zhong20,GravLetter20}}
\label{tab:2ndOrder_Examples}
\end{table}
Our second-order model equations cover the first-order models considered in \cite{lu2019nonparametric, Tang2019, Zhong20} as special cases, see Appendix \ref{sec:1stOrder}, which is why we choose second-order models as the main focus of this work.  Furthermore, the dynamical characteristics produced by second-order models are much richer and can model more complicated collective motions and emergent behavior of the agents. Note that our second-order model in \eqref{eq:2ndOrder}, even when written as a first-order system in more variables, is a strict generalization of the previous first-order analysis.

\subsection*{Rate of convergence notation}
For the system \eqref{eq:2ndOrder}, depending on the number of variables, we will have a different rate of convergence as the dimension of the function(s) we are learning changes. In order to present a unified theorem, we adopt the following notation. Let $\bm{\mathcal{V}}$ denote the set of features in the range of $\topwdT$ that are arguments of ${\bm{\phi}}^E$ or ${\bm{\phi}}^A$ across each of the $(k,k')$ pairs.
For many collective dynamical systems $\bm{\mathcal{V}} = \{r\}$. In the case where the system has both ${\bm{\phi}}^{E}$ and ${\bm{\phi}}^{A}$ but $\bm{\mathcal{V}} = \{r\}$, it is easy to see from Theorem \ref{main:convrate} below that we still only pay the $1$-dimensional rate.
Analogously, we define $\bm{\mathcal{V}}_{\xi}$ to be the set of all features in the range of $\topwdT$ that are arguments  involved in the ${\bm{\phi}}^{\xi}$ part of \eqref{eq:2ndOrder} across all of the $(k,k')$ pairs. This notation is used to frame the convergence rate theorem on the $\xi$ variable as well.

\section{Preliminaries and notation} \label{sec:notation}
We vectorize the models given in \eqref{eq:2ndOrder} in order to give them a more compact description.  Letting $\mbf{v}_i(t) = \dot\mbf{x}_i(t)$, we take the following notations
\[
\mbf{X}(t) = \begin{bmatrix} \mbf{x}_1(t) \\ \vdots \\ \mbf{x}_N(t) \end{bmatrix} \in \mathbb{R}^{Nd}, \quad \mbf{V}(t) := \begin{bmatrix} \mbf{v}_1(t) \\ \vdots \\ \mbf{v}_N(t) \end{bmatrix} \in \mathbb{R}^{Nd}, \quad \mbf{\Xi}(t) := \begin{bmatrix} \xi_1(t) \\ \vdots \\ \xi_N(t) \end{bmatrix} \in \mathbb{R}^N.
\]
We introduce a weighted norm to measure the system variables, denoted $\norm{\cdot}_{\mathcal{S}(\cdot)}$ and given by,
\begin{equation}\label{eq:sys_var_norm}
\norm{\mbf{Z}}_{\mathcal{S}}^2 := \sum_{i = 1}^N\frac{1}{N_{\mathpzc{k}_i}}\norm{\mbf{z}_i}^2,
\end{equation}
for $\mbf{Z} = \begin{bmatrix} \mbf{z}_1^\top, & \ldots, & \mbf{z}_N^\top \end{bmatrix}^\top$ with each $\mbf{z}_i \in \real^{d}$ or $\real$. Here $\norm{\cdot}$ is the same norm used in the construction of pairwise distance data for the interaction  kernels. In the subsequent equations we {drop the explicit dependence on} $t$ for simplicity.  The weight factor, $\frac{1}{N_{\mathpzc{k}_i}}$, is introduced so that agents of different types are considered equally and we learn well even in the case that the number of agents in the classes is highly non-uniform.
With thsee vectorized notations, the model in \eqref{eq:2ndOrder} becomes,
\[
\begin{dcases}
\vec{m} \circ \ddot\mbf{X} &= \rhsf^{\text{nc}, \dot\bx}(\mbf{X}, \mbf{V}, \mbf{\Xi}) + \mathbf{f}^{\bm{\phi}^E}(\mbf{X}, \mbf{V}, \mbf{\Xi}) + \mathbf{f}^{\bm{\phi}^A}(\mbf{X}, \mbf{V}, \mbf{\Xi}) \\
\dot\mbf{\Xi} &= \rhsf^{\text{nc}, \xi}(\mbf{X}, \mbf{V}, \mbf{\Xi}) + \mathbf{f}^{{\bm{\phi}}^{\xi}}(\mbf{X}, \mbf{V}, \mbf{\Xi}).
\end{dcases}
\]
Here $\vec{m} = \begin{bmatrix} m_1, & \ldots, & m_N \end{bmatrix}^\top \in \mathbb{R}^N$, $\circ$ is the Hadamard product, and we use boldface fonts to denote the vectorized form of our estimators (with some once-for-all-fixed ordering of the pairs $(k,k')_{k,k'=1,\dots,K}$):
\begin{equation} \label{not:vectorized}
\bm{\phi}^E = \{\phi^{E}_{k k'}\}_{k, k' = 1}^{K}, \quad \bm{\phi}^A = \{\phi^{A}_{k k'}\}_{k, k' = 1}^{K}, \quad \bm{\phi}^{\xi} = \{\phi^{\xi}_{k k'}\}_{k, k' = 1}^{K}.
\end{equation}
We also use the shorthand,
\begin{equation} \label{not:shorthand}
\bm{\phi}^{EA} = \bm{\phi}^{E} \oplus \bm{\phi}^{A},
\end{equation}
to denote the element of the direct sum of the function spaces containing $\bm{\phi}^{E}, \bm{\phi}^{A}$. This notation will be used throughout on the energy and alignment ($EA$ for short) portion of the system in order to simplify the notation.

We denote by $\rhsf^{\text{nc}, \dot\bx}$, the vectorized notation for the non-collective force defined as follows,
$\rhsf^{\text{nc}, \dot\bx}(\mbf{X}, \mbf{V}, \mbf{\Xi}) := \begin{bmatrix} \force^{\dot\bx}(\mbf{x}_1, \dot\mbf{x}_1, \xi_1), \hdots, \force^{\dot\bx}(\mbf{x}_N, \dot\mbf{x}_N, \xi_N)  \end{bmatrix}^T \in \mathbb{R}^{Nd}$,
and
$$\mathbf{f}^{\bm{\phi}^E} := \bigg[
\sum_{i'=1}^N \frac{1}{N_{\mathpzc{k}_{i'}}} \phi^{E}_{\mathpzc{k}_1 \mathpzc{k}_{i'}}(r_{1i'}, \topwd^E_{1i'})(\mbf{x}_{i'} - \mbf{x}_1),
\hdots,
\sum_{i'=1}^N \frac{1}{N_{\mathpzc{k}_{i'}}} \phi^{E}_{\mathpzc{k}_N \mathpzc{k}_{i'}}(r_{Ni'}, \topwd^E_{Ni'})(\mbf{x}_{i'} - \mbf{x}_N)
\bigg]^T \in \mathbb{R}^{Nd}.
$$
We omit the analogous definitions for $\mathbf{f}^{\bm{\phi}^A}$ and $\mathbf{f}^{\bm{\phi}^{\xi}}$.

\subsection{Trajectory Performance Measurement}  \label{subsec:trajperf}
We will also consider another measurement to assess the learning performance of the estimated kernels in terms of trajectory error. We compare the observed trajectories to the estimated trajectories evolved from the same initial conditions but with the estimated interaction kernels.   Let $\bY(t) = [ \mbf{X}^\top(t),  \mbf{V}^\top(t),  \mbf{\Xi}^\top(t) ]^\top$ be the trajectory from dynamics generated by the true/unknown interaction kernels with initial condition, $\bY(0)$; and $\hat\bY(t) = [ \hat\mbf{X}^\top(t), \hat\mbf{V}^\top(t), \hat\mbf{\Xi}^\top(t)]^\top$ be the trajectory from dynamics generated by the estimated interaction kernels learned from observation of $\{\bY(t_l)\}_{l = 1}^L$ with the same initial condition, $\bY(0)$ (i.e., $\hat\bY(0) = \bY(0)$).  We define a norm for the difference between $\bY$ and $\hat\bY$ at time $t_l$:

\begin{equation}\label{eq:y_norm}
\norm{\bY(t_l) - \hat\bY(t_l)}_{\mathcal{Y}}^2 = \norm{\mbf{X}(t_l) - \hat\mbf{X}(t_l)}_{\mathcal{S}}^2 + \norm{\mbf{V}(t_l) - \hat\mbf{V}(t_l)}_{\mathcal{S}}^2 + \norm{\mbf{\Xi}(t_l) - \hat\mbf{\Xi}(t_l)}_{\mathcal{S}}^2\,,
\end{equation}
and a corresponding norm on the trajectory $\bY_{[0, T]} = \{\bY(t_l)\}_{l = 1}^L$ ($0 = t_1 < \dots < t_L = T$):
\begin{equation}\label{eq:traj_norm}
\norm{\bY_{[0, T]} - \hat\bY_{[0, T]}}_{\text{traj}} = \max_{l = 1, \ldots, L} \norm{\bY(t_l) - \hat\bY(t_l)}_{\mathcal{Y}}\,.
\end{equation}
We also consider a relative version, invariant under changes of units of measure:
\[
\norm{\bY_{[0, T]} - \hat\bY_{[0, T]}}_{\text{traj}^*} = \frac{\norm{\bY_{[0, T]} - \hat\bY_{[0, T]}}_{\text{traj}}}{\norm{\bY_{[0, T]}}_{\text{traj}}}\,.
\]
Lastly, we report errors between $\mbf{X}_{[0, T]}$ and $\hat\mbf{X}_{[0, T]}$,
\[
\norm{\mbf{X}_{[0, T]} - \hat\mbf{X}_{[0, T]}}_{\mathcal{S}^*} = \frac{\max_{l = 1,\ldots,L}\{\norm{\mbf{X}(t_l) - \hat\mbf{X}(t_l)}_{\mathcal{S}}\}}{\max_{l = 1, \ldots,L}\{\norm{\mbf{X}(t_l)}_{\mathcal{S}}\}}.
\]
Similar re-scaled norms are used for the difference between $\mbf{V}_{[0, T]}$ and $\hat\mbf{V}_{[0, T]}$, and for the difference between $\mbf{\Xi}_{[0, T]}$ and $\hat\mbf{\Xi}_{[0, T]}$.

\subsection{Function spaces} \label{sec:functionspaces}
We begin by describing some basic ideas about measures and function spaces.
Consider a compact or precompact set $\mathcal{U} \subset \mathbb{R}^{p}$ for some $p$, then we define the infinity norm as,
$$\|h\|_{\infty}=\operatorname{ess} \sup _{x \in \mathcal{U}}|h(x)|.$$
Further define, $L^{\infty}(\mathcal{U})$ as the space of real valued functions defined on $\mathcal{U}$ with finite $\infty$-norm.
A key function space we need to consider is,
$$C_{c}^{k, \alpha}(\mathcal{\mathcal{U}}) \text { for } k \in \mathbb{N}, 0<\alpha \leq 1,$$ defined as the space of compactly supported, $k$ times continuously differentiable functions with a $k$-th derivative that is H\"older continuous of order $\alpha$.
We can then consider vectorizations of these spaces as
$$\boldsymbol{L}^{\infty}(\mathcal{U}):=\bigoplus_{k, k^{\prime}=1,1}^{K, K} L^{\infty}(\mathcal{U}),$$ which has the vectorized infinity norm given by
$$\|\boldsymbol{f}\|_{\infty}=\max _{k, k^{'}}\left\|f_{k k^{\prime}}\right\|_{\infty}, \forall \boldsymbol{f} \in \boldsymbol{L}^{\infty}(\mathcal{U}).$$
Similarly, we consider direct sums of measures, with corresponding vectorized function spaces. This is done explicitly in section \ref{sec:PerfMeasures} where we define a weighted $L^2$ space (under a particular measure with an associated norm that induces a weighting) and consider vectorized versions of it.

In order to develop a theoretical foundation, and in line with the literature, we make assumptions on the class of functions that can arise as interaction kernels in the model \eqref{eq:2ndOrder}. As the agents get farther and farther apart, they eventually should have no influence on each other. This is an approximation to the vanishing, or rapidly decreasing to zero, nature of pairwise interaction as distance increases that is observed in many physical models. We thus assume a maximum interaction radius for the interaction kernels which represents the maximal distance at which one agent can influence another. Similar assumptions will be made on the feature maps.

More precisely, for each pair $(k,k')$ we consider the following spaces,
\begin{align}
\mathcal{K}_{k k'}^{\{E,A,\xi\}} := L^{\infty}([R_{k k'}^{\min}, R_{k k'}^{\max}] \times \mathbb{S}^{{\{E,A,\xi\}}}_{k k'})\,,
\end{align}
for $k, k' = 1, \ldots, K$ where recall that the ${\{E,A,\xi\}}$ notation means that there is an admissible space (or more generally an expression) for each element of the set ${\{E,A,\xi\}}$.
Here, $R_{k\idxcl'}^{\min},R_{k\idxcl'}^{\max}$ are the minimum or maximum, respectively, possible interaction radii for agents in $C_{k'}$ influencing agents in $C_{k}$.
Similarly, $\mathbb{S}^{E}_{k k'}, \mathbb{S}^{A}_{k k'}, \mathbb{S}^{\xi}_{k k'}$ are compact sets in $\real^{p_{kk'}^E}, \real^{p_{kk'}^A}, \real^{p_{kk'}^{\xi}}$ which contain the ranges of the feature maps, $\topwd^E_{kk'}, \topwd^A_{kk'}$ and $\topwd^{\xi}_{kk'}$.

We can also define the vectorizations of these spaces, which we will look at subsets of when we define the admissible spaces below, given by
$$
\bm{\mathcal{K}}^{{\{E,A,\xi\}}} :=\bigoplus_{k,k'=1,1}^{K,K} \mathcal{K}_{k k'}^{\{E,A,\xi\}}.
$$
In order to provide uniform bounds, we introduce the following sets:
\begin{equation} \label{eq:compactsets}
\mathbf{S}^E := \prod_{k,k'}\mathbb{S}_{kk'}^E \qquad
\mathbf{S}^A := \prod_{k,k'}\mathbb{S}_{kk'}^A \qquad
\mathbf{S}^{\xi} := \prod_{k,k'}\mathbb{S}_{kk'}^{\xi}
\end{equation}

Next, we introduce notation to bound the interaction radii on the pairwise distances and pairwise velocities.
\begin{remark}
Here we note an important distinction. There is the range of the norms of pairwise interactions generated by the dynamics,  and the underlying support of the interaction kernels themselves. These two notions of interaction radius are distinct, we will comment on the estimation and subtleties of both in the numerical algorithm section.
\end{remark}
We let
\begin{align} \label{eq:Rdef}
\mathbf{R} &:= \prod_{k,k'}[R_{kk'}^{\min},R_{kk'}^{\max}]\\
R &:= \max_{k,k'}R_{kk'}^{\max}.
\end{align}
Notice that a uniform support for all interaction kernels on the pairwise distance variable $r$ is $[0,R]$.

We denote the distribution of the initial conditions by $\bm{\mu^{\bY}}$. This measure is unknown and is the source of randomness in our system. It reflects that we will observe trajectories which start at different initial conditions, but that evolve from the same dynamical system, which allows for learnability. For the numerical experiments we will choose our initial conditions to be sampled uniformly over a system dependent range.

Due to the form of \eqref{eq:2ndOrder}, and the norms defined below, we consider the following dynamics induced ranges of the variables. Note that the first supremum is taken over the initial conditions, each of which generate different solutions $\mbf{x}_i(t)$ which are used in the second supremum.
\begin{align}
R_{\dot{x}} &:= \sup_{\bY(0) \sim \bm{\mu^{\bY}}} \sup_{t \in [0,T]} \max_{i,i'} \|\dot{\mbf{x}}_i(t) - \dot{\mbf{x}}_{i'}(t) \| \label{eq:Rdyndef} \\
R_{\xi}   &:= \sup_{\bY(0) \sim \bm{\mu^{\bY}}} \sup_{t \in [0,T]} \max_{i,i'} \|\xi_i(t) - \xi_{i'}(t) \|
\label{eq:Rdyndef:xi}
\end{align}
We assume in this work that both of these quantities $R_{\dot{x}}, R_{\xi}$ are finite. This will be easily satisfied if the measures $\bm{\mu^{\dot{x}}}, \bm{\mu}^{\xi}$ (specifying the distribution of the initial conditions on the velocities and the environment variable) are compactly supported. This follows by the assumptions on the interaction kernels below and that we only consider finite final time $T$.

In order for the second-order systems given by \eqref{eq:2ndOrder} to be well-posed, we assume that the interaction kernels lie in \textit{admissible sets}. For each of the kernels, let $\mathcal{U}_{kk'}:=[R_{k k'}^{\min}, R_{k k'}^{\max}] \times \mathbb{S}^{{\{E,A,\xi\}}}_{k k'}$ and define

\begin{align} \label{eq:E-Admissibility}
\boldsymbol{\mathcal{K}}_{S_{{\{E,A,\xi\}}}}^{{\{E,A,\xi\}}}&:=\Big\{\left(\phi^{{\{E,A,\xi\}}}_{k k'}\right)_{k,k' = 1,1}^{K,K}: \text{{ for all }} k,k'=1,\dots,K\,, \\
&\qquad\phi^{{\{E,A,\xi\}}}_{k k'} \in C_c^{0,1}\left(\mathcal{U}_{kk'}\right),
 \norm{\phi^{{\{E,A,\xi\}}}_{k k'}}_{\infty} + \operatorname{Lip}\left[\phi^{{\{E,A,\xi\}}}_{k k'}\right] \leq S_{{\{E,A,\xi\}}}\Big\}\,. \nonumber
\end{align}
When estimating the $EA$ part of the system, we will consider the direct sum admissible space, for $S_{EA}\geq \max\{S_E,S_A\}$,
\begin{equation} \label{eq:dsum-admissibility}
\boldsymbol{\mathcal{K}}_{S_{EA}}^{EA} := \boldsymbol{\mathcal{K}}_{S_{E}}^{E} \oplus \boldsymbol{\mathcal{K}}_{S_{A}}^{A}
\end{equation}

The admissibility assumptions allow us to establish properties such as existence and uniqueness of solutions to \eqref{eq:2ndOrder} as well as to have control on the trajectory errors in finite time $[0, T]$. It further allows us to show regularity and absolute continuity with respect to Lebesgue measure of the appropriate performance measures defined in section \ref{sec:PerfMeasures}.

In the learning approach, we will consider hypothesis spaces that we will search in order to estimate the various interaction kernels.
The hypothesis spaces corresponding to $\{\phi_{kk'}^{{\{E,A,\xi\}}}\}$ are denoted as $\{\mathcal{H}_{kk'}^{{\{E,A,\xi\}}}\}$ and we vectorize them as,
\begin{align} \label{not:hypspaces}
\mbf{\mathcal{H}}^{{\{E,A,\xi\}}} := \bigoplus_{k,k'=1,1}^{K,K}\mathcal{H}_{kk'}^{{\{E,A,\xi\}}}.
\end{align}
Analogous to our simplified notation for $\bm{\phi}^{EA}, \bm{\varphi}^{EA}$ described in \eqref{not:shorthand}, we define the direct sum of the hypothesis spaces as,
\begin{equation} \label{not:bhypspaceEA}
\mbf{\mathcal{H}}^{EA} := \mbf{\mathcal{H}}^E\oplus \mbf{\mathcal{H}}^A
\end{equation}
We will consider specific hypothesis spaces during the learning theory and numerical algorithm sections.

\section{Inference problem and learning approach}\label{sec:learn_framework}
In this section, we first introduce the problem of inferring the interaction kernels from observations of trajectory data and give a brief review and generalization of the learning approach proposed in the works \cite{lu2019nonparametric} and \cite{Zhong20}.

\subsection{Problem setting}
 Our observation data is given by $\{\bY^{(m)}(t_l), \smash{\dot\bY^{(m)}(t_l)}\}_{m=1,l= 1}^{M,L}$ for $0 = t_1 < t_2 < \cdots < t_L = T$. Here $\dot\bY=[\mbf{y}_1^\top,\cdots, \mbf{y}_N^\top]$ and $ \mbf{y}_i = \begin{bmatrix}\mbf{x}_i^\top(t), \dot\mbf{x}_i^\top(t), \xi_i(t)\end{bmatrix}^\top$.  For simplicity, we only consider equidistant observation points: $t_{l} - t_{l - 1} = h$ for $l = 2, \ldots, L$.  The proposed algorithm with slight modifications also works for non-equispaced points. The $M$ sets of discrete trajectory data are generated by the system (2.1) with the unknown set of interaction kernels, i.e. $\bm{\phi}^E, \bm{\phi}^A, \bm{\phi}^{\xi}$, whose  initial conditions $\bY^{(m)}(0)$ are drawn i.i.d from $\bm{\mu^{\bY}}$, a probability measure defined on the space $\real^{N(2d+2)}$. The goal is to infer the unknown interaction kernels directly from data.

\subsection{Loss functionals}  \label{sec:lossfunctionals}
Given observations, the references \cite{lu2019nonparametric, Tang2019, Zhong20} proposed an empirical error functional, recalling the notational convention \eqref{not:shorthand},
\begin{eqnarray}
\mbf{\mathcal{E}}_{M}^{EA}(\bvpEA)&=&\frac{1}{LM}\sum_{l=1, m=1}^{L,M}
\Big\Vert\ddot{\mbf{X}}^{(m)}(t_l) -  \rhsf^{\text{nc}, \dot\bx}(\bX^{(m)}(t_l),\bV^{(m)}(t_l),\bXi^{(m)}(t_l)) \label{eq:Error_Func_Empirical:EA} \\
&-& \mathbf{f}^{\bm{\varphi}^{E}}(\bX^{(m)}(t_l),\bV^{(m)}(t_l),\bXi^{(m)}(t_l)) - \mathbf{f}^{\bm{\varphi}^{A}}(\bX^{(m)}(t_l),\bV^{(m)}(t_l),\bXi^{(m)}(t_l))\Big\Vert^2_{\mathcal{S}},  \nonumber\\
\mbf{\mathcal{E}}_{M}^{\xi}(\bm{\varphi}^{\xi})&=& \frac{1}{LM} \sum_{l=1, m=1}^{L,M} \Big\Vert\dot{\bm{\Xi}}(t_l) - \rhsf^{\text{nc}, \xi}(\bX^{(m)}(t_l),\bV^{(m)}(t_l),\bXi^{(m)}(t_l)) - \label{eq:Error_Func_Empirical:xi} \\
& &\mathbf{f}^{\bm{\varphi}^{\xi}}(\bX^{(m)}(t_l),\bV^{(m)}(t_l),\bXi^{(m)}(t_l))\Big\Vert^2_{\mathcal{S}}\,.\nonumber
\end{eqnarray}
The estimators of interaction kernels are defined as the minimizers of  the error functionals $\mbf{\mathcal{E}}_{M}^{EA}$ and $\mbf{\mathcal{E}}_{M}^{\xi}$ over $\mbf{\mathcal{H}}^{EA}$ and $\mbf{\mathcal{H}}^{\xi}$ respectively, i.e.
\begin{equation} \label{not:minimizers}
\widehat{\bm{\phi}}^{EA} = \argmin{\bm{\varphi}^{EA} \in \mbf{\mathcal{H}}^{EA}}\mbf{\mathcal{E}}_{M}^{EA}(\bm{\varphi}^{EA})\,,
\qquad
{\widehat{\bm{\phi}}}^{\xi} = \argmin{{\bm{\varphi}}^{\xi} \in \mbf{\mathcal{H}}^{\xi}}\mbf{\mathcal{E}}_{M}^{\xi}({\bm{\varphi}}^{\xi}).
\end{equation}
For the learning theory, we will consider the following error functionals.
On the energy and alignment portion, we consider,
\begin{align}
\mbf{\mathcal{E}}_{\infty}^{EA}(\bvpEA)&:=\mathbb{E}_{\bm{\mu^{\bY}}}\frac{1}{L}\sum_{l=1}^{L}\Big[
\Big\Vert\ddot{\mbf{X}}(t_l) -  \rhsf^{\text{nc}, \dot\bx}(\bX(t_l),\bV(t_l),\bXi(t_l)) \nonumber \\
&- \mathbf{f}^{\bm{\varphi}^{E}}(\bX(t_l),\bV(t_l),\bXi(t_l)) - \mathbf{f}^{\bm{\varphi}^{A}}(\bX(t_l),\bV(t_l),\bXi(t_l))\Big\Vert^2_{\mathcal{S}}
\Big].  \label{eq:Error_Func:EA}
\end{align}
Similarly, on the $\xi$ portion, we consider,
\begin{equation}
\mbf{\mathcal{E}}_{\infty}^{\xi}(\bm{\varphi}^{\xi})= \mathbb{E}_{\bm{\mu^{\bY}}}\frac{1}{L} \sum_{l=1}^L \Big[
\Big\Vert\dot{\bm{\Xi}}(t_l) - \rhsf^{\text{nc}, \xi}(\bX(t_l),\bV(t_l),\bXi(t_l)) - \mathbf{f}^{\bm{\varphi}^{\xi}}(\bX(t_l),\bV(t_l),\bXi(t_l))\Big\Vert^2_{\mathcal{S}}
\Big]. \label{eq:Error_Func:xi}
\end{equation}
We can relate these error functionals to the natural empirical error functionals introduced at the start of the section as follows. By the Strong Law of Large Numbers we have that,
\[
\mbf{\mathcal{E}}_{\infty}^{EA}(\bvpEA)= \lim_{M\to \infty} \mbf{\mathcal{E}}_{M}^{EA}(\bvpEA) \quad \text{and} \quad
\mbf{\mathcal{E}}_{\infty}^{\xi}(\bm{\varphi}^{\xi}) =\lim_{M\to \infty} \mbf{\mathcal{E}}_{M}^{\xi}(\bm{\varphi}^{\xi}).
\]
The relationship between the theoretical and empirical error functionals will play a key role in the learning theory.

\subsection{Overview of theoretical contributions} \label{sec:theory_overview}
The papers \cite{Zhong20, lu2019nonparametric,Tang2019}  have applied this learning approach to a variety of systems and the extensive numerical simulations demonstrate the effectiveness of the approach. However, theoretical guarantees of the proposed approach for second order systems had not been fully developed and will be the main focus of this paper. Our theory contains the first-order theory in \cite{Tang2019} as a special case, as discussed in Appendix \ref{sec:1stOrder}.  We focus on the regime where $L$ is fixed but $M \rightarrow \infty$. We provide a learning theory that answers the fundamental questions:
\begin{itemize}
\item \textbf{Quantitative description of estimator errors}. We will introduce measures to describe how close the estimators are to the true interaction kernels, that lead to novel dynamics-adapted norms. See section \ref{sec:LT}.
\item \textbf{Identifiability of kernels}. We will establish the existence and uniqueness of the estimators as well as relate the solvability of our inverse problem to a fundamental coercivity property.  See section \ref{sec:coercivity}.
\item  \textbf{Consistency and optimal convergence rate of the  estimators}. We will prove theorems on strong consistency and optimal minimax rates of convergence of the estimators, which exploit the separability of the learning on the energy and alignment from the learning on the environment variable. See section \ref{sec:main_results}.
\item \textbf{Trajectory Prediction} We prove a theorem that describes the performance of the estimated dynamics using the estimated kernels compared to the true dynamics. Our result demonstrates how the expected supremum error (over the entire time interval) of our trajectories is controlled by the norm of the difference between the true and estimated kernels, further justifying our choice of norms and estimation procedure. See section \ref{sec:trajectory}.
\end{itemize}

\subsection{Hypothesis Space and Algorithm}  \label{sec:hypspaceandalg}
First, we choose finite-dimensional subspaces of $\mathcal{H}^{E}_{k\idxcl'}$, i.e., $\tilde{\mathcal{H}}_{k\idxcl'}^E \subset \mathcal{H}^{E}_{k\idxcl'}$, whose basis functions are piece-wise polynomials of varying degrees (other type of basis functions are also possible, e.g., clamped B-splines as shown in \cite{lu2019nonparametric}); similarly for $\tilde{\mathcal{H}}_{k\idxcl'}^A \subset \mathcal{H}^{A}$.  Hence, each test function, $\varphi^{E}_{k k'}, \varphi^{A}_{k k'}$, can be expressed in terms of the linear combination of the basis functions as follows
\begin{align*}
\varphi^{E}_{k\idxcl'}(r, \topwd^E) &= \sum_{\eta_{k\idxcl'}^{E} = 1}^{n_{k\idxcl'}^{E}} \alpha_{k, k', \eta_{k\idxcl'}^{E}}^{E}\psi^{\bx}_{k, k', \eta_{k\idxcl'}^{E}}(r, \topwd^E), \\
\varphi^{A}_{k\idxcl'}(r, \dot{r}, \topwd^A) &= \sum_{\eta_{k\idxcl'}^{A} = 1}^{n_{k\idxcl'}^{A}} \alpha_{k, k', \eta_{k\idxcl'}^{A}}^{A}\psi^{\bx}_{k, k', \eta_{k\idxcl'}^{A}}(r, \dot{r}, \topwd^A),.
\end{align*}
Substituting this linear combination back into \eqref{eq:Error_Func_Empirical:EA}, we obtain a system of linear equations,
\[
A^{EA}\vec{\alpha}^{EA} = \vec{b}^{EA}.
\]
Here, $\vec{\alpha}^{EA} \in \real^{n^{EA}} = \begin{bmatrix} (\vec{\alpha}^E)^\top & (\vec{\alpha}^A)^\top \end{bmatrix}^\top$ with $\vec{\alpha}^E$ and $\vec{\alpha}^A$ being the collection of $\alpha_{k, k', \eta_{k\idxcl'}^{E}}^{E}$ or $\alpha_{k, k', \eta_{k\idxcl'}^{A}}^{A}$ respectively.  Moreover, $A^{EA} \in \real^{n^{EA} \times n^{EA}}$ and $\vec{b}^{EA} \in \real^{n^{EA}}$.  See Sec. \ref{sec:algorithm_numeric} for full details.

The overhead memory storage needed is $MLN(d(5 + n^{EA} + n^{\xi}) + 3)$, with $MLN(4d + 2)$ needed for trajectory data, $MLNd(n^{EA} + n^{\xi}$ (here $n^{EA} = n^E + n^A$, the sum of the number of basis functions on $E$ and $A$) for learning matrices, and $MLN(d + 1)$ for right hand side vectors.  Hence if $M \gg \mathcal{O}(1)$, we can consider parallelization in $m$ in order to reduce the overhead memory, ending up with $M_{\text{per core}}\big(LN(d(5 + n^{EA} + n^{\xi}) + 3)\big)$ with $M_{\text{per core}} = \frac{M}{n_{\text{cores}}}$.  The final storage of $A$ and $\vec{b}$ only needs $n^{EA}(n^{EA} + 1) + n^{\xi}(n^{\xi} + 1)$.

The computational cost of for solving the system $A^{EA}\vec{\alpha}^{EA} = \vec{b}^{EA}$ is $MLN^2 + MLd(n^{EA})^2 + (n^{EA})^2\log(n^{EA})$, with $MLN^2$ for computing pairwise data, $MLd(n^{EA})^2$ for constructing the learning matrix and right hand side vector, and $(n^{EA})^2\log(n^{EA})$ for solving the linear system.  In the case of choosing the optimal number of basis functions, i.e., $n^{EA}_{\text{opt}} = M^{\frac{|V|}{2s + |V|}}$, we end up with the total computational cost, of the order $M^{1 + \frac{2|V|}{2s + |V|}}$, which is slightly super linear in $M$ , but less than quadratic in the common situation when $2s+|V|\ge 2|V|$.

Similar analysis on solving $A^{\xi}\vec{\alpha}^{\xi} = \vec{b}^{\xi}$ also shows that the computational cost is slightly super linear in $M$.

\section{Learning theory}\label{sec:LT}
Given estimators, how to best measure the estimation error? This is the first question to address in order to have a full understanding of the performance. In earlier works, a set of probability measures that are adapted to the dynamical system and learning setting are introduced to describe how close  the estimators  are to the true interaction kernels.  We will generalize these ideas and measures to our learning problem.

\subsection{Probability measures} \label{sec:PerfMeasures}
The variables we are working with have a natural distribution on the space of pairwise distances and features. Together with the fact that our functions have as arguments the variables ($r,\dot{r},\topwd^E,\topwd^A,\topwd^{\xi}$) which are functions of the state space, we are led to consider probability measures which account for the distribution of the data, while respecting the interaction structure of the system. For further intuition into these measures, see \cite{Tang2019,lu2019nonparametric,BFHM17}.
For each interacting pair $(k,k')$, we introduce the following  probability measures,

\begin{equation}\label{meas:fullEA}
\begin{dcases}
\rho_T^{EA, k, k'}(r, \topwd^E, \dot{r}, \topwd^A)
&:= \mathbb{E}_{\bY_0 \sim \bm{\mu^{\bY}}}\frac{1}{TN_{k k'}}\int_{t = 0}^T\sum_{\substack{i \in C_{k}, i' \in C_{k'} \\ i \neq i'}}\delta_{i i', t}(r, \topwd^E, \dot{r}, \topwd^A)\, dt \\
\rho_T^{EA, L, k,k'}(r, \topwd^E, \dot{r}, \topwd^A) &:= \mathbb{E}_{\bY_0 \sim \bm{\mu^{\bY}}}\frac{1}{LN_{k\idxcl'}}\sum_{l = 1}^L\sum_{\substack{i \in C_{k}, i' \in C_{k'} \\ i \neq i'}}\delta_{i i', t_l}(r, \topwd^E, \dot{r}, \topwd^A) \\
\rho_T^{EA, L, M, k, k'}(r, \topwd^E, \dot{r}, \topwd^A)
&:= \frac{1}{MLN_{k k'}}\sum_{l, m = 1}^{L, M}\sum_{\substack{i \in C_{k}, i' \in C_{k'} \\ i \neq i'}}\delta_{i i', t_l, m}(r, \topwd^E, \dot{r}, \topwd^A)
\end{dcases}
\end{equation}
where $N_{kk'} = N_kN_{k'}$ for $k\neq k'$ and $N_{k k^{\prime}}={N_{k} \choose 2}$ for $k=k'$, and the dirac measures are defined as
\[
\begin{dcases}
\delta_{i i', t}(r, \topwd^E, \dot{r}, \topwd^A) &:= \delta_{r_{i i'}(t), \topwd^E_{i i'}(t), \dot{r}_{i i'}(t), \topwd^A_{i i'}(t)}(r, \topwd^E, \dot{r}, \topwd^A) \\
\delta_{i i', t, m}(r, \topwd^E, \dot{r}, \topwd^A) &:= \delta_{r_{i i'}^{(m)}(t), \topwd^{E, (m)}_{i i'}(t), \dot{r}_{i i'}^{(m)}(t), \topwd^{A, (m)}_{i i'}(t)}(r, \topwd^E, \dot{r}, \topwd^A).
\end{dcases}
\] We use a superscript $(m)$ to denote that the variable is calculated from the data from that $m^{\text{th}}$ trajectory. For example, $\topwd^{E, (m)}_{i i'}(t)$ denotes the feature maps between agents $i$ and $i'$ at time $t$ along the $m^{\text{th}}$ trajectory.

The measure $\rho_T^{EA, L, k, k'}$ is  the discrete counterpart of $\rho_T^{EA, k, k'}$ at the observation time instances. In practice, we can use $\rho_T^{EA, L, M, k, k'}$ to approximate $\rho_T^{EA, L, k, k'}$ since it can be computed from observational data and will converge to $\rho_T^{EA, L, k, k'}$ as $M \rightarrow \infty$.

We also consider the marginal distributions
\begin{equation}\label{meas:rhoE}
\rho_T^{E, k, k'}(r, \topwd^E)       := \int_{\dot{r}}\int_{\topwd^A} \rho_T^{EA, k, k'} \, d\topwd^A \, d\dot{r}\quad,\quad
\rho_T^{A, k, k'}(r, \dot{r}, \topwd^A)       := \int_{\topwd^E} \rho_T^{EA, k, k'} \, d\topwd^E
\end{equation}
and $\rho_T^{E, L, k, k'}(r, \topwd^E) $, $\rho_T^{E, L, M, k, k'}(r, \topwd^E)$, $\rho_T^{A, L, k, k'}(r, \dot{r}, \topwd^A)$,  $\rho_T^{A, L, M, k, k'}(r, \dot{r}, \topwd^A)$ defined analogously as above.
The empirical measures, $\rho_T^{E, L, M, k, k'}, \rho_T^{A, L, M, k ,k'}$, are the ones used in the actual algorithm to quantify the learning performances of the estimators $\widehat{\intkernel}^{E}_{k k'}$ and $\widehat{\intkernel}^{A}_{k k'}$ respectively.  They are also crucial in discussing the separability of $\widehat{\intkernel}^{E}_{k k'}$ and $\widehat{\intkernel}^{A}_{k k'}$.

For ease of notation, we introduce the following measures to handle the heterogeneity of the system, and which are used to describe error over all of the pairs $(k,k')$.
\begin{equation} \label{meas:vectorized}
\bm{\rho}_T^{EA,L}=\bigoplus_{k, k^{\prime}=1,1}^{K, K} \rho_{T}^{EA,L, k k^{\prime}}, \quad
\bm{\rho}_T^{EA}=\bigoplus_{k, k^{\prime}=1,1}^{K, K} \rho_{T}^{EA,k k^{\prime}}, \quad \LtwoB\left(\bm{\rho}_T^{EA,L}\right)=\bigoplus_{k, k^{\prime}=1,1}^{K, K} L^2\left(\rho_{T}^{EA,L, k k^{\prime}}\right)
\end{equation}

Similar definitions apply for measures related to learning the $\xi$-based interaction kernels, see Appendix \ref{sec:pm_xi}.
We discuss some key properties of the measures in Appendix \ref{sec:app:analytical}.

\subsection{Learning performance}
We now discuss the performance measures for the estimated interaction kernels. We have already treated the trajectory estimation error in section \ref{subsec:trajperf}.
We use weighted $L^2$-norms (with mild abuse of notation, we omit the weight from the notation) based on the dynamics-adapted measures introduced above, and with analogous definitions when discrete over $L$:
\begin{eqnarray}
\norm{\widehat{\intkernel}^{E}_{k k'} - \phi^{E}_{k k'}}_{L^2(\rho_T^{E, k, k'})}^2 &:=& \int_{(r,\topwd^E)}(\widehat{\intkernel}^{E}_{k k'}(r, \topwd^E) - \phi^{E}_{k k'}(r, \topwd^E))^2r^2 \, d\rho_T^{E, k, k'} \nonumber\\
\norm{\widehat{\intkernel}^{A}_{k k'} - \phi^{A}_{k k'}}_{L^2(\rho_T^{A, k, k'})}^2 &:=& \int_{(r,\dot{r},\topwd^A)}(\widehat{\intkernel}^{E}_{k k'}(r, \dot{r}, \topwd^A) - \phi^{E}_{k k'}(r, \dot{r}, \topwd^A))^2\dot{r}^2 \, d\rho_T^{A, k, k'} \nonumber \\
\norm{\widehat{\intkernel}_{k k'}^{EA}  - \phi_{k k'}^{EA}}_{L^2(\rho_T^{EA, k, k'})}^2
&:=& \int_{r,\topwd^E,\dot{r},\topwd^A} \Big[(\widehat{\intkernel}^{E}_{k k'}(r, \topwd^E) - \phi^{E}_{k k'}(r, \topwd^E))r \nonumber \\
&  & \qquad + (\widehat{\intkernel}^{E}_{k k'}(r, \dot{r}, \topwd^A) - \phi^{E}_{k k'}(r, \dot{r}, \topwd^A))\dot{r}\Big]^2 \, d\rho_T^{EA, k, k'}\,.\label{e:L2NormDefs}
\end{eqnarray}

Our learning theory focuses on minimizing the difference between $\widehat{\intkernel}^{E}_{k k'} \oplus \widehat{\intkernel}^{A}_{k k'}$ and $\phi^{E}_{k k'} \oplus \phi^{A}_{k k'}$ in the joint norm given by \eqref{e:L2NormDefs}.  As long as the joint norm is small, our estimators  produce faithful approximations of the right hand side function of the original system and  trajectories. However, it does not necessarily imply that both $\widehat{\intkernel}^{E}_{k k'} - \phi^{E}_{k k'}$'s and $\widehat{\intkernel}^{A}_{k k'} - \phi^{A}_{k k'}$'s are small in their corresponding energy- and alignment-based norms, since  the joint norm is a weaker norm. It would be interesting to study if there is any equivalence between these two norms, but the problem appears to be quite delicate.  The theoretical investigation is still ongoing.

Now, we have all the tools needed to establish a theoretical framework: dynamics induced probability measures, performance measurements in appropriate norms, and loss functionals. These will allow us to discuss the convergence properties of our estimators.  Full details of the numerical algorithm are given in Appendix \ref{sec:algorithm_numeric}.
\subsection*{Notational summary}
A summary of the learning theory notation introduced in sections \ref{sec:notation}, \ref{sec:learn_framework}, and the notation above, is given below in table \ref{tab:2ndOrderlearning_def}.
\begin{table}[h]
\centering
\small{
\small{\begin{tabular}{ c | c | c}
\hline
Notation                    & Definition & Ref \\
\hline\hline
$M$   & number of trajectories & Sec. \ref{sec:intro} \\
\hline
$L$ & number of times in $[0,T]$ for each trajectory & Sec. \ref{sec:ModelDesc} \\
\hline
$\bY(t)$ & full state space vector containing $\mbf{X}(t),\mbf{V}(t),\mbf{\Xi}(t)$& Sec. \ref{sec:ModelDesc}  \\
\hline
${\{E,A,\xi\}}$  & wildcard, means the notation applies for all $3$ variables &  Sec. \ref{sec:ModelDesc} \\
\hline
 $\| \mbf{X} \|_{\mathcal{S}}$   &    $\sum_{i = 1}^N\frac{1}{{N_{\mathpzc{k}_i}}}\norm{\mbf{x}_i}^2$ &  \eqref{eq:sys_var_norm}\\
\hline
$\bm{\mu^{\bY}}$ & distribution on the initial conditions $\bY(0)$ & Sec. \ref{sec:functionspaces}\\
\hline
${\bm{\phi}}^{{\{E,A,\xi\}}}=(\phi_{kk'}^{{\{E,A,\xi\}}})_{k,k'}$ &  vectorized true $E,A,\xi$ interaction kernels & \eqref{not:vectorized} \\

\hline
${\bm{\varphi}}^{{\{E,A,\xi\}}}=(\varphi_{kk'}^{{\{E,A,\xi\}}})_{k,k'}$ &  ${\bm{\varphi}}^{{\{E,A,\xi\}}} \in \mbf{\mathcal{H}}^{{\{E,A,\xi\}}}$ with $\varphi_{kk'}^{{\{E,A,\xi\}}} \in \mathcal{H}_{kk'}^{{\{E,A,\xi\}}}$ & \eqref{not:vectorized} \\
 \hline
$EA$ & shorthand denoting energy and alignment part of system & \eqref{not:shorthand} \\
\hline
$\bm{\phi}^{EA}$  & represents the joint function ${\bm{\phi}}^{E}\oplus{\bm{\phi}}^{A}\in \mbf{\mathcal{H}}^{EA}$ & \eqref{not:shorthand} \\
 \hline
$ \norm{\bY}_{\mathcal{Y}}^2$ & $\norm{\mbf{X}}_{\mathcal{S}}^2 + \norm{\mbf{V}}_{\mathcal{S}}^2 + \norm{\mbf{\Xi}}_{\mathcal{S}}^2.$ & \eqref{eq:y_norm}\\
\hline
$\textbf{S}^{{\{E,A,\xi\}}}$ &  $\prod_{k,k'}\mathbb{S}_{kk'}^{{\{E,A,\xi\}}}$ & \eqref{eq:compactsets}\\
\hline
$\mathbf{R}$  & $ \prod_{k,k'}[R_{kk'}^{\min},R_{kk'}^{\max}]$ & \eqref{eq:Rdef} \\
\hline
$R$ &  $\max_{k,k'}R_{kk'}^{\max}$ & \eqref{eq:Rdef}\\
\hline
$\boldsymbol{\mathcal{K}}_{S_{{\{E,A,\xi\}}}}^{{\{E,A,\xi\}}}$ & admissible spaces for the $E,A,\xi$ kernels & \eqref{eq:E-Admissibility} \\
\hline
$\mathcal{H}_{kk'}^{{\{E,A,\xi\}}}$ &  the hypothesis spaces for $\phi_{kk'}^{{\{E,A,\xi\}}}$ & \eqref{not:hypspaces} \\
\hline
$\mbf{\mathcal{H}}^{{\{E,A,\xi\}}}=\oplus_{kk'} \mathcal{H}_{kk'}^{{\{E,A,\xi\}}}$ &  the hypothesis spaces for ${\bm{\phi}}^{{\{E,A,\xi\}}}$ & \eqref{not:hypspaces} \\
\hline
$\mbf{\mathcal{H}}^{EA}$  & direct sum of hypothesis spaces $\mbf{\mathcal{H}}^E\oplus \mbf{\mathcal{H}}^A$ & \eqref{not:bhypspaceEA}\\
\hline
$\mbf{\mathcal{E}}_{M}^{EA}(\cdot),\mbf{\mathcal{E}}_{M}^{\xi}(\cdot)$ &  empirical $EA$ error functional, $\xi$ error functional & Sec. \ref{sec:lossfunctionals} \\
\hline
$\widehat{\bm{\phi}}_{M}^{EA}:=\widehat{\bm{\phi}}_{L,M,\mbf{\mathcal{H}}^{EA}}^{EA}$ & $\mathrm{argmin}_{\bm{\varphi}^{EA} \in \mbf{\mathcal{H}}^{EA}} \bm{\mathcal{E}}_{M}^{EA}(\bm{\varphi}^{EA})$ & \eqref{not:minimizers} \\
\hline
$\widehat{\bm{\phi}}_{M}^{\xi}:=\widehat{\bm{\phi}}_{L,M,\mbf{\mathcal{H}}^{\xi}}^{\xi}$ & $\mathrm{argmin}_{\bm{\varphi}^{\xi} \in \mbf{\mathcal{H}}^{\xi}} \bm{\mathcal{E}}_{M}^{\xi}(\bm{\varphi}^{\xi})$ & \eqref{not:minimizers} \\
\hline
$\{\psi_{kk',p}^{{\{E,A,\xi\}}}\}_{p=1}^{n_{kk'}^{{\{E,A,\xi\}}}}$ & basis for  $\mathcal{H}_{kk'}^{{\{E,A,\xi\}}}$ & Sec. \ref{sec:hypspaceandalg}\\
\hline
$\bm{\rho}_T^{EA}$, $\bm{\rho}_T^{\xi}$  & measure for $EA$, $\xi$ with continuous time, infinite trajectories & \eqref{meas:vectorized}, \ref{sec:pm_xi}  \\
\hline
$\bm{\rho}_T^{EA,L}$, $\bm{\rho}_T^{\xi, L}$ & measure for $EA$, $\xi$ discrete in time, infinite trajectories & \eqref{meas:vectorized}, \ref{sec:pm_xi} \\
\hline
$\LtwoB\left(\bm{\rho}_T^{EA,L}\right)$ & $\bigoplus_{k, k^{\prime}=1,1}^{K, K} L^2\left(\rho_{T}^{EA,L, k k^{\prime}}\right)$ & \eqref{meas:vectorized} \\
\hline
$c_{\mbf{\mathcal{H}}^{EA}}$, $c_{\mbf{\mathcal{H}}^{\xi}}$ & coercivity constant on the $\mbf{\mathcal{H}}^{EA}$, $\mbf{\mathcal{H}}^{\xi}$ hypothesis spaces &  Def.\ref{def_coercivity_2ndorder}  \\
\hline
$\mbf{\mathcal{H}}_M^{EA}$,$\mbf{\mathcal{H}}_M^{\xi}$ & hypothesis spaces on $EA,\xi$ depending on $M$ & Sec. \ref{sec:main_results}\\
\hline
$A^{EA}$, $A^{\xi}$  & Learning matrices for the inverse problem & Sec. \ref{sec:numerics} \\
\hline
$\mathcal{N}(\mbf{\mathcal{H}},\delta)$ &  $\delta$-covering number, under the $\infty$-norm, of a set $\mbf{\mathcal{H}}$ &  \cite{Wellner}\\
\hline
\end{tabular}}
}
\caption{Notation used throughout the paper}
\label{tab:2ndOrderlearning_def}
\end{table}

\section{Identifiability of kernels from data} \label{sec:coercivity}
In this section we introduce a technical condition on the dynamical system that relates to the solvability of the inverse problem and plays a key role in the learning theory. We establish theorems in two directions:
\begin{enumerate}
\item Showing how identifiability of the kernels can be derived from the coercivity condition by relating the coercivity constant to the singular values of the learning matrices associated to our inverse problem, for both finitely and infinitely many trajectories.
\item Establishing the coercivity condition for a wide class of dynamical systems of the form \eqref{eq:2ndOrder}, under mild assumptions on the distribution $\bm{\mu^{\bY}}$ of the initial conditions. From our numerical experiments, we expect that coercivity holds even more generally.
\end{enumerate}

We will make the following assumptions on the hypothesis spaces used in the learning approach for the remainder of the paper:
\begin{assumption}\label{compactnesscondition}
 $\mbf{\mathcal{H}}^{EA}$ is a compact convex subset of $\boldsymbol{\mathcal{K}}_{S_{EA}}^{EA}:= \boldsymbol{\mathcal{K}}_{S_E}^E \oplus \boldsymbol{\mathcal{K}}_{S_A}^A$ (see \ref{eq:dsum-admissibility}) which implies that the infinity norm of all elements in $\mbf{\mathcal{H}}^{EA}$ is bounded above by a constant $S_{EA}\geq \max\{{S_E, S_A}\}$.
\end{assumption}
\begin{assumption}\label{compactnessconditionxi}
 $\mbf{\mathcal{H}}^{\xi}$ is a compact convex subset of $\boldsymbol{\mathcal{K}}_{S_{\xi}}^{\xi}$ (see \ref{eq:E-Admissibility}) and is bounded above by $S_{0}\geq S_{\xi}$.
\end{assumption}
It is easy to see that $\mbf{\mathcal{H}}^{EA}$ can be naturally embedded as a compact subset of $\LtwoB(\bm{\rho}_T^{EA,L})$ and that $\mbf{\mathcal{H}}^{\xi}$ can be naturally embedded as a compact subset of $\LtwoB(\bm{\rho}_T^{\xi,L})$ (recall these measures are defined in section \ref{sec:PerfMeasures}). Assumptions \ref{compactnesscondition}, \ref{compactnessconditionxi} ensure the existence of minimizers to the loss functionals $\mbf{\mathcal{E}}_{M}^{EA}, \mbf{\mathcal{E}}_{M}^{\xi}$ defined in \eqref{eq:Error_Func_Empirical:EA} and \eqref{eq:Error_Func_Empirical:xi}, which will be proven in Appendix \ref{sec:app:learntechnical}.

In order to ensure learnability we introduce a coercivity condition, with terminology coming from the Lax-Milgram theorem.
In the second-order case we will have two coercivity conditions, one for the energy and alignment, and the other for the $\xi$ variable. These conditions serve the same purpose in both cases. They first ensure that the minimizers to the error functionals are unique, and second that when the expected error functional is small, then the distance from the estimator to the true kernels is small in the appropriate $\bm\rho_T$ norm.

Due to their connection to the error functional and the learnability of the kernels, coercivity plays an important role in the theorems of Section \ref{sec:main_results}.

\begin{definition}[Coercivity condition] \label{def_coercivity_2ndorder}
For the dynamical system \eqref{eq:2ndOrder} observed at time instants $0=t_1<t_2<\dots <t_L=T$ and with initial condition distributed $\bm{\mu^{\bY}}$ on $\real^{(2d+1)N}$, it satisfies the coercivity condition on the hypothesis space $\mbf{\mathcal{H}}^{EA}$ with constant $c_{\mbf{\mathcal{H}}^{EA}}$ if
\begin{align}\label{2ndorder:gencoer}
 c_{\mbf{\mathcal{H}}^{EA}}  := \inf_{ \bvpEA \in \mbf{\mathcal{H}}^{EA} \backslash \{\mbf 0  \}} \,\frac{ \frac{1}{L}\sum_{l=1}^{L}\mathbb{E}_{\bm{\mu^{\bY}}} \bigg[ \big\|  \mathbf{f}_{\bvpEA}(\bX(t_l),\bV(t_l),\bXi(t_l)) \big\|_{\mathcal{S}}^2\bigg] }{\|\bvpEA\|_{\LtwoB(\bm{\rho}_T^{EA,L})}^2} >0 .
\end{align}
An analogous definition holds for continuous observations on the time interval $[0,T]$, by replacing the sum over observations at discrete times with an integral over $[0,T]$.
Similarly, the system satisfies the coercivity condition on the hypothesis space $\mbf{\mathcal{H}}^{\xi}$ with constant $c_{\mbf{\mathcal{H}}^{\xi}}$ if
\begin{align}\label{2ndorder:gencoer:xi}
 c_{\mbf{\mathcal{H}}^{\xi}}  := \inf_{ \bm{\varphi}^{\xi} \in \mbf{\mathcal{H}}^{\xi} \backslash \{\mbf 0  \}} \,\frac{ \frac{1}{L}\sum_{l=1}^{L}\mathbb{E}_{\bm{\mu^{\bY}}} \bigg[ \big\|  \mathbf{f}_{\bm{\varphi}^{\xi}}(\bX(t_l),\bV(t_l),\bXi(t_l)) \big\|_{\mathcal{S}}^2\bigg] }{\| \bm{\varphi}^{\xi} \|_{\LtwoB(\bm{\rho}_T^{\xi,L})}^2} >0.
\end{align}
An analogous definition holds for continuous observations on the time interval $[0,T]$, by replacing the sum over observations at discrete times with an integral over $[0,T]$.
\end{definition}

 In the following, we prove the coercivity condition on  general compact sets of $\LtwoB([0,R],\bm{\rho}_T^{EA,L})$ under suitable hypotheses. Our result is independent of $N$, which implies that the finite sample bounds of Theorem \ref{t:2ndordersystem:thm_optRate} can be dimension free -- in that the coercivity constant has no dependence on $N$. This result implies that coercivity may be a fundamental property of the dynamical system, including in the mean field regime ($N\to\infty$).

\subsection{Identifiability from coercivity} \label{sec:Conditioning}

By choosing the hypothesis space to be compact and convex, we are able to show that the error functional has a unique minimizer. However, many possible bases exist that could potentially yield good performance (in terms of the error functional and $\LtwoB$ error to the true kernel).
We want to choose a basis such that the regression matrix, $A^{EA}$
defined in Appendix \ref{sec:algorithm_numeric}, is well-conditioned -- which will ensure that the inverse problem can be solved and thus an estimator can be learned that will have good performance (in terms of the error functional and $\LtwoB(\bm{\rho}_T^{EA,L})$ error to the true kernel).
In the proposition below we establish two results in this direction. The key for both results is that the basis is chosen to be orthonormal in $\LtwoB(\bm{\rho}_T^{EA,L})$, versus the naive choice of basis in the underlying direct sum of $\bm{L}^{\infty}$ spaces that the interaction kernels live in. The first result is theoretical and shows that, under appropriate assumptions on the basis, the minimal singular value of the expected regression matrix (denoted $A_{\infty}^{EA}$) equals the coercivity constant. While the second result is critical for the practical implementation as it lower bounds the minimal singular value of the regression matrix by the coercivity constant with high probability. In both cases, the numerical performance is affected by the size of the coercivity constant of $\mbf{\mathcal{H}}^{EA}$ and if the hypothesis space is well-chosen, then the coercivity constant will be sufficiently large and the regression matrix will be well-conditioned.

To ease the notation, we introduce the bilinear functional $\dbinnerp {\cdot, \cdot}$ on $\mbf{\mathcal{H}}^{EA} \times \mbf{\mathcal{H}}^{EA}$, defined by
\begin{align}
\label{eq:bilinearFn}
 \dbinnerp {\bm{\varphi}^{EA}_1, \bm{\varphi}^{EA}_2}:= \nonumber
 \\ & \hspace{-3cm} \frac{1}{L}\sum_{l,i=1}^{L,N}\frac{1}{N_{\mathpzc{k}_{i}}} \mathbb{E}_{\bm{\mu^{\bY}}}\bigg[\bigg\langle\sum_{i'=1}^{N} \frac{1}{N_{\mathpzc{k}_{i'}}} \Big( \varphi_{1,\mathpzc{k}_i\mathpzc{k}_{i'}}^E(r_{ii'}(t),\topwd^E_{ii'})\mbf{r}_{ii'}(t) +
\varphi_{1,\mathpzc{k}_i\mathpzc{k}_{i'}}^A(r_{ii'}(t),\topwd^A_{ii'})\mbf{\dot{r}}_{ii'}(t) \Big)
, \nonumber \\
& \hspace{-2cm} \sum_{i'=1}^{N}\frac{1}{N_{\mathpzc{k}_{i'}}}\Big( \varphi_{2,\mathpzc{k}_i\mathpzc{k}_{i'}}^E(r_{ii'}(t),\topwd^E_{ii'})\mbf{r}_{ii'}(t) +
\varphi_{2,\mathpzc{k}_i\mathpzc{k}_{i'}}^A(r_{ii'}(t),\topwd^A_{ii'})\mbf{\dot{r}}_{ii'}(t) \Big) \bigg\rangle \bigg]
\end{align} for any
$\bm{\varphi}^{EA}_1=(\combf{\varphi_{1,kk'}^E}{\varphi_{1,kk'}^A})_{k,k'=1,1}^{K,K} \in \mbf{\mathcal{H}}^{EA}$, and
$\bm{\varphi}^{EA}_2=(\combf{\varphi_{2,kk'}^E}{\varphi_{2,kk'}^A})_{k,k'=1,1}^{K,K} \in \mbf{\mathcal{H}}^{EA}$.
For every pair $(k,k')$ let $(\combf{\psi_{kk',i}^E}{\psi_{kk',i}^A})_{i=1}^{n_{kk'}}$
be a basis of
$$\mathcal{H}_{{kk'}}^{EA} \subset L^\infty([0,R]\times \mathbb{S}_{kk'}^E) \oplus L^{\infty}([0,R]\times \mathbb{S}_{kk'}^A) $$
satisfying the orthonomality and boundedness conditions
\begin{align}\label{onb}
\langle \psi_{kk',p}^{EA},\psi_{kk',p'}^{EA} \rangle_{L^2{(\rho_T^{EA,L,kk'})}}&=\delta_{p,p'}\quad,\quad
\|\psi_{kk',p}^{EA}\|_{\infty} &\leq S_{EA}.
\end{align}
We note that multivariable basis functions arise naturally in this setting due to the model. A tensor product basis of splines or piecewise polynomials can be used, as one explicit example. The $n_{kk'}$ notation allows multivariable functions, different choices for the number of basis functions across pairs $(k,k')$, and a different number of basis functions within a pair with respect to the underlying coordinates of the tensor product.

By convention, we use the lexicographic ordering to order within pairs $(k,k')$ (with order $r,\topwd^E, \topwd^A$), and then across pairs (with the lexicographic ordering on pairs of integers).
Set $\mathbf{n}=\sum_{k,k'}n_{kk'}=dim(\mbf{\mathcal{H}}^{EA})$; then for any function ${\bm{\varphi}}^{EA} \in \mbf{\mathcal{H}}^{EA}$, we can write $${\bm{\varphi}}^{EA}=\sum_{p=1}^{\mathbf{n}} a_{p}\mbf{\psi}_{p}^{EA}.$$
Under the setting above, we have the following relationship between the coercivity constant and the minimal singular value of the empirical and expected learning matrix:

\begin{proposition}
\label{2ndordersystem:coercivityconstant}
Consider the matrices $$A_{\infty}^{EA}=\big(\dbinnerp{\mbf{\psi}_{p}^{EA},\mbf{\psi}_{p'}^{EA}}\big)_{p, p'} \in \mathbb{R}^{\mathbf{n} \times \mathbf{n}}\,,\quad A_{\infty}^{\xi}=\big(\dbinnerp{\mbf{\psi}_{p}^{\xi},\mbf{\psi}_{p'}^{\xi}}\big)_{p, p'} \in \mathbb{R}^{\mathbf{n}_{\xi} \times \mathbf{n}_{\xi}}\,,$$ and choose the hypothesis spaces as $\mbf{\mathcal{H}}^{EA} = \mathrm{span}\{ \combf{\mbf{\psi}_{p}^E}{\mbf{\psi}_{p}^A}\}_{p=1}^{\mathbf{n}}$ and $\mbf{\mathcal{H}}^{\xi} = \mathrm{span}\{ \mbf{\psi}_{p}^{\xi}\}_{p=1}^{\mathbf{n}_{\xi}}$. Then the coercivity constants for $\mbf{\mathcal{H}}^{EA}, \mbf{\mathcal{H}}^{\xi}$ are the smallest singular value of $A_{\infty}^{EA}$, $A_{\infty}^{\xi}$, respectively:
\begin{align}
\sigma_{\min}(A_{\infty}^{EA}) =c_{\mbf{\mathcal{H}}^{EA}}\, \qquad  \sigma_{\min}(A_{\infty}^{\xi}) =c_{\mbf{\mathcal{H}}^{\xi}}
\label{2ndorder:gencoer:result}
\end{align}
with $c_{\mbf{\mathcal{H}}^{EA}}$ defined in \eqref{2ndorder:gencoer} and $c_{\mbf{\mathcal{H}}^{\xi}}$ defined in \eqref{2ndorder:gencoer:xi}.
Additionally, for large $M$, the smallest singular value of $A^{EA}$ satisfies the inequality
\begin{align*}
\sigma_{\min}(A^{EA}) \geq  0.8 c_{\mbf{\mathcal{H}}^{EA}}
\end{align*}
with probability at least $1-2\mathbf{n}\exp\bigg(-\frac{c_{\mbf{\mathcal{H}}^{EA}}^2M}{100\mathbf{n}^2c_1^2+\frac{20}3\cdot c_1\cdot c_{\mbf{\mathcal{H}}^{EA}}\cdot\mathbf{n}} \bigg)$ with $c_1=2K^4\max\{R,R_{\dot{x}}\}^2S_{EA}^2+2$.
Similarly, for large $M$, the smallest singular value of $A_{M}^{\xi}$ satisfies the inequality
\begin{align*}
\sigma_{\min}(A_{M}^{\xi}) \geq  0.8 c_{\mbf{\mathcal{H}}^{\xi}}
\end{align*}
with probability at least $1-2\mathbf{n}_{\xi}\exp\bigg(-\frac{c_{\mbf{\mathcal{H}}^{\xi}}^2M}{100n^2c_1^2+\frac{20}{3}c_2\cdot c_{\mbf{\mathcal{H}}^{\xi}}\cdot n} \bigg)$ with $c_2=2K^4R_{\xi}^2S_{\xi}^2+2$.
\label{p:sigmaminAL}
Therefore, with high probability, a system and its associated hypothesis space satisfying the coercivity condition and $M$ sufficiently large, the inverse problem will be solvable with condition number controlled by the coercivity constant.
\end{proposition}

\begin{proof}
We prove the result in the $EA$ case, the proof of the results about the $\xi$ part of the system is analogous. The orthonormality of the component functions given in \eqref{onb}, implies that
$\langle \mbf{\psi}_{p}^{EA}, \mbf{\psi}_{p'}^{EA}\rangle_{\LtwoB(\bm{\rho}_T^{EA,L})}=\delta_{pp'}$.
Expand $\bm{\varphi}^{EA} \in \mbf{\mathcal{H}}^{EA}$ in this basis as $\bm{\varphi}^{EA} = \sum_{p=1}^{\mathbf{n}} a_p\mbf{\psi}_{p}^{EA}.$ Let the vector $v = (a_1, \ldots, a_{\mathbf{n}}) \in \mathbb{R}^{\mathbf{n}}$, and notice that
\begin{align*}
\frac{1}{L}\sum_{l=1}^{L}\mathbb{E}_{\bm{\mu^{\bY}}} \bigg[ \big\|  \mathbf{f}_{\bvpEA}(\bX(t_l),\bV(t_l),\bXi(t_l)) \big\|_{\mathcal{S}}^2\bigg]  =\dbinnerp {\sum_{p=1}^{\mathbf{n}} a_{p}\mbf{\psi}_{p}^{EA}, \sum_{p=1}^{\mathbf{n}} a_{p} \mbf{\psi}_{p}^{EA}} \\
= v^T A_{\infty}^{EA}v \geq \sigma_{\min}(A_{\infty}^{EA})\|v\|^2
= \sigma_{\min}(A_{\infty}^{EA}) \big\| \bm{\varphi}^{EA}\|_{\LtwoB(\bm{\rho}_T^{EA,L})}^2
\end{align*}
This lower bound is achieved by the singular vector corresponding to the singular value $\sigma_{\min}(A_{\infty}^{EA})$, so that by definition \eqref{2ndorder:gencoer} we have that $\sigma_{\min}(A_{\infty}^{EA}) =c_{\mbf{\mathcal{H}}^{EA}}$. \\

For the second statement, we consider the learning matrix $A^{EA}$ (defined in section \ref{sec:algorithm_numeric}), which we can also write as $A_M^{EA}$ to emphasize the dependence on $M$ as needed built, from the $M$ observed trajectories. By construction, for each $m$, $A_{\infty}^{EA}=\mathbb{E}_{\bm{\mu^{\bY}}}[A^{EA,(m)}]$ and $\lim_{M\rightarrow\infty} A_{M}^{EA}=A_{\infty}^{EA}$ by the Strong Law of Large Numbers.
Next we will derive some important properties of the learning matrix that will allow us to apply the matrix Bernstein inequality (see \cite{tropp2015introduction}, Theorem 6.1.1, Corollary 6.1.2). Note that we will use the notation from this reference.
First we note an elementary matrix analysis result (see \cite{bhatia97} Problem III.6.13); for any two square matrices $A,B$,
$$
\max_j |\sigma_j(A) - \sigma_j(B)| \leq \|A-B \|.
$$
All norms in this proof are the spectral norm, unless otherwise specified.
Thus if we get a concentration inequality of the form $\mathbb{P}_{\bm{\mu^{\bY}}}\{\|A_{\infty}^{EA} - A_{M}^{EA}  \| \geq t \}$ we will get the desired result relating the minimal singular values of $A_{M}^{EA}$ to $c_{\mbf{\mathcal{H}}^{EA}}$. First, notice that $\mathbb{E}_{\bm{\mu^{\bY}}}[A_{M}^{EA}]= A_{\infty}^{EA}$.
Additionally, using the definition of the regression matrix, and our assumptions on the kernels and the dynamics, we can bound every entry by $c_1 = 2K^4\max\{R, R_{\dot{x}} \}^2S_{EA}^2 + 2$. This immediately implies the bound
$$
\|A_{M}^{EA} - A_{\infty}^{EA}\| \leq 2\mathbf{n}c_1.
$$
Next, we upper bound the matrix variance statistic (in our case $Z = \sum_{k=1}^M S_k$ where $S_k = \frac{1}{M}A^{EA,(k)}$), defined as
$$
v(Z) = \max\{ \| \mathbb{E}[(Z - \mathbb{E}Z)(Z - \mathbb{E}Z)^* \|, \| \mathbb{E}[(Z - \mathbb{E}Z)^{*}(Z - \mathbb{E}Z) \| \}.
$$
Using a similar analysis to bound each entry of the matrices, we can arrive at the result that $v(Z) \leq 2\mathbf{n}^2c_1^2$. Now, we apply the matrix Bernstein inequality to see that
$$
\mathbb{P}_{\bm{\mu^{\bY}}}\Big\{\|A_{M}^{EA} - A_{\infty}^{EA} \| \geq t  \Big\} \leq 2\mathbf{n}
\exp \bigg( -\frac{c_{\mbf{\mathcal{H}}^{EA}}^2M}{100\mathbf{n}^2c_1^2 + \frac{20c_1c_{\mbf{\mathcal{H}}^{EA}}}{3}\mathbf{n} }  \bigg).
$$
Note that the $M$ in the numerator comes because $A_{M}^{EA}$ has a factor of $\frac{1}{M}$ on it.
Lastly, choose $t = \frac{c_{\mbf{\mathcal{H}}^{EA}}}{5}$, which together with the results above yield the desired inequality.
\end{proof}

From Proposition \ref{p:sigmaminAL} we see that, for each hypothesis space $\mathcal{H}_{kk'}$, it is important to choose a basis that is well-conditioned in $\LtwoB(\bm{\rho}_T^{EA,L}), \LtwoB(\bm{\rho}_T^{\xi, L})$,  instead of in the corresponding $\bm{L}^{\infty}$ spaces.  If not, the learning matrices, defined in Appendix \ref{sec:algorithm_numeric}, $A^{EA}, A^{\xi}$ may be ill-conditioned or even singular. This would lead to fundamental numerical challenges in solving for the kernels.
In order to mitigate these issues, one can use piecewise polynomials on a partition of the support of the empirical measure and/or use the pseudo-inverse with an adaptive tolerance.

 \subsection{Discussions on the coercivity condition}

The coercivity condition is key to the identifiability of the kernels from data. It is determined by the distribution of the solution to the agent system and introduces constraints on the hypothesis space. For the second-order system, it is therefore related to the distribution $\bm{\mu^{\bY}}$ of the initial conditions, the true interaction kernels, and the non-collective force. The coercivity condition has been studied for first-order systems in \cite{lu2019nonparametric,Tang2019, li2019identifiability}. Below, we give a brief review.

For homogeneous systems, \cite{lu2019nonparametric,Tang2019} showed that the coercivity condition holds true on any compact subset of the corresponding $L^2$ space for the case of $L=1$. This result has been generalized to cover heterogeneous systems in \cite{Tang2019} and a few examples of the stochastic homogeneous system including linear systems and nonlinear systems with stationary distributions for general $L$ in \cite{li2019identifiability}.

In this paper, we shall employ a similar idea as for first-order systems and extend the result to second-order systems. One key in the proof is  to show the positiveness of integral operators that arise in the expectation in Eq. \eqref{2ndorder:gencoer}. We focus on a representative model of second-order homogeneous systems,
\begin{equation}
\begin{dcases}
m_i\ddot\mbf{x}_i &= \force^{\dot\bx}(\mbf{x}_i, \dot\mbf{x}_i, \xi_i) + \sum_{i'=1}^N \frac{1}{N} (\phi^{E}(\norm{\mbf{x}_{i'} - \mbf{x}_i})(\mbf{x}_{i'} - \mbf{x}_i) + \phi^{A}(\norm{\mbf{x}_{i'} - \mbf{x}_i})(\dot\mbf{x}_{i'} - \dot\mbf{x}_i) \\
\dot\xi_i &= \force^{\xi}(\mbf{x}_i, \dot\mbf{x}_i, \xi_i) + \sum_{i'=1}^N \frac{1}{N} \phi^{\xi}(\norm{\mbf{x}_{i'} - \mbf{x}_i})(\xi_{i'} - \xi_i)
\end{dcases}\,
\end{equation}
which includes the first-order systems considered in  \cite{BFHM17,lu2019nonparametric,Tang2019} as special cases and various second-order system examples in \cite{lu2019nonparametric, Zhong20} as specific applications. We shall prove the coercivity condition holds true for the case $L=1$:
\begin{theorem}\label{2ndordersingle:coercivity}
Consider the system \eqref{eq:2ndOrderCoercivity} at time $t_1=0$ with the initial distribution $\mu_0^{\bY}=\begin{bmatrix}\mu_0^{\mbf{X}}\\ \mu_0^{\dot\mbf{X}} \\ \mu_0^{\mbf{\Xi}} \end{bmatrix}$, where ${ \mu}_0^{\mbf{X}}$  is exchangeable Gaussian with $\mathrm{cov}(\mbf{x}_i(t_1))-\mathrm{cov}(\mbf{x}_i(t_1),\mbf{x}_j(t_1))=\lambda I_d$\,  for a constant $\lambda>0$,  ${ \mu}_0^{\dot\mbf{X}}, { \mu}_0^{\mbf{\Xi}}$  are exchangeable with finite second moment, and they are independent of ${ \mu}_0^{\mbf{X}}$.   Then

\begin{align*}
\mathbb{E}_{\mu_0^{\bY}} \|\mathbf{f}_{\combf{\varphi^E}{\varphi^A}} (\mbf{X}(0), \mbf{V}(0))\|^2_{\mathcal{S}} &\geq c_{1,N, \mathcal{H}^{EA}}\vertiii{\combf{\varphi^E}{\varphi^A}}_{L^2(\rho_T^{EA,1})},\\
\mathbb{E}_{\mu_0^{\bY}} \|\mathbf{f}_{\varphi^{\xi}} (\mbf{X}(0), \mbf{\Xi}(0))\|^2_{\mathcal{S}} &\geq c_{1,N, \mathcal{H}^{\xi}}\vertiii{\varphi^{\xi}}_{L^2(\rho_T^{\xi,1})},
\end{align*} where
\begin{itemize}

\item $c_{1,N, \mathcal{H}^{EA}}\geq (\frac{N-1}{2N^2}+ \frac{(N-1)(N-2)}{2N^2}c ),
c=\min\bigg\{c_{\mathcal{H}^{EA}}^{E},  c_{\mathcal{H}^{EA}}^{A}c_{\mu_0^{\dot\mbf{X}}}\bigg\}$, where $c_{{\bm \mu}_0^{\dot\mbf{X}}}=1-\frac{\mathbb{E} \langle \dot \mbf{x}_i(0),\dot \mbf{x}_{i'}(0)\rangle}{\mathbb{E} \|\dot \mbf{x}_i(0)\|^2}$ ($i\neq i'$) and $c_{\mathcal{H}^{EA}}^{E}$ and $c_{\mathcal{H}^{EA}}^{A}$ are non-negative constants independent of $N$, and are strictly positive for compact  $\mathcal{H}^{EA}$ of $L^2(\rho^{EA,1}_T)$.

\item $c_{1,N, \mathcal{H}^{\xi}}\geq (\frac{N-1}{N^2}+ \frac{(N-1)(N-2)}{N^2}c ),
c=  c_{\mathcal{H}^{\xi}}c_{\mu_0^{\mbf{\Xi}}} $ with $c_{{ \mu}_0^{\mbf{\Xi}}}=1-\frac{\mathbb{E} \langle  \xi_i(0), \xi_{i'}(0)\rangle}{\mathbb{E} \| \xi_i(0)\|^2}$ ($i\neq i'$) and  $c_{\mathcal{H}^{\xi}}$ is  a non-negative constant independent of $N$, which is strictly positive for compact  $\mathcal{H}^{\xi}$ of $L^2(\rho^{\xi,1}_T)$.

\end{itemize}
 \end{theorem}

This exhibits a particular case that even the coercivity constant is independent of the number of agents $N$. Consequently, the estimated errors of our estimators are independent of $N$, and therefore not only is the convergence rate of our estimators independent of the dimension $(2d+1)N$ of the phase space, but even the constants in front of the rate term are independent of $N$.  Our results extend those for first-order systems from \cite{lu2019nonparametric, Tang2019, Zhong20}. The empirical numerical experiments on second-order systems, already conducted in \cite{Zhong20}, support that the coercivity condition is satisfied by large classes of second-order systems, and is ``generally'' satisfied for general $L$  on relevant hypothesis spaces, with a constant independent of the number of agents $N$ thanks to the exchangeability of the distribution of the initial conditions, and of the agents at any time $t$.  The proof of the result above is given in Appendix \ref{sec:app:Coercivity}.

\section{Consistency and optimal convergence rate of estimators} \label{sec:main_results}

The final preparatory results for our main theorems combine concentration with a union bound. Here we control the probability that the supremum of the difference between the expected and empirical normalized errors over the whole hypothesis space is large.

\subsection{Concentration}

Our first main result is a concentration estimate that relates the coercivity condition to an appropriate bias-variance tradeoff in our setting.
Let $\mathcal{N}(\mbf{\mathcal{H}},\delta)$ be the $\delta$-covering number, with respect to  the $\infty$-norm, of the set $\mbf{\mathcal{H}}$.

\begin{theorem}[Concentration]\label{2ndorder:maintheorem}
Suppose that ${\bm{\phi}}^{\{E,A,\xi\}} \in \mbf{\mathcal{K}}_{S_{\{E,A,\xi\}}}^{\{E,A,\xi\}}$. Consider a convex, compact (with respect to the $\infty$-norm) hypothesis spaces
$$\mbf{\mathcal{H}}_M^{EA} \subset \mbf{L^\infty}(\mathbf{R} \times \mathbf{S}^E) \oplus \mbf{L^\infty}(\mathbf{R} \times \mathbf{S}^A),\quad \mbf{\mathcal{H}}_M^{\xi} \subset \mbf{L^\infty}(\mathbf{R} \times \mathbf{S}^{\xi})$$
bounded above by $S_0 \geq \max\{S_E, S_A,S_{\xi}\}$ respectively. Additionally, assume that the coercivity condition \eqref{2ndorder:gencoer} holds on $\mbf{\mathcal{H}}_M^{EA}$ and condition \eqref{2ndorder:gencoer:xi} on $\mbf{\mathcal{H}}_M^{\xi}$.

Then for all $\epsilon >0$, with  probability (with respect to $\bm{\mu^{\bY}}$) at least $1-\delta$, we have the estimates
\begin{equation}
\begin{aligned}
\label{2ndorder:mainestimate1}
c_{\mbf{\mathcal{H}}_M^{EA}}\| \widehat{\bm{\phi}}^{EA}_{M} - \bpEA  \|^2_{\LtwoB(\bm{\rho}_T^{EA,L})} &\leq   2\inf_{\bvpEA \in \mbf{\mathcal{H}}_M^{EA}}\|\bvpEA-\bpEA\|^2_{\LtwoB(\bm{\rho}_T^{EA,L})} +2\epsilon,\\
c_{\mbf{\mathcal{H}}_M^{\xi}}\| \widehat{{\bm{\phi}}}_{M}^{\xi} - \bm{\phi}^{\xi}  \|^2_{\LtwoB(\bm{\rho}_T^{\xi, L})} &\leq   2\inf_{\bm{\varphi}^{\xi} \in \mbf{\mathcal{H}}_M^{\xi}}\|\bm{\varphi}^{\xi}-\bm{\phi}^{\xi}\|^2_{\LtwoB(\bm{\rho}_T^{\xi, L})} + 2\epsilon\,,
\end{aligned}
\end{equation}
provided that, for the first bound to hold,
$$M \geq \frac{1152S_{0}^2\max\{R, R_{\dot{x}} \}^2K^4}{\epsilon c_{\mbf{\mathcal{H}}_M^{EA}}}\bigg(\log \Big(\mathcal{N}\Big(\mbf{\mathcal{H}}^{EA}_M,\frac{\epsilon}{48S_{0}\max\{R, R_{\dot{x}} \}^2K^4} \Big)\Big)+\log\Big(\frac{1}{\delta}\Big)  \bigg)\,,$$
and similarly for the second inequality, using $\mbf{\mathcal{H}}_M^{\xi}$.
\end{theorem}

\begin{proof} [\textbf{of Theorem \rm{\ref{2ndorder:maintheorem}}}]
We start out by setting $\alpha=\frac{1}{6}$ in Proposition \ref{learning:fullconcentrationEA}, it is easy to see that this chosen value yields the tightest bound in the argument below. To ease the notation we let
$\widehat {\bm{\phi}}_{L,M,\mbf{\mathcal{H}}^{EA}}^{EA} = \combf{\widehat {\bm{\phi}}_{L,M,\mbf{\mathcal{H}}^{EA}}^E}{\widehat {\bm{\phi}}_{L,M,\mbf{\mathcal{H}}^{EA}}^A}$ and similarly for $\widehat {\bm{\phi}}_{L,\infty,\mbf{\mathcal{H}}^{EA}}^{EA}$.
From the Proposition, we have that
 $$\sup_{ \bvpEA \in \mbf{\mathcal{H}}^{EA}}  \frac{\mathcal{D}_{\infty}(\bvpEA)-\mathcal{D}_{M}(\bvpEA)}{\mathcal{D}_{\infty}(\bvpEA)+\epsilon}  <\frac{1}{2},$$
 holds true with probability
\begin{equation} \label{eq:usefulProbBound}
\mathcal{P}\geq 1-\mathcal{N}\bigg(\mbf{\mathcal{H}}^{EA}, \frac{\epsilon}{48S_{EA}\max\{R, R_{\dot{x}}\}^2K^4}\bigg) \exp \bigg(-\frac{c_{\mbf{\mathcal{H}}^{EA}}M\epsilon}{1152S_{EA}^2\max\{R, R_{\dot{x}}\}^2K^6}\bigg).
\end{equation}
This immediately implies, by choosing $\bvpEA=\widehat {\bm{\phi}}_{L,M,\mbf{\mathcal{H}}^{EA}}^{EA}$ and reorganizing, that with probability $\mathcal{P}$
$$  \mathcal{D}_{\infty}(\widehat {\bm{\phi}}_{L,M,\mbf{\mathcal{H}}^{EA}}^{EA})< 2\mathcal{D}_{M}(\widehat {\bm{\phi}}_{L,M,\mbf{\mathcal{H}}^{EA}}^{EA})+\epsilon\,.$$
By definition of $\widehat {\bm{\phi}}_{L,M,\mbf{\mathcal{H}}^{EA}}^{EA}$ as the minimizer of the empirical error functional $\bm{\mathcal{E}}_{M}^{EA}$, we see that
$$\mathcal{D}_{M}(\widehat {\bm{\phi}}_{L,M,\mbf{\mathcal{H}}^{EA}}^{EA})=\bm{\mathcal{E}}_{M}^{EA}(\widehat {\bm{\phi}}_{L,M,\mbf{\mathcal{H}}^{EA}}^{EA})-\bm{\mathcal{E}}_{M}^{EA}(\widehat {\bm{\phi}}_{L,\infty,\mbf{\mathcal{H}}^{EA}}^{EA}) \leq 0,$$
and combining this result with equation (\ref{eq_minH}) from Proposition \ref{2ndordersystem:convexity}, we have
\begin{equation} \label{eq:concentrationEAintermediate}
c_{\mbf{\mathcal{H}}^{EA}} \|\widehat {\bm{\phi}}_{L,M,\mbf{\mathcal{H}}^{EA}}^{EA}-\widehat {\bm{\phi}}_{L,\infty,\mbf{\mathcal{H}}^{EA}}^{EA}\|_{\LtwoB( \bm{\rho}_T^{EA,L})}^2
 \leq \mathcal{D}_{\infty}(\widehat {\bm{\phi}}_{L,M,\mbf{\mathcal{H}}^{EA}}^{EA} )<\epsilon,
\end{equation}
with the final inequality holding with probability $\mathcal{P}$.
Now we bound the $\LtwoB$ error of the empirical estimator to the true interaction kernel, so that with probability $\mathcal{P}$
 \begin{align*}
 \|\widehat {\bm{\phi}}_{L,M,\mbf{\mathcal{H}}^{EA}}^{EA} - {\bm{\phi}}^{EA}\|_{\LtwoB(\bm{\rho}_T^{EA,L})}^2
 & \leq 2 \|\widehat {\bm{\phi}}_{L,M,\mbf{\mathcal{H}}^{EA}}^{EA}-\widehat {\bm{\phi}}_{L,\infty,\mbf{\mathcal{H}}^{EA}}^{EA}\|_{\LtwoB(\bm{\rho}_T^{EA,L})}^2
+2\|\widehat {\bm{\phi}}_{L,\infty,\mbf{\mathcal{H}}^{EA}}^{EA} - {\bm{\phi}}^{EA}\|_{\LtwoB(\bm{\rho}_T^{EA,L})}^2 \\
&   \leq  \frac{2}{c_{\mbf{\mathcal{H}}^{EA}}}\bigg( \epsilon +\inf_{\bvpEA \in \mbf{\mathcal{H}}^{EA}} K^2\| \bvpEA -\bpEA \|_{\LtwoB(\bm{\rho}_T^{EA,L})}^2\bigg) \\
&  \leq  \frac{2}{c_{\mbf{\mathcal{H}}^{EA}}}\bigg( \epsilon +\inf_{\bvpEA \in \mbf{\mathcal{H}}^{EA}} K^2\max\{R, R_{\dot{x}} \}^2\| \bvpEA -\bpEA \|_{\infty}^2\bigg).
\end{align*}
The first inequality follow from the coercivity condition (\ref{2ndorder:gencoer}) and the definition of
$\widehat{\bm{\phi}}^{EA}_{\infty}$.
The second follows by the definition of the norms.
Now for a chosen $0<\delta<1$, let $$1-\mathcal{N}\bigg(\mbf{\mathcal{H}}^{EA}, \frac{\epsilon}{48S_{EA}\max\{R,R_{\dot{x}} \}^2K^4 }\bigg) \exp \bigg(-\frac{c_{\mbf{\mathcal{H}}^{EA}}M\epsilon}{1152S_{EA}^2\max\{R,R_{\dot{x}} \}^2K^6} \bigg) \geq 1-\delta$$ and solve for $M$. The proof for the $\xi$ part of the system result is similar.
\end{proof}

\subsection{Consistency}  \label{sec:consistency}
In the regime where $M \to \infty$, we will choose an increasing sequence of hypothesis spaces, each satisfying the conditions of Theorem \ref{2ndorder:maintheorem}. By our assumptions on the kernels, we can also choose the sequence of $\mbf{\mathcal{H}}_M^{EA}$'s such that the approximation error goes to $0$ as $M\to \infty$.
This enables us to control the infimum on the right hand side of (\ref{2ndorder:mainestimate1}). From here we can apply Theorem \ref{2ndorder:maintheorem} on each $M$ to prove the consistency of our estimators with respect to the $\LtwoB(\bm{\rho}_T^{EA,L})$ norm and derive the following consistency theorem.

\begin{theorem}[Strong Consistency]\label{main:consistency}  Suppose that
$$\{\mbf{\mathcal{H}}_M^{EA}\}_{M=1}^{\infty} \subset \mbf{L^\infty}(\mathbf{R} \times \mathbf{S}^E) \oplus \mbf{L^\infty}(\mathbf{R} \times \mathbf{S}^A)$$ is a family of compact and convex subsets such that the approximation error goes to zero,
$$\inf_{\bvpEA\in \mbf{\mathcal{H}}_M^{EA}}\|\bvpEA-\bpEA\|_{\infty} \xrightarrow{M\rightarrow \infty} 0.$$
Further suppose that the coercivity condition holds on $\bigcup_{M}\mbf{\mathcal{H}}_M^{EA}$, and that $ \bigcup_{M}\mbf{\mathcal{H}}_M^{EA}$ is compact in  $\mbf{L^\infty}(\mathbf{R} \times \mathbf{S}^E) \oplus \mbf{L^\infty}(\mathbf{R} \times \mathbf{S}^A)$. Then the estimator is strongly consistent with respect to the $\LtwoB(\bm{\rho}_T^{EA,L})$ norm:
$$
\lim_{M\rightarrow \infty}\|\widehat{\bm{\phi}}^{EA}_{M} - \bpEA\|_{\mbf{L}^2(\bm{\rho}_T^{EA,L})} =0 \text{ with probability one.}
$$
An analogous consistency result holds for the estimator in the $\xi$ variable.

\end{theorem}
These two results together provide a consistency result on the full estimation of the triple $(\widehat{{\bm{\phi}}^{\xi}}, \widehat{{\bm{\phi}}^{E}}, \widehat{{\bm{\phi}}^{A}})$ and thus consistency of our estimation procedure on the full system (\ref{eq:2ndOrder}).  \\

\begin{proof}[\textbf{of Theorem \rm{\ref{main:consistency}} }]
To simplify the notation, we use the same conventions as the proof of Theorem \ref{2ndorder:maintheorem} and let
$\mathcal{D}_{\infty} = \mathcal{D}_{L, \infty, \mbf{\mathcal{H}}^{EA}_M}$. By definition of the coercivity constant in (\ref{2ndorder:gencoer}), we have the inequality $c_{\cup_{M}{\mbf{\mathcal{H}}^{EA}_M}} \leq c_{\mbf{\mathcal{H}}^{EA}_M} $.
From an analogous argument used to arrive at equation (\ref{eq:concentrationEAintermediate}) in the proof of Theorem \ref{2ndorder:maintheorem}, we obtain that

\begin{align}\label{importantinequality}
c_{\cup_{M}{\mbf{\mathcal{H}}^{EA}_M}}\|\widehat{\bm{\phi}}^{EA}_M - \bpEA\|_{\LtwoB(\bm{\rho}_T^{EA,L})}^2
&\leq \mathcal{D}_{\infty}(\widehat{\bm{\phi}}^{EA}_M)+\mbf{\mathcal{E}}_{\infty}^{EA}(\widehat{\bm{\phi}}^{EA}_{\infty}).
\end{align}
Let $\epsilon > 0$, the inequality \eqref{importantinequality} gives us that
\begin{align*}
P_{\bm{\mu^{\bY}}} \{c_{\cup_{M}{\mbf{\mathcal{H}}^{EA}_M}} \|\widehat{\bm{\phi}}^{EA}_M - \bpEA\|_{\mbf{L}^2(\bm{\rho}_T^{EA,L})}^2 \geq \epsilon \}
&\leq P_{\bm{\mu^{\bY}}} \{\mathcal{D}_{\infty}(\widehat{\bm{\phi}}^{EA}_M)+\mbf{\mathcal{E}}_{\infty}^{EA}(\widehat{\bm{\phi}}^{EA}_{\infty}) \geq \epsilon \}\\
&\leq P_{\bm{\mu^{\bY}}}\bigg\{\mathcal{D}_{\infty}(\widehat{\bm{\phi}}^{EA}_M)\geq \frac{\epsilon}{2}\bigg\}+
P_{\bm{\mu^{\bY}}}\bigg\{\mbf{\mathcal{E}}_{\infty}^{EA}(\widehat{\bm{\phi}}^{EA}_{\infty}) \geq \frac{\epsilon}{2}\bigg\}.
\end{align*}
We now bound the two terms in the above expression separately. For the first term, the proof of Theorem  \ref{2ndorder:maintheorem} shows that
\begin{align*}
P_{\bm{\mu^{\bY}}}\{\mbf{\mathcal{D}}_{\infty}(\widehat{\bm{\phi}}^{EA}_M)
\geq \frac{\epsilon}{2}\}
&\leq \mathcal{N}\bigg(\mbf{\mathcal{H}}^{EA}_M,\frac{\epsilon}{C_1}\bigg)\exp\bigg(-\frac{c_{\mbf{\mathcal{H}}^{EA}_M}M\epsilon}{C_2}\bigg)\\
&\leq \mathcal{N}\bigg(\cup_{M}{\mbf{\mathcal{H}}^{EA}_M},\frac{\epsilon}{C_1}\bigg)\exp\bigg(-\frac{c_{\cup_{M}\mbf{\mathcal{H}}^{EA}_M}M\epsilon}{C_2}\bigg)
\end{align*}
where $C_1 = 96S_{EA}^2\max\{R, R_{\dot{x}} \}^2K^4$, $C_2 = 2304S_{EA}^2\max\{R, R_{\dot{x}} \}^2K^4$, and
$\mathcal{N}(\cup_{M}{\mbf{\mathcal{H}}^{EA}_M},\frac{\epsilon}{C_1})$ is finite because of the compactness assumption on $\cup_{M}\mbf{\mathcal{H}}^{EA}_M$.

Summing this bound in $M$ we get that,
\begin{align*}
\sum_{M=1}^{\infty} P_{\bm{\mu^{\bY}}}\{\mathcal{D}_{\infty}(\widehat{\bm{\phi}}^{EA}_M) \geq \frac{\epsilon}{2}\}&\leq \mathcal{N}\bigg(\cup_{M}{\mbf{\mathcal{H}}^{EA}_M},\frac{\epsilon}{C_1}\bigg)\sum_{M=1}^{\infty}\exp\bigg(-\frac{c_{\cup_M\mbf{\mathcal{H}}^{EA}_M}M\epsilon}{C_2}\bigg) < \infty.
\end{align*}
For the second term, the bound  \eqref{2ndordersystem:expectationerrorfunctional} yields that
$$\mbf{\mathcal{E}}_{\infty}^{EA}(\widehat{\bm{\phi}}^{EA}_{\infty}) \leq 4K^4S_{EA}\max\{R, R_{\dot{x}}\}^2 \inf_{\bvpEA \in \mbf{\mathcal{H}}^{EA}_M} \|\bvpEA-\bpEA\|_{\infty}\xrightarrow{M\rightarrow \infty} 0.$$
Since $\epsilon$ is fixed, the above result, together with our assumption on the sequence of hypothesis spaces, implies that $P_{\bm{\mu^{\bY}}}\Big\{\mbf{\mathcal{E}}_{\infty}^{EA}(\widehat{\bm{\phi}}^{EA}_{\infty}) \geq \frac{\epsilon}{2}\Big\}=0$ for $M$ sufficiently large. So we have
$\sum_{M=1}^{\infty}P_{\bm{\mu^{\bY}}}\{\mbf{\mathcal{E}}_{\infty}^{EA}(\widehat{\bm{\phi}}^{EA}_{\infty}) \geq \frac{\epsilon}{2}\} < \infty$.
The finiteness of the two sums above implies, by the first Borel-Cantelli Lemma, that
$$P_{\bm{\mu^{\bY}}}\big\{\limsup_{{M\rightarrow \infty}}\{ c_{\cup_{M}{\mbf{\mathcal{H}}^{EA}_M}} \|\widehat{\bm{\phi}}^{EA}_M-\bpEA \|_{\LtwoB(\bm{\rho}_T^{EA,L})}^2  \geq \epsilon\}\big\}=0,$$

As $\epsilon$ was arbitrary, we have the desired strong consistency of the estimator. An exactly analogous argument gives the result on the $\xi$ part of the system.
\end{proof}

\subsection{Rate of convergence} \label{sec:rateofconv}

Given data collected from $M$ trajectories, we would like to choose the best hypothesis space to maximize the accuracy of the estimators. Theorem \ref{2ndorder:maintheorem} highlights the classical bias-variance tradeoff in our setting. On the one hand, we would like the hypothesis space $\mbf{\mathcal{H}}^{EA}_M$ to be large so that the bias $$\inf_{\bm{\varphi}^{EA} \in \mbf{\mathcal{H}}^{EA}_M}\|\bm{\varphi}^{EA} -\bm{\phi}^{EA} \|^2_{\LtwoB(\bm{\rho}_T^{EA,L})}\,,\text{ or }\,\,\inf_{\bm{\varphi}^{EA} \in \mbf{\mathcal{H}}^{EA}}\|\bm{\varphi}^{EA}-\bm{\phi}^{EA}\|^2_{\infty}\,,$$ is small. Simultaneously,  we would like $\mbf{\mathcal{H}}^{EA}_M$  to be small enough so that the covering number
$\mathcal{N}(\mbf{\mathcal{H}}^{EA}_M,\epsilon) $
is small. Just as in nonparametric regression, our rate of convergence depends on a regularity condition of the true interaction kernels and properties of the hypothesis space, as is demonstrated in the following theorem. We establish the optimal (up to a log factor) min-max rate of convergence by choosing a hypothesis space in a sample size dependent manner.

\subsection*{Comments}
Parts (a) and (c) of the theorem concern an approximation theory type rate of convergence where $M$ plays no role in the choice of hypothesis space, whereas parts (b) and (d) of the theorem present a minimax rate of convergence that chooses an adaptive hypothesis space depending on $M$ to achieve the optimal rate.

The splitting of the convergence result, between $EA$ and $\xi$ parts, emphasizes a common theme of the paper:
we leverage that the system can be decoupled for the learning process to improve the rate of convergence and performance in learning the estimators, but analytically we study it as the full coupled system, see the trajectory prediction result in Theorem \ref{secondordersystem:TrajDiff_NEW}.

A final important comment is that even though the dimension of the space in which we measure the error may be large, namely the dimension of $\bm{\rho}_T^{EA,L}$ (which can be calculated as $\dim(\bm{\rho}_T^{EA,L}) = \sum_{kk'}p_{(k,k')}^E + \sum_{(k,k')}p_{(k,k')}^A$), we exploit the structure of the system in such a way that our convergence rate only depends on the number of unique variables across all of the estimators. This number is given by the number of variables in $\bm{\mathcal{V}}$, denoted as $|\bm{\mathcal{V}}|$, for the $EA$ portion and by $|\bm{\mathcal{V}}_{\xi}|$ for the $\xi$ portion of the system.
We note that we are not predicting the number of variables nor their form, they are assumed known, this applies to both the pairwise interaction variables and the feature maps.

\begin{theorem}[Rate of Convergence] \label{main:convrate}
Let $\widehat{\bm{\phi}}^{EA}:=\combf{\widehat{\bm{\phi}}_M^E}{\widehat{\bm{\phi}}_M^A}$
denote the minimizer of the empirical error functional $\bm{\mathcal{E}}_{M}^{EA}$ (defined in
\eqref{eq:Error_Func_Empirical:EA}) over the hypothesis space $\mbf{\mathcal{H}}^{EA}_M$.  \\
(a)  Let the hypothesis space be chosen as the direct sum of the admissible spaces, namely $\mbf{\mathcal{H}}^{EA} = \combf{\mbf{\mathcal{K}}_{S_E}^E}{\mbf{\mathcal{K}}_{S_A}^A},$ and assume that the coercivity condition \eqref{2ndorder:gencoer} holds true on it.

Then, there exists a constant $C$ depending only on $K,S_{EA},R, R_{\dot{x}}$ such that
 \[
 \mathbb{E}_{\bm{\mu^{\bY}}}\Big[\| \widehat{\bm{\phi}}^{EA}_M-\bm{\phi}^{EA}\|^2_{\LtwoB(\bm{\rho}_T^{EA,L})} \Big]\leq \frac{C}{c_{\mbf{\mathcal{H}}^{EA}}} M^{-\frac{1}{|\bm{\mathcal{V}}|+1}}.
 \]
(b) Assume that $\{\mbf{\mathcal{L}}_n\}_{n=1}^{\infty}$  is a sequence of finite-dimensional linear subspaces of  $\mbf{L}^{\infty}(\mathbf{R}\times \mathbf{S}^E) \oplus \mbf{L}^{\infty}(\mathbf{R} \times \mathbf{S}^A)$  satisfying the dimension and approximation constraints
 \begin{align}\label{assumptions}
\text{dim}(\mbf{\mathcal{L}}_n) \leq c_0K^2n^{|\bm{\mathcal{V}}|}\,,\quad \inf_{\bm{\varphi}^{EA} \in \mbf{\mathcal{L}}_n}\|\bm{\varphi}^{EA}-\bm{\phi}^{EA}\|_{\infty}  \leq c_1 n^{-s},
\end{align}
for some fixed constants $c_0,c_1$ representing dimension-independent approximation characteristics of the linear subspaces, and $s>0$ related to the regularity of the kernels. The value $n$ can be thought of as the number of basis functions along each of the $|\bm{\mathcal{V}}|$ axes. Suppose the coercivity condition holds true on the set $\cup_n\mbf{\mathcal{L}}_n$.  Define $\mbf{\mathcal{B}}_n$ to be the closed ball centered at the origin of radius $(c_1+S_{EA})$ in $\mbf{\mathcal{L}}_{n}$.  Let $c_M^{EA} := c_{\bigcup_M \mbf{\mathcal{H}}^{EA}_M}$.
If we choose the hypothesis space as $\mbf{\mathcal{H}}_M=\mbf{\mathcal{B}}_{(\frac{M}{\log M})^{\frac{1}{2s+|\bm{\mathcal{V}}|}}}$, then there exists a constant $C$ depending on $K,R, R_{\dot{x}}, S_{EA},c_0, c_1, s$ such that we achieve the convergence rate,
\begin{align}\label{rate}
 \mathbb{E}_{\bm{\mu^{\bY}}}\Big[\|
\widehat{\bm{\phi}}^{EA}_{M} -\bm{\phi}^{EA} \|^2_{\LtwoB(\bm{\rho}_T^{EA,L})} \Big] \leq \frac{C}{c_M^{EA}} \left(\frac{\log M}{M}\right)^{\frac{2s}{2s+|\bm{\mathcal{V}}|}}\,.
 \end{align}
(c) under the corresponding assumptions as in (a), there exists a constant $C$ depending only on $K,S_{\xi},R$ such that
 \[
 \mathbb{E}_{\bm{\mu^{\bY}}}\Big[\|\widehat{\bm{\phi}}_M^{\xi}-{\bm{\phi}}^{\xi}\|^2_{\LtwoB(\bm{\rho}_T^{\xi, L})} \Big]\leq \frac{C}{c_{\mbf{\mathcal{H}}^{\xi}}} M^{-\frac{1}{|\bm{\mathcal{V}}_{\xi}|+1}}.
 \]
(d) under the corresponding assumptions as in (b), there exists a constant $C$ depending only on $K,R, S_{\xi},c_0, c_1, s$ such that, and for $c_M^{\xi}:=c_{\bigcup_M \mbf{\mathcal{H}}^{\xi}_M}$,
\begin{align}\label{rate:xi}
 \mathbb{E}_{\bm{\mu^{\bY}}}\Big[\|
\widehat{\bm{\phi}}^{\xi}_{M}
 -{\bm{\phi}}^{\xi} \|^2_{\LtwoB(\bm{\rho}_T^{\xi, L})} \Big] \leq \frac{C}{c_M^{\xi}} \left(\frac{\log M}{M}\right)^{\frac{2s}{2s+|\bm{\mathcal{V}}_{\xi}|}}\,. \end{align}
 \label{t:2ndordersystem:thm_optRate}
\end{theorem}
We in fact prove bounds not only in expectation, but also with high probability, for every fixed large-enough $M$, as the proof will show. \\

\begin{proof}[\textbf{of Theorem \rm{\ref{main:convrate}}}]
For part (a), let $\mbf{\mathcal{H}}=\combf{\mbf{\mathcal{K}}_{S_E}^E}{\mbf{\mathcal{K}}_{S_A}^A}$.
Standard results on covering numbers of function spaces (see
theorem 2.7.1 of \cite{Wellner}) give us that the covering number of $\bm{\mathcal{H}}$ satisfies
$$
\mathcal{N}(\mbf{\mathcal{H}}, \epsilon, \|\cdot \|_{\infty}) \leq C_{\mbf{\mathcal{H}}}\exp\bigg(K^2\bigg(\frac{1}{\epsilon}\bigg)^{|\bm{\mathcal{V}}|}\bigg)
$$
for some absolute constant $C_{\mbf{\mathcal{H}}}$ depending only on $\mbf{\mathcal{H}}$ and $|\bm{\mathcal{V}}|$. By assumption on the hypothesis space, we have that
$$
\inf_{\bm{\varphi}^{EA} \in \mbf{\mathcal{H}}}\|\bm{\varphi}^{EA}-\bm{\phi}^{EA}\|^2_{\infty}=0.
$$
From this, the concentration estimate \eqref{2ndorder:mainestimate1} together with the covering number bound imply that,
\begin{align}
P_{\bm{\mu^{\bY}}}\{\|
\widehat{\bm{\phi}}^{EA}_{L,M,\mbf{\mathcal{H}}}
- \bm{\phi}^{EA} \|_{\LtwoB(\bm{\rho}_T^{EA,L})}^2
> \epsilon \}
&\leq  \mathcal{N}(\mbf{\mathcal{H}}, C_1 \epsilon, \|\cdot \|_{\infty})\exp(-C_2M\epsilon) \nonumber \\
& \leq C_{\mbf{\mathcal{H}}} \exp( K^2(C_1\epsilon)^{-|\bm{\mathcal{V}}|}  - C_2M\epsilon) \label{eq:maintheorem_intResult}
\end{align}
where $C_1 = \frac{c_{\mbf{\mathcal{H}}}}{48S_{EA}\max\{R, R_{\dot{x}} \}^2K^4}$ and $C_2  = \frac{c_{\mbf{\mathcal{H}}}}{1152S_{EA}^2\max\{R, R_{\dot{x}} \}^2K^4}$.
Next, define the function
$$g(\epsilon) := K^2(C_1\epsilon)^{-|\bm{\mathcal{V}}|} - \frac{C_2M\epsilon}{2},$$ which we will minimize to achieve the desired probability bound.

By direct calculation, $g(\epsilon)=0$ if we choose $\epsilon=\epsilon_{M}=(\frac{C_3}{M})^{\frac{1}{|\bm{\mathcal{V}}|+1}}$, where $C_3=\Big(\frac{2K^2}{C_2C_1^{|\bm{\mathcal{V}}|}}\Big)^{\frac{1}{|\bm{\mathcal{V}}|+1}}$.  It is then an easy computation to see that the derivative of $g(\epsilon)$ is $\leq 0$ for all $\epsilon \geq \epsilon_M$. Thus, we can put this result into the bound \eqref{eq:maintheorem_intResult} to arrive at the probability bound,
\begin{equation} \label{eq:maintheorem_intResult2}
P_{\bm{\mu^{\bY}}}\{\|
\widehat{\bm{\phi}}^{EA}_{L,M,\mbf{\mathcal{H}}} -
\bm{\phi}^{EA} \|_{\LtwoB(\bm{\rho}_T^{EA,L})}^2
> \epsilon \} \leq  \begin{cases} \exp\Big(\frac{-C_2M\epsilon}{2}\Big), &  \epsilon \geq \epsilon_M \\ 1,  &  \epsilon \leq \epsilon_M \end{cases}
\end{equation}
Integrating this bound over $\epsilon\in(0,+\infty)$ and using the elementary inequality $e^{-x}\leq x+1$ for all $x\geq 0$, we get that
$$
\int_0^{\infty} P_{\bm{\mu^{\bY}}}\{\|
\widehat{\bm{\phi}}^{EA}_{L,M,\mbf{\mathcal{H}}} -
\bm{\phi}^{EA} \|_{\LtwoB(\bm{\rho}_T^{EA,L})}^2
> \epsilon \}  d\epsilon \leq \Big(\frac{C_4}{M}\Big)^{\frac{1}{|\bm{\mathcal{V}}|+1}} + O\Big(\frac{1}{M}\Big)
$$
Now, bringing the coercivity part from \eqref{2ndorder:mainestimate1} back in, we achieve the rate,
$$
\mathbb{E}_{\bm{\mu^{\bY}}}[\|
\widehat{\bm{\phi}}^{EA}_{L,M,\mbf{\mathcal{H}}}-
\bm{\phi}^{EA}\|_{\LtwoB(\bm{\rho}_T^{EA,L})}^2]  \le \frac{C_4}{c_{\mbf{\mathcal{H}}}} M^{-\frac{1}{|\bm{\mathcal{V}}|+1}}, $$ where $C_4$ is an absolute constant that only depends on $K, S_{EA}, R, R_{\dot{x}}$.

For part (b), we note the following basic result on the covering number of $\mbf{\mathcal{B}}_n$  by $\epsilon$-balls (see \cite[Proposition 5]{CS02}),
$$
\mathcal{N}(\bm{\mathcal{B}}_n, \epsilon, \|\cdot \|_{\infty}) \leq \bigg(\frac{4(c_1+S_{EA})}{\epsilon}\bigg)^{c_0K^2n^{|\bm{\mathcal{V}}|}}.
$$
Using \eqref{2ndorder:mainestimate1}, and the approximation assumption, we bound the probability as
\begin{equation}
 \begin{aligned}
P_{\bm{\mu^{\bY}}}\{\|
\widehat{\bm{\phi}}^{EA}_{L,M,\bm{\mathcal{B}}_n}
- \bm{\phi}^{EA}\|_{\LtwoB(\bm{\rho}_T^{EA,L})}^2
> \epsilon + c_2n^{-2s} \} \\
& \hspace{-7cm} = P_{\bm{\mu^{\bY}}}\{\|
\widehat{\bm{\phi}}^{EA}_{L,M,\bm{\mathcal{B}}_n} - \bm{\phi}^{EA}\|_{\LtwoB(\bm{\rho}_T^{EA,L})}^2
> t'n^{-2s} + c_2n^{-2s} \} \\
& \hspace{-7cm} = P_{\bm{\mu^{\bY}}}\{\|
\widehat{\bm{\phi}}^{EA}_{L,M,\bm{\mathcal{B}}_n} - \bm{\phi}^{EA} \|_{\LtwoB(\bm{\rho}_T^{EA,L})}^2
> tn^{-2s} \} \\
& \hspace{-7cm} \leq \mathcal{N}\Big(\bm{\mathcal{B}}_n, c_3'tn^{-2s} , \|\cdot \|_{\infty}\Big)\exp(-c_4Mtn^{-2s}) \\
& \hspace{-7cm} \leq \Big(\frac{c_3}{tn^{-2s}}\Big)^{c_0K^2n^{|\bm{\mathcal{V}}|}} \exp(-c_4Mtn^{-2s})\\
& \hspace{-7cm} \leq \exp(c_0K^2n^{|\bm{\mathcal{V}}|}\log(c_3)+c_0K^2n^{|\bm{\mathcal{V}}|}|\log(tn^{-2s})|-c_4Mtn^{-2s}),
\label{est2}
\end{aligned}
\end{equation}
where $c_2=\frac{1}{c_{\cup_n\mbf{\mathcal{L}}_n}}c_1$, $ c_3' = \frac{c_{\cup_n\mbf{\mathcal{L}}_n}}{48(S_{EA}+c_1)\max\{R, R_{\dot{x}}\}^2K^4}$, $c_3=\frac{192 (S_{EA}+c_1)^2\max\{R, R_{\dot{x}}\}^2K^4}{c_{\cup_n\mbf{\mathcal{L}}_n}}$, and $c_4=\frac{c_{\cup_n\mbf{\mathcal{L}}_n}}{1152(S_{EA}+c_1)^2\max\{R, R_{\dot{x}}\}^2K^4}$ are absolute constants independent of $M$.
Define
$$
g(n) :=c_0n^{|\bm{\mathcal{V}}|}K^2\log(c_3)+c_0n^{|\bm{\mathcal{V}}|}K^2|\log(tn^{-2s})|-\frac{c_4}{2}Mtn^{-2s}.
$$
To find the optimal $n$ in terms of $M$, we minimize $g$ in $n$. By taking a derivative, and solving the corresponding equation, we see that the optimal $n$ is
$$n_* = O\bigg(\Big(\frac{M}{\log M}\Big)^{\frac{1}{2s+|\bm{\mathcal{V}}|}}\bigg),$$ with constant independent of $M$ and only depending on $c_3, c_4, c_2$.
For convenience we will choose $n_*= \lfloor(\frac{M}{\log M})^{\frac{1}{2s+|\bm{\mathcal{V}}|}}\rfloor$. Now let $\epsilon_M = (\frac{M}{\log M})^{\frac{2s}{2s+|\bm{\mathcal{V}}|}}$ and  consider
$$h(\epsilon) = c_0n_*K^2\log(c_3) + c_0n_*K^2|\log(\epsilon)| - \frac{c_4}{2}M\epsilon.$$
As before, let $\epsilon = tn_*^{-2s} = t\epsilon_M$ and consider $h(t\epsilon_M)$. It is easy to see that $\lim_{t \rightarrow 0^+} h(t\epsilon_M) = \infty$ and $\lim_{t \rightarrow \infty}h(t \epsilon_M)=-\infty$. Together with the continuity of $h$, these facts imply that there exists a constant $c_5$, depending on $K,c_0, c_2,c_3,c_4$ such that $h(c_5\epsilon_M) = 0$.
We further need that $h'(\epsilon) \leq 0$ for all $\epsilon \geq c_5\epsilon_M$. By taking the derivative of $h$, setting it $\leq 0$, we find that this condition eventually holds by basic calculus on $h$.
Therefore, if needed to satisfy the derivative condition, we can enlarge the constant $c_5$ to a constant $c_6$ (independent of $M$) such that $h(\epsilon) \leq 0$ and $ h'(\epsilon) \leq 0$ for all $\epsilon \geq c_6\epsilon_M$.
These results imply the probability bound,
\begin{equation} \label{eq:maintheorem_intResult2}
P_{\bm{\mu^{\bY}}}\{\|
\widehat{\bm{\phi}}^{EA}_{L,M,\mbf{\mathcal{B}}_{n_*}}
- \bm{\phi}^{EA} \|_{\LtwoB(\bm{\rho}_T^{EA,L})}^2
> \epsilon \} \leq  \begin{cases} \exp\Big(\frac{-c_4M\epsilon}{2}\Big), &  \epsilon \geq c_6\epsilon_M \\ 1,  &  \epsilon \leq c_6\epsilon_M \end{cases}
\end{equation}
Integrating this bound over $\epsilon\in(0,+\infty)$ and using the elementary inequality $e^{-x}\leq x+1$ for all $x\geq 0$, we get that
$$
\int_0^{\infty} P_{\bm{\mu^{\bY}}}\{\|
\widehat{\bm{\phi}}^{EA}_{L,M,\mbf{\mathcal{B}}_{n_*}} - \bm{\phi}^{EA}\|_{\LtwoB(\bm{\rho}_T^{EA,L})}^2
> \epsilon \}  d\epsilon \leq C_1\Big(\frac{\log M}{M}\Big)^{\frac{2s}{2s+|\bm{\mathcal{V}}|}},
$$
where $C_1$ is a constant depending on $c_0,c_1,s, K, S_{EA}, R, R_{\dot{x}}$.
Now with $\mbf{\mathcal{H}}^{EA}_M = \mbf{\mathcal{B}}_{n_*}$ and using \eqref{2ndorder:mainestimate1}, we have shown the convergence rate,
$$
\mathbb{E}_{\bm{\mu^{\bY}}}[\|
\widehat{\bm{\phi}}^{EA}_{L,M,\mbf{\mathcal{H}}^{EA}_M}
- \bm{\phi}^{EA} \|_{\LtwoB(\bm{\rho}_T^{EA,L})}^2]  \le \frac{c_7}{c_{\bigcup_M \mbf{\mathcal{H}}^{EA}_M}} \Big(\frac{M}{\log M}\Big)^{-\frac{2s}{2s+|\bm{\mathcal{V}}|}},
$$ where $c_7$ is an absolute constant that only depends on $s,K,c_0, c_1, S_{EA}, R, R_{\dot{x}}$.

\end{proof}

In both theorems, the convergence rates $\frac{2s}{2s+|\bm{\mathcal{V}}|}$ and $\frac{2s}{2s+|\bm{\mathcal{V}}_{\xi}|}$ coincide with the minimax rate of convergence for nonparametric regression in the corresponding dimension -- up to the logarithmic factor. It is possible that this logarithmic factor could be removed (see techniques in Chapter 11-15 of \cite{Gyorfi06}), but with considerable additional complexity of the proof.
Achieving the same rate of convergence as if we had observed the noisy values of the interaction kernels directly, rather than through the dynamics, is a major strength of our approach. The strong consistency results show the asymptotic optimality of our method, and for wide classes of systems the assumptions of the theorems apply.
Specifically, for part (b) of the theorems, the dimension and approximation conditions can be explicitly achieved by piecewise polynomials or splines appropriately adapted to the regularity of the kernel. In the conditions of theorem \ref{t:2ndordersystem:thm_optRate}, $n$ can be the number of partitions along each axis of the variables in $\bm{\mathcal{V}}$. Then, using multivariate splines or piecewise polynomials we will have a fixed constant $c_0$ (corresponding to the number of parameters to estimate for each function) times $Kn^{\bm{\mathcal{V}}}$ as the dimension of the linear space. Furthermore, by standard approximation theory results, see \cite{Schumaker} (Chapters 12,13), \cite{DeVore1980},\cite{DeBoor1983}, for $s$ the regularity of the interaction kernels we achieve the desired approximation condition with piecewise polynomials of degree $\lfloor s\rfloor$.
In our admissible spaces we have $s=1$, note that the theorems are stronger if we have a kernel of higher regularity.

We next briefly examine the convergence rate on a few systems of fundamental interest.
Recall that in table \ref{tab:2ndOrder_Examples} we have as the final two columns the values $|\bm{\mathcal{V}}|, |\bm{\mathcal{V}}_{\xi}|$. These correspond directly to the rate of convergence of each of the system under our learning approach.

Some specific highlights:
\begin{itemize}
\item For Anticipation Dynamics (AD), even though we are learning both an energy and alignment kernel, because there are only $2$ unique variables shared across both of them we learn at the $2$-dimensional rate.
\item For the Synchronized Oscillator we achieve the $2$-dimensional optimal learning rate on each of the $EA$ and $\xi$ portions (rather than a $4$-dimensional rate) due to the decoupled nature of the system, similarly we only pay the $1$-dimensional rate twice for the Phototaxis system. This is a key reason for splitting our learning theory between $EA$- and $\xi$-interaction kernels and accounting for shared and non-shared variables: this enables us to substantially improve the performance guarantees in actual applications.
\item Due to the design of the measures, norms and the associated learning algorithm, even in the heterogeneous case for celestial mechanics and predator-swarm, we only pay the $1$-dimensional learning rate, although the constants are of course affected by the heterogeneity and the algorithm requires a larger learning matrix.
\item The rates of convergence of our estimators for all previously-studied first-order systems (see \cite{lu2019nonparametric,Tang2019,Zhong20}) can be derived from Theorem \ref{main:convrate}.
\end{itemize}

One downside of the results above is the lack of dependence on $L$ as it seems natural that finer time samples in each trajectory should improve the results. Indeed, the numerical experiments of \cite{Zhong20, lu2019nonparametric, Tang2019} demonstrate that more data in $L$ may indeed be helpful to improve the performance. One technique used in \cite{Zhong20} for very long trajectory data (large $L$, medium to small $M$) is to split each trajectory into larger $M$ with smaller $L$ in each.
In this way, one can explicitly show the desired convergence rate in a form agreeing with the above theorems, we do not believe that it leads to significantly different performance compared to using the original data, only that we can easily transform the data into a form that the theorems apply to and with no loss of performance.  Explicit dependence on $N$ is not the objective of this work, see \cite{BFHM17}, but further study of the mean-field regime is of interest to the authors and work is ongoing.

\section{Performance of trajectory prediction} \label{sec:trajectory}
Once estimators $\widehat{\bm{\varphi}}^{EA},\widehat{\bm{\varphi}}^{\xi}$ are obtained, a natural question is the accuracy of the evolved trajectory based on these estimated kernels. The next theorem shows that the error in prediction is (i) bounded trajectory-wise by a continuous-time version of the error functional, and (ii) bounded on average by the $\LtwoB(\bm{\rho}_T^{EA}), \LtwoB(\bm{\rho}_T^{\xi})$, respectively, error of the estimator. This further validates the effectiveness of our error functional and $\LtwoB(\bm{\rho}_T)$-metric to assess the quality of the estimator.
In particular, this emphasizes that although the system is a coupled system of ODE's, our decoupled learning procedure with our choice of norm will lead to control of the expected supremum error as long as we minimize the $\LtwoB(\bm{\rho}_T^{EA}), \LtwoB(\bm{\rho}_T^{\xi})$ norms in obtaining our estimators.

\begin{theorem}\label{secondordersystem:TrajDiff_NEW}
Suppose that $\widehat{\bm{\phi}}^E \in \boldsymbol{\mathcal{K}}_{S_E}^E$, $\widehat{\bm{\phi}}^A \in \boldsymbol{\mathcal{K}}_{S_A}^A$ and $\widehat{\bm{\phi}}^{\xi} \in \boldsymbol{\mathcal{K}}_{S_{\xi}}^{\xi}$.
Denote by $\widehat{\bY}(t)$ and $\bY(t)$
the solutions of the systems with kernels $\widehat{\bm{\phi}}^E=(\widehat\phi_{kk'}^{E})_{k,k'=1}^{K,K}, \widehat{\bm{\phi}}^A =(\widehat\phi_{kk'}^{A})_{k,k'=1}^{K,K}$, and $\widehat{\bm{\phi}}^{\xi} =(\widehat\phi_{kk'}^{\xi})_{k,k'=1}^{K,K}$ and ${\bm{\phi}}^E, {\bm{\phi}}^A,{\bm{\phi}}^{\xi}$ respectively, both with the same initial condition. Then
\begin{align*}
\sup_{t\in[0,T]}\Vert \widehat{\bY}(t)- \bY(t)\Vert_{\mathcal{Y}}^2 &\leq g(T) \bigg[2T^2 \int_{p=0}^t \int_{s=0}^p \| \ddot{\mbf{X}} - \rhsf^{\text{nc}, \dot\bx}(\bX,\bV,\bXi) - \mathbf{f}^{\widehat{\bm{\phi}}^{EA}}(\bX,\bV,\bXi) \|_{\mathcal{S}}^2 \text{ ds dp} \\
&+ 2T  \int_{s=0}^t  \| \ddot{\mbf{X}} - \rhsf^{\text{nc}, \dot\bx}(\bX,\bV,\bXi) - \mathbf{f}^{\widehat{\bm{\phi}}^{EA}}(\bX,\bV,\bXi) \|_{\mathcal{S}}^2 ds \\
&+ 2T\int_{s=0}^t \| \dot{\mbf{\Xi}} - \rhsf^{\text{nc}, \xi}(\bX,\bV,\bXi) - \mathbf{f}^{\widehat{{\bm{\phi}}}_{\xi}}(\bX,\bV,\bXi) \|_{\mathcal{S}}^2 ds \bigg]
\end{align*}
where $g(T) =1+(1+B_1T)T\exp(A_1T+T^2/2)$. The constants are $A_1 = 2T(8KP + \mathcal{L} +8QK + \mathcal{L}^{\xi})$ and $B_1 = 2T^2(8KP + \mathcal{L})$, with any unspecified constants made precise in the proof and only depending on the Lipschitz constants of the noncollective forces and the feature maps, as well as the values $S^E, S^A, S^{\xi}$ coming from the admissible spaces.
It is bounded on average, with respect to the initial distribution $\bm{\mu^{\bY}}$, by
\begin{align} \label{trajPred:KeyExpectation}
\mathbb{E}_{\bm{\mu^{\bY}}}[\sup_{t\in[0,T]} \|\widehat{\bY}(t)- \bY(t)\|_{\mathcal{Y}}^2] &\leq g(T)\bigg((T^2K^2 + TK^2)\Vert \widehat{\bm{\phi}}^{EA} - \bm{\phi}^{EA} \Vert_{\LtwoB(\bm{\rho}_T^{EA})}^2 \nonumber \\
&+ TK^2\Vert \widehat{\bm{\phi}}^{\xi} - \bm{\phi}^{\xi} \Vert_{\LtwoB(\bm{\rho}_T^{\xi})}^2 \bigg)
\end{align}
with the measures $\bm{\rho}_T^{EA},\bm{\rho}_T^{\xi}$ defined in (\ref{meas:fullEA}, \ref{meas:rhoxibar}). Expression (\ref{trajPred:KeyExpectation}) shows that by minimizing the right hand side, we can control the expected $\mathcal{Y}$-supremum error of the estimated trajectories.

\end{theorem}
We postpone the somewhat lengthy proof to Appendix \ref{s:trajectoryerrorproof}.

\section{Applications} \label{sec:numerics}
Our learning theory, as well as measures, norms, functionals etc. can be applied to study all the examples considered in the works \cite{lu2019nonparametric, Tang2019, Zhong20}. These examples, particularly those of \cite{Zhong20}, can thus be considered as applications of the theoretical results as well as of the algorithm in section \ref{sec:algorithm_numeric}.

We choose to study two new dynamics, which are not considered in \cite{lu2019nonparametric, Zhong20} since they exhibit some unique features of our generalized model.  In particular, we choose them due to their special form of having both energy-based and alignment-based interactions.  These are the flocking with external potential (FwEP) model in \cite{ST2020} and the anticipation dynamics (AD) model in \cite{shu2019anticipation}.

Table \ref{tab:app_params} shows the value of learning parameters for these dynamics.
\begin{table}[H]
\centering
\begin{tabular}{ c | c | c | c | c | c | c }
\hline
$M_{\rho}$ & $L$   & $T_f$ & $T$ & $\mu^{\mbf{x}}$             & $\mu^{\dot\mbf{x}}$          & Num. of learning trials\\
\hline
$2000$     & $500$ & $10$  & $5$ & $\text{Unif.}([0, 5]^2)$ & $\text{Unif.}([0, 5]^2)$ & $10$\\
\hline
\end{tabular}
\caption{Values of parameters for the learning.}
\label{tab:app_params}
\end{table}

The setup of the learning experiment is as follows.  We use $M_{\rho}$ different initial conditions to evolve the dynamics\footnote{We use the built-in MATLAB integrating routine, $ode15s$, with relative tolerance at $10^{-8}$ and absolute tolerance at $10^{-11}$.} from $0$ to $T$ to generate a good approximation to $\rho_T^{L, EA}, \rho_T^{L, E}$ and $\rho_T^{L, A}$.  Then we use another set of $M$ ($M = 500$ for FwEP and $M = 750$ for AD) initial conditions to generate training data to learn the corresponding $\phi^{E}$ and $\phi^{A}$ from the empirical distributions, $\rho_T^{L, M, EA}$'s, etc.  We report the relative learning errors calculated via \eqref{e:L2NormDefs} for $\widehat{\intkernel}^{E} \oplus \widehat{\intkernel}^{A}$, \eqref{e:L2NormDefs} for $\widehat{\intkernel}^{E}$, and \eqref{e:L2NormDefs} for $\widehat{\intkernel}^{A}$, along with pictorial comparison of those interaction kernels as well as a visualization on the pairwise data which is used to learn the estimated kernels.  Then we evolve the dynamics either from the training set of $M$ initial conditions or another set of $M$ randomly chosen initial conditions with $\phi^{E} \oplus \phi^{A}$ and $\widehat{\intkernel}^{E} \oplus \widehat{\intkernel}^{A}$ from $0$ to $T_f > T$, and report the trajectory errors calculated using \eqref{eq:traj_norm} on $\mbf{y}$ (the whole system), and for $\mbf{x}$ (the position) and $\mbf{v}$ (the velocity).  Again, pictorial comparison of the trajectories are also shown.  We report the trajectory errors over $[0, T]$ and $[T, T_f]$.  The learning results are shown in the following sections.

\subsection{Learning results for flocking with external potential}
We consider the FwEP model for its simplicity and clustering behavior in both position and velocity (hence flocking occurs).  The dynamics of the FwEP model is given as follows,
\begin{equation}
\ddot\mbf{x}_i = \frac{1}{N}\sum_{i' = 1, i' \neq i}^N a(\dot\mbf{x}_{i'} - \dot\mbf{x}_i) \nonumber + \frac{1}{N}\sum_{i' = 1, i' \neq i}^N\phi(\norm{\mbf{x}_{i'} - \mbf{x}_i})(\dot\mbf{x}_{i'} - \dot\mbf{x}_i).\label{eq:FwEP_model}
\end{equation}
Here $a > 0$ is a constant representing an attraction force, and $\phi = \frac{1}{(1 + r^2)^\beta}$\footnote{This choice of interaction for alignment is not mandatory.  It is a good comparison between this FwEP model with the Cucker-Smale mode with this choice of interaction functions.} with $\beta = \frac{1}{2}$.  To fit into our learning regime, we take, $m_i = 1$, $K = 1$, no $\xi_i$, no non-collective force, and
\[
\phi^{E} = a \quad \text{and} \quad \phi^{A} = \frac{1}{(1 + r^2)^\beta}.
\]
For the FwEP model, we use the space of $1^{st}$ degree piece-wise polynomials with dimension $n^E = 122$ for learning $\phi^{E}$; and for $\phi^{A}$, we use the same space.  First, consider the comparison of energy-based interactions shown in Fig. \ref{fig:FwEP_phiE}.
\begin{figure}[H]
  \centering
  \includegraphics[width=.8\linewidth]{./Figures/FwEP_phiE}
  \caption{$\phi^{E}$ vs. $\widehat{\intkernel}^{E}$, Err: $3.9 \cdot 10^{-6} \pm 6.6 \cdot 10^{-7}$.  The lines shown in blue are the estimated interaction kernels, and the lines shown in black are the true interaction kernels. The colored areas shown in the background are the learned distributions of pairwise distance data.}
  \label{fig:FwEP_phiE}
\end{figure}
Fig. \ref{fig:FwEP_phiE} shows that our learning performance on constant functions using piecewise linear polynomials shows promising results.  However, we still have trouble learning the behavior of the interaction at $r = 0$, part of it due to the weight of $\vec{0}$ in the model, and the other part of it being lack of available data towards $r = 0$.  Next, we show the comparison of alignment-base interactions shown in Fig. \ref{fig:FwEP_phiA_results} with distribution of the pairwise data.
\begin{figure}[H]
\centering
\begin{subfigure}{.5\textwidth}
  \centering
  \includegraphics[width=.9\linewidth]{./Figures/FwEP_phiA}
  \caption{$\phi^{A}$ vs. $\widehat{\intkernel}^{A}$, Err: $9.1 \cdot 10^{-3} \pm 2.5 \cdot 10^{-4}$.}
  \label{fig:FwEP_phiA}
\end{subfigure}%
\begin{subfigure}{.5\textwidth}
  \centering
  \includegraphics[width=.9\linewidth]{./Figures/FwEP_rhoA}
  \caption{$\rho_{T}^{A,L}$ vs. $\rho_{T}^{A,L, M}$.}
  \label{fig:FwEP_rhoA}
\end{subfigure}
\caption{The lines shown in blue are the estimated interaction kernels, and the lines shown in black are the true interaction kernels.
The colored areas shown in the background are the learned distributions of pairwise distance data.}
\label{fig:FwEP_phiA_results}
\end{figure}
Again, in Fig. \ref{fig:FwEP_phiA}, it shows a faithful approximation from our estimated kernels.  The $\widehat{\intkernel}^{E} \oplus \widehat{\intkernel}^{A}$ error is: $5.8 \cdot 10^{-3} \pm 1.6 \cdot 10^{-4}$.  The comparison of trajectories are shown in Fig. \ref{fig:FwEP_trajs}.
\begin{figure}[H]
  \centering
  \includegraphics[width=.8\linewidth]{./Figures/FwEP_trajs}
  \caption{Trajectory Comparison.}
  \label{fig:FwEP_trajs}
\end{figure}
Fig. \ref{fig:FwEP_trajs} shows little visual difference between the learned and observed trajectories.  A more quantitative description of the trajectory errors are shown in table \ref{tab:FwEP_traj_err}.
\begin{table}[H]
\centering
\begin{tabular}{ c | c | c }
\hline
                                   & $[0, T]$                                   & $[T, T_f]$ \\
\hline
$\text{mean}_{\text{IC}}$ on $\mbf{x}$ & $7.2 \cdot 10^{-4} \pm 1.9 \cdot 10^{-5}$  & $6.7 \cdot 10^{-4} \pm 1.8 \cdot 10^{-5}$ \\
$\text{mean}_{\text{IC}}$ on $\mbf{v}$ & $1.15 \cdot 10^{-3} \pm 3.1 \cdot 10^{-5}$ & $1.5 \cdot 10^{-3} \pm 4.1 \cdot 10^{-3}$ \\
$\text{mean}_{\text{IC}}$ on $\mbf{y}$ & $6.2 \cdot 10^{-6} \pm 1.7 \cdot 10^{-7}$  & $2.22 \cdot 10^{-6} \pm 6.8 \cdot 10^{-8}$ \\
\hline
$\text{std}_{\text{IC}}$ on $\mbf{x}$  & $1.28 \cdot 10^{-4} \pm 4.2 \cdot 10^{-6}$ & $1.22 \cdot 10^{-4} \pm 4.0 \cdot 10^{-6}$ \\
$\text{std}_{\text{IC}}$ on $\mbf{v}$  & $2.20 \cdot 10^{-4} \pm 7.0 \cdot 10^{-6}$ & $2.5 \cdot 10^{-4} \pm 1.0 \cdot 10^{-5}$ \\
$\text{std}_{\text{IC}}$ on $\mbf{y}$  & $1.52 \cdot 10^{-6} \pm 5.9 \cdot 10^{-8}$ & $6.0 \cdot 10^{-7} \pm 2.6 \cdot 10^{-8}$ \\
\hline
\hline
$\text{mean}_{\text{IC}}$ on $\mbf{x}$ & $7.2 \cdot 10^{-4} \pm 1.7 \cdot 10^{-5}$  & $6.7 \cdot 10^{-4} \pm 1.6 \cdot 10^{-5}$ \\
$\text{mean}_{\text{IC}}$ on $\mbf{v}$ & $1.15 \cdot 10^{-3} \pm 2.7 \cdot 10^{-5}$ & $1.46 \cdot 10^{-3} \pm 3.3 \cdot 10^{-5}$ \\
$\text{mean}_{\text{IC}}$ on $\mbf{y}$ & $6.2 \cdot 10^{-6} \pm 1.6 \cdot 10^{-7}$  & $2.22 \cdot 10^{-6} \pm 5.4 \cdot 10^{-8}$ \\
\hline
$\text{std}_{\text{IC}}$ on $\mbf{x}$  & $1.30 \cdot 10^{-4} \pm 5.8 \cdot 10^{-6}$ & $1.24 \cdot 10^{-4} \pm 5.3 \cdot 10^{-6}$ \\
$\text{std}_{\text{IC}}$ on $\mbf{v}$  & $2.25 \cdot 10^{-4} \pm 9.9 \cdot 10^{-6}$ & $2.56 \cdot 10^{-4} \pm 9.1 \cdot 10^{-6}$ \\
$\text{std}_{\text{IC}}$ on $\mbf{y}$  & $1.6 \cdot 10^{-6} \pm 7.6 \cdot 10^{-8}$  & $6.2 \cdot 10^{-7} \pm 2.4 \cdot 10^{-8}$ \\
\hline
\end{tabular}
\caption{Trajectory Errors.  The first three rows of mean trajectory errors are from the training set of initial conditions.  The next three rows are standard deviation of the trajectory errors from the training set of initial conditions.  The following three rows are mean trajectory errors from a new set of initial conditions.  Finally, the last three rows report the standard deviation of the trajectory errors from a new set of initial conditions.}
\label{tab:FwEP_traj_err}
\end{table}
We are maintaining a relative four-digit accuracy in estimating the position, and a relative three-digit accuracy in estimating the velocity of the agents in the system.  Although we are able to reconstruct $\phi^{E}$ with a $6$-digit accuracy, we are not able to do the same for $\phi^{A}$.  The error in $\widehat{\intkernel}^{E} \oplus \widehat{\intkernel}^{A}$ reflects this discrepancy by considering the two functions together.
\subsection{Learning results for anticipation dynamics with $U(r) = \frac{r^p}{p}$}\label{sec:example_AD}
The energy-based interactions are constants in the FwEP models, if we want to consider more complicated models, i.e., interactions depending on pairwise distance and more, the AD models are suitable candidates.  The dynamics of the AD model is given as follows,
\begin{align}
\ddot\mbf{x}_i &= \frac{1}{N}\sum_{i' = 1, i' \neq i}^N \frac{\tau U'(\norm{\mbf{x}_{i'} - \mbf{x}_i})}{\norm{\mbf{x}_{i'} - \mbf{x}_i}}(\dot\mbf{x}_{i'} - \dot\mbf{x}_i) \nonumber \\
& \quad + \frac{1}{N}\sum_{i' = 1, i' \neq i}^N\Big\{\frac{-\tau U'(\norm{\mbf{x}_{i'} - \mbf{x}_i})(\mbf{x}_{i'} - \mbf{x}_i)\cdot(\dot\mbf{x}_{i'} - \dot\mbf{x}_i)}{\norm{\mbf{x}_{i'} - \mbf{x}_i}^3} \nonumber \\
& \quad + \frac{\tau U''(\norm{\mbf{x}_{i'} - \mbf{x}_i})(\mbf{x}_{i'} - \mbf{x}_i)\cdot(\dot\mbf{x}_{i'} - \dot\mbf{x}_i)}{\norm{\mbf{x}_{i'} - \mbf{x}_i}^2} + \frac{U'(\norm{\mbf{x}_{i'} - \mbf{x}_i})}{\norm{\mbf{x}_{i'} - \mbf{x}_i}}\Big\}(\mbf{x}_{i'} - \mbf{x}_i). \label{eq:AD_model}
\end{align}
Compared to the original model in \cite{shu2019anticipation}, we take $\tau_{i, i'} = 0$.  In order to fit the model into our learning regime, we take
\[
\phi^{A}(r) = \frac{\tau U'(r)}{r} \quad \text{and} \quad \phi^{E}(r, s) = \frac{-\tau U'(r)s}{r^3} + \frac{\tau U''(r)s}{r^2} + \frac{U'(r)}{r}.
\]
Here we have no $\xi_i$, $K = 1$, $m_i = 1$, and
\[
s^E_{i, i'} = s^A_{i, i'} = (\mbf{x}_{i'} - \mbf{x}_i)\cdot(\dot\mbf{x}_{i'} - \dot\mbf{x}_i).
\]
We also use $\tau = 0.1$.

It is shown in \cite{shu2019anticipation} that if $U''$ is bounded when $r \rightarrow \infty$ with $U(0) = U'(0) = 0$, then unconditionally flocking would occur.  We take $U(r) = \frac{r^p}{p}$ for $1 < p \le 2$, then the system would show unconditional flocking.  We choose $p = 1.5$ for our learning trials\footnote{$p = 2$ induces constant forces on the dynamics.}.  We use a tensor grid of $1^{st}$ degree piece-wise standard polynomials with $n^E = 28^2$ for learning $\phi^{E}(r, s)$, then a set of $1^{st}$ degree piecewise standard polynomials with $n^A = 138$ for learning $\phi^{A}(r)$.  For the energy-based interactions we have the following results.
\begin{figure}[H]
\centering
\begin{subfigure}{.5\textwidth}
  \centering
  \includegraphics[width=.9\linewidth]{./Figures/AD01_phiE}
  \caption{$U(r) = \frac{r^{1.5}}{1.5}$: $\phi^{E}$ vs. $\widehat{\intkernel}^{E}$, Err: $6 \cdot 10^{-1} \pm 2.6 \cdot 10^{-1}$.}
  \label{fig:AD01_phiE}
\end{subfigure}%
\begin{subfigure}{.5\textwidth}
  \centering
  \includegraphics[width=.9\linewidth]{./Figures/AD01_rhoE}
  \caption{$U(r) = \frac{r^{1.5}}{1.5}$: $\rho_T^{E,L}$ vs. $\rho_T^{E,L, M}$.}
  \label{fig:AD01_rhoE}
\end{subfigure}
\caption{The lines shown in blue are the estimated interaction kernels, and the lines shown in black are the true interaction kernels.
The colored areas shown in the background are the learned distributions of pairwise distance data.}
\label{fig:AD01_phiE_results}
\end{figure}

As is shown in Fig. \ref{fig:AD01_rhoE}, the concentration of pairwise distance data is away from $0$, making the estimation of the behavior of $\phi^{E}(r, s)$ at $r$ close to $0$ extremely difficult, meanwhile, since $\phi^{E}$ is also weighted by the pairwise difference, $\mbf{x}_{i'} - \mbf{x}_i$, and at $r_{i, i'}$ close to $0$, the information is also lost.  Next, we present the alignment-based interaction kernels.

\begin{figure}[H]
\centering
\begin{subfigure}{.5\textwidth}
  \centering
  \includegraphics[width=.9\linewidth]{./Figures/AD01_phiA}
  \caption{$U(r) = \frac{r^{1.5}}{1.5}$: $\phi^{A}$ vs. $\widehat{\intkernel}^{A}$. Err: $1.7 \cdot 10^{-1} \cdot 3.9 \cdot 10^{-2}$.}
  \label{fig:AD01_phiA}
\end{subfigure}%
\begin{subfigure}{.5\textwidth}
  \centering
  \includegraphics[width=.9\linewidth]{./Figures/AD01_rhoA}
  \caption{$U(r) = \frac{r^{1.5}}{1.5}$: $\rho_T^{A,L}$ vs. $\rho_T^{A,L, M}$.}
  \label{fig:AD01_rhoA}
\end{subfigure}
\caption{The lines shown in blue are the estimated interaction kernels, and the lines shown in black are the true interaction kernels.
The colored areas shown in the background are the learned distributions of pairwise distance data.}
\label{fig:AD01_phiA_results}
\end{figure}
We have less trouble estimating the behavior of $\phi^{A}$ at $r = 0$.  We have trouble estimating $\phi^{A}$ at the other end of the spectrum of $r$, since the agents have aligned their velocities, hence the weight $\mbf{v}_{i'} -  \mbf{v}_i$ is close to a zero vector.  The overall learning performance for estimating $\phi^{A}$ is better compared to estimating $\phi^{E}$.  The $\widehat{\intkernel}^{E} \oplus \widehat{\intkernel}^{A}$ error is: $6 \cdot 10^{-1} \pm 3.0 \cdot 10^{-1}$.  The comparison of trajectories between the true kernels (LHS) and the estimators (RHS) is shown below.
\begin{figure}[H]
  \centering
  \includegraphics[width=.8\linewidth]{./Figures/AD01_trajs}
  \caption{$U(r) = \frac{r^{1.5}}{1.5}$: Trajectory Comparison.}
  \label{fig:AD01_traj}
\end{figure}
Visually, there is no difference between the true dynamics and the estimated dynamics.  We offer more quantitative insight into the difference between the two in table \ref{tab:AD01_traj_err}.
\begin{table}[H]
\centering
\begin{tabular}{ c | c | c }
\hline
                                   & $[0, T]$                                   & $[T, T_f]$ \\
\hline
$\text{mean}_{\text{IC}}$ on $\mbf{x}$ & $2.22 \cdot 10^{-3} \pm 8.5 \cdot 10^{-5}$ & $2.4 \cdot 10^{-3} \pm 1.0 \cdot 10^{-4}$ \\
$\text{mean}_{\text{IC}}$ on $\mbf{v}$ & $7.8 \cdot 10^{-3} \pm 3.9 \cdot 10^{-4}$  & $1.78 \cdot 10^{-2} \pm 8.9 \cdot 10^{-4}$ \\
$\text{mean}_{\text{IC}}$ on $\mbf{y}$ & $1.91 \cdot 10^{-5} \pm 9.0 \cdot 10^{-7}$ & $6.2 \cdot 10^{-6} \pm 3.6 \cdot 10^{-7}$ \\
\hline
$\text{std}_{\text{IC}}$ on $\mbf{x}$  & $2.8 \cdot 10^{-4} \pm 1.2 \cdot 10^{-5}$  & $3.4 \cdot 10^{-4} \pm 1.5 \cdot 10^{-5}$ \\
$\text{std}_{\text{IC}}$ on $\mbf{v}$  & $1.14 \cdot 10^{-3} \pm 7.8 \cdot 10^{-5}$ & $2.7 \cdot 10^{-3} \pm 1.5 \cdot 10^{-4}$ \\
$\text{std}_{\text{IC}}$ on $\mbf{y}$  & $4.7 \cdot 10^{-6} \pm 3.0 \cdot 10^{-7}$  & $6.2 \cdot 10^{-6} \pm 3.6 \cdot 10^{-7}$ \\
\hline
\hline
$\text{mean}_{\text{IC}}$ on $\mbf{x}$ & $2.22 \cdot 10^{-3} \pm 8.4 \cdot 10^{-5}$ & $2.4 \cdot 10^{-3} \pm 1.0 \cdot 10^{-4}$ \\
$\text{mean}_{\text{IC}}$ on $\mbf{v}$ & $7.8 \cdot 10^{-3} \pm 3.8 \cdot 10^{-4}$  & $1.78 \cdot 10^{-2} \pm 8.7 \cdot 10^{-4}$ \\
$\text{mean}_{\text{IC}}$ on $\mbf{y}$ & $1.91 \cdot 10^{-5} \pm 8.6 \cdot 10^{-7}$ & $2.4 \cdot 10^{-5} \pm 1.1 \cdot 10^{-6}$ \\
\hline
$\text{std}_{\text{IC}}$ on $\mbf{x}$  & $3.0 \cdot 10^{-4} \pm 2.4 \cdot 10^{-5}$  & $3.4 \cdot 10^{-4} \pm 1.3 \cdot 10^{-5}$ \\
$\text{std}_{\text{IC}}$ on $\mbf{v}$  & $1.15 \cdot 10^{-3} \pm 6.8 \cdot 10^{-5}$ & $2.7 \cdot 10^{-3} \pm 1.5 \cdot 10^{-4}$ \\
$\text{std}_{\text{IC}}$ on $\mbf{y}$  & $4.7 \cdot 10^{-6} \pm 2.6 \cdot 10^{-7}$  & $6.2 \cdot 10^{-6} \pm 3.1 \cdot 10^{-7}$ \\
\hline
\end{tabular}
\caption{$U(r) = \frac{r^{1.5}}{1.5}$: Trajectory Errors.}

\label{tab:AD01_traj_err}
\end{table}
We maintain a $3$-digit relative accuracy in estimating the position/velocity of the agents, even though for the interaction kernels, we are only able to maintain a $1$-digit relative accuracy.
\section{Conclusion and further directions} \label{sec:conclusion}
We have described a second-order model of interacting agents that incorporates multiple agent types, an environment, external forces, and multivariable interaction kernels. The inference procedure described exploits the structure of the system to achieve a learning rate that only depends on the dimension of the interaction kernels, which is much smaller than the full ambient dimension $(2d+1)N$. %
Our estimators are strongly consistent, and in fact have learning rates that are min-max optimal within the nonparametric class, under mild assumptions on the interaction kernels and the system. We described how one can relate the expected supremum error of the trajectories for the system driven by the estimated interaction kernels to the difference between the true interaction kernels and the estimated ones -- this result gives strong support to the use of our weighted $L^2$ norms as the correct way to measure performance and derive estimators. A detailed discussion of the full numerical algorithm, including the inverse problem derived from data and a coercivity condition to ensure learnability, along with complex examples, were presented and we showed how the formulation presented covers a very wide range of systems coming from many disciplines.

There are various ways that one could build on this work to handle different systems and for many of these further directions, the theoretical framework, techniques, and theorems presented here would be directly useful. In particular, one could consider second-order stochastic systems or a similar system but on a manifold, more complex environments, having more unknowns within the model beyond just the interaction kernels (say estimating the non-collective forces as well),
identifying the best feature maps to model the data, and considering semiparametric problems where there are hidden parameters within the interaction kernels or other parts of the model that we wish to estimate along with the interaction kernels.
The generality of the model and its broad coverage of models across the sciences, together with the scalability and performance of the algorithm, could inspire new models -- both explicit equations and nonparametric estimators learned from data -- which are theoretically justified and highly practical.

\acks{MM is grateful for discussions with Fei Lu and Yannis Kevrekidis, and for partial support from NSF-1837991, NSF-1913243, NSF-1934979, NSF-Simons-2031985, AFOSR-FA9550-17-1-0280 and FA9550-20-1-0288, ARO W911NF-18-C-0082, and to the Simons Foundation for the Simons Fellowship for the year '20-'21; ST for support from an AMS Simons travel grant; JM for support from NIH - T32GM11999. Please direct correspondence to any of the first three authors.

All authors jointly designed research and wrote the manuscript; JM and ST derived theoretical results; MZ developed algorithms and applications; JM and MZ analyzed data.}