{'fraction_non_alphanumeric': 0.04065181890830757, 'fraction_numerical': 0.011407989066861239, 'mean_word_length': 4.2997790326540635, 'pattern_counts': {'":': 0, '<': 18, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 20, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Learning Theory for Inferring Interaction Kernels in Second-Order Interacting Agent Systems

Modeling the complex interactions of systems of particles or agents is a fundamental scientific and mathematical problem that is studied in diverse fields, ranging from physics and biology, to economics and machine learning. In this work, we describe a very general second-order, heterogeneous, multivariable, interacting agent model, with an environment, that encompasses a wide variety of known systems. We describe an inference framework that uses nonparametric regression and approximation theory based techniques to efficiently derive estimators of the interaction kernels which drive these dynamical systems. We develop a complete learning theory which establishes strong consistency and optimal nonparametric min-max rates of convergence for the estimators, as well as provably accurate predicted trajectories. The estimators exploit the structure of the equations in order to overcome the curse of dimensionality and we describe a fundamental coercivity condition on the inverse problem which ensures that the kernels can be learned and relates to the minimal singular value of the learning matrix. The numerical algorithm presented to build the estimators is parallelizable, performs well on high-dimensional problems, and is demonstrated on complex dynamical systems.


Introduction
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 2-5 includes a very large number of relevant systems and allows for their modeling. Clustering of opinions [45,27,10,56] is a simple first-order case that exhibits clustering. Flocking of birds [32,29,28] provides a simple example of a second-order system that exhibits an emergent shared velocity of all agents. Milling of fish [25,1,2,24] 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 ξ, was studied by [73,59,58,57]. There are also models that include both energy and alignment interaction kernels, a particular case of this is the anticipation dynamics model from [71], 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 [83] and is further studied in the upcoming [54]. All of these models fit into our framework and we have presented detailed studies of these, and others in this work as well as [83,51,52]. These dynamics exhibit a wide range of emergent behaviors, and as shown in [79,74,29,38,23,5,56], 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 ( [4,68,46]), modeling dynamical systems ( [37,6,47]), governing equations ( [20,80]), biology ( [21]), fluid mechanics ( [63,41]), many-body problems in quantum systems ( [19]), mean-field games ( [67]), meteorology ( [40]), and dynamical systems ( [26,17,81,11]). 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 [12,69]. 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 ( [16,66,15]), a LASSO-type penalty ( [42,43]), and sparse Bayesian regression ( [82]). Other approaches consider multiscale methods, statistical mechanics, or force-based models, see [3,8]. Deep learning has also been applied to learn dynamical systems, for ODEs see [62,65] and for PDEs see [60,61,50], as well as the references therein.
The majority of the earliest work in inferring interaction kernels in systems of the type (1.1), (2.2) 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 [53,44]. 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. [78,14,49,18,64]. The question of identifiability of the parameter emerges, see e.g. [34,55]. Our work is closely related to this viewpoint but our parameter is now infinite-dimensional, with identifiability discusses in section 6.
Our learning approach is based on exploiting the structure of collective dynamical systems and nonparametric estimation techniques (see [30,76,39,36,9]). We focus here on second-order models and the form of the equations, generalizing the first order models (see discussion in Appendix F), derived from Newton's second law: for i = 1, . . . , N m iẍi (t) = Fẋ(x i ,ẋ i ) Here, m i is the mass of the i th agent, x i is its position, Fẋ is a non-collective force, and φ E : R + → R, φ A : R + → R 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 [13], 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 [52] and learning theory on these first-order models was developed in [51,48]. Further extension of the model and algorithm to more complicated second-order systems, with particular emphasis on emergent collective behaviors, was discussed in [83]. A big data application to real celestial motion ephemerides is developed and discussed in [54]. In this work, we provide a rigorous learning theory covering the models presented in [83], as well as the second-order models introduced in [52]. We consider generalizations of the models in [83], to include models with higher-dimensional interaction kernels, that do not depend only on pairwise distances. Compared to the theories studied in [51,48], our theory focuses on second-order models with interaction kernels of the form φ E (r)r + φ A (r)ṙ (with r andṙ representing norms of differences of positions and, respectively, velocities of pairs of agents); additionally, we discuss the identifiability and separability of φ E and φ 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 (2.2), to learn the underlying interaction functions. This variational approach was initially developed in [13,52] and further studied and extended in [51,83]. 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 (2.2) to move the inference task to just the unknown functions (φ E , φ A , φ ξ ) 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 4 and for the full theorem see section 7.2, which for the model (1.1) yields, where ρ EA,L T is a dynamics-adapted measure on pairwise distances -and we use a weighted L 2 space detailed in section 5, particularly (5.3). 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 |V| in the interaction kernels. See section 7.3 for the full theorem, (see section 4 for relevant definitions) which is given by: In the case of model (1.1), |V| = 1, as in the results for first-order systems [13,52,51]. 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 µ Y throughout this work, see (section 7).
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 6. 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 [83], 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 9.
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 2-5. If one wishes to quickly jump to the theoretical sections, and then refer back to the definitions as needed, we have provided table 1,3 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 6-8) 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 4.3, with full details in the corresponding sections. Lastly, we consider applications in section 9, as well as have many additional proofs and details in appendices A-H.

Comparison with existing work
Our learning approach discovers the governing structure of a particular subset of dynamical systems of the form, from observations {Y (t l ),Ẏ (t l )} L l=1 , by implicitly inferring the right hand side, F φ EA ,φ ξ . The main difficulties in establishing an effective theory of learning F φ EA ,φ ξ are the curse of dimensionality caused by the size of Y , which is often N (2d + 1), where N is the number of agents, d the dimension of physical space; and the dependence of the observation data, for example Y (t l+1 ) depends on Y (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 [51], but many of the approaches above successfully address this. Our work, however, exploits the interactingagent 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 (1.1). 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 φ E , φ 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 φ E 's considered in [52,51]. The current theoretical framework is not able to conclusively show that φ E and φ A can be learned separately; however we demonstrate in various numerical experiments that by learning φ E and φ A jointly, we still achieve strong performance. Finally, we note that the first-order theory developed in [51] is a special case of our second-order theory, see details in appendix F.

Model description
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, Here U is a potential energy depending on pairwise distance. From the Lagrange equation, d dt ∂ẋ i L = ∂ x i L, we obtain a simple second-order collective dynamics model Here, φ E (r) = U (r) r represents an energy-based interaction between agents. For example, taking U (r) = N Gm 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 k=1 , with N k being the number of agents of type k. They interact according to the system of ODEs , · · · , K} are the indices of the agent type of the agents i and i respectively; the interaction kernels, φ E kk , φ A kk , φ ξ kk , 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 s E ii (t), s A ii (t), s ξ ii (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 φ E ; everything else is assumed given. Table 1 gives a detailed explanation for the definition of the variables used in (2.2). We note that in what follows the notation, {E, A, ξ}, attached to any expression means that there is one of those maps, functions, etc. for each element in the set {E, A, ξ}. It is a convenient way to avoid excessive repetition of similar definitions.

Variable
Definition i, i index of agent, from 1, . . . , N m i mass of agent i x i (t),ẋ i (t),ẍ i (t) ∈ R d position/velocity/acceleration vector of agent i at time t ξ i ,ξ i auxiliary variable, and its derivative · Euclidean norm in R d K number of types of agents k, k index of type of agents, from 1, . . . , K N k number of agents in type k k i type index of agent i C k set of indices for type k agent, subset of {1, . . . , N } φ kk (or φ kk ) interaction kernel: influence of agent k (or k ) on agent k Fẋ, F ξ non-collective forces onẍ i andξ i respectively φ E , φ A , φ ξ energy, alignment, and environment-based interaction kernels respectively The specific instances of the feature map V together with corresponding projections π {E,A,ξ} kk 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  Table 2: Summary of the models studied in this work and in [52,51,83,54] function V is Lipschitz and known and so are all the π {E,A,ξ} kk 's. The Lipschitz assumption is necessary for us to control the trajectory error and ensure the well-posedness of the system. The function 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, ξ 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 s {E,A,ξ} (k,k ) 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, s The models encompassed by the form (2.2) are quite diverse. For a concrete example, please see section 9.2. We summarize many examples in table 2 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, V, V ξ used within the EA or ξ portions of the system. The number of these unique variables specifies the dimension in the minimax convergence rate, see section 7.3.
Our second-order model equations cover the first-order models considered in [52,51,83] as special cases, see Appendix F, 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 (2.2), even when written as a first-order system in more variables, is a strict generalization of the previous first-order analysis.

Rate of convergence notation
For the system (2.2), 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 V denote the set of features in the range of V that are arguments of φ E or φ A across each of the (k, k ) pairs. For many collective dynamical systems V = {r}. In the case where the system has both φ E and φ A but V = {r}, it is easy to see from Theorem 9 below that we still only pay the 1dimensional rate. Analogously, we define V ξ to be the set of all features in the range of V that are arguments involved in the φ ξ part of (2.2) across all of the (k, k ) pairs. This notation is used to frame the convergence rate theorem on the ξ variable as well.

Preliminaries and notation
We vectorize the models given in (2.2) in order to give them a more compact description.
Letting v i (t) =ẋ i (t), we take the following notations We introduce a weighted norm to measure the system variables, denoted · S(·) and given by, Here · 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, 1 , 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 (2.2) becomes, Here m = m 1 , . . . , m N T ∈ R N , • 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,...,K ): We also use the shorthand, to denote the element of the direct sum of the function spaces containing φ E , φ 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.

Trajectory Performance Measurement
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 T be the trajectory from dynamics generated by the true/unknown interaction kernels with initial condition, Y (0); and T be the trajectory from dynamics generated by the estimated interaction kernels learned from observation of {Y (t l )} L l=1 with the same initial condition, Y (0) (i.e.,Ŷ (0) = Y (0)). We define a norm for the difference between Y andŶ at time t l : and a corresponding norm on the trajectory Y [0,T ] = {Y (t l )} L l=1 (0 = t 1 < · · · < t L = T ): (3.5) We also consider a relative version, invariant under changes of units of measure: Lastly, we report errors between X [0,T ] andX [0,T ] ,

Function spaces
We begin by describing some basic ideas about measures and function spaces. Consider a compact or precompact set U ⊂ R p for some p, then we define the infinity norm as, Further define, L ∞ (U) as the space of real valued functions defined on U with finite ∞-norm.
A key function space we need to consider is, C k,α c (U) for k ∈ N, 0 < α ≤ 1, defined as the space of compactly supported, k times continuously differentiable functions with a k-th derivative that is Hölder continuous of order α. We can then consider vectorizations of these spaces as which has the vectorized infinity norm given by Similarly, we consider direct sums of measures, with corresponding vectorized function spaces. This is done explicitly in section 5.1 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 (2.2). 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, for k, k = 1, . . . , K where recall that the {E, A, ξ} notation means that there is an admissible space (or more generally an expression) for each element of the set {E, A, ξ}.
Here, R min kk , R max kk are the minimum or maximum, respectively, possible interaction radii for agents in C k influencing agents in C k . Similarly, S E kk , S A kk , S ξ kk are compact sets in R p E kk , R p A kk , R p ξ kk which contain the ranges of the feature maps, s E kk , s A kk and s ξ 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 In order to provide uniform bounds, we introduce the following sets: Next, we introduce notation to bound the interaction radii on the pairwise distances and pairwise velocities.

Remark 1
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.

We let
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 µ Y . 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 (2.2), 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 x i (t) which are used in the second supremum.
We assume in this work that both of these quantities Rẋ, R ξ are finite. This will be easily satisfied if the measures µẋ, µ ξ (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 (2.2) to be well-posed, we assume that the interaction kernels lie in admissible sets. : for all k, k = 1, . . . , K , When estimating the EA part of the system, we will consider the direct sum admissible space, for S EA ≥ max{S E , S A }, The admissibility assumptions allow us to establish properties such as existence and uniqueness of solutions to (2.2) 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 5.1.
In the learning approach, we will consider hypothesis spaces that we will search in order to estimate the various interaction kernels.
(3.14) Analogous to our simplified notation for φ EA , ϕ EA described in (3.3), we define the direct sum of the hypothesis spaces as, We will consider specific hypothesis spaces during the learning theory and numerical algorithm sections.

Inference problem and learning approach
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 [52] and [83].

Problem setting
Our observation data is given by For simplicity, we only consider equidistant observation points: t l − t l−1 = h for l = 2, . . . , 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. φ E , φ A , φ ξ , whose initial conditions Y (m) (0) are drawn i.i.d from µ Y , a probability measure defined on the space R N (2d+2) . The goal is to infer the unknown interaction kernels directly from data.

Loss functionals
Given observations, the references [52,51,83] proposed an empirical error functional, recalling the notational convention (3.3), The estimators of interaction kernels are defined as the minimizers of the error functionals E EA M and E ξ M over H EA and H ξ respectively, i.e.
For the learning theory, we will consider the following error functionals. On the energy and alignment portion, we consider, Similarly, on the ξ portion, we consider, (4.5) 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, The relationship between the theoretical and empirical error functionals will play a key role in the learning theory.

Overview of theoretical contributions
The papers [83,52,51] 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 [51] as a special case, as discussed in Appendix F. We focus on the regime where L is fixed but M → ∞. We provide a learning theory that answers the fundamental questions: • 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 5.
• 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 6.
• 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 7.
• 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 8.

Hypothesis Space and Algorithm
First, we choose finite-dimensional subspaces of H E kk , i.e.,H E kk ⊂ H E kk , 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 [52]); similarly forH A kk ⊂ H A . Hence, each test function, ϕ E kk , ϕ A kk , can be expressed in terms of the linear combination of the basis functions as follows Substituting this linear combination back into (4.1), we obtain a system of linear equations, The overhead memory storage needed is M LN (d(5 + n EA + n ξ ) + 3), with M LN (4d + 2) needed for trajectory data, M LN d(n EA + n ξ (here n EA = n E + n A , the sum of the number of basis functions on E and A) for learning matrices, and M LN (d + 1) for right hand side vectors. Hence if M O(1), we can consider parallelization in m in order to reduce the overhead memory, ending up with M per core LN (d(5+n EA +n ξ )+3) with M per core = M ncores . The final storage of A and b only needs n EA (n EA + 1) + n ξ (n ξ + 1).
The computational cost of for solving the system 2s+|V | , which is slightly super linear in M , but less than quadratic in the common situation when 2s + |V | ≥ 2|V |.
Similar analysis on solving A ξ α ξ = b ξ also shows that the computational cost is slightly super linear in M .

Learning theory
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.

Probability measures
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,ṙ, s E , s A , s ξ ) 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 [51,52,13]. For each interacting pair (k, k ), we introduce the following probability measures, for k = k , and the dirac measures are defined as We use a superscript (m) to denote that the variable is calculated from the data from that m th trajectory. For example, s E,(m) ii (t) denotes the feature maps between agents i and i at time t along the m th trajectory.
The measure ρ EA,L,k,k T is the discrete counterpart of ρ EA,k,k T at the observation time instances. In practice, we can use ρ EA,L,M,k,k T to approximate ρ EA,L,k,k T since it can be computed from observational data and will converge to ρ EA,L,k,k T as M → ∞. We also consider the marginal distributions (r,ṙ, s A ) defined analogously as above. The empirical measures, ρ E,L,M,k,k T , ρ A,L,M,k,k T , are the ones used in the actual algorithm to quantify the learning performances of the estimators φ E kk and φ A kk respectively. They are also crucial in discussing the separability of φ E kk and φ A kk . 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 ).
3) Similar definitions apply for measures related to learning the ξ-based interaction kernels, see Appendix G. We discuss some key properties of the measures in Appendix D.

