Learning interacting particle systems from unlabeled data

Viska Wei; Fei Lu

stat.MLApr 2, 2026

Learning interacting particle systems from unlabeled data

Authors: Viska Wei, Fei Lu

Organizations: Department of Applied Mathematics and Statistics and Department of Physics and Astronomy, Johns Hopkins University, Baltimore, MD 21218, USA. · Department of Mathematics, Johns Hopkins University, Baltimore, MD 21218, USA.

Abstract

Learning the potentials of interacting particle systems is a fundamental task across various scientific disciplines. A major challenge is that unlabeled data collected at discrete time points lack trajectory information due to limitations in data collection methods or privacy constraints. We address this challenge by introducing a trajectory-free self-test loss function that leverages the weak-form stochastic evolution equation of the empirical distribution. The loss function is quadratic in potentials, supporting parametric and nonparametric regression algorithms for robust estimation that scale to large, high-dimensional systems with big data. Systematic numerical tests show that our method outperforms baseline methods that regress on trajectories recovered via label matching, tolerating large observation time steps. We establish the convergence of parametric estimators as the sample size increases, providing a theoretical foundation for the proposed approach.

Figures & tables

\mathcal{L}(\theta)=\frac{1}{B}\sum_{\ell,m\in\mathcal{I}_{batch}}\left[\frac{1}{2}J_{diss}^{(\ell,m)}\Delta t-\frac{\sigma^{2}}{2}J_{diff}^{(\ell,m)}\Delta t+\delta E_{f}^{(\ell,m)}\right],

Algorithm 2 Neural Network Estimation via Trajectory-Free Loss

Figure 1: Workflow of both estimation algorithms using the self-test loss function. Left: Least squares regression expands

V

and

\Phi

in prescribed basis functions (or RBFs), and solves the minimizer via least squares with regularization. Right: Neural network regression parameterizes

V_{\theta}

and

\Phi_{\theta}

by neural networks, computes derivatives via automatic differentiation, and minimizes the loss using stochastic gradient descent.

Method	Labeled	#Params	Complexity
Labeled MLE	✓	$K$	$O\bigl(MLd\,NK(N{+}K)+K^{3}\bigr)$
Self-Test LSE	—	$K$	$O\bigl(MLd\,NK(N{+}K)+K^{3}\bigr)$
Sinkhorn MLE	—	$K$	$O\bigl(MLd[NK(N{+}K)+T_{S}dN^{2}]{+}K^{3}\bigr)$
Self-Test NN	—	$P$	$O\bigl(EMLd(NP_{V}+N^{2}P_{\Phi})\bigr)$

Table 1: Computational complexity of all estimation methods. Here,

M

is the number of trajectories,

L

is the number of time snapshots,

N

is the number of particles,

d

is the dimension,

K=K_{V}{+}K_{\Phi}

is the basis size,

P_{V},P_{\Phi}

are the network parameters,

E

is the number of epochs, and

T_{S}

is the number of Sinkhorn iterations.

Figure 2:

M

-scaling on the reference model under Riemann-sum (left pair) and trapezoidal (right pair) time integrations. Both with data generated with

\delta t=10^{-4}

and integrated with the various

\Delta t

values. The Riemann-sum has error bound

O(\Delta t+M^{-1/2})

, so the four

\Delta t

values track the

O(M^{-1/2})

rate (green line) until they saturate at an

O(\Delta t)

floor (left pair), while the trapezoidal rule has error bound

O((\Delta t)^{2}+M^{-1/2})

, so all four

\Delta t

values track the

O(M^{-1/2})

rate without saturation (right pair). The mean and standard deviation are computed over 10 trials per point; detailed numbers in Table 6 .

Figure 3: Convergence for discrete-time model (i.e.,

\delta t=\Delta t

, zero gap). The intrinsic

O(\Delta t)

bias of the discrete-time model (see Section 2.3 ) dominates the error when

\Delta t

is large, and the convergence in

M

is only visible at

\Delta t\leq 10^{-3}

(and

\Delta t=10^{-4}

for

\Phi

) where the bias is small enough to allow the

O(M^{-1/2})

statistical error to dominate. See also in Table 6 .

Figure 4: Condition numbers scaling with

N

for the normal matrix and its diagonal blocks (

\kappa_{*}

\kappa_{VV}

\kappa_{\Phi\Phi}

). Slopes (italic numbers at

N{=}100

) are log-log regression rates from Table 7 . Dashed/dotted green lines show the

O(N)

theoretical rates. Left:

\kappa_{*}

increases in

N

, with a rate near

O(N)

when

d=10

. Center:

\kappa_{VV}

N

-independent. Right:

\kappa_{\Phi\Phi}

has a rate near

O(N)

when

d=10

due to distance concentration.

$d{=}2$ , $N{=}10$ , $M{=}2{,}000$ , $T{=}1$ , $\sigma{=}1$ , $\delta t{=}10^{-4}$
	Labeled MLE		Self-Test LSE		Sinkhorn MLE		Self-Test NN
$\Delta t$	$\nabla V$	$\nabla\Phi$	$\nabla V$	$\nabla\Phi$	$\nabla V$	$\nabla\Phi$	$\nabla V$	$\nabla\Phi$
$10^{-4}$	$0.80{\pm}0.55$	$0.35{\pm}0.19$	$1.35{\pm}0.79$	$1.24{\pm}0.29$	$1.34{\pm}0.80$	$0.82{\pm}0.18$	—	—
$10^{-3}$	$1.41{\pm}0.45$	$5.18{\pm}0.16$	$0.67{\pm}0.37$	$0.74{\pm}0.15$	$3.13{\pm}0.58$	$7.17{\pm}0.15$	$0.91{\pm}0.17$	$1.74{\pm}0.13$
$10^{-2}$	$10.44{\pm}0.41$	$37.49{\pm}0.28$	$0.80{\pm}0.43$	$1.10{\pm}0.17$	$21.64{\pm}0.34$	$47.38{\pm}0.25$	$1.25{\pm}0.33$	$2.47{\pm}0.22$
$10^{-1}$	$11.30{\pm}0.58$	$89.38{\pm}0.18$	$6.84{\pm}0.58$	$6.93{\pm}0.72$	$45.37{\pm}0.56$	$95.78{\pm}0.23$	$3.15{\pm}0.68$	$5.80{\pm}0.42$

Table 2: Method comparison across observation intervals

\Delta t

for the reference model (relative gradient errors, %). Oracle: true basis, mean

\pm

std over 10 trials. NN: MLP

[64{,}64{,}64]

, mean

\pm

std over 10 trials. Green = best per

\Delta t

Figure 5: Non-radial potential recovery (

d{=}2

). Percentages are relative

L^{2}(\rho)

errors.

Model	$\Phi$ type	Violation / stress point
Smoothness test	Smoothed indicator $\mathbf{1}_{[a,b]}^{\epsilon}$	$C^{4}_{b}$ smoothness (Assumption 4.1 (i))
Conditioning test	$\gamma/(r+1)$ inverse	Coercivity — near-singular Gram matrix
Singularity test	$r^{-12}-r^{-6}$ singularity	$C^{4}_{b}$ + gradient blowup
Smooth control	$D(1-e^{-a(r-r_{0})})^{2}$	None (well-behaved baseline)

Table 3: Boundary test models and the assumptions they violate or stress.

Figure 6: Typical estimators of

\Phi^{\prime}(r)

for the four boundary test models (

d{=}2

M{=}2{,}000

\sigma{=}1

\delta t{=}10^{-4}

\Delta t{=}10^{-2}

). Gray shading indicates the approximate exploration measure

\rho_{\Phi}

. Smoothness: Smoothed indicator with sharp transitions at

r=0.5

and

r=1

. Conditioning: Slow-decaying

1/(r{+}1)

interaction that challenges conditioning. Singularity: Lennard-Jones

r^{-12}{-}r^{-6}

singularity; gradient diverges as

r\to 0

. Smooth control: Morse potential with equilibrium at

r_{0}=0.8

$d{=}2$ , $N{=}10$ , $M_{\text{total}}{=}20{,}000$ , $M{=}2{,}000$ , $T{=}1$ , $\sigma{=}1$ , $\delta t{=}10^{-4}$
		Labeled MLE		Self-Test		Sinkhorn		Self-Test NN
Model	$\Delta t$	$\nabla V$	$\nabla\Phi$	$\nabla V$	$\nabla\Phi$	$\nabla V$	$\nabla\Phi$	$\nabla V$	$\nabla\Phi$
Smoothness	$10^{-3}$	$6.70{\pm}0.62$	$12.98{\pm}0.19$	${\color[rgb]{0,0.4,0}\bm{1.07{\pm}0.52}}$	${\color[rgb]{0,0.4,0}\bm{2.01{\pm}0.20}}$	$6.30{\pm}0.53$	$14.42{\pm}0.19$	$1.76{\pm}0.54$	$3.32{\pm}0.29$
	$10^{-2}$	$21.05{\pm}0.55$	$52.63{\pm}0.18$	${\color[rgb]{0,0.4,0}\bm{0.85{\pm}0.47}}$	${\color[rgb]{0,0.4,0}\bm{2.55{\pm}0.29}}$	$18.23{\pm}0.45$	$59.31{\pm}0.16$	$1.41{\pm}0.44$	$4.60{\pm}0.86$
	$10^{-1}$	$3.57{\pm}0.32$	$92.15{\pm}0.22$	$9.07{\pm}1.33$	$10.01{\pm}0.74$	$31.32{\pm}0.20$	$95.89{\pm}0.12$	${\color[rgb]{0,0.4,0}\bm{3.64{\pm}1.42}}$	${\color[rgb]{0,0.4,0}\bm{8.77{\pm}0.88}}$
Conditioning	$10^{-3}$	$1.15{\pm}0.59$	$11.97{\pm}9.08$	${\color[rgb]{0,0.4,0}\bm{1.46{\pm}0.73}}$	$31.02{\pm}9.63$	$13.58{\pm}1.06$	$271{\pm}11$	$3.83{\pm}0.38$	${\color[rgb]{0,0.4,0}\bm{26.7{\pm}8.8}}$

Table 4: Relative errors (%) as

\Delta t

varies in boundary tests. Parametric estimators: mean

\pm

std in 10 trials (

M{=}2{,}000

subsampled from

M_{\text{total}}{=}20{,}000

). Green = best among unlabeled methods.

$d{=}2$ , $\Delta t{=}10^{-2}$ , $M{=}2000$ , oracle basis
Model	$N{=}5$	$N{=}10$	$N{=}20$	$N{=}50$	$N{=}100$	Rate in $N$
Reference	37	53	80	123	155	$0.5$
Smoothness	42	46	88	208	396	$0.8$
Conditioning	81K	146K	262K	564K	962K	$0.8$
Singularity	10.8K	12.1K	11.8K	12.1K	11.8K	$0.0$ (constant)
Smooth ctrl	7.9K	14.2K	25.4K	53.5K	89.5K	$0.8$

Table 5: Condition number

\kappa(\mathbf{A}_{*})

across all models. Rate in

N

: log-log regression

\partial\!\log\kappa/\partial\!\log N

; Proposition B.3 predicts slope

{\leq}\,1

Appendix figures & tables3 assets

Supplementary material from the paper’s appendix.

Appendix

$d{=}2$ , $N{=}10$ , $T{=}1$ , $\sigma{=}1$ , oracle basis, $\delta t{=}10^{-4}$
	$\Delta t{=}10^{-4}$		$10^{-3}$		$10^{-2}$		$10^{-1}$
$M$	$\nabla V$	$\nabla\Phi$	$\nabla V$	$\nabla\Phi$	$\nabla V$	$\nabla\Phi$	$\nabla V$	$\nabla\Phi$
20	$6.6{\pm}2.2$	$2.8{\pm}1.5$	$6.4{\pm}2.9$	$3.1{\pm}1.0$	$5.9{\pm}3.2$	$3.5{\pm}1.9$	$12.3{\pm}6.7$	$5.7{\pm}3.4$
50	$4.1{\pm}1.9$	$2.2{\pm}1.4$	$4.2{\pm}1.9$	$1.7{\pm}0.8$	$5.1{\pm}1.7$	$2.2{\pm}0.7$	$5.5{\pm}2.4$	$3.8{\pm}2.9$
200	$2.2{\pm}1.1$	$1.2{\pm}0.6$	$2.2{\pm}1.2$	$1.2{\pm}0.6$	$2.1{\pm}0.9$	$1.4{\pm}0.6$	$4.2{\pm}2.0$	$2.4{\pm}0.8$
1000	$1.4{\pm}0.6$	$0.7{\pm}0.3$	$0.8{\pm}0.5$	$0.9{\pm}0.2$	$1.0{\pm}0.5$	$0.7{\pm}0.2$	$2.0{\pm}1.0$	$1.0{\pm}0.6$

Appendix

Table 6:

M

-scaling of self-test LSE: relative

L^{2}(\rho)

errors (%) for the reference model (

d{=}2

N{=}10

, 10 trials per cell). (a)-(b) : Riemann sum and Trapezoidal quadrature for time integration. The observation time step

\Delta t\in\{10^{-3},10^{-2},10^{-1}\}

, and the data is generated with time step fixed at

\delta t=10^{-4}

. (c) : Zero-gap regime (

\delta t=\Delta t

). The slope is estimated by log-log regression of error vs.

M

. Those slopes matching the theoretical

O(M^{-1/2})

rate are highlighted in green . See Figure 2 for the corresponding plots.

$\Delta t{=}10^{-2}$ , $M{=}2000$ , $T{=}1$ , $\sigma{=}1$ , oracle basis
	$d=2$			$d=5$			$d=10$
$N$	$\kappa_{*}$	$\kappa_{VV}$	$\kappa_{\Phi\Phi}$	$\kappa_{*}$	$\kappa_{VV}$	$\kappa_{\Phi\Phi}$	$\kappa_{*}$	$\kappa_{VV}$	$\kappa_{\Phi\Phi}$
5	$37$	$19.7$	$12.8$	$42$	$35.4$	$6.2$	$170$	$156$	$3.0$
10	$53$	$20.2$	$6.8$	$44$	$32.4$	$4.7$	$232$	$166$	$6.9$
20	$80$	$20.5$	$3.7$	$50$	$31.4$	$4.5$	$533$	$171$	$14$
50	$123$	$20.8$	$1.9$	$69$	$30.7$	$6.1$	$1629$	$177$	$31$

Appendix

Table 7: Condition numbers of the normal matrix and its diagonal blocks in self-test LSE for the reference model. Rate in

N

is estimated via log-log regression

\partial\log\kappa/\partial\log N

; Proposition B.3 predicts that

\kappa_{*}=O(N)

, i.e. a rate

\,1

Figure 7: Typical estimators of

V^{\prime}

and

\Phi^{\prime}

for the radial potentials when

d=2

. Top left: Smoothness test. Top right: Conditioning test. Bottom left: Singularity test (LJ). Bottom right: Smooth control (Morse). Gray shading indicates the empirical density of observed distances.

System Identification

System identification infers mathematical models and parameters of dynamical systems from input-output, trajectory, sensor, or video data.

See all 157 papers

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
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
Learning Latent Energy-Based Models via Interacting Particle Langevin Dynamics
Oct 14, 2025Joanna Marks, Tim Y. J. Wang, O. Deniz AkyildizParameter Estimation Latent Variable Models
Department of Mathematics, Imperial College London

Sep 10, 2026J/K move · Enter open · S save

Learning Interaction Kernels from Collective Steady States

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

Open paper AskPDF arXiv

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.

Abstract

System Identification

Explore similar work

Learning Interaction Kernels from Collective Steady States

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

Learning Latent Energy-Based Models via Interacting Particle Langevin Dynamics