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

Read PDF
Ask AI
Cite

Copy citation

DOINot available for this paper

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)c)Candidate ModelComplexity
F1F_{1}x˙=FϕE\dot{x}=F_{\phi^{E}}1
F2F_{2}x˙=Ff+FϕE\dot{x}=F_{f}+F_{\phi^{E}}2
S1S_{1}x¨=FϕE\ddot{x}=F_{\phi^{E}}3
S2S_{2}x¨=FϕA\ddot{x}=F_{\phi^{A}}3
S3S_{3}x¨=Ff+FϕE\ddot{x}=F_{f}+F_{\phi^{E}}4
S4S_{4}x¨=Ff+FϕA\ddot{x}=F_{f}+F_{\phi^{A}}4
Table 1: Hierarchy of candidate models used for model selection. In the Candidate column, FF denote a first-order framework, while SS denote the second-order framework. The notation utilized for vector fields in RdN\mathbb{R}^{dN} is as described in Sections 2.1.1 and Appendix A.1 (e.g. equation 7 ). We also suppress x˙i=vi\dot{x}_{i}=v_{i} for second-order systems for notational simplicity.
Data setFunction
Training ( MtrainM_{\mathrm{train}} samples)Obtain estimators for mechanistic forces for each model
Validation ( MvalM_{\mathrm{val}} samples)Compare candidate models and tune hyperparameters to select candidate model
Testing ( MtestM_{\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 MM is partitioned between the three classifications.
MetricDataTime horizonFunction
EreconrelE_{\mathrm{recon}}^{\mathrm{rel}}MtrainM_{\mathrm{train}}[0,T][0,T]Obtain admissible frameworks (trajectory data)
Eres,ϵrelE_{\mathrm{res},\epsilon}^{\mathrm{rel}}MtrainM_{\mathrm{train}}[0,T][0,T]Obtain admissible frameworks (residual data)
Eˉtraj\bar{E}_{\mathrm{traj}}MvalM_{\mathrm{val}}[0,Tf][0,T_{f}]Primary validation metric (trajectories)
EresE_{\mathrm{res}}MtrainM_{\mathrm{train}}[0,T][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.
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 ff and its recovered mean estimates f^ˉ\bar{\hat{f}} . Mean estimates are obtained over Tr=10T_{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=20M=20 trajectories with N=10N=10 oscillators in dimension d=1d=1 , observed over the interval [0,T]=[0,0.5][0,T]=[0,0.5] with L=51L=51 time points.
MetricNon-parametricSemi-Parametric
EresrelE^{\mathrm{rel}}_{\mathrm{res}}9.133 297 54\times10−039.133\,297\,54\text{\times}{10}^{-03} ±\pm 1.069 999 32\times10−031.069\,999\,32\text{\times}{10}^{-03}9.213 417 49\times10−039.213\,417\,49\text{\times}{10}^{-03} ±\pm 1.074 338 81\times10−031.074\,338\,81\text{\times}{10}^{-03}
EϕrelE^{\mathrm{rel}}_{\phi}1.151 949 00\times10−021.151\,949\,00\text{\times}{10}^{-02} ±\pm 6.567 552 90\times10−046.567\,552\,90\text{\times}{10}^{-04}1.150 350 07\times10−021.150\,350\,07\text{\times}{10}^{-02} ±\pm 6.505 249 25\times10−046.505\,249\,25\text{\times}{10}^{-04}
EfrelE^{\mathrm{rel}}_{f}1.992 020 20\times10−021.992\,020\,20\text{\times}{10}^{-02} ±\pm 5.681 235 00\times10−035.681\,235\,00\text{\times}{10}^{-03}1.712 688 22\times10−141.712\,688\,22\text{\times}{10}^{-14} ±\pm 1.184 873 90\times10−141.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 Tr=10T_{r}=10 independent learning trials. For additional details on algorithm parameters, we refer the reader to the caption of Figure 1 .
Figure 2: (Kuramoto model) Comparison of the true phase trajectories θi\theta_{i} (blue) and the estimated phase trajectories θ^i\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 MtrainM_{\mathrm{train}} training trajectories, whereas the bottom rows displays trajectories generated from MtestM_{\mathrm{test}} previously unseen testing initial conditions. Note that trajectories are plotted over the full time horizon [0,Tf][0,T_{f}] , with the end of the training interval time TT indicated on the horizontal axis. Recall that models are trained using observations collected only over the interval [0,T]=[0,0.5][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 setStatisticFully non-parametricSemi-parametric
MtrainM_{\mathrm{train}}Eˉrecon\bar{E}_{\mathrm{recon}}1.486 115 21\times10−031.486\,115\,21\text{\times}{10}^{-03} ±\pm 1.290 008 01\times10−041.290\,008\,01\text{\times}{10}^{-04}1.515 174 72\times10−031.515\,174\,72\text{\times}{10}^{-03} ±\pm 1.326 100 95\times10−041.326\,100\,95\text{\times}{10}^{-04}
σErecon\sigma_{E_{\mathrm{recon}}}7.059 050 63\times10−047.059\,050\,63\text{\times}{10}^{-04} ±\pm 1.149 350 14\times10−041.149\,350\,14\text{\times}{10}^{-04}7.258 467 47\times10−047.258\,467\,47\text{\times}{10}^{-04} ±\pm 1.074 815 46\times10−041.074\,815\,46\text{\times}{10}^{-04}
Eˉpred\bar{E}_{\mathrm{pred}}5.457 046 59\times10−035.457\,046\,59\text{\times}{10}^{-03} ±\pm 1.454 816 20\times10−031.454\,816\,20\text{\times}{10}^{-03}5.459 184 78\times10−035.459\,184\,78\text{\times}{10}^{-03} ±\pm 1.409 017 06\times10−031.409\,017\,06\text{\times}{10}^{-03}
σEpred\sigma_{E_{\mathrm{pred}}}5.665 476 87\times10−035.665\,476\,87\text{\times}{10}^{-03} ±\pm 2.895 039 32\times10−032.895\,039\,32\text{\times}{10}^{-03}6.020 934 52\times10−036.020\,934\,52\text{\times}{10}^{-03} ±\pm 2.947 784 88\times10−032.947\,784\,88\text{\times}{10}^{-03}
MtestM_{\mathrm{test}}Eˉrecon\bar{E}_{\mathrm{recon}}1.635 224 03\times10−031.635\,224\,03\text{\times}{10}^{-03} ±\pm 1.994 618 41\times10−041.994\,618\,41\text{\times}{10}^{-04}1.598 611 49\times10−031.598\,611\,49\text{\times}{10}^{-03} ±\pm 1.293 745 35\times10−041.293\,745\,35\text{\times}{10}^{-04}
σErecon\sigma_{E_{\mathrm{recon}}}9.230 112 42\times10−049.230\,112\,42\text{\times}{10}^{-04} ±\pm 2.374 169 98\times10−042.374\,169\,98\text{\times}{10}^{-04}9.094 150 34\times10−049.094\,150\,34\text{\times}{10}^{-04} ±\pm 2.159 352 77\times10−042.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 Eˉrecon\bar{E}_{\mathrm{recon}} , Eˉpred\bar{E}_{\mathrm{pred}} , σErecon\sigma_{E_{\mathrm{recon}}} , and σEpred\sigma_{E_{\mathrm{pred}}} are computed according to ( 63 ) - ( 64 ). The table reports the mean ±\pm one standard deviation of these quantities computed on Tr=10T_{r}=10 independent learning trials. The first section corresponds to trajectories generated from MtrainM_{\mathrm{train}} training initial conditions, while the second section corresponds to MtestM_{\mathrm{test}} unseen predicted initial conditions sampled from μ0\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 .
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 ff with the corresponding mean learned estimates f^ˉ\bar{\hat{f}} . As f:R2→R2f:\mathbb{R}^{2}\to\mathbb{R}^{2} , the top and bottom rows correspond to the first and second components of ff , 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 Tr=10T_{r}=10 independent learning trials. The background density in Subfigure 3(a) (right vertical axis) corresponds to the empirical pairwise distance distribution ρ^R\hat{\rho}_{R} , while the density visualized in Subfigure 3(b) (left column) corresponds to the empirical velocity distribution ρ^V\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=20M=20 trajectories with N=10N=10 agents in dimension d=2d=2 , observed over [0,T]=[0.0,0.5][0,T]=[0.0,0.5] with L=51L=51 time points.
MetricNon-parametricSemi-Parametric
EresrelE^{\mathrm{rel}}_{\mathrm{res}}3.250 636 15\times10−023.250\,636\,15\text{\times}{10}^{-02} ±\pm 3.804 678 17\times10−043.804\,678\,17\text{\times}{10}^{-04}1.464 225 90\times10−021.464\,225\,90\text{\times}{10}^{-02} ±\pm 7.416 366 49\times10−047.416\,366\,49\text{\times}{10}^{-04}
EϕrelE^{\mathrm{rel}}_{\phi}1.199 870 60\times10−011.199\,870\,60\text{\times}{10}^{-01} ±\pm 1.993 494 02\times10−031.993\,494\,02\text{\times}{10}^{-03}1.197 529 15\times10−011.197\,529\,15\text{\times}{10}^{-01} ±\pm 1.841 751 39\times10−031.841\,751\,39\text{\times}{10}^{-03}
EfrelE^{\mathrm{rel}}_{f}2.974 167 42\times10−022.974\,167\,42\text{\times}{10}^{-02} ±\pm 3.449 852 01\times10−043.449\,852\,01\text{\times}{10}^{-04}6.480 458 97\times10−046.480\,458\,97\text{\times}{10}^{-04} ±\pm 1.344 640 77\times10−041.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 Tr=10T_{r}=10 independent learning trials. For additional details on algorithm parameters, we refer the reader to the caption of Figure 3 .
Figure 4: (SPP model) Comparison of the true trajectories xx (left column) and the estimated trajectories x^\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 MtrainM_{\mathrm{train}} training initial-conditions, whereas the bottom row displays trajectories generated from MtestM_{\mathrm{test}} previously unseen testing initial conditions. Note that trajectories are plotted over the full time horizon [0,Tf][0,T_{f}] , with the color gradient representing the temporal evolution of the dynamics, from yellow at t=0t=0 to orange at t=Tft=T_{f} . The black segment indicates the training interval [0,T][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][0,T]=[0.0,0.5] , with testing evaluated beyond the fitting horizon [0.0,0.5][0.0,0.5] until Tf=2.0T_{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 setStatisticFully non-parametricSemi-parametric
MtrainM_{\mathrm{train}}Eˉrecon\bar{E}_{\mathrm{recon}}4.291 656 12\times10−034.291\,656\,12\text{\times}{10}^{-03} ±\pm 7.329 969 06\times10−057.329\,969\,06\text{\times}{10}^{-05}2.762 636 52\times10−032.762\,636\,52\text{\times}{10}^{-03} ±\pm 8.642 627 33\times10−058.642\,627\,33\text{\times}{10}^{-05}
σErecon\sigma_{E_{\mathrm{recon}}}1.018 780 56\times10−031.018\,780\,56\text{\times}{10}^{-03} ±\pm 1.047 847 99\times10−041.047\,847\,99\text{\times}{10}^{-04}9.363 603 47\times10−049.363\,603\,47\text{\times}{10}^{-04} ±\pm 1.232 918 49\times10−041.232\,918\,49\text{\times}{10}^{-04}
Eˉpred\bar{E}_{\mathrm{pred}}4.969 235 82\times10−024.969\,235\,82\text{\times}{10}^{-02} ±\pm 2.224 500 17\times10−032.224\,500\,17\text{\times}{10}^{-03}2.484 601 58\times10−022.484\,601\,58\text{\times}{10}^{-02} ±\pm 1.481 024 68\times10−031.481\,024\,68\text{\times}{10}^{-03}
σEpred\sigma_{E_{\mathrm{pred}}}9.785 216 05\times10−039.785\,216\,05\text{\times}{10}^{-03} ±\pm 9.877 813 09\times10−049.877\,813\,09\text{\times}{10}^{-04}9.042 348 73\times10−039.042\,348\,73\text{\times}{10}^{-03} ±\pm 1.314 571 31\times10−031.314\,571\,31\text{\times}{10}^{-03}
MtestM_{\mathrm{test}}Eˉrecon\bar{E}_{\mathrm{recon}}4.256 561 62\times10−034.256\,561\,62\text{\times}{10}^{-03} ±\pm 1.291 396 91\times10−041.291\,396\,91\text{\times}{10}^{-04}2.767 011 67\times10−032.767\,011\,67\text{\times}{10}^{-03} ±\pm 1.018 333 24\times10−041.018\,333\,24\text{\times}{10}^{-04}
σErecon\sigma_{E_{\mathrm{recon}}}1.133 254 76\times10−031.133\,254\,76\text{\times}{10}^{-03} ±\pm 1.405 147 45\times10−041.405\,147\,45\text{\times}{10}^{-04}8.268 659 24\times10−048.268\,659\,24\text{\times}{10}^{-04} ±\pm 9.900 607 36\times10−059.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 Eˉrecon\bar{E}_{\mathrm{recon}} , Eˉpred\bar{E}_{\mathrm{pred}} , σErecon\sigma_{E_{\mathrm{recon}}} , and σEpred\sigma_{E_{\mathrm{pred}}} are computed according to ( 63 ) - ( 64 ). The table reports the mean ±\pm one standard deviation of these quantities computed over Tr=10T_{r}=10 independent learning trials. The first section corresponds to trajectories generated from MtrainM_{\mathrm{train}} training initial conditions, while the second section corresponds to MtestM_{\mathrm{test}} unseen initial conditions sampled from μ0\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 3 .
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 ff . Here ff is assumed of the form f(v;p)=p1v+p2∣∣v∣∣1vf(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 MtrainM_{\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,Tf][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 ff , 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 .
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 ff and its recovered mean estimates f^ˉ\bar{\hat{f}} . As f:R2→R2f:\mathbb{R}^{2}\to\mathbb{R}^{2} , the top and bottom rows correspond to first and second components of ff , 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 Tr=10T_{r}=10 independent learning trials. The background density in Subfigure 6(a) (right vertical axis) corresponds to the empirical pairwise distance distribution ρ^R\hat{\rho}_{R} , while the density visualized in Subfigure 6(b) (left column) corresponds to the empirical velocity distribution ρ^V\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=20M=20 trajectories with N=10N=10 agents in dimension d=2d=2 , observed over [0,T]=[0,0.5][0,T]=[0,0.5] with L=51L=51 time points.
MetricNon-parametricSemi-Parametric
EresrelE^{\mathrm{rel}}_{\mathrm{res}}6.726 628 52\times10−056.726\,628\,52\text{\times}{10}^{-05} ±\pm 6.398 876 10\times10−066.398\,876\,10\text{\times}{10}^{-06}7.291 303 63\times10−057.291\,303\,63\text{\times}{10}^{-05} ±\pm 7.334 376 82\times10−067.334\,376\,82\text{\times}{10}^{-06}
EϕrelE^{\mathrm{rel}}_{\phi}1.985 512 11\times10−041.985\,512\,11\text{\times}{10}^{-04} ±\pm 2.292 855 28\times10−052.292\,855\,28\text{\times}{10}^{-05}1.980 978 68\times10−041.980\,978\,68\text{\times}{10}^{-04} ±\pm 2.288 845 60\times10−052.288\,845\,60\text{\times}{10}^{-05}
EfrelE^{\mathrm{rel}}_{f}5.533 440 49\times10−055.533\,440\,49\text{\times}{10}^{-05} ±\pm 8.754 940 29\times10−068.754\,940\,29\text{\times}{10}^{-06}3.738 354 94\times10−063.738\,354\,94\text{\times}{10}^{-06} ±\pm 2.418 218 40\times10−062.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 Tr=10T_{r}=10 independent learning trials. For additional details on algorithm parameters, we refer the reader to the caption of Figure 6 .
Figure 7: (Phototaxis model) Comparison of true trajectories xx (left column) and the estimated trajectories x^\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 MtrainM_{\mathrm{train}} training initial-conditions, whereas the bottom row displays trajectories generated from MtestM_{\mathrm{test}} previously unseen testing initial conditions. Note that trajectories are plotted over the full time horizon [0,Tf][0,T_{f}] , with the color gradient representing the temporal evolution of the dynamics, from yellow at t=0t=0 to orange at t=Tft=T_{f} . The black segment indicates the training interval [0,T][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][0,T]=[0,0.5] , with testing evaluated beyond the fitting horizon [0.0,0.5][0.0,0.5] until Tf=2.0T_{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 setStatisticNon-parametricSemi-Parametric
MtrainM_{\mathrm{train}}Eˉrecon\bar{E}_{\mathrm{recon}}4.027 552 51\times10−064.027\,552\,51\text{\times}{10}^{-06} ±\pm 6.549 230 71\times10−076.549\,230\,71\text{\times}{10}^{-07}4.496 428 52\times10−064.496\,428\,52\text{\times}{10}^{-06} ±\pm 7.364 670 99\times10−077.364\,670\,99\text{\times}{10}^{-07}
σErecon\sigma_{E_{\mathrm{recon}}}1.009 259 16\times10−061.009\,259\,16\text{\times}{10}^{-06} ±\pm 3.055 190 48\times10−073.055\,190\,48\text{\times}{10}^{-07}1.268 869 62\times10−061.268\,869\,62\text{\times}{10}^{-06} ±\pm 3.899 608 97\times10−073.899\,608\,97\text{\times}{10}^{-07}
Eˉpred\bar{E}_{\mathrm{pred}}1.367 294 68\times10−051.367\,294\,68\text{\times}{10}^{-05} ±\pm 2.262 020 00\times10−062.262\,020\,00\text{\times}{10}^{-06}1.178 536 84\times10−051.178\,536\,84\text{\times}{10}^{-05} ±\pm 1.560 367 45\times10−061.560\,367\,45\text{\times}{10}^{-06}
σEpred\sigma_{E_{\mathrm{pred}}}2.768 822 80\times10−062.768\,822\,80\text{\times}{10}^{-06} ±\pm 7.388 936 98\times10−077.388\,936\,98\text{\times}{10}^{-07}3.210 520 55\times10−063.210\,520\,55\text{\times}{10}^{-06} ±\pm 6.862 167 86\times10−076.862\,167\,86\text{\times}{10}^{-07}
MtestM_{\mathrm{test}}Eˉrecon\bar{E}_{\mathrm{recon}}5.067 807 84\times10−065.067\,807\,84\text{\times}{10}^{-06} ±\pm 6.953 979 44\times10−076.953\,979\,44\text{\times}{10}^{-07}4.523 697 17\times10−064.523\,697\,17\text{\times}{10}^{-06} ±\pm 5.193 499 62\times10−075.193\,499\,62\text{\times}{10}^{-07}
σErecon\sigma_{E_{\mathrm{recon}}}1.163 612 56\times10−061.163\,612\,56\text{\times}{10}^{-06} ±\pm 3.301 187 82\times10−073.301\,187\,82\text{\times}{10}^{-07}1.012 524 41\times10−061.012\,524\,41\text{\times}{10}^{-06} ±\pm 3.038 581 17\times10−073.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 Eˉrecon\bar{E}_{\mathrm{recon}} , Eˉpred\bar{E}_{\mathrm{pred}} , σErecon\sigma_{E_{\mathrm{recon}}} , and σEpred\sigma_{E_{\mathrm{pred}}} are computed according to ( 63 ) - ( 64 ). The table reports the mean ±\pm one standard deviation of these quantities computed over Tr=10T_{r}=10 independent learning trials. The first section corresponds to trajectories generated from MtrainM_{\mathrm{train}} training initial conditions, while the second section corresponds to MtestM_{\mathrm{test}} unseen initial conditions sampled from μ0\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 6 .
Data setStatisticNon-parametricSemi-parametric
KuramotoEˉrecon\bar{E}_{\mathrm{recon}}8.639 226 09\times10−048.639\,226\,09\text{\times}{10}^{-04} ±\pm 1.436 067 53\times10−041.436\,067\,53\text{\times}{10}^{-04}7.933 406 60\times10−047.933\,406\,60\text{\times}{10}^{-04} ±\pm 9.072 559 42\times10−059.072\,559\,42\text{\times}{10}^{-05}
σErecon\sigma_{E_{\mathrm{recon}}}2.033 986 40\times10−042.033\,986\,40\text{\times}{10}^{-04} ±\pm 2.378 269 91\times10−052.378\,269\,91\text{\times}{10}^{-05}2.029 332 07\times10−042.029\,332\,07\text{\times}{10}^{-04} ±\pm 2.688 566 76\times10−052.688\,566\,76\text{\times}{10}^{-05}
Eˉpred\bar{E}_{\mathrm{pred}}3.360 057 73\times10−033.360\,057\,73\text{\times}{10}^{-03} ±\pm 1.260 603 89\times10−031.260\,603\,89\text{\times}{10}^{-03}2.721 694 99\times10−032.721\,694\,99\text{\times}{10}^{-03} ±\pm 9.891 757 78\times10−049.891\,757\,78\text{\times}{10}^{-04}
σEpred\sigma_{E_{\mathrm{pred}}}1.502 954 18\times10−031.502\,954\,18\text{\times}{10}^{-03} ±\pm 7.714 561 34\times10−047.714\,561\,34\text{\times}{10}^{-04}1.372 298 20\times10−031.372\,298\,20\text{\times}{10}^{-03} ±\pm 6.667 325 58\times10−046.667\,325\,58\text{\times}{10}^{-04}
SPPEˉrecon\bar{E}_{\mathrm{recon}}7.886 463 53\times10−037.886\,463\,53\text{\times}{10}^{-03} ±\pm 2.475 329 39\times10−042.475\,329\,39\text{\times}{10}^{-04}0.006 884 921 442 267 940.006\,884\,921\,442\,267\,94 ±\pm 0.000 253 715 904 779 536 20.000\,253\,715\,904\,779\,536\,2
σErecon\sigma_{E_{\mathrm{recon}}}1.737 702 05\times10−031.737\,702\,05\text{\times}{10}^{-03} ±\pm 4.804 516 38\times10−044.804\,516\,38\text{\times}{10}^{-04}0.002 113 167 912 756 090 80.002\,113\,167\,912\,756\,090\,8 ±\pm 0.000 543 637 710 861 935 60.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 Eˉrecon\bar{E}_{\mathrm{recon}} , Eˉpred\bar{E}_{\mathrm{pred}} , σErecon\sigma_{E_{\mathrm{recon}}} , and σEpred\sigma_{E_{\mathrm{pred}}} are computed according to ( 63 ) - ( 64 ). The table reports the mean ±\pm one standard deviation of these quantities computed Tr=10T_{r}=10 independent learning trials. Predicted trajectories are made with Mtest4NM_{\mathrm{test}}^{4N} untrained initial conditions generated for 4N4N agents. In numerical simulations, Mtest4N=MtestM_{\mathrm{test}}^{4N}=M_{\mathrm{test}} . Further information regarding learning and model parameters may be found in Appendix E .
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 Nnew=4NN_{\mathrm{new}}=4N agents, when N=10N=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 Cr=1,lr=0.5,Ca=0.5,la=2,α=1.6C_{r}=1,l_{r}=0.5,C_{a}=0.5,l_{a}=2,\alpha=1.6 , and β=0.5\beta=0.5 . The top row shows the trajectory predictions over the training time interval [0,T]=[0,4][0,T]=[0,4] for N=10N=10 agents. The bottom row shows the trajectory predictions for a scaled system with Nnew=100N_{\mathrm{new}}=100 agents, evolved over a time horizon five times longer than the training interval, with Tf=20T_{f}=20 . Further information regarding the learning parameters, system parameters, and initial conditions may be found in Appendix E .
Figure 10: Mean mechanistic relative feature recovery errors (log scale) as functions of temporal resolution LL and independent replicate size MM computed over Tr=10T_{r}=10 independent learning trials. The relatively large errors observed for M=1M=1 replicate, even when increasing the number of time samples LL , 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 ff (Figure 10(b) ), suggests that the number of replicates MM 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 .
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 ff . The estimated empirical measures ρ^R\hat{\rho}_{R} and ρ^V\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 ζ=0.01\zeta=0.01 and ζ=0.1\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][0,T] , and testing (predictions) are provided on (T,Tf](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 .
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 ζ=0.01\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 ff 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ϕE\phi^{E}ϕA\phi^{A}ff
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.
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 ϕ^E\hat{\phi}^{E} , an alignment interaction kernel ϕ^A\hat{\phi}^{A} , and an environmental force f^\hat{f} . From observations of trajectories, the algorithm correctly identify the system as alignment-dominated. Estimation is performed over Tr=10T_{r}=10 independent learning trials, and the inferred alignment kernel (Figure 13(b) ) has weighted L2L^{2} norm ∥ϕ^A∥L2(ρ^R)=\|\hat{\phi}^{A}\|_{L^{2}(\hat{\rho}_{R})}= 0.742 1530.742\,153 ±\pm 0.003 4680.003\,468 , which is several orders of magnitude larger than the inferred energy kernel ∥ϕ^E∥L2(ρ^R)=\|\hat{\phi}^{E}\|_{L^{2}(\hat{\rho}_{R})}= 1.9985\times10−041.9985\text{\times}{10}^{-04} ±\pm 1.0221\times10−041.0221\text{\times}{10}^{-04} (Figure 13(a) ) and environmental force ∥f^∥L2(ρ^V)=\|\hat{f}\|_{L^{2}(\hat{\rho}_{V})}= 0.000 5280.000\,528 ±\pm 0.000 0700.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 .
CandidateOrderEresrelϵE_{\mathrm{res}}^{\mathrm{rel}_{\epsilon}}EreconrelE_{\mathrm{recon}}^{\mathrm{rel}}Eˉtraj(Mval)\bar{E}_{\mathrm{traj}}(M_{\mathrm{val}})ComplexityEresE_{\mathrm{res}}RankSelected
S2S_{2}20.0006260.0000040.00025930.0005091Yes
S4S_{4}20.0007080.0000040.00043740.0005752No
S5S_{5}20.0006240.0000040.00025550.0005073No
S6S_{6}20.0007170.0000040.00046460.0005834No
F2F_{2}10.5937260.0032660.73206120.4824665No
S3S_{3}20.5937260.0032660.73206140.4824666No
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ϕE\phi^{E}ϕA\phi^{A}ff
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.
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 ϕ^E\hat{\phi}^{E} , an alignment interaction kernel ϕ^A\hat{\phi}^{A} , and an environmental force f^\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 Tr=10T_{r}=10 independent learning trials, and the inferred learned alignment kernel (Figure 14(b) ) has weighted L2L^{2} norm ∥ϕ^A∥L2(ρ^R)=\|\hat{\phi}^{A}\|_{L^{2}(\hat{\rho}_{R})}= 0.931 0350.931\,035 ±\pm 0.002 0160.002\,016 and the estimated environmental force has weighted L2L^{2} norm ∥f^∥L2(ρ^V)=\|\hat{f}\|_{L^{2}(\hat{\rho}_{V})}= 0.600 5040.600\,504 ±\pm 0.006 2110.006\,211 (Figure 14(c) ), both of which are several orders of magnitude larger than the inferred energy kernel ∥ϕ^E∥L2(ρ^R)=\|\hat{\phi}^{E}\|_{L^{2}(\hat{\rho}_{R})}= 7.7845\times10−047.7845\text{\times}{10}^{-04} ±\pm 6.9402\times10−046.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 .
CandidateOrderEresrelϵE_{\mathrm{res}}^{\mathrm{rel}_{\epsilon}}EreconrelE_{\mathrm{recon}}^{\mathrm{rel}}Eˉtraj(Mval)\bar{E}_{\mathrm{traj}}(M_{\mathrm{val}})ComplexityEresE_{\mathrm{res}}RankSelected
S4S_{4}20.0000670.0000020.00001640.0000401Yes
S6S_{6}20.0000680.0000020.00001760.0000412No
F2F_{2}10.2164970.0090560.10212420.1303413No
S3S_{3}20.2164970.0090560.10212440.1303414No
S2S_{2}20.2241090.0156540.27298130.1349245No
S5S_{5}20.2241090.0156540.27298150.1349246No
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ϕE\phi^{E}ϕA\phi^{A}ff
2✓✓✓
Table 15: Mechanisms of the SPPCS system. The symbol “-" indicates that the mechanism is not present, while “✓" indicates that the mechanism is present.
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 ϕ^E\hat{\phi}^{E} , an alignment interaction kernel ϕ^A\hat{\phi}^{A} , and an environmental force f^\hat{f} . From observations of trajectories, the algorithm correctly identify the system as being governed by all three mechanisms. Estimation is performed over Tr=10T_{r}=10 independent learning trials, and the inferred alignment kernel has weighted L2L^{2} norm ∥ϕ^A∥L2(ρ^R)=\|\hat{\phi}^{A}\|_{L^{2}(\hat{\rho}_{R})}= 0.344 5570.344\,557 ±\pm 0.031 9420.031\,942 (Figure 15(b) ), the inferred environmental force has weighted L2L^{2} norm ∥f^∥L2(ρ^V)=\|\hat{f}\|_{L^{2}(\hat{\rho}_{V})}= 3.062 9653.062\,965 ±\pm 0.047 2180.047\,218 (Figure 15(c) ), and the inferred energy kernel has weighted L2L^{2} norm ∥ϕ^E∥L2(ρR)=\|\hat{\phi}^{E}\|_{L^{2}(\rho_{R})}= 6.5314\times10−016.5314\text{\times}{10}^{-01} ±\pm 6.0227\times10−026.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 .
CandidateOrderEresrelϵE_{\mathrm{res}}^{\mathrm{rel_{\epsilon}}}EreconrelE_{\mathrm{recon}}^{\mathrm{rel}}Eˉtraj(Mval)\bar{E}_{\mathrm{traj}}(M_{\mathrm{val}})ComplexityEresE_{\mathrm{res}}RankSelected
S6S_{6}20.0631760.0006080.07223160.1052551Yes
S4S_{4}20.1273280.0034800.35443040.2121382No
F2F_{2}10.2370080.0038000.39650520.3948733No
S3S_{3}20.2370080.0038000.39650540.3948734No
S2S_{2}2N/AN/AN/A3N/A5No
S5S_{5}2N/AN/AN/A5N/A6No
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ϕE\phi^{E}ϕA\phi^{A}ff
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.
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 ϕ^E\hat{\phi}^{E} and an environmental force f^\hat{f} from observations of trajectories, and correctly identifies a first-order model dominated by interaction forces. Estimation is performed over Tr=10T_{r}=10 independent learning trials, the inferred interaction kernel has weighted L2L^{2} ∥ϕ^E∥L2(ρ^R)=\|\hat{\phi}^{E}\|_{L^{2}(\hat{\rho}_{R})}= 3.5858\times10−013.5858\text{\times}{10}^{-01} ±\pm 1.7518\times10−021.7518\text{\times}{10}^{-02} (Figure 16(a) ), while the inferred environmental force has weighted L2L^{2} norm ∥f^∥L2(ρ^V)=\|\hat{f}\|_{L^{2}(\hat{\rho}_{V})}= 0.0000\times10000.0000\text{\times}{10}^{00} ±\pm 0.0000\times10000.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 .
CandidateOrderEresrelϵE_{\mathrm{res}}^{\mathrm{rel_{\epsilon}}}EreconrelE_{\mathrm{recon}}^{\mathrm{rel}}Eˉtraj(Mval)\bar{E}_{\mathrm{traj}}(M_{\mathrm{val}})ComplexityEresE_{\mathrm{res}}RankSelected
F1F_{1}10.1868860.0062270.01643010.0052281Yes
F2F_{2}10.1865800.0063510.01586720.0052202No
S2S_{2}20.9304540.0207640.09966230.0370023No
S5S_{5}20.9275060.0420130.14710150.0368854No
S1S_{1}20.9831270.0348450.33578530.0390975No
S3S_{3}2N/AN/AN/A4N/A6No
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

Figure 17: (Kuramoto model) Linear B-spline basis functions used to construct hypothesis spaces. Subfigure 17(a) plots the interaction kernel basis functions (ψk)k=1nϕ(\psi_{k})_{k=1}^{n_{\phi}} in HϕH_{\phi} with nϕ=10n_{\phi}=10 . Subfigure 17(b) plots the environmental force basis functions (hk)k=1nf(h_{k})_{k=1}^{n_{f}} in HfH_{f} with nf=10n_{f}=10 .
Figure 18: (Kuramoto model) Comparison of the empirical measures obtained from the MtrainM_{\mathrm{train}} training data and the reference probability measures from Mρ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.
Figure 19: (Phototaxis model) Linear B-spline basis functions used to construct hypothesis spaces. Subfigure 19(a) plots the interaction kernel basis functions (ψk)k=1nϕ(\psi_{k})_{k=1}^{n_{\phi}} in HϕH_{\phi} with nϕ=8n_{\phi}=8 . Subfigure 19(b) plots the environmental force basis functions (hk)k=1nf(h_{k})_{k=1}^{n_{f}} in HfH_{f} , where the velocity space is discretized using nb=10n_{b}=10 basis functions in each dimension, resulting in nf=10×10=100n_{f}=10\times 10=100 tensor product basis functions for Hf\mathcal{H}_{f} .
Figure 20: (Phototaxis model) Comparison of the empirical measures obtained from the MtrainM_{\mathrm{train}} training data and the reference measure from MρM_{\rho} samples, as discussed in Section 3.1 . See Table 22 for parameters utilized. Subfigure 20(a) plots ρ^R\hat{\rho}_{R} versus ρR\rho_{R} , and Subfigure 20(b) plots ρ^V\hat{\rho}_{V} versus ρV\rho_{V} . The close agreement suggests that the training trajectory data provides an accurate estimate for the regions explored by the dynamics.
Figure 21: (SPP model) Linear B-spline basis functions used to construct the hypothesis spaces. Subfigure 21(a) plots the interaction kernel basis functions (ψk)k=1nϕ(\psi_{k})_{k=1}^{n_{\phi}} in HϕH_{\phi} with nϕ=8n_{\phi}=8 . Subfigure 19(b) plots the environmental force basis functions (hk)k=1nf(h_{k})_{k=1}^{n_{f}} in HfH_{f} . The velocity space is discretized using nb=10n_{b}=10 basis functions in each dimension, resulting in nf=10×10=100n_{f}=10\times 10=100 tensor product basis functions for Hf\mathcal{H}_{f} .
Figure 22: (SPP) Comparison of the empirical measures obtained from the MtrainM_{\mathrm{train}} training data and the reference measure from MρM_{\rho} samples, as discussed in Section 3.1 . See Table 22 for parameters utilized. Subfigure 22(a) plots ρ^R\hat{\rho}_{R} versus ρR\rho_{R} , and Subfigure 22(b) plots ρ^V\hat{\rho}_{V} versus ρV\rho_{V} . The close agreement suggests that the training trajectory data provides an accurate estimate for the regions explored by the dynamics.
MLEϕrelE_{\phi}^{\mathrm{rel}}EfrelE_{f}^{\mathrm{rel}}
1511.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}
11011.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}
12012.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}
14011.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}
18012.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}
4511.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 MM (number of replicates) and LL (number of time samples). Each entry reports the mean ±\pm sample standard deviation over Tr=10T_{r}=10 independent learning trials.
MLEϕrelE_{\phi}^{\mathrm{rel}}EfrelE_{f}^{\mathrm{rel}}
1516.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}
11014.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}
12011.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}
14015.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}
18018.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}
4511.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 MM (number of replicates) and LL (number of time samples). Each entry reports the mean ±\pm sample standard deviation over Tr=10T_{r}=10 independent learning trials.
MLEϕrelE_{\phi}^{\mathrm{rel}}EfrelE_{f}^{\mathrm{rel}}
1511.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}
11011.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}
12015.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}
14016.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}
18018.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}
4511.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 MM (number of replicates) and LL (number of time samples). Each entry reports the mean ±\pm sample standard deviation over Tr=10T_{r}=10 independent learning trials.
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 ff under the same noise levels. For the SPP and phototaxis systems, the recovery of ff 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: ρ^R\hat{\rho}_{R} for interaction kernels in all systems except Kuramoto, where ρ^Δθ\hat{\rho}_{\Delta\theta} is used; ρ^V\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).
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,Tf][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 t0t_{0} and orange corresponding to TfT_{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.
Figure 25: (Initial condition distribution, SPP system) Comparison of the variational learning algorithm under two different initial condition distributions μ0\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 ff . 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.
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 L2L^{2} error, and the shaded regions indicate one standard deviation across the corresponding configurations. Subfigure 26(a) shows the interaction kernel error, EϕrelE_{\phi}^{\mathrm{rel}} , as the number of kernel basis functions, nϕn_{\phi} , varies. Subfigure 26(b) shows the environmental-force error, EfrelE_{f}^{\mathrm{rel}} , as the number of environmental force basis functions, nfn_{f} , varies. In both cases, the lowest average recovery error is achieved with approximately ten basis functions.
ParameterKuramotoPhototaxisSPPSPPCS2DSPPCSOP
profile--clumpsmilling-clumps-
System order ss1222221
State dimension dd1222221
Number of agents NN10101010101010
training trajectories MtrainM_{\mathrm{train}}20205050205020
testing trajectories MtestM_{\mathrm{test}}10101010101010
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ρM_{\rho} data set is used to construct reference probability densities.
ParameterKuramotoPhototaxisSPPCS2DSPPCSOP
Interaction typeEnergyAlignmentEnergyAlignmentEnergy ++ AlignmentEnergy
Radial basis functions nϕn_{\phi}10810101050
Energy basis functions nϕEn_{\phi^{E}}–10–101050
Alignment basis functions nϕAn_{\phi^{A}}–10–101050
Environment basis parameter nbn_{b}101010101010
Environment basis size per state nfd=(nb)dn_{f^{d}}=(n_{b})^{d}1010010010010010
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.
SystemInitial condition u0u_{0}
Kuramotoθ0(m,i)∼U(0,2π)\theta_{0}^{(m,i)}{\sim}\mathcal{U}(0,2\pi)
Phototaxis[x0(m,i)∼U([0,2]2),  v0(m,i)∼μvPhottaxis]\Big[x_{0}^{(m,i)}{\sim}\mathcal{U}([0,2]^{2}),\;v_{0}^{(m,i)}\sim\mu_{v}^{\mathrm{Phottaxis}}\Big]
SPP[x0(m,i)∼μxSPP,  v0(m,i)∼μvSPP]\Big[x_{0}^{(m,i)}\sim\mu_{x}^{\mathrm{SPP}},\;v_{0}^{(m,i)}\sim\mu_{v}^{\mathrm{SPP}}\Big]
CS2D[ x0(m,i)∼U(0.1,10),v0(m,i)∼U(0.1,5)]\Big[\,x_{0}^{(m,i)}\sim\mathcal{U}(0.1,10),v_{0}^{(m,i)}\sim\mathcal{U}(0.1,5)\Big]
SPPCS[x0(m,i)∼μxSPP,  v0(m,i)∼μvSPP]\Big[x_{0}^{(m,i)}\sim\mu_{x}^{\mathrm{SPP}},\;v_{0}^{(m,i)}\sim\mu_{v}^{\mathrm{SPP}}\Big]
OPx0(m,i)∼U(0,10)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, u0=θ0u_{0}=\theta_{0} . For second-order systems, u0=(x0,v0)u_{0}=(x_{0},v_{0}) . μvPhototaxis\mu_{v}^{\mathrm{Phototaxis}} denotes the perturbed tensor-grid measure described in E.3 The spatial measures μxSPP\mu_{x}^{\mathrm{SPP}} and μxSPPCS\mu_{x}^{\mathrm{SPPCS}} are cluster-based sampling distributions. The velocity measures μvSPP\mu_{v}^{\mathrm{SPP}} and μvSPPCS\mu_{v}^{\mathrm{SPPCS}} are perturbed tensor-grid distributions on [−22,22]2[-2\sqrt{2},2\sqrt{2}]^{2} .
SystemComponentTermcoefftruecoeff_{\mathrm{true}}E(coeff^)±σ(coeff^)\mathbb{E}(\hat{coeff})\pm\sigma(\hat{coeff})
Kuramoto11θ\theta0.010.010.009\,999\,999\,999\,999\,9$\pm$2.087\,722\,866\,931\,242\,4\text{\times}{10}^{-16}
Phototaxis1111I0U∞eℓ(1)=−0.12I_{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}
v1v_{1}I0=1.0I_{0}=1.01.000\,000\,319\,453\,960\,8$\pm$3.538\,386\,027\,474\,474\text{\times}{10}^{-06}
v2v_{2}00-2.967\,810\,805\,171\,508\text{\times}{10}^{-07}$\pm$2.676\,347\,740\,402\,615\text{\times}{10}^{-06}
2211I0U∞eℓ(2)=0.12I_{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}
v1v_{1}00-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 Tr=10T_{r}=10 independent trials. Coefficients corresponding to non-present mechanisms are expected to be close to zero.
ParameterDefinition
Dynamical system
NNNumber of agents in the system
ddDimension of the state space of each agent
xi(t)∈Rdx_{i}(t)\in\mathbb{R}^{d}State of the ii -th agent at time tt
x(t)∈RdNx(t)\in\mathbb{R}^{dN}Full stacked state of the system at time tt
f:Rd→Rdf:\mathbb{R}^{d}\to\mathbb{R}^{d}Environmental force function (assumed identical for all agents)
Appendix
Table 26: Parameters and their definitions.
More in this topic

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

Explore similar work

CardsList

  1. Learning Interaction Kernels from Collective Steady States

    Sep 10, 2026Baoli Hao, Mauro Maggioni, Ming ZhongSystem Identification

    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

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

    Mar 16, 2026Giacomo Albi, Alessandro Alla, Elisa CalzolaSystem Identification

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

  3. Swarm-Inspired Generation of Collective Behaviors in Graph Dynamical Systems

    Jun 23, 2026Ji Chen, Song Chen, Chengzhang Gong +2Dynamical SystemsGraph Representation Learning

    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