Simultaneous inference of environmental and interaction forces in collective dynamics

Nipuni de Silva; Ming Zhong; James M. Greene

cs.LGAug 25, 2026

Simultaneous inference of environmental and interaction forces in collective dynamics

Authors: Nipuni de Silva, Ming Zhong, James M. Greene

Organizations: Department of Mathematics, Clarkson University, Potsdam, NY, United States · Department of Mathematics, University of Houston, Houston, TX, United States

Abstract

Collective dynamics arise in a wide range of physical, biological, and engineering applications. Examples include cell migration, swarm robotics, social dynamics, and animal behavior. A defining characteristic of these systems is the emergence of large-scale coordination from local interactions among agents; a fundamental question is thus to understand the local interactions that give rise to the observed emergent dynamics. We are interested in methods for learning interactions generally, which can describe a wide class of physical systems exhibiting collective dynamics defined by an interaction kernel, without a priori assumptions on the analytical form of this kernel (i.e. it is nonparametric). The advantage of this kernel-based approach is that it incorporates the underlying physics of the model (i.e. collective dynamics), which more general equation-learning approaches may ignore, potentially limiting their effectiveness for model accuracy and predictions. In this work, we extend existing variational learning approaches to collective systems with both interaction kernels and environmental/intra-agent forces. The proposed framework simultaneously infers the interaction kernel non-parametrically while learning the environmental force using either semi-parametric or fully nonparametric representations. The methodology is validated on several benchmark models exhibiting synchronization, alignment, attraction-repulsion, and external environmental forces. We also introduce a model-selection procedure based on our nonparametric learning framework to identify models that optimally explain a given set of trajectory observations. By exploiting the feature-identification capability of the learned models, the proposed procedure can distinguish among different collective dynamics frameworks and recover mechanistic interaction mechanisms directly from trajectory data.

Figures & tables

Candidate ( $c)$	Candidate Model	Complexity
$F_{1}$	$\dot{x}=F_{\phi^{E}}$	1
$F_{2}$	$\dot{x}=F_{f}+F_{\phi^{E}}$	2
$S_{1}$	$\ddot{x}=F_{\phi^{E}}$	3
$S_{2}$	$\ddot{x}=F_{\phi^{A}}$	3
$S_{3}$	$\ddot{x}=F_{f}+F_{\phi^{E}}$	4
$S_{4}$	$\ddot{x}=F_{f}+F_{\phi^{A}}$	4

Table 1: Hierarchy of candidate models used for model selection. In the Candidate column,

F

denote a first-order framework, while

S

denote the second-order framework. The notation utilized for vector fields in

\mathbb{R}^{dN}

is as described in Sections 2.1.1 and Appendix A.1 (e.g. equation 7 ). We also suppress

\dot{x}_{i}=v_{i}

for second-order systems for notational simplicity.

Data set	Function
Training ( $M_{\mathrm{train}}$ samples)	Obtain estimators for mechanistic forces for each model
Validation ( $M_{\mathrm{val}}$ samples)	Compare candidate models and tune hyperparameters to select candidate model
Testing ( $M_{\mathrm{test}}$ samples)	Measure final evaluation of models and report performance

Table 2: Description of training, validation, and testing data sets for model selection algorithm. Note that the full set of initial condition replicates

M

is partitioned between the three classifications.

Metric	Data	Time horizon	Function
$E_{\mathrm{recon}}^{\mathrm{rel}}$	$M_{\mathrm{train}}$	$[0,T]$	Obtain admissible frameworks (trajectory data)
$E_{\mathrm{res},\epsilon}^{\mathrm{rel}}$	$M_{\mathrm{train}}$	$[0,T]$	Obtain admissible frameworks (residual data)
$\bar{E}_{\mathrm{traj}}$	$M_{\mathrm{val}}$	$[0,T_{f}]$	Primary validation metric (trajectories)
$E_{\mathrm{res}}$	$M_{\mathrm{train}}$	$[0,T]$	Residual consistency

Table 3: Summary of the evaluation metrics, the datasets on which they are computed, their associated time horizons, and their intended purposes.

Panel of Figure 1 — Figure 1: (Kuramoto model) Comparison of the mechanistic feature recovery ability of the fully non-parametric (NP) and semi-parametric (SP) estimation algorithms for the Kuramoto model ( 78 ). In Subfigure 1(a) we visualize the true interaction kernel $\phi$ (blue) together with the mean recovered interaction kernels $\bar{\hat{\phi}}$ obtained via the NP (yellow) and SP (orange) approaches. Subfigure 1(b) provides the corresponding true environmental force $f$ and its recovered mean estimates $\bar{\hat{f}}$ . Mean estimates are obtained over $T_{r}=10$ independent learning trials, and the corresponding standard deviations are also plotted. The background density in Subfigure 1(a) (right vertical axis) corresponds to the empirical distribution $\hat{\rho}_{\Delta\theta}$ of pairwise phase differences, while in Subfigure 1(b) this corresponds to the empirical distribution $\hat{\rho}_{\theta}$ of oscillator phases observed in the training data. The corresponding relative feature recovery errors, computed via equation ( 52 ), are reported in Table 4 . Both formulations recover the interaction kernel with comparable accuracy, while the semi-parametric formulation achieves nearly exact recovery of the environmental force. The training data consists of $M=20$ trajectories with $N=10$ oscillators in dimension $d=1$ , observed over the interval $[0,T]=[0,0.5]$ with $L=51$ time points.

Metric	Non-parametric	Semi-Parametric
$E^{\mathrm{rel}}_{\mathrm{res}}$	$9.133\,297\,54\text{\times}{10}^{-03}$ $\pm$ $1.069\,999\,32\text{\times}{10}^{-03}$	$9.213\,417\,49\text{\times}{10}^{-03}$ $\pm$ $1.074\,338\,81\text{\times}{10}^{-03}$
$E^{\mathrm{rel}}_{\phi}$	$1.151\,949\,00\text{\times}{10}^{-02}$ $\pm$ $6.567\,552\,90\text{\times}{10}^{-04}$	$1.150\,350\,07\text{\times}{10}^{-02}$ $\pm$ $6.505\,249\,25\text{\times}{10}^{-04}$
$E^{\mathrm{rel}}_{f}$	$1.992\,020\,20\text{\times}{10}^{-02}$ $\pm$ $5.681\,235\,00\text{\times}{10}^{-03}$	$1.712\,688\,22\text{\times}{10}^{-14}$ $\pm$ $1.184\,873\,90\text{\times}{10}^{-14}$

Table 4: (Kuramoto model) Comparison of relative residual error ( 55 ) and feature recovery metrics ( 52 ) for the fully non-parametric and semi-parametric estimation algorithms for the Kuramoto model ( 78 ). The table reports the mean error

\pm

one standard deviation over

T_{r}=10

independent learning trials. For additional details on algorithm parameters, we refer the reader to the caption of Figure 1 .

Panel of Figure 2 — Figure 2: (Kuramoto model) Comparison of the true phase trajectories $\theta_{i}$ (blue) and the estimated phase trajectories $\hat{\theta}_{i}$ (yellow) for both the semi-parametric (left column) and non-parametric (right column) methods. For trajectory errors, see Table 5 . The top row displays trajectories initialized from the $M_{\mathrm{train}}$ training trajectories, whereas the bottom rows displays trajectories generated from $M_{\mathrm{test}}$ previously unseen testing initial conditions. Note that trajectories are plotted over the full time horizon $[0,T_{f}]$ , with the end of the training interval time $T$ indicated on the horizontal axis. Recall that models are trained using observations collected only over the interval $[0,T]=[0,0.5]$ , and are evaluated beyond the fitting horizon to assess both reconstruction and predictive performance. Both method variants closely reproduce the true dynamics on the training trajectories and generalize well to unseen initial conditions, with no visual differences between the semi-parametric and non-parametric models. For additional details on algorithm parameters, we refer the reader to the caption of Figure 1 .