Learning performance
We now discuss the performance measures for the estimated interaction kernels. We have already treated the trajectory estimation error in section 3.1. We use weighted L 2 -norms (with mild abuse of notation, we omit the weight from the notation) based on the dynamicsadapted measures introduced above, and with analogous definitions when discrete over L: Our learning theory focuses on minimizing the difference between φ E kk ⊕ φ A kk and φ E kk ⊕ φ A kk in the joint norm given by (5.4). 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 φ E kk − φ E kk 's and φ A kk − φ A kk '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 H.
wildcard, means the notation applies for all 3 variables shorthand denoting energy and alignment part of system   Table 3: Notation used throughout the paper

Notational summary
A summary of the learning theory notation introduced in sections 3, 4, and the notation above, is given below in table 3.

Identifiability of kernels from data
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: 1. 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.
2. Establishing the coercivity condition for a wide class of dynamical systems of the form (2.2), under mild assumptions on the distribution µ Y of the initial conditions. From our numerical experiments, we expect that coercivity holds even more generally.
We will make the following assumptions on the hypothesis spaces used in the learning approach for the remainder of the paper: Assumption 3 H ξ is a compact convex subset of K ξ S ξ (see 3.12) and is bounded above by It is easy to see that H EA can be naturally embedded as a compact subset of L 2 (ρ EA,L T ) and that H ξ can be naturally embedded as a compact subset of L 2 (ρ ξ,L T ) (recall these measures are defined in section 5.1). Assumptions 2, 3 ensure the existence of minimizers to the loss functionals E EA M , E ξ M defined in (4.1) and (4.2), which will be proven in Appendix B. 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 ξ 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 ρ 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 7.
Definition 4 (Coercivity condition) For the dynamical system (2.2) observed at time instants 0 = t 1 < t 2 < · · · < t L = T and with initial condition distributed µ Y on R (2d+1)N , it satisfies the coercivity condition on the hypothesis space 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 ].
In the following, we prove the coercivity condition on general compact sets of L 2 ([0, R], ρ EA,L T ) under suitable hypotheses. Our result is independent of N , which implies that the finite sample bounds of Theorem 9 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 → ∞).

Identifiability from coercivity
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 L 2 error to the true kernel). We want to choose a basis such that the regression matrix, A EA defined in Appendix H, 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 L 2 (ρ EA,L T ) 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 L 2 (ρ EA,L T ), versus the naive choice of basis in the underlying direct sum of L ∞ 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 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 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 ·, · on H EA × H EA , defined by satisfying the orthonomality and boundedness conditions 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, s E , s A ), and then across pairs (with the lexicographic ordering on pairs of integers). Set n = k,k n kk = dim(H EA ); then for any function ϕ EA ∈ H EA , we can write 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:

Proposition 5 Consider the matrices
and choose the hypothesis spaces as . Then the coercivity constants for H EA , H ξ are the smallest singular value of A EA ∞ , A ξ ∞ , respectively: with c H EA defined in (6.1) and c H ξ defined in (6.2). Additionally, for large M , the smallest singular value of A EA satisfies the inequality Similarly, for large M , the smallest singular value of A ξ M satisfies the inequality 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.
Proof We prove the result in the EA case, the proof of the results about the ξ part of the system is analogous. The orthonormality of the component functions given in (6.4), implies that ψ EA p , ψ EA p L 2 (ρ EA,L T ) = δ pp . Expand ϕ EA ∈ H EA in this basis as ϕ EA = n p=1 a p ψ EA p . Let the vector v = (a 1 , . . . , a n ) ∈ R n , and notice that This lower bound is achieved by the singular vector corresponding to the singular value σ min (A EA ∞ ), so that by definition (6.1) we have that σ min (A EA ∞ ) = c H EA .
For the second statement, we consider the learning matrix A EA (defined in section H), which we can also write as A EA M to emphasize the dependence on M as needed built, from the M observed trajectories. By construction, for each m, ∞ 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 [75], 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 [7] Problem III.6.13); for any two square matrices A, B, All norms in this proof are the spectral norm, unless otherwise specified. Thus if we get a concentration inequality of the form ≥ t} we will get the desired result relating the minimal singular values of 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ẋ} 2 S 2 EA + 2. This immediately implies the bound Next, we upper bound the matrix variance statistic (in our case Using a similar analysis to bound each entry of the matrices, we can arrive at the result that v(Z) ≤ 2n 2 c 2 1 . Now, we apply the matrix Bernstein inequality to see that Note that the M in the numerator comes because A EA M has a factor of 1 M on it. Lastly, choose t = c H EA 5 , which together with the results above yield the desired inequality.
From Proposition 5 we see that, for each hypothesis space H kk , it is important to choose a basis that is well-conditioned in L 2 (ρ EA,L T ), L 2 (ρ ξ,L T ), instead of in the corresponding L ∞ spaces. If not, the learning matrices, defined in Appendix H, A EA , A ξ 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.

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 µ Y of the initial conditions, the true interaction kernels, and the non-collective force. The coercivity condition has been studied for first-order systems in [52,51,48]. Below, we give a brief review.
For homogeneous systems, [52,51] 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 [51] and a few examples of the stochastic homogeneous system including linear systems and nonlinear systems with stationary distributions for general L in [48].
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. (6.1). We focus on a representative model of second-order homogeneous systems, 6) which includes the first-order systems considered in [13,52,51] as special cases and various second-order system examples in [52,83] as specific applications. We shall prove the coercivity condition holds true for the case L = 1: for a constant λ > 0, µẊ 0 , µ Ξ 0 are exchangeable with finite second moment, and they are independent of µ X 0 . Then Learning governing laws in second-order agent systems and c E H EA and c A H EA are non-negative constants independent of N , and are strictly positive for compact H EA of L 2 (ρ EA,1 T ).
and c H ξ is a non-negative constant independent of N , which is strictly positive for compact H ξ of L 2 (ρ ξ,1 T ). 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 [52,51,83]. The empirical numerical experiments on second-order systems, already conducted in [83], 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 C.

Consistency and optimal convergence rate of estimators
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.

Concentration
Our first main result is a concentration estimate that relates the coercivity condition to an appropriate bias-variance tradeoff in our setting. Let N (H, δ) be the δ-covering number, with respect to the ∞-norm, of the set H.
bounded above by S 0 ≥ max{S E , S A , S ξ } respectively. Additionally, assume that the coercivity condition (6.1) holds on H EA M and condition (6.2) on H ξ M .
Then for all > 0, with probability (with respect to µ Y ) at least 1 − δ, we have the estimates holds true with probability EA max{R, Rẋ} 2 K 6 .
with the final inequality holding with probability P. Now we bound the L 2 error of the empirical estimator to the true interaction kernel, so that with probability P The first inequality follow from the coercivity condition (6.1) and the definition of φ EA ∞ . The second follows by the definition of the norms. Now for a chosen 0 < δ < 1, let and solve for M . The proof for the ξ part of the system result is similar.

Consistency
In the regime where M → ∞, we will choose an increasing sequence of hypothesis spaces, each satisfying the conditions of Theorem 7. By our assumptions on the kernels, we can also choose the sequence of H EA M 's such that the approximation error goes to 0 as M → ∞. This enables us to control the infimum on the right hand side of (7.1). From here we can apply Theorem 7 on each M to prove the consistency of our estimators with respect to the L 2 (ρ EA,L T ) norm and derive the following consistency theorem.
is a family of compact and convex subsets such that the approximation error goes to zero, An analogous consistency result holds for the estimator in the ξ variable.
These two results together provide a consistency result on the full estimation of the triple ( φ ξ , φ E , φ A ) and thus consistency of our estimation procedure on the full system (2.2).
Let > 0, the inequality (7.4) gives us that We now bound the two terms in the above expression separately. For the first term, the proof of Theorem 7 shows that For the second term, the bound (B.3) yields that Since is fixed, the above result, together with our assumption on the sequence of hypothesis spaces, implies that The finiteness of the two sums above implies, by the first Borel-Cantelli Lemma, that As was arbitrary, we have the desired strong consistency of the estimator. An exactly analogous argument gives the result on the ξ part of the system.

Rate of convergence
Given data collected from M trajectories, we would like to choose the best hypothesis space to maximize the accuracy of the estimators. Theorem 7 highlights the classical bias-variance tradeoff in our setting. On the one hand, we would like the hypothesis space H EA M to be large so that the bias is small. Simultaneously, we would like H EA M to be small enough so that the covering number N (H EA M , ) 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.

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 ξ 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 10.
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 ρ EA,L T (which can be calculated as d(ρ EA,L T ) = kk p E (k,k ) + (k,k ) p A (k,k ) ), 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 V, denoted as |V|, for the EA portion and by |V ξ | for the ξ 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.
and assume that the coercivity condition (6.1) holds true on it. Then, there exists a constant C depending only on K, S EA , R, Rẋ such that , then there exists a constant C depending on K, R, Rẋ, S EA , c 0 , c 1 , s such that we achieve the convergence rate, (7.6) (c) under the corresponding assumptions as in (a), there exists a constant C depending only on K, S ξ , R such that (7.7) 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.
From this, the concentration estimate (7.1) together with the covering number bound imply that, which we will minimize to achieve the desired probability bound.
It is then an easy computation to see that the derivative of g( ) is ≤ 0 for all ≥ M . Thus, we can put this result into the bound (7.8) to arrive at the probability bound, Integrating this bound over ∈ (0, +∞) and using the elementary inequality e −x ≤ x + 1 for all x ≥ 0, we get that Now, bringing the coercivity part from (7.1) back in, we achieve the rate, where C 4 is an absolute constant that only depends on K, S EA , R, Rẋ.
For part (b), we note the following basic result on the covering number of B n by -balls (see [31,Proposition 5]), Using (7.1), and the approximation assumption, we bound the probability as 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 As before, let = tn −2s * = t M and consider h(t M ). It is easy to see that lim t→0 + h(t M ) = ∞ and lim t→∞ h(t M ) = −∞. 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 M ) = 0. We further need that h ( ) ≤ 0 for all ≥ c 5 M . By taking the derivative of h, setting it ≤ 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( ) ≤ 0 and h ( ) ≤ 0 for all ≥ c 6 M . These results imply the probability bound, Integrating this bound over ∈ (0, +∞) and using the elementary inequality e −x ≤ x + 1 for all x ≥ 0, we get that where C 1 is a constant depending on c 0 , c 1 , s, K, S EA , R, Rẋ. Now with H EA M = B n * and using (7.1), we have shown the convergence rate, where c 7 is an absolute constant that only depends on s, K, c 0 , c 1 , S EA , R, Rẋ.
In both theorems, the convergence rates 2s 2s+|V| and 2s 2s+|V ξ | 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 [39]), 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 9, n can be the number of partitions along each axis of the variables in 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 V as the dimension of the linear space. Furthermore, by standard approximation theory results, see [70] (Chapters 12,13), [33], [35], for s the regularity of the interaction kernels we achieve the desired approximation condition with piecewise polynomials of degree s . 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 2 we have as the final two columns the values |V|, |V ξ |. These correspond directly to the rate of convergence of each of the system under our learning approach. Some specific highlights: • 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.
• For the Synchronized Oscillator we achieve the 2-dimensional optimal learning rate on each of the EA and ξ 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 ξ-interaction kernels and accounting for shared and non-shared variables: this enables us to substantially improve the performance guarantees in actual applications.
• 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.
• The rates of convergence of our estimators for all previously-studied first-order systems (see [52,51,83]) can be derived from Theorem 9.
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 [83,52,51] demonstrate that more data in L may indeed be helpful to improve the performance. One technique used in [83] 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 [13], but further study of the mean-field regime is of interest to the authors and work is ongoing.

Performance of trajectory prediction
Once estimators ϕ EA , ϕ ξ 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 L 2 (ρ EA T ), L 2 (ρ ξ T ), respectively, error of the estimator. This further validates the effectiveness of our error functional and L 2 (ρ 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 L 2 (ρ EA T ), L 2 (ρ ξ T ) norms in obtaining our estimators.
and φ E , φ A , φ ξ respectively, both with the same initial condition. Then The constants are A 1 = 2T (8KP + L + 8QK + L ξ ) and B 1 = 2T 2 (8KP + 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 ξ coming from the admissible spaces. It is bounded on average, with respect to the initial distribution µ Y , by with the measures ρ EA T , ρ ξ T defined in (5.1, G.1). Expression (8.1) shows that by minimizing the right hand side, we can control the expected Y-supremum error of the estimated trajectories.
We postpone the somewhat lengthy proof to Appendix A.

Applications
Our learning theory, as well as measures, norms, functionals etc. can be applied to study all the examples considered in the works [52,51,83]. These examples, particularly those of [83], can thus be considered as applications of the theoretical results as well as of the algorithm in section H.
We choose to study two new dynamics, which are not considered in [52,83] 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 [72] and the anticipation dynamics (AD) model in [71]. Table 4 shows the value of learning parameters for these dynamics.  Table 4: Values of parameters for the learning.
The setup of the learning experiment is as follows. We use M ρ different initial conditions to evolve the dynamics 1 from 0 to T to generate a good approximation to ρ L,EA T , ρ L,E T and ρ L,A T . 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 φ E and φ A from the empirical distributions, ρ L,M,EA T 's, etc. We report the relative learning errors calculated via (5.4) for φ E ⊕ φ A , (5.4) for φ E , and (5.4) for φ 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 φ E ⊕ φ A and φ E ⊕ φ A from 0 to T f > T , and report the trajectory errors calculated using (3.5) on y (the whole system), and for x (the position) and 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.
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 φ E with a 6-digit accuracy, we are not able to do the same for φ A . The error in φ E ⊕ φ A reflects this discrepancy by considering the two functions together. 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, Compared to the original model in [71], we take τ i,i = 0. In order to fit the model into our learning regime, we take Here we have no ξ i , K = 1, m i = 1, and We also use τ = 0.1. It is shown in [71] that if U is bounded when r → ∞ with U (0) = U (0) = 0, then unconditionally flocking would occur. We take U (r) = r p p for 1 < p ≤ 2, then the system would show unconditional flocking. We choose p = 1.5 for our learning trials 3 . We use a tensor grid of 1 st degree piece-wise standard polynomials with n E = 28 2 for learning φ E (r, s), then a set of 1 st degree piecewise standard polynomials with n A = 138 for learning φ A (r). For the energy-based interactions we have the following results.
(b) U (r) = r 1.  Figure 4: 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.
As is shown in Fig. 4b, the concentration of pairwise distance data is away from 0, making the estimation of the behavior of φ E (r, s) at r close to 0 extremely difficult, meanwhile, since φ E is also weighted by the pairwise difference, x i − x i , and at r i,i close to 0, the information is also lost. Next, we present the alignment-based interaction kernels.
(b) U (r) = r 1. We have less trouble estimating the behavior of φ A at r = 0. We have trouble estimating φ A at the other end of the spectrum of r, since the agents have aligned their velocities, hence the weight v i − v i is close to a zero vector. The overall learning performance for estimating φ A is better compared to estimating φ E . The φ E ⊕ φ A error is: 6 · 10 −1 ± 3.0 · 10 −1 . The comparison of trajectories between the true kernels (LHS) and the estimators (RHS) is shown below. 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 6.
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.

Conclusion and further directions
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.

Acknowledgments and contributions
Similarly, let F [ϕ ξ ](x, ξ, s ξ ) := ϕ ξ (||x||, s ξ )ξ. Now we have by assumption that all trajectories start from the same initial conditions on both the position and velocity, which implies that Y (0) = Y (0) and˙ Y (0) =Ẏ (0). For every t ∈ [0, T ] we have that, by the fundamental theorem of calculus and the triangle inequality, Here we have introduced the term, which can be expressed explicitly as, Note that in I, j is the index of the type among the {1, . . . K} and i indexes within each type C j . This holds similarly in later expressions I 1 , I 2 . For the third term of (A.1), we exploit the Lipschitz property of the non-collective force: First we introduce the convenient notations of ii . Now we break up I using the triangle inequality and get that I ≤ I 1 + I 2 where So using the Lipschitz property of F [ φ EA jj ] we get that, since then, By the assumptions on the feature maps, we have that Combining these bounds we see that, , s A (j,j ) ] + 1) 2 , and then let P =SJ and we get by Young's inequality that, and performing a similar analysis we get that So gathering terms, we can reexpress (A.1) as Performing an analogous analysis on S , with some additional effort, one can get the following result on the phase variable Similarly, we have that, Gathering the bounds (A.7, A.8, A.9), we have that where we denote the last three lines by a(t) and notice that this is a nondecreasing function in t. We also denote A 1 = 2T (8KP + F + 8QK + F ξ ) and B 1 = 2T 2 (8KP + F ). Now use theorem 30, which is in [22] and is originally in Bainov and Simeonov. With this notation, we can rewrite the above bound as And so in the notation of Theorem 30 we have So that, So that we can immediately conclude the first assertion of the theorem, Lastly, we can use the results of section B.1 to get the key result on the expected supremum error. We take expectation on each of the three terms of a(T ) and normalize them so they are in the form of the results of B.1.
We similarly get that, and can get an analogous bound for the remaining term of a(T ). These bounds together lead to which implies the desired result.

B.1 Continuity of the error functionals
For any t ∈ [0, T ], consider the two random variables, These will be used in various places throughout the technical proofs and easily relate to the natural error functionals defined in (4.4), (4.5) as, We begin by establishing basic continuity results for our error functionals over the hypothesis space. The specific structure of the governing equations plays a critical role in the analysis.

Alignment and energy based kernels
Proposition 11 For ϕ EA , φ EA ∈ H EA the true and empirical error functionals are bounded as follows, Recall the definitions of R, Rẋ in equations (3.8), and (3.10).
Proof Using Jensen's inequality, and ρ t T = K,K k,k =1,1 ρ t,kk T . Therefore, we have that Taking the expectation with respect to µ Y on each side of (B.7) we get the first inequality. The second inequality follows by noticing that,

Environment interaction kernels
Here we show an analogous result to the alignment and energy result above. The techniques are similar and the result serves an identical purpose in the theory. Recall the definition of R ξ in (3.11).
Lemma 13 For all ϕ EA ∈ H EA , define the defect function L EA M (ϕ EA ) as Then, given two functions ϕ EA 1 , ϕ EA 2 ∈ H EA , the defect function is bounded by A similar lemma can be immediately deduced on the ξ variable .
Lemma 14 For all ϕ ξ ∈ H ξ , define the defect function L ξ M (ϕ ξ ) as Then, given two functions ϕ ξ 1 , ϕ ξ 2 ∈ H ξ , the defect function is bounded by almost surely with respect to µ Y .

B.2 Uniqueness of minimizers over a compact convex space
Recall the energy and alignment bilinear functional ·, · EA , previously defined in equation ∈ H EA . The S-inner product is the inner product induced by the · S norm by the polarization identity, which holds as we are working in an L 2 space, so the parallelogram law holds. Then our coercivity condition (6.1). can be written in terms of this bilinear functional as: for all ϕ EA ∈ H EA then for all ϕ EA ∈ H EA , the difference of the error functional at this element of H EA and the minimizer is lower bounded as, Thus, the minimizer of E EA ∞ over H EA is unique in L 2 (ρ EA,L T ).
Proof For ϕ EA ∈ H EA , and to ease the notation let φ EA := φ EA L,∞,H EA , we have Using that X, X − Y, Y = X − Y, X + Y , which holds by bilinearity and the definition of the form, we get that By the coercivity condition, the first term in (B.15) is at least as large as We are left to show the second term in (B.15) is nonnegative. Since H EA is convex, for all t ∈ [0, 1], t(ϕ EA ) + (1 − t)( φ EA ) ∈ H EA . By definition of φ EA as an argmin, which means, using a decomposition analogous to the one above in (B.15), that By the results of section B.1, we have (Lipschitz) continuity of the bilinear functional ·, · over H EA × H EA . Next, take the lim t→0 + of (B.16) and by the dominated convergence theorem (which holds due to the boundedness and continuity assumptions on the kernels) we pass the limit through the expectations in ·, · . This gives that (B.15) is greater than 0, giving the desired result on the uniqueness of the minimizer.

Proposition 16
Let the minimizer of the error functional be denoted φ L,∞,H ξ := arg min then for all ϕ ξ ∈ H ξ , the difference of the error functional at this element of H ξ and the minimizer is lower bounded as, Thus, the minimizer of E ξ ∞ over H ξ is unique in L 2 (ρ ξ,L T ).

B.3 Uniform estimates on defect functions
We start this section by introducing normalized errors of the estimators. Denote the minimizer of E EA ∞ (·) over H EA by For any ϕ EA ∈ H EA , define the normalized errors as These quantities capture the difference between the expected/empirical errors of the estimator and the function in the hypothesis space minimizing the expected error functional. We begin by proving a lemma that assumes the distance between the expected and empirical normalized errors are small for a given estimator. We then show that we have similar control on these distances for all points in a neighborhood of this particular one. This control enables us to apply a covering argument in the main proposition of this section due to the compactness of the hypothesis space.

Remark 17
Exactly analogous definitions hold for the ξ variable and we will simply state the results in that case.
Lemma 18 For all > 0 and 0 < α < 1, if the function ϕ EA Proof To ease the notation, write φ EA := φ EA L,∞,H EA , and using definition (B.11), we have that By definition we have that D ∞ (ϕ EA 2 ) ≥ 0 implying that, For the second term, by equation (B.4) and the assumption that α < 1, we obtain that We conclude that

B.4 Concentration
Proposition 20 For all > 0, 0 < α < 1, ϕ EA ∈ H EA , the following concentration bound holds Proof Consider the random variable Θ (with randomness coming from the random initial condition distributed µ Y ), and to ease the notation let φ EA := φ EA L,∞,H EA , The coercivity condition given in Definition (4), Proposition 15 and (B.3) allow us to bound the variance, denoted σ 2 , of Θ as follows.
By applying equation (B.8) from the proof of Proposition 11, we have that Θ ≤ 8S 2 EA max{R, Rẋ} 2 K 4 almost surely. We then apply the one-sided Bernstein inequality to Θ and recalling the definitions (B.1) together with the definitions of the normalized errors in (B.20, B.21), we get that: . Now we provide a lower bound for the exponent to simplify the dependencies. We show that, , or equivalently, By the estimate (B.22), since 0 < α ≤ 1, and 0 < c L,N,H EA < K 2 it is sufficient to show that This follows from Young's inequality as 2D ∞ (ϕ EA ) + 2 ≤ (D ∞ (ϕ EA ) + ) 2 , and together these results give the desired bound of the proposition.
We can easily derive the desired supremum bound by a covering argument. The estimation of the covering numbers involved will play a critical role in the main theorems and will be done in a dimension dependent way in order to get optimal minimax rates. Proposition 21 In the notation of Proposition 20, Now, by Proposition 20, for each i, By definition, H EA ⊆ i D i , so that Finally, we state the results for the ξ variable, the proofs are analogous. The advantage of splitting the theorems will become apparent. Specifically, it allows us to control the covering numbers on the EA and ξ hypothesis spacees separately, enabling us to get faster rates than if we viewed the task as estimating all functions simultaneously. This is possible due to the fundamentally decoupled nature of the dynamical system.

Proposition 23
In the notation of Proposition 22, denotes the covering number of H ξ with radius

Appendix C. Further verification of coercivity condition
In this appendix, we study the coercivity condition for the second-order system of the form: We prove the coercivity condition for the system (C.1) in the case of L = 1.
When the system does not have a ξ variable, the system (C.1) is related to the anticipation dynamics studied in [71]. When φ A ≡ 0 or φ E ≡ 0, the system is called energy-based or alignment-based respectively, and has found application in various disciplines including opinion dynamics, particle dynamics, fish-milling dynamics, Cucker-Smale flocking dynamics and phototaxix dynamics. We refer the reader to [52,51,83] where extensive numerical experiments were conducted to demonstrate the effectiveness of the proposed learning approach on the aforementioned dynamics.
For conciseness, we only present the proof of coercivity for learning of φ E and φ A . A similar argument can be conducted to prove the coercivity for learning φ ξ . Our aguments also work for special cases when φ A ≡ 0 or φ E ≡ 0. Therefore we obtain strict generalization of the coercivity results in [52,51] which are only for first-order energy-based systems. We may go further to analyze the coercivity condition for heterogeneous systems. Compared to the homogeneous system, the coercivity condition would impose constraints on the angle between components of subspaces defined on the directed sum of measures, suggesting that the learning task is more difficult in the case K ≥ 2.
Theorem 24 Consider the system (C.1) at time for a constant λ > 0, µẊ 0 , µ Ξ 0 are exchangeable with finite second moment, and they are independent of µ X 0 . Then (1) where we have and c E H EA and c A H EA are non-negative constants and are positive for compact H EA of L 2 (ρ EA,1 T ) and independent of N .
is a non-negative constant and is positive for compact H ξ of L 2 (ρ ξ,1 T ) and independent of N .

Proof
The proof of part (1) follows from the definition of measures, norms and the properties of the initial distributions. For part (2), we have By the property of µ Y 0 , we have for all (i, j, k), where X i s are exchangeable Gaussian random vectors with cov(X 1 ) − cov(X 1 , X 2 ) = λI d and Y i s are exchangeable random vectors who have the same distribution with the initial velocities of agentsẋ i s. From the Lemma 10 in [51] and Lemma 25 below, where the constants c E H EA , c A H EA ≥ 0 are always positive and independent of N for compact H EA , and we used the fact Therefore, For part (3), the proof follows a similar path as for part (2).
From Theorem 24, we see a particular case when the coercivity constant c 1,N,H EA is positive uniformly in N if c µẊ 0 > 0. In fact, many distributions on R dN with non-i.i.d R d components make the constant c µẊ 0 positive. For example, the components ofẊ are exchangeable Gaussian but not i.i.d, and d ≥ 2. In this particular case, coercivity is a property also of the system in the limit as N → ∞, satisfying the mean-field equations. As a result, the estimation error of our estimators is independent of N .
The proof of Theorem 24 used the following lemma, whose proof is the same with the proof of Lemma 10 in [51]. To be self-contained, we listed the statement here.
Proposition 29 Suppose the distribution µ Y of the initial condition is compactly supported. Then for each (k, k ), the support of the measures ρ EA,kk T , ρ ξ,kk T (and therefore ρ EA,L,kk T , ρ ξ,L,kk T ) is also compact.
Proof The compact support of the variables r,ṙ, ξ and the feature maps follows by the global well-posedness of the system in finite time, together with the Lipschitz assumptions on the non-collective forces. This compact support over a fixed, finite time is what is claimed in Proposition 29.
The main point is that by making reasonable assumptions on the non-collective forces, feature maps, interaction kernels, and the interval of time, together with the assumption that our agents' initial conditions cannot be arbitrarily far apart, we can derive that the pairwise distance, velocity and ξ will be controlled. Thus, the measures in section 5.1, if given enough trajectories, will be well approximated by the discretized version using the numerical approach described in section 9. Meaning that if we have a reasonable number of trajectories, we can look at the set of pairwise distances, velocities, etc. that these agents explore and bin them to set the support of the interaction kernels. Explicit values for the constants claimed in the proposition depend on the properties of the non-collective forces, the support and sup-norm of the interaction kernels, the interval T , and the number of agents.

Appendix E. Background results
In this section, for the convenience of the reader, we gather a few of the technical tools used in the analysis of the system. These are fundamental results necessary for developing the trajectory prediction, the measure support, and the existence and uniqueness results. We also include some of the necessary results on covering numbers of function spaces used for the learning theory.
The first theorem we present is an iterated Grönwall type result that allows us to analyze the trajectory error of the full system Y (t).

Appendix F. Additional comments on first-order models and theory
Our second-order model formulation covers the first-order equations of [51,52] as a special case. When Fẋ( ≡ 0 for all k, k = 1, . . . , K, and m i 1, (2.2) becomes, It extends the first-order models considered in [52,51,83] by adding non-collective forces, F x , F ξ , multi-dimensional interaction kernels, φ E k,k , φ ξ k,k , and auxiliary variables, ξ i . The first-order theory considered in [52,51] was focused on the learnability of functions of the form, φ E (r)r, which is a special case of our second-order theory, where we study the functions of the form φ E (r)r + φ A (r)ṙ, with φ A (r) ≡ 0.
is used in the actual algorithm. We consider direct sums of the measures for the phase variable for ease of notation. Lastly, for the ξ-based interaction kernels, i.e., φ ξ kk versus φ ξ kk , we consider the following norm,

Appendix H. Numerical algorithm
In this section, we will detail the construction of the linear systems to learn α EA and α ξ . We start from the procedure of solving for α EA . First, we build the basis functions for the finite dimensional subspace H EA ⊂ H EA .

Remark 31
The support of the unknown interaction kernels is not assumed to be known. We build our finite dimensional subspaces, H EA for example, based on the empirical observation data. For the support-detection capability of our estimators, see the examples of opinion dynamics in [52,83].
For the energy-based learning matrix, Ψ E,(m) , we use a lexicographical order on (k, k ) for k, k = 1, . . . , K. We define n E k,k ,prev = (k 1 ,k 2 )<(k,k ) n E k 1 ,k 2 ; if (k, k ) = (1, 1), we take n E 1,1,prev = 0. Then for η E kk = 1, . . . , n E k,k , Ψ E,(m) is given as follows, , s E i,i (t l ))(x i (t l )−x i (t l )), i ∈ C k , and for l = 1, . . . , L. Similar process of construction is done for Ψ A,(m) . Then we define, Then, α EA = ( α E ) T ( α A ) T T , is obtained by solving Then, we assemble Similar assembly from α A is done for φ A kk . In the case of using finite difference approximation to approximate the second derivatives of x i , we end up with where ζ = O( T L ) when a first-order finite difference scheme is used. Next for α ξ , we build the basis functions for the finite dimensional subspace H ξ kk ⊂ H ξ kk . We utilize the tensor grid of basis functions, i.e., tensor product  Thus, α ξ is obtained by solving A ξ α ξ = b ξ .
Then, we assemble