Data set	Statistic	Fully non-parametric	Semi-parametric
$M_{\mathrm{train}}$	$\bar{E}_{\mathrm{recon}}$	$1.486\,115\,21\text{\times}{10}^{-03}$ $\pm$ $1.290\,008\,01\text{\times}{10}^{-04}$	$1.515\,174\,72\text{\times}{10}^{-03}$ $\pm$ $1.326\,100\,95\text{\times}{10}^{-04}$
	$\sigma_{E_{\mathrm{recon}}}$	$7.059\,050\,63\text{\times}{10}^{-04}$ $\pm$ $1.149\,350\,14\text{\times}{10}^{-04}$	$7.258\,467\,47\text{\times}{10}^{-04}$ $\pm$ $1.074\,815\,46\text{\times}{10}^{-04}$
	$\bar{E}_{\mathrm{pred}}$	$5.457\,046\,59\text{\times}{10}^{-03}$ $\pm$ $1.454\,816\,20\text{\times}{10}^{-03}$	$5.459\,184\,78\text{\times}{10}^{-03}$ $\pm$ $1.409\,017\,06\text{\times}{10}^{-03}$
	$\sigma_{E_{\mathrm{pred}}}$	$5.665\,476\,87\text{\times}{10}^{-03}$ $\pm$ $2.895\,039\,32\text{\times}{10}^{-03}$	$6.020\,934\,52\text{\times}{10}^{-03}$ $\pm$ $2.947\,784\,88\text{\times}{10}^{-03}$
$M_{\mathrm{test}}$	$\bar{E}_{\mathrm{recon}}$	$1.635\,224\,03\text{\times}{10}^{-03}$ $\pm$ $1.994\,618\,41\text{\times}{10}^{-04}$	$1.598\,611\,49\text{\times}{10}^{-03}$ $\pm$ $1.293\,745\,35\text{\times}{10}^{-04}$
	$\sigma_{E_{\mathrm{recon}}}$	$9.230\,112\,42\text{\times}{10}^{-04}$ $\pm$ $2.374\,169\,98\text{\times}{10}^{-04}$	$9.094\,150\,34\text{\times}{10}^{-04}$ $\pm$ $2.159\,352\,77\text{\times}{10}^{-04}$

Table 5: (Kuramoto model) Comparison of trajectory-level errors for the fully non-parametric and semi-parametric estimation algorithms for Kuramoto model ( 78 ). For each learning trial, the trajectory statistics

\bar{E}_{\mathrm{recon}}

\bar{E}_{\mathrm{pred}}

\sigma_{E_{\mathrm{recon}}}

, and

\sigma_{E_{\mathrm{pred}}}

are computed according to ( 63 ) - ( 64 ). The table reports the mean

\pm

one standard deviation of these quantities computed on

T_{r}=10

independent learning trials. The first section corresponds to trajectories generated from

M_{\mathrm{train}}

training initial conditions, while the second section corresponds to

M_{\mathrm{test}}

unseen predicted initial conditions sampled from

\mu_{0}

. Note that across all trajectory measures, error estimates are both small and similar for both methods. For additional details on algorithm parameters, we refer the reader to the caption of Figure 1 .

Panel of Figure 3 — Figure 3: (SPP model) Comparison of mechanistic feature recovery ability of the fully non-parametric (NP) and semi-parametric (SP) estimation algorithms for the SPP model ( 79 ). In Subfigure 3(a) we visualize the true interaction kernel $\phi$ (blue) with the mean recovered kernels $\bar{\hat{\phi}}$ obtained via the NP (yellow) and SP (orange) formulations. Subfigure 3(b) provides the true environmental force $f$ with the corresponding mean learned estimates $\bar{\hat{f}}$ . As $f:\mathbb{R}^{2}\to\mathbb{R}^{2}$ , the top and bottom rows correspond to the first and second components of $f$ , respectively. Within each row, the true force (left column), the fully non-parametric mean estimate (middle column), and the semi-parametric mean estimate (right column) are shown, together with one standard deviation, which are obtained over $T_{r}=10$ independent learning trials. The background density in Subfigure 3(a) (right vertical axis) corresponds to the empirical pairwise distance distribution $\hat{\rho}_{R}$ , while the density visualized in Subfigure 3(b) (left column) corresponds to the empirical velocity distribution $\hat{\rho}_{V}$ induced by the training data. The corresponding relative feature recovery errors are reported in Table 6 . The training data consists of $M=20$ trajectories with $N=10$ agents in dimension $d=2$ , observed over $[0,T]=[0.0,0.5]$ with $L=51$ time points.

Metric	Non-parametric	Semi-Parametric
$E^{\mathrm{rel}}_{\mathrm{res}}$	$3.250\,636\,15\text{\times}{10}^{-02}$ $\pm$ $3.804\,678\,17\text{\times}{10}^{-04}$	$1.464\,225\,90\text{\times}{10}^{-02}$ $\pm$ $7.416\,366\,49\text{\times}{10}^{-04}$
$E^{\mathrm{rel}}_{\phi}$	$1.199\,870\,60\text{\times}{10}^{-01}$ $\pm$ $1.993\,494\,02\text{\times}{10}^{-03}$	$1.197\,529\,15\text{\times}{10}^{-01}$ $\pm$ $1.841\,751\,39\text{\times}{10}^{-03}$
$E^{\mathrm{rel}}_{f}$	$2.974\,167\,42\text{\times}{10}^{-02}$ $\pm$ $3.449\,852\,01\text{\times}{10}^{-04}$	$6.480\,458\,97\text{\times}{10}^{-04}$ $\pm$ $1.344\,640\,77\text{\times}{10}^{-04}$

Table 6: (SPP model) Comparison of relative residual error ( 55 ) and feature recovery metrics ( 52 ) for the fully non-parametric and semi-parametric estimation algorithms for the SPP model ( 79 ). The table reports the mean error

\pm

one standard deviation of these quantities computed over

T_{r}=10

independent learning trials. For additional details on algorithm parameters, we refer the reader to the caption of Figure 3 .

Panel of Figure 4 — Figure 4: (SPP model) Comparison of the true trajectories $x$ (left column) and the estimated trajectories $\hat{x}$ for both the semi-parametric (middle column) and fully non-parametric (right column) methods in the phase plane. For trajectory errors, see Table 7 . The top row displays trajectories initialized from the $M_{\mathrm{train}}$ training initial-conditions, whereas the bottom row displays trajectories generated from $M_{\mathrm{test}}$ previously unseen testing initial conditions. Note that trajectories are plotted over the full time horizon $[0,T_{f}]$ , with the color gradient representing the temporal evolution of the dynamics, from yellow at $t=0$ to orange at $t=T_{f}$ . The black segment indicates the training interval $[0,T]$ , which is utilized for model fitting; recall that no fitting is performed in the testing data in Figure 4(b) . Here $[0,T]=[0.0,0.5]$ , with testing evaluated beyond the fitting horizon $[0.0,0.5]$ until $T_{f}=2.0$ to assess both reconstruction and prediction in Figure 4(b) . Both method variants closely reproduce the true dynamics on the training trajectories and generalize well to unseen initial conditions. For additional details on algorithm parameters, we refer the reader to the caption of Figure 3 .

Data set	Statistic	Fully non-parametric	Semi-parametric
$M_{\mathrm{train}}$	$\bar{E}_{\mathrm{recon}}$	$4.291\,656\,12\text{\times}{10}^{-03}$ $\pm$ $7.329\,969\,06\text{\times}{10}^{-05}$	$2.762\,636\,52\text{\times}{10}^{-03}$ $\pm$ $8.642\,627\,33\text{\times}{10}^{-05}$
	$\sigma_{E_{\mathrm{recon}}}$	$1.018\,780\,56\text{\times}{10}^{-03}$ $\pm$ $1.047\,847\,99\text{\times}{10}^{-04}$	$9.363\,603\,47\text{\times}{10}^{-04}$ $\pm$ $1.232\,918\,49\text{\times}{10}^{-04}$
	$\bar{E}_{\mathrm{pred}}$	$4.969\,235\,82\text{\times}{10}^{-02}$ $\pm$ $2.224\,500\,17\text{\times}{10}^{-03}$	$2.484\,601\,58\text{\times}{10}^{-02}$ $\pm$ $1.481\,024\,68\text{\times}{10}^{-03}$
	$\sigma_{E_{\mathrm{pred}}}$	$9.785\,216\,05\text{\times}{10}^{-03}$ $\pm$ $9.877\,813\,09\text{\times}{10}^{-04}$	$9.042\,348\,73\text{\times}{10}^{-03}$ $\pm$ $1.314\,571\,31\text{\times}{10}^{-03}$
$M_{\mathrm{test}}$	$\bar{E}_{\mathrm{recon}}$	$4.256\,561\,62\text{\times}{10}^{-03}$ $\pm$ $1.291\,396\,91\text{\times}{10}^{-04}$	$2.767\,011\,67\text{\times}{10}^{-03}$ $\pm$ $1.018\,333\,24\text{\times}{10}^{-04}$
	$\sigma_{E_{\mathrm{recon}}}$	$1.133\,254\,76\text{\times}{10}^{-03}$ $\pm$ $1.405\,147\,45\text{\times}{10}^{-04}$	$8.268\,659\,24\text{\times}{10}^{-04}$ $\pm$ $9.900\,607\,36\text{\times}{10}^{-05}$

Table 7: (SPP model) Comparison of trajectory-level errors for the fully non-parametric and semi-parametric estimation algorithms for the SPP model ( 79 ). For each learning trial, the trajectory statistics

\bar{E}_{\mathrm{recon}}

\bar{E}_{\mathrm{pred}}

\sigma_{E_{\mathrm{recon}}}

, and

\sigma_{E_{\mathrm{pred}}}

are computed according to ( 63 ) - ( 64 ). The table reports the mean

\pm

one standard deviation of these quantities computed over

T_{r}=10

independent learning trials. The first section corresponds to trajectories generated from

M_{\mathrm{train}}

training initial conditions, while the second section corresponds to

M_{\mathrm{test}}

unseen initial conditions sampled from

\mu_{0}

Panel of Figure 5 — Figure 5: (Misspecified parametric model, SPP) Results for the SPP system ( 79 ) when the semi-parametric algorithm is supplied with an incorrect functional form for the environmental force $f$ . Here $f$ is assumed of the form $f(v;p)=p_{1}v+p_{2}||v||_{1}v$ , while the data is generated from the model ( 79 ). Note that the algorithm is unable to recover the correct mechanistic features and dynamics. Subfigure 5(a) plots the resulting phase portraits for one initial condition from the $M_{\mathrm{train}}$ training data; observe that even on the training trajectories, the learned dynamics deviate noticeably from the true dynamics within the predictive interval $[T,T_{f}]$ . Subfigure 5(b) compares the true interaction kernel $\phi$ (blue solid line) with the learned kernel $\hat{\phi}$ obtained from the misspecified parametric model (red dashed line). Subfigure 5(c) compares the true (left column) and estimated (right column) environmental forces $f$ , with the upper and lower panels corresponding to the first and second components, respectively. For additional details on algorithm parameters, we refer the reader to the caption of Figure 3 .

Panel of Figure 6 — Figure 6: (Phototaxis model) Comparison of mechanistic feature recovery ability of the fully non-parametric (NP) and semi-parametric (SP) estimation algorithms for the phototaxis model ( 82 ). In Subfigure 6(a) we visualize the true interaction kernel $\phi$ (blue) together with the mean recovered kernels $\bar{\hat{\phi}}$ obtained via the NP (yellow) and P (orange) approaches. Subfigure 6(b) provides the corresponding true environmental force $f$ and its recovered mean estimates $\bar{\hat{f}}$ . As $f:\mathbb{R}^{2}\to\mathbb{R}^{2}$ , the top and bottom rows correspond to first and second components of $f$ , respectively. Within each row, the true force (left column), the fully non-parametric mean estimate (middle column), and semi-parametric mean estimate (right column) are shown, together with one standard deviation, which are obtained over $T_{r}=10$ independent learning trials. The background density in Subfigure 6(a) (right vertical axis) corresponds to the empirical pairwise distance distribution $\hat{\rho}_{R}$ , while the density visualized in Subfigure 6(b) (left column) corresponds to the empirical velocity distribution $\hat{\rho}_{V}$ induced by the training data. The corresponding relative feature recovery errors are reported in Table 8 . The training data consists of $M=20$ trajectories with $N=10$ agents in dimension $d=2$ , observed over $[0,T]=[0,0.5]$ with $L=51$ time points.

Metric	Non-parametric	Semi-Parametric
$E^{\mathrm{rel}}_{\mathrm{res}}$	$6.726\,628\,52\text{\times}{10}^{-05}$ $\pm$ $6.398\,876\,10\text{\times}{10}^{-06}$	$7.291\,303\,63\text{\times}{10}^{-05}$ $\pm$ $7.334\,376\,82\text{\times}{10}^{-06}$
$E^{\mathrm{rel}}_{\phi}$	$1.985\,512\,11\text{\times}{10}^{-04}$ $\pm$ $2.292\,855\,28\text{\times}{10}^{-05}$	$1.980\,978\,68\text{\times}{10}^{-04}$ $\pm$ $2.288\,845\,60\text{\times}{10}^{-05}$
$E^{\mathrm{rel}}_{f}$	$5.533\,440\,49\text{\times}{10}^{-05}$ $\pm$ $8.754\,940\,29\text{\times}{10}^{-06}$	$3.738\,354\,94\text{\times}{10}^{-06}$ $\pm$ $2.418\,218\,40\text{\times}{10}^{-06}$

Table 8: (Phototaxis model) Comparison of relative residual error ( 55 ) and feature recovery metrics ( 52 ) for the fully non-parametric and semi-parametric estimation algorithms for the phototaxis model ( 82 ). The table reports the mean error

\pm

one standard deviation of these quantities computed over

T_{r}=10

independent learning trials. For additional details on algorithm parameters, we refer the reader to the caption of Figure 6 .

Panel of Figure 7 — Figure 7: (Phototaxis model) Comparison of true trajectories $x$ (left column) and the estimated trajectories $\hat{x}$ for both the semi-parametric (middle column) and full non-parametric (right column) methods in the phase plane. For trajectory errors, see Table 9 . The top row displays trajectories initialized from the $M_{\mathrm{train}}$ training initial-conditions, whereas the bottom row displays trajectories generated from $M_{\mathrm{test}}$ previously unseen testing initial conditions. Note that trajectories are plotted over the full time horizon $[0,T_{f}]$ , with the color gradient representing the temporal evolution of the dynamics, from yellow at $t=0$ to orange at $t=T_{f}$ . The black segment indicates the training interval $[0,T]$ , which is utilized for model fitting; recall that no fitting is performed in the testing data in Figure 4(b) . Here $[0,T]=[0,0.5]$ , with testing evaluated beyond the fitting horizon $[0.0,0.5]$ until $T_{f}=2.0$ to assess both reconstruction and prediction in Figure 7(b) . Both method variants closely reproduce the true dynamics on the training trajectories and generalize well to unseen initial conditions. For additional details on algorithm parameters, we refer the reader to the caption of Figure 6 .

Data set	Statistic	Non-parametric	Semi-Parametric
$M_{\mathrm{train}}$	$\bar{E}_{\mathrm{recon}}$	$4.027\,552\,51\text{\times}{10}^{-06}$ $\pm$ $6.549\,230\,71\text{\times}{10}^{-07}$	$4.496\,428\,52\text{\times}{10}^{-06}$ $\pm$ $7.364\,670\,99\text{\times}{10}^{-07}$
	$\sigma_{E_{\mathrm{recon}}}$	$1.009\,259\,16\text{\times}{10}^{-06}$ $\pm$ $3.055\,190\,48\text{\times}{10}^{-07}$	$1.268\,869\,62\text{\times}{10}^{-06}$ $\pm$ $3.899\,608\,97\text{\times}{10}^{-07}$
	$\bar{E}_{\mathrm{pred}}$	$1.367\,294\,68\text{\times}{10}^{-05}$ $\pm$ $2.262\,020\,00\text{\times}{10}^{-06}$	$1.178\,536\,84\text{\times}{10}^{-05}$ $\pm$ $1.560\,367\,45\text{\times}{10}^{-06}$
	$\sigma_{E_{\mathrm{pred}}}$	$2.768\,822\,80\text{\times}{10}^{-06}$ $\pm$ $7.388\,936\,98\text{\times}{10}^{-07}$	$3.210\,520\,55\text{\times}{10}^{-06}$ $\pm$ $6.862\,167\,86\text{\times}{10}^{-07}$
$M_{\mathrm{test}}$	$\bar{E}_{\mathrm{recon}}$	$5.067\,807\,84\text{\times}{10}^{-06}$ $\pm$ $6.953\,979\,44\text{\times}{10}^{-07}$	$4.523\,697\,17\text{\times}{10}^{-06}$ $\pm$ $5.193\,499\,62\text{\times}{10}^{-07}$
	$\sigma_{E_{\mathrm{recon}}}$	$1.163\,612\,56\text{\times}{10}^{-06}$ $\pm$ $3.301\,187\,82\text{\times}{10}^{-07}$	$1.012\,524\,41\text{\times}{10}^{-06}$ $\pm$ $3.038\,581\,17\text{\times}{10}^{-07}$

Table 9: (Phototaxis model) Comparison of trajectory-level errors for the fully non-parametric and semi-parametric estimation algorithms for the phototaxis model ( 82 ). For each learning trial, the trajectory statistics

\bar{E}_{\mathrm{recon}}

\bar{E}_{\mathrm{pred}}

\sigma_{E_{\mathrm{recon}}}

, and

\sigma_{E_{\mathrm{pred}}}

are computed according to ( 63 ) - ( 64 ). The table reports the mean

\pm

one standard deviation of these quantities computed over

T_{r}=10

independent learning trials. The first section corresponds to trajectories generated from

M_{\mathrm{train}}

training initial conditions, while the second section corresponds to

M_{\mathrm{test}}

unseen initial conditions sampled from

\mu_{0}

Data set	Statistic	Non-parametric	Semi-parametric
Kuramoto	$\bar{E}_{\mathrm{recon}}$	$8.639\,226\,09\text{\times}{10}^{-04}$ $\pm$ $1.436\,067\,53\text{\times}{10}^{-04}$	$7.933\,406\,60\text{\times}{10}^{-04}$ $\pm$ $9.072\,559\,42\text{\times}{10}^{-05}$
	$\sigma_{E_{\mathrm{recon}}}$	$2.033\,986\,40\text{\times}{10}^{-04}$ $\pm$ $2.378\,269\,91\text{\times}{10}^{-05}$	$2.029\,332\,07\text{\times}{10}^{-04}$ $\pm$ $2.688\,566\,76\text{\times}{10}^{-05}$
	$\bar{E}_{\mathrm{pred}}$	$3.360\,057\,73\text{\times}{10}^{-03}$ $\pm$ $1.260\,603\,89\text{\times}{10}^{-03}$	$2.721\,694\,99\text{\times}{10}^{-03}$ $\pm$ $9.891\,757\,78\text{\times}{10}^{-04}$
	$\sigma_{E_{\mathrm{pred}}}$	$1.502\,954\,18\text{\times}{10}^{-03}$ $\pm$ $7.714\,561\,34\text{\times}{10}^{-04}$	$1.372\,298\,20\text{\times}{10}^{-03}$ $\pm$ $6.667\,325\,58\text{\times}{10}^{-04}$
SPP	$\bar{E}_{\mathrm{recon}}$	$7.886\,463\,53\text{\times}{10}^{-03}$ $\pm$ $2.475\,329\,39\text{\times}{10}^{-04}$	$0.006\,884\,921\,442\,267\,94$ $\pm$ $0.000\,253\,715\,904\,779\,536\,2$
	$\sigma_{E_{\mathrm{recon}}}$	$1.737\,702\,05\text{\times}{10}^{-03}$ $\pm$ $4.804\,516\,38\text{\times}{10}^{-04}$	$0.002\,113\,167\,912\,756\,090\,8$ $\pm$ $0.000\,543\,637\,710\,861\,935\,6$

Table 10: (Scaling number of agents) Comparison of trajectory level errors for the fully non-parametric and semi-parametric estimation algorithms for the Kuramoto ( 78 ), SPP ( 79 ), and phototaxis ( 82 ) models. For each learning trial, the trajectory statistics

\bar{E}_{\mathrm{recon}}

\bar{E}_{\mathrm{pred}}

\sigma_{E_{\mathrm{recon}}}

, and

\sigma_{E_{\mathrm{pred}}}

are computed according to ( 63 ) - ( 64 ). The table reports the mean

\pm

one standard deviation of these quantities computed

T_{r}=10

independent learning trials. Predicted trajectories are made with

M_{\mathrm{test}}^{4N}

untrained initial conditions generated for

4N

agents. In numerical simulations,

M_{\mathrm{test}}^{4N}=M_{\mathrm{test}}

. Further information regarding learning and model parameters may be found in Appendix E .

Panel of Figure 8 — Figure 8: (Scaling number of agents) The true and predicted trajectories for the Kuramoto (Figure 8(a) ), SPP (Figure 8(b) ), and phototaxis (Figure 8(c) ) systems obtained via the fully non-parametric method, when obtaining by the estimators by training on a smaller number of agents. Trajectory predictions are provided on $N_{\mathrm{new}}=4N$ agents, when $N=10$ agents were utilized in the inference. Trajectory errors are reported in Table 10 in the semi-parametric column. Further information regarding learning and model parameters may be found in Appendix E .

Figure 9: (Scaling number of agents: SPP milling model) The true and predicted trajectories for the SPP milling system obtained using the learned model, with milling parameters $C_{r}=1,l_{r}=0.5,C_{a}=0.5,l_{a}=2,\alpha=1.6$ , and $\beta=0.5$ . The top row shows the trajectory predictions over the training time interval $[0,T]=[0,4]$ for $N=10$ agents. The bottom row shows the trajectory predictions for a scaled system with $N_{\mathrm{new}}=100$ agents, evolved over a time horizon five times longer than the training interval, with $T_{f}=20$ . Further information regarding the learning parameters, system parameters, and initial conditions may be found in Appendix E .

Panel of Figure 10 — Figure 10: Mean mechanistic relative feature recovery errors (log scale) as functions of temporal resolution $L$ and independent replicate size $M$ computed over $T_{r}=10$ independent learning trials. The relatively large errors observed for $M=1$ replicate, even when increasing the number of time samples $L$ , suggests that a single trajectory does not contain sufficient information to accurately recover the underlying dynamics. The above figure, for both interaction kernel $\phi$ (Figure 10(a) ) and environmental force $f$ (Figure 10(b) ), suggests that the number of replicates $M$ is the dominant sampling parameter which controls the accuracy of the estimation algorithm. Further information regarding learning and model parameters may be found in Appendix E , as well as in the figure captions of Section 5.2 .

Panel of Figure 11 — Figure 11: (Noise robustness, phototaxis model) Mechanistic feature recovery under multiplicative observational noise with noise levels $\zeta$ . Subfigure 11(a) plots the estimated interaction kernels $\phi$ for different noise levels, while Subfigure 11(d) plots the corresponding estimated environmental forces $f$ . The estimated empirical measures $\hat{\rho}_{R}$ and $\hat{\rho}_{V}$ induced by the training data are also included in each subfigure. Subfigures 11(b) and 11(c) compare the true, noisy, and learned trajectories for $\zeta=0.01$ and $\zeta=0.1$ , respectively. True (blue), noisy (grey), and estimated (yellow to orange, varying by time) trajectory data are all plotted; recall that training occurs on the noisy (grey) data. Training occurs on $[0,T]$ , and testing (predictions) are provided on $(T,T_{f}]$ . We observe that as $\zeta$ increases, feature recovery gradually deteriorates, primarily in regions with limited trajectory data (small measures), while the learned trajectories remain close to the true trajectories from low to moderate noise levels. Further information regarding learning and model parameters may be found in Appendix E , as well as in the figure captions of Section 5.2 .

Panel of Figure 12 — Figure 12: (Velocity estimation with noisy data, Kuramoto model) Comparison of filtering methods for trajectory data with multiplicative observational noise in the Kuramoto model ( 78 ). Here the noise level is fixed at $\zeta=0.01$ (see ( 68 ) and ( 69 )). Subfigures 12(a) and 12(b) provide plots for the noisy observations of state $\theta$ and their derivatives $\dot{\theta}$ , respectively, for a representative trajectory, with true trajectories (black), noisy observations (blue), a sliding window (SW) filter (orange), and a Kalman filter with the Rauch–Tung–Striebel (KF-RTS) smoother (green) included. Subfigure 12(c) plots the estimated interaction kernels $\phi$ , while Subfigure 12(d) plot the recovered environmental forces $f$ obtained from the true, noisy, and filtered observations. The estimated empirical measures $\hat{\rho}_{\Delta\theta}$ and $\hat{\rho}_{\theta}$ induced by the training data are also included in each subfigure. Subfigure 12(e) compares the true, noisy, smoothed, and learned trajectories for a representative training trajectory (i.e. a single agent). Further information regarding learning and model parameters may be found in Appendix E , as well as in the caption of Figure 1 .

Order	$\phi^{E}$	$\phi^{A}$	$f$
2	-	✓	-

Table 11: Mechanisms of the Cucker-Smale system ( 83 ). The symbol “-" indicates that the mechanism is not present, while “✓" indicates that the mechanism is present.

Panel of Figure 13 — Figure 13: (Model selection, Cucker-Smale) Mechanisms present in the Cucker-Smale system ( 83 ) within the generalized framework ( 72 ). The model selection framework estimates three quantities: an energy-based interaction kernel $\hat{\phi}^{E}$ , an alignment interaction kernel $\hat{\phi}^{A}$ , and an environmental force $\hat{f}$ . From observations of trajectories, the algorithm correctly identify the system as alignment-dominated. Estimation is performed over $T_{r}=10$ independent learning trials, and the inferred alignment kernel (Figure 13(b) ) has weighted $L^{2}$ norm $\|\hat{\phi}^{A}\|_{L^{2}(\hat{\rho}_{R})}=$ $0.742\,153$ $\pm$ $0.003\,468$ , which is several orders of magnitude larger than the inferred energy kernel $\|\hat{\phi}^{E}\|_{L^{2}(\hat{\rho}_{R})}=$ $1.9985\text{\times}{10}^{-04}$ $\pm$ $1.0221\text{\times}{10}^{-04}$ (Figure 13(a) ) and environmental force $\|\hat{f}\|_{L^{2}(\hat{\rho}_{V})}=$ $0.000\,528$ $\pm$ $0.000\,070$ (Figure 13(c) ). Note the scale of the vertical axes. Further information regarding learning and model parameters may be found in Appendix E .

Candidate	Order	$E_{\mathrm{res}}^{\mathrm{rel}_{\epsilon}}$	$E_{\mathrm{recon}}^{\mathrm{rel}}$	$\bar{E}_{\mathrm{traj}}(M_{\mathrm{val}})$	Complexity	$E_{\mathrm{res}}$	Rank	Selected
$S_{2}$	2	0.000626	0.000004	0.000259	3	0.000509	1	Yes
$S_{4}$	2	0.000708	0.000004	0.000437	4	0.000575	2	No
$S_{5}$	2	0.000624	0.000004	0.000255	5	0.000507	3	No
$S_{6}$	2	0.000717	0.000004	0.000464	6	0.000583	4	No
$F_{2}$	1	0.593726	0.003266	0.732061	2	0.482466	5	No
$S_{3}$	2	0.593726	0.003266	0.732061	4	0.482466	6	No

Table 12: (Model selection, Cucker-Smale) Summary of the candidate models considered during model selection. For each candidate, we report the system order, framework compatibility diagnostics, validation trajectory error, model complexity, residual error, ranking, and whether the candidate was ultimately selected based on the observed trajectory data. In this case, the framework correctly identifies the Cucker-Smale flocking model ( 83 ).

Order	$\phi^{E}$	$\phi^{A}$	$f$
2	-	✓	✓

Table 13: Mechanisms of the phototaxis system ( 82 ). The symbol “-" indicates that the mechanism is not present, while “✓" indicates that the mechanism is present.

Panel of Figure 14 — Figure 14: (Model selection, phototaxis) Mechanisms present in the phototaxis system ( 82 ) within the generalized framework ( 72 ). The model selection framework estimates three quantities: an energy-based interaction kernel $\hat{\phi}^{E}$ , an alignment interaction kernel $\hat{\phi}^{A}$ , and an environmental force $\hat{f}$ . From observation of trajectories, the algorithm correctly identify the system as being governed by alignment and environmental force interactions. Estimation is performed over $T_{r}=10$ independent learning trials, and the inferred learned alignment kernel (Figure 14(b) ) has weighted $L^{2}$ norm $\|\hat{\phi}^{A}\|_{L^{2}(\hat{\rho}_{R})}=$ $0.931\,035$ $\pm$ $0.002\,016$ and the estimated environmental force has weighted $L^{2}$ norm $\|\hat{f}\|_{L^{2}(\hat{\rho}_{V})}=$ $0.600\,504$ $\pm$ $0.006\,211$ (Figure 14(c) ), both of which are several orders of magnitude larger than the inferred energy kernel $\|\hat{\phi}^{E}\|_{L^{2}(\hat{\rho}_{R})}=$ $7.7845\text{\times}{10}^{-04}$ $\pm$ $6.9402\text{\times}{10}^{-04}$ (Figure 14(a) ). Note the scale of the vertical axes. These results indicate that the model-selection framework correctly identifies the active mechanisms while suppressing the inactive energy interaction. Further information regarding learning and model parameters may be found in Appendix E .

Candidate	Order	$E_{\mathrm{res}}^{\mathrm{rel}_{\epsilon}}$	$E_{\mathrm{recon}}^{\mathrm{rel}}$	$\bar{E}_{\mathrm{traj}}(M_{\mathrm{val}})$	Complexity	$E_{\mathrm{res}}$	Rank	Selected
$S_{4}$	2	0.000067	0.000002	0.000016	4	0.000040	1	Yes
$S_{6}$	2	0.000068	0.000002	0.000017	6	0.000041	2	No
$F_{2}$	1	0.216497	0.009056	0.102124	2	0.130341	3	No
$S_{3}$	2	0.216497	0.009056	0.102124	4	0.130341	4	No
$S_{2}$	2	0.224109	0.015654	0.272981	3	0.134924	5	No
$S_{5}$	2	0.224109	0.015654	0.272981	5	0.134924	6	No

Table 14: (Model selection, phototaxis) Summary of the candidate models considered during model selection. For each candidate, we report the system order, framework compatibility diagnostics, validation trajectory error, model complexity, residual error, ranking, and whether the candidate was ultimately selected based on the observed trajectory data. In this case, the framework correctly identifies the phototaxis model ( 82 ).

Order	$\phi^{E}$	$\phi^{A}$	$f$
2	✓	✓	✓

Table 15: Mechanisms of the SPPCS system. The symbol “-" indicates that the mechanism is not present, while “✓" indicates that the mechanism is present.

Panel of Figure 15 — Figure 15: (Model selection, SPPCS) Mechanisms present in the SPPCS system within the generalized framework ( 72 ). The model selection framework estimates three quantities: an energy-based interaction kernel $\hat{\phi}^{E}$ , an alignment interaction kernel $\hat{\phi}^{A}$ , and an environmental force $\hat{f}$ . From observations of trajectories, the algorithm correctly identify the system as being governed by all three mechanisms. Estimation is performed over $T_{r}=10$ independent learning trials, and the inferred alignment kernel has weighted $L^{2}$ norm $\|\hat{\phi}^{A}\|_{L^{2}(\hat{\rho}_{R})}=$ $0.344\,557$ $\pm$ $0.031\,942$ (Figure 15(b) ), the inferred environmental force has weighted $L^{2}$ norm $\|\hat{f}\|_{L^{2}(\hat{\rho}_{V})}=$ $3.062\,965$ $\pm$ $0.047\,218$ (Figure 15(c) ), and the inferred energy kernel has weighted $L^{2}$ norm $\|\hat{\phi}^{E}\|_{L^{2}(\rho_{R})}=$ $6.5314\text{\times}{10}^{-01}$ $\pm$ $6.0227\text{\times}{10}^{-02}$ (Figure 15(a) ). Note the scale of the vertical axes. Further information regarding learning and model parameters may be found in Appendix E .

Candidate	Order	$E_{\mathrm{res}}^{\mathrm{rel_{\epsilon}}}$	$E_{\mathrm{recon}}^{\mathrm{rel}}$	$\bar{E}_{\mathrm{traj}}(M_{\mathrm{val}})$	Complexity	$E_{\mathrm{res}}$	Rank	Selected
$S_{6}$	2	0.063176	0.000608	0.072231	6	0.105255	1	Yes
$S_{4}$	2	0.127328	0.003480	0.354430	4	0.212138	2	No
$F_{2}$	1	0.237008	0.003800	0.396505	2	0.394873	3	No
$S_{3}$	2	0.237008	0.003800	0.396505	4	0.394873	4	No
$S_{2}$	2	N/A	N/A	N/A	3	N/A	5	No
$S_{5}$	2	N/A	N/A	N/A	5	N/A	6	No

Table 16: (Model selection, SPPCS) Summary of the candidate models considered during model selection. For each candidate, we report the system order, framework compatibility diagnostics, validation trajectory error, model complexity, residual error, ranking, and whether the candidate was ultimately selected. In this case, the framework correctly identifies the SPPCS model with all three mechanisms active. N/A values indicate errors in trajectory reconstruction due to numerical instability of obtained estimators.

Order	$\phi^{E}$	$\phi^{A}$	$f$
1	✓	–	–

Table 17: Mechanisms of the phototaxis system ( 84 ). The symbol “-" indicates that the mechanism is not present, while “✓" indicates that the mechanism is present.

Panel of Figure 16 — Figure 16: (Model selection, opinion dynamics) Mechanisms present for the opinion dynamics system ( 84 ) within the generalized first-order framework ( 71 ). The model selection framework estimates an energy interaction kernel $\hat{\phi}^{E}$ and an environmental force $\hat{f}$ from observations of trajectories, and correctly identifies a first-order model dominated by interaction forces. Estimation is performed over $T_{r}=10$ independent learning trials, the inferred interaction kernel has weighted $L^{2}$ $\|\hat{\phi}^{E}\|_{L^{2}(\hat{\rho}_{R})}=$ $3.5858\text{\times}{10}^{-01}$ $\pm$ $1.7518\text{\times}{10}^{-02}$ (Figure 16(a) ), while the inferred environmental force has weighted $L^{2}$ norm $\|\hat{f}\|_{L^{2}(\hat{\rho}_{V})}=$ $0.0000\text{\times}{10}^{00}$ $\pm$ $0.0000\text{\times}{10}^{00}$ (Figure 16(b) ). Note the scale of the vertical axes. Further information regarding learning and model parameters may be found in Appendix E .

Candidate	Order	$E_{\mathrm{res}}^{\mathrm{rel_{\epsilon}}}$	$E_{\mathrm{recon}}^{\mathrm{rel}}$	$\bar{E}_{\mathrm{traj}}(M_{\mathrm{val}})$	Complexity	$E_{\mathrm{res}}$	Rank	Selected
$F_{1}$	1	0.186886	0.006227	0.016430	1	0.005228	1	Yes
$F_{2}$	1	0.186580	0.006351	0.015867	2	0.005220	2	No
$S_{2}$	2	0.930454	0.020764	0.099662	3	0.037002	3	No
$S_{5}$	2	0.927506	0.042013	0.147101	5	0.036885	4	No
$S_{1}$	2	0.983127	0.034845	0.335785	3	0.039097	5	No
$S_{3}$	2	N/A	N/A	N/A	4	N/A	6	No

Table 18: (Model selection, opinion dynamics) Summary of the candidate models considered during model selection. For each candidate, we report the system order, framework compatibility diagnostics, validation trajectory error, model complexity, residual error, ranking, and whether the candidate was ultimately selected. In this case, the framework correctly identifies the opinion dynamics model described by interactions only. N/A values indicate errors in trajectory reconstruction due to numerical instability of obtained estimators.

Appendix figures & tables18 assets

Supplementary material from the paper’s appendix.

Appendix

Panel of Figure 17 — Figure 17: (Kuramoto model) Linear B-spline basis functions used to construct hypothesis spaces. Subfigure 17(a) plots the interaction kernel basis functions $(\psi_{k})_{k=1}^{n_{\phi}}$ in $H_{\phi}$ with $n_{\phi}=10$ . Subfigure 17(b) plots the environmental force basis functions $(h_{k})_{k=1}^{n_{f}}$ in $H_{f}$ with $n_{f}=10$ .

Panel of Figure 18 — Figure 18: (Kuramoto model) Comparison of the empirical measures obtained from the $M_{\mathrm{train}}$ training data and the reference probability measures from $M_{\rho}$ samples, as discussed in Section 3.1 . See Table 22 for parameters utilized. Subfigure 18(a) plots $\hat{\rho}_{\Delta\theta}$ versus $\rho_{\Delta\theta}$ , and Subfigure 18(b) plots $\hat{\rho}_{\theta}$ versus $\rho_{\theta}$ . The close agreement suggests that the training trajectory data provides an accurate estimate for the regions explored by the dynamics.

Panel of Figure 19 — Figure 19: (Phototaxis model) Linear B-spline basis functions used to construct hypothesis spaces. Subfigure 19(a) plots the interaction kernel basis functions $(\psi_{k})_{k=1}^{n_{\phi}}$ in $H_{\phi}$ with $n_{\phi}=8$ . Subfigure 19(b) plots the environmental force basis functions $(h_{k})_{k=1}^{n_{f}}$ in $H_{f}$ , where the velocity space is discretized using $n_{b}=10$ basis functions in each dimension, resulting in $n_{f}=10\times 10=100$ tensor product basis functions for $\mathcal{H}_{f}$ .

Panel of Figure 20 — Figure 20: (Phototaxis model) Comparison of the empirical measures obtained from the $M_{\mathrm{train}}$ training data and the reference measure from $M_{\rho}$ samples, as discussed in Section 3.1 . See Table 22 for parameters utilized. Subfigure 20(a) plots $\hat{\rho}_{R}$ versus $\rho_{R}$ , and Subfigure 20(b) plots $\hat{\rho}_{V}$ versus $\rho_{V}$ . The close agreement suggests that the training trajectory data provides an accurate estimate for the regions explored by the dynamics.

Panel of Figure 21 — Figure 21: (SPP model) Linear B-spline basis functions used to construct the hypothesis spaces. Subfigure 21(a) plots the interaction kernel basis functions $(\psi_{k})_{k=1}^{n_{\phi}}$ in $H_{\phi}$ with $n_{\phi}=8$ . Subfigure 19(b) plots the environmental force basis functions $(h_{k})_{k=1}^{n_{f}}$ in $H_{f}$ . The velocity space is discretized using $n_{b}=10$ basis functions in each dimension, resulting in $n_{f}=10\times 10=100$ tensor product basis functions for $\mathcal{H}_{f}$ .

Panel of Figure 22 — Figure 22: (SPP) Comparison of the empirical measures obtained from the $M_{\mathrm{train}}$ training data and the reference measure from $M_{\rho}$ samples, as discussed in Section 3.1 . See Table 22 for parameters utilized. Subfigure 22(a) plots $\hat{\rho}_{R}$ versus $\rho_{R}$ , and Subfigure 22(b) plots $\hat{\rho}_{V}$ versus $\rho_{V}$ . The close agreement suggests that the training trajectory data provides an accurate estimate for the regions explored by the dynamics.

M	L	$E_{\phi}^{\mathrm{rel}}$	$E_{f}^{\mathrm{rel}}$
1	51	$1.3881\text{\times}{10}^{00}$\pm$2.2499\text{\times}{10}^{00}$	$3.6814\text{\times}{10}^{01}$\pm$6.4836\text{\times}{10}^{01}$
1	101	$1.8696\text{\times}{10}^{00}$\pm$2.0156\text{\times}{10}^{00}$	$4.8435\text{\times}{10}^{01}$\pm$5.6091\text{\times}{10}^{01}$
1	201	$2.8650\text{\times}{10}^{00}$\pm$3.9872\text{\times}{10}^{00}$	$6.4122\text{\times}{10}^{01}$\pm$8.9776\text{\times}{10}^{01}$
1	401	$1.7684\text{\times}{10}^{00}$\pm$2.3126\text{\times}{10}^{00}$	$4.2537\text{\times}{10}^{01}$\pm$5.9298\text{\times}{10}^{01}$
1	801	$2.6533\text{\times}{10}^{00}$\pm$3.9945\text{\times}{10}^{00}$	$7.0362\text{\times}{10}^{01}$\pm$1.1683\text{\times}{10}^{02}$
4	51	$1.2803\text{\times}{10}^{-02}$\pm$2.1729\text{\times}{10}^{-03}$	$8.1026\text{\times}{10}^{-02}$\pm$6.6889\text{\times}{10}^{-02}$

Appendix

Table 19: (Kuramoto model) Kuramoto relative interaction kernel and environmental force errors for different

M

(number of replicates) and

L

(number of time samples). Each entry reports the mean

\pm

sample standard deviation over

T_{r}=10

independent learning trials.

M	L	$E_{\phi}^{\mathrm{rel}}$	$E_{f}^{\mathrm{rel}}$
1	51	$6.1893\text{\times}{10}^{-03}$\pm$6.2809\text{\times}{10}^{-03}$	$6.1180\text{\times}{10}^{-02}$\pm$1.3356\text{\times}{10}^{-02}$
1	101	$4.4905\text{\times}{10}^{-03}$\pm$4.3087\text{\times}{10}^{-03}$	$1.1115\text{\times}{10}^{-01}$\pm$8.8252\text{\times}{10}^{-02}$
1	201	$1.3141\text{\times}{10}^{-02}$\pm$2.9193\text{\times}{10}^{-02}$	$8.5239\text{\times}{10}^{-02}$\pm$3.5008\text{\times}{10}^{-02}$
1	401	$5.4909\text{\times}{10}^{-03}$\pm$3.4038\text{\times}{10}^{-03}$	$8.5608\text{\times}{10}^{-02}$\pm$3.6644\text{\times}{10}^{-02}$
1	801	$8.6348\text{\times}{10}^{-03}$\pm$1.2248\text{\times}{10}^{-02}$	$7.7069\text{\times}{10}^{-02}$\pm$1.3662\text{\times}{10}^{-02}$
4	51	$1.9382\text{\times}{10}^{-04}$\pm$3.1394\text{\times}{10}^{-05}$	$1.2214\text{\times}{10}^{-02}$\pm$9.3570\text{\times}{10}^{-03}$

Appendix

Table 20: (Phototaxis model) Phototaxis model relative interaction kernel and environmental force errors for different

M

(number of replicates) and

L

(number of time samples). Each entry reports the mean

\pm

sample standard deviation over

T_{r}=10

independent learning trials.

M	L	$E_{\phi}^{\mathrm{rel}}$	$E_{f}^{\mathrm{rel}}$
1	51	$1.1381\text{\times}{10}^{00}$\pm$1.2043\text{\times}{10}^{00}$	$2.7564\text{\times}{10}^{00}$\pm$3.3192\text{\times}{10}^{00}$
1	101	$1.4278\text{\times}{10}^{00}$\pm$1.7454\text{\times}{10}^{00}$	$1.3453\text{\times}{10}^{00}$\pm$1.0674\text{\times}{10}^{00}$
1	201	$5.3246\text{\times}{10}^{-01}$\pm$5.6657\text{\times}{10}^{-01}$	$2.1987\text{\times}{10}^{00}$\pm$2.1506\text{\times}{10}^{00}$
1	401	$6.1445\text{\times}{10}^{-01}$\pm$5.5824\text{\times}{10}^{-01}$	$1.0916\text{\times}{10}^{00}$\pm$1.5755\text{\times}{10}^{00}$
1	801	$8.0262\text{\times}{10}^{-01}$\pm$3.6660\text{\times}{10}^{-01}$	$1.1428\text{\times}{10}^{00}$\pm$5.8255\text{\times}{10}^{-01}$
4	51	$1.3722\text{\times}{10}^{-01}$\pm$3.5001\text{\times}{10}^{-02}$	$1.9606\text{\times}{10}^{00}$\pm$4.5055\text{\times}{10}^{00}$

Appendix

Table 21: (SPP model) SPP model relative interaction kernel and environmental force errors for different

M

(number of replicates) and

L

(number of time samples). Each entry reports the mean

\pm

sample standard deviation over

T_{r}=10

independent learning trials.

Panel of Figure 23 — Figure 23: (Noise robustness, mechanism recovery) Each subfigure illustrates feature learning under multiplicative observational noise added to the training data at different noise levels $\zeta$ , as described in Section 3.5.4 . Subfigures 23(a) – 23(c) show the learned interaction kernels $\phi$ for increasing values of $\zeta$ . Subfigures 23(d) – 23(f) show the corresponding learned environmental force $f$ under the same noise levels. For the SPP and phototaxis systems, the recovery of $f$ is shown component-wise, where the upper and lower rows correspond to the first and second components, respectively.The background shading indicates the empirical training measure associated with each feature: $\hat{\rho}_{R}$ for interaction kernels in all systems except Kuramoto, where $\hat{\rho}_{\Delta\theta}$ is used; $\hat{\rho}_{V}$ for the environmental force in second-order systems; and $\hat{\rho}_{\theta}$ for the environmental variable in the Kuramoto system. We notice that, as the noise level $\zeta$ increases, the learned features gradually deviate from the true features, yet overall recovery is gnerally accurate and estimated features are minimally affected, primarily in regions with small induced measure (i.e. a small amount of trajectory data).

Panel of Figure 24 — Figure 24: (Noise robustness, trajectory recovery and prediction) Each row represents the Kuramoto, Phototaxis, and SPP systems, respectively. The subfigures show one of the observed trajectories before and after being perturbed by each noise level $\zeta$ described in Section 3.5.4 , over the full time horizon $[0,T_{f}]$ . For the Kuramoto system, the solid blue line represents the true trajectory, the dashed line represents the predicted trajectory learned from the noisy observations, and the semitransparent gray dots denote the noisy trajectory used as the training data (together with noisy observations of the derivative). For the SPP and phototaxis systems, the solid blue lines represent the true trajectories, the semitransparent gray dotted lines represent the noisy trajectories used as the training data (together with noisy observations of the derivative), and the color-gradient lines represent the predicted trajectories learned from the noisy observations, with yellow corresponding to $t_{0}$ and orange corresponding to $T_{f}$ .We observe that the learned trajectories remain close to the true trajectories at lower noise levels. As the noise level increases, gradual deviations in the predicted dynamics become apparent, particularly for the SPP system. However, the overall trajectory recovery remains stable, indicating a degree of robustness of the learned models to observational noise.

Panel of Figure 25 — Figure 25: (Initial condition distribution, SPP system) Comparison of the variational learning algorithm under two different initial condition distributions $\mu_{0}$ . Profile 1 (blue) is trained using a narrower initial-condition distribution, while Profile 2 (orange) is trained using a distribution with larger variance. The true kernel is shown in black. Subfigure 25(a) compares the recovered interaction kernel $\phi$ . Both profiles recover the overall kernel accurately; however, Profile 2 fails to capture the narrow, steep peak near the origin. Subfigure 25(b) compares the recovered environmental force $f$ . Profile 1 provides a poor approximation, whereas Profile 2 accurately recovers the true environmental force, demonstrating that a broader initial condition distribution improves recovery of the environmental dynamics.

Panel of Figure 26 — Figure 26: Effect of the number of B-spline basis functions on the recovery accuracy of the interaction kernel and environmental force for the SPP system. In each panel, one basis dimension is varied while the recovery error is averaged over all tested values of the other basis dimension. The solid curves represent the mean relative weighted $L^{2}$ error, and the shaded regions indicate one standard deviation across the corresponding configurations. Subfigure 26(a) shows the interaction kernel error, $E_{\phi}^{\mathrm{rel}}$ , as the number of kernel basis functions, $n_{\phi}$ , varies. Subfigure 26(b) shows the environmental-force error, $E_{f}^{\mathrm{rel}}$ , as the number of environmental force basis functions, $n_{f}$ , varies. In both cases, the lowest average recovery error is achieved with approximately ten basis functions.

Parameter	Kuramoto	Phototaxis	SPP	SPP	CS2D	SPPCS	OP
profile	-	-	clumps	milling	-	clumps	-
System order $s$	1	2	2	2	2	2	1
State dimension $d$	1	2	2	2	2	2	1
Number of agents $N$	10	10	10	10	10	10	10
training trajectories $M_{\mathrm{train}}$	20	20	50	50	20	50	20
testing trajectories $M_{\mathrm{test}}$	10	10	10	10	10	10	10

Appendix

Table 22: Data-generation parameters used in the numerical experiments. Training data are used to estimate the interaction kernels and environmental forces, testing data are used for trajectory validation, and the

M_{\rho}

data set is used to construct reference probability densities.

Parameter	Kuramoto	Phototaxis	SPP	CS2D	SPPCS	OP
Interaction type	Energy	Alignment	Energy	Alignment	Energy $+$ Alignment	Energy
Radial basis functions $n_{\phi}$	10	8	10	10	10	50
Energy basis functions $n_{\phi^{E}}$	–	10	–	10	10	50
Alignment basis functions $n_{\phi^{A}}$	–	10	–	10	10	50
Environment basis parameter $n_{b}$	10	10	10	10	10	10
Environment basis size per state $n_{f^{d}}=(n_{b})^{d}$	10	100	100	100	100	10

Appendix

Table 23: Algorithm parameters utilized in the variational least-squares framework. The interaction kernels are approximated using radial B-spline bases, while the environment force is represented using a tensor product of B-splines. For first-order systems the environment variable is the state, whereas for second-order systems it is generally the velocity.

System	Initial condition $u_{0}$
Kuramoto	$\theta_{0}^{(m,i)}{\sim}\mathcal{U}(0,2\pi)$
Phototaxis	$\Big[x_{0}^{(m,i)}{\sim}\mathcal{U}([0,2]^{2}),\;v_{0}^{(m,i)}\sim\mu_{v}^{\mathrm{Phottaxis}}\Big]$
SPP	$\Big[x_{0}^{(m,i)}\sim\mu_{x}^{\mathrm{SPP}},\;v_{0}^{(m,i)}\sim\mu_{v}^{\mathrm{SPP}}\Big]$
CS2D	$\Big[\,x_{0}^{(m,i)}\sim\mathcal{U}(0.1,10),v_{0}^{(m,i)}\sim\mathcal{U}(0.1,5)\Big]$
SPPCS	$\Big[x_{0}^{(m,i)}\sim\mu_{x}^{\mathrm{SPP}},\;v_{0}^{(m,i)}\sim\mu_{v}^{\mathrm{SPP}}\Big]$
OP	$x_{0}^{(m,i)}{\sim}\mathcal{U}(0,10)$

Appendix

Table 24: System-specific initial conditions used in numerical experiments. For the first-order Kuramoto system,

u_{0}=\theta_{0}

. For second-order systems,

u_{0}=(x_{0},v_{0})

\mu_{v}^{\mathrm{Phototaxis}}

denotes the perturbed tensor-grid measure described in E.3 The spatial measures

\mu_{x}^{\mathrm{SPP}}

and

\mu_{x}^{\mathrm{SPPCS}}

are cluster-based sampling distributions. The velocity measures

\mu_{v}^{\mathrm{SPP}}

and

\mu_{v}^{\mathrm{SPPCS}}

are perturbed tensor-grid distributions on

[-2\sqrt{2},2\sqrt{2}]^{2}

System	Component	Term	$coeff_{\mathrm{true}}$	$\mathbb{E}(\hat{coeff})\pm\sigma(\hat{coeff})$
Kuramoto	$1$	$\theta$	$0.01$	$0.009\,999\,999\,999\,999\,9$\pm$2.087\,722\,866\,931\,242\,4\text{\times}{10}^{-16}$
Phototaxis	$1$	$1$	$I_{0}U_{\infty}e_{\ell}^{(1)}=-\frac{0.1}{\sqrt{2}}$	$-0.070\,710\,681\,492\,353\,88$\pm$1.423\,730\,942\,286\,253\,9\text{\times}{10}^{-07}$
		$v_{1}$	$I_{0}=1.0$	$1.000\,000\,319\,453\,960\,8$\pm$3.538\,386\,027\,474\,474\text{\times}{10}^{-06}$
		$v_{2}$	$0$	$-2.967\,810\,805\,171\,508\text{\times}{10}^{-07}$\pm$2.676\,347\,740\,402\,615\text{\times}{10}^{-06}$
	$2$	$1$	$I_{0}U_{\infty}e_{\ell}^{(2)}=\frac{0.1}{\sqrt{2}}$	$0.070\,710\,663\,656\,462\,13$\pm$5.863\,941\,042\,682\,641\text{\times}{10}^{-08}$
		$v_{1}$	$0$	$-4.931\,084\,525\,950\,616\text{\times}{10}^{-07}$\pm$2.709\,497\,747\,258\,749\text{\times}{10}^{-06}$

Appendix

Table 25: Recovery of the parametric coefficients for the candidate basis functions. For each coefficient, we report the mean and standard deviation over

T_{r}=10

independent trials. Coefficients corresponding to non-present mechanisms are expected to be close to zero.

Parameter	Definition
Dynamical system
$N$	Number of agents in the system
$d$	Dimension of the state space of each agent
$x_{i}(t)\in\mathbb{R}^{d}$	State of the $i$ -th agent at time $t$
$x(t)\in\mathbb{R}^{dN}$	Full stacked state of the system at time $t$
$f:\mathbb{R}^{d}\to\mathbb{R}^{d}$	Environmental force function (assumed identical for all agents)

Appendix

Table 26: Parameters and their definitions.

Emergent Behavior in Multi-Agent Systems

Emergent behavior in multi-agent systems studies how local agent interactions produce collective patterns such as cooperation, norms, and instability.

See all 172 papers

Sep 10, 2026

Learning Interaction Kernels from Collective Steady States

We propose a learning procedure for system identification in interacting particle systems from single-snapshot observations of collective behaviors, unlike existing approaches that rely on observations of trajectories. This setting leads to a fundamentally ill-posed inverse problem, which we solve by using a regularization strategy based on the empirical distribution of observed configurations, drawn from different, unobserved initial conditions. We test our learning procedure on a variety of representative models with steady-state and quasi-stationary patterns, where collective behaviors encode implicit information about the interaction mechanisms, demonstrating that our approach enables stable and accurate recovery of the underlying interaction laws, leading to faithful reproduction of the collective behavior, and in many cases even of the dynamics leading up to it.

Baoli Hao, Mauro Maggioni, Ming Zhong

Department of Applied Mathematics, Illinois Institute of Technology, Chicago, 60616, IL, USA · Department of Mathematics, Department of Applied Mathematics and Statistics, Johns Hopkins University, Baltimore, 21218, MD, USA · Department of Mathematics, University of Houston, Houston, 77204, TX, USA

System Identification

Ask PDF

Mar 16, 2026

Discovery of interaction and diffusion kernels in particle-to-mean-field multi-agent systems

We propose a data-driven framework to learn interaction kernels in stochastic multi-agent systems. Our approach aims at identifying the functional form of nonlocal interaction and diffusion terms directly from trajectory data, without any a priori knowledge of the underlying interaction structure. Starting from a discrete stochastic binary-interaction model, we formulate the inverse problem as a sequence of sparse regression tasks in structured finite-dimensional spaces spanned by compactly supported basis functions, such as piecewise linear polynomials. In particular, we assume that pairwise interactions between agents are not directly observed and that only limited trajectory data are available. To address these challenges, we propose two complementary identification strategies. The first based on random-batch sampling, which compensates for latent interactions while preserving the statistical structure of the full dynamics in expectation. The second based on a mean-field approximation, where the empirical particle density reconstructed from the data defines a continuous nonlocal regression problem. Numerical experiments demonstrate the effectiveness and robustness of the proposed framework, showing accurate reconstruction of both interaction and diffusion kernels even from partially observed. The method is validated on benchmark models, including bounded-confidence and attraction-repulsion dynamics, where the two proposed strategies achieve comparable levels of accuracy.

Giacomo Albi, Alessandro Alla, Elisa Calzola

Department of Computer Science, University of Verona, Verona, Italy · Department of Mathematics, La Sapienza University of Rome, Rome, Italy

System Identification

Ask PDF

Jun 23, 2026

Swarm-Inspired Generation of Collective Behaviors in Graph Dynamical Systems

Collective behavior arises when locally interacting units produce coordinated global organization, from synchronization in dynamical systems to task-relevant information flow on graphs. The central challenge is not only to explain how collective behavior emerges, but to design local interaction rules that can produce desired global organization and generalize across graphs, dynamics and tasks.To address this challenge, we introduce the Swarm-Inspired Emergent Synchronizer (SIES), a graph-dynamical framework that learns generalizable local-interaction laws for controllable collective organization. Each node is an agent-like dynamical unit with a state and task cue, and signed source-target-conditioned attention acts as an adaptive coupling term inside an explicit evolution model. Therefore, SIES combines an explicit dynamical engine with local agent intelligence, similar to biological swarms. For synchronization control, SIES learns a generalizable coupling operator that produces prescribed synchronization patterns for CDSs across untrained network scales, target phase relations, and intrinsic node dynamics without retraining. The learned operator also reaches gait-related modes faster than three oscillator baselines and generalizes synchronization-driven locomotion to simulated multi-legged robots of different scales and a physical hexapod after leg disablement. For graph representation learning, SIES applies the same signed interaction principle to message passing and achieves the highest performance among the compared methods on heterophilous node-classification benchmarks. Together, these results position SIES as a generalizable and learnable graph-dynamical interaction framework with promise for synchronization control, adaptive robot coordination, and heterophilous graph representation learning.

Ji Chen, Song Chen, Chengzhang Gong +2

College of Control Science and Engineering, Zhejiang University, Hangzhou, 310000, Zhejiang, China. · Department of Mathematics, National University of Singapore, 10 Lower Kent Ridge Road, 119076, Singapore. · Huzhou Institute of Zhejiang University, Huzhou, 313000, Zhejiang, China. +3

Dynamical Systems Graph Representation Learning Multi-Robot Coordination

Ask PDF

Abstract

Emergent Behavior in Multi-Agent Systems

Explore similar work

Learning Interaction Kernels from Collective Steady States

Discovery of interaction and diffusion kernels in particle-to-mean-field multi-agent systems

Swarm-Inspired Generation of Collective Behaviors in Graph Dynamical Systems