Learning Multi-type heterogeneous interacting particle systems

We propose a framework for the joint inference of network topology, multi-type interaction kernels, and latent type assignments in heterogeneous interacting particle systems from multi-trajectory data. This learning task is a challenging non-convex mixed-integer optimization problem, which we address through a novel three-stage approach. First, we leverage shared structure across agent interactions to recover a low-rank embedding of the system parameters via matrix sensing. Second, we identify discrete interaction types by clustering within the learned embedding. Third, we recover the network weight matrix and kernel coefficients through matrix factorization and a post-processing refinement. We provide theoretical guarantees with estimation error bounds under a Restricted Isometry Property (RIP) assumption and establish conditions for the exact recovery of interaction types based on cluster separability. Numerical experiments on synthetic datasets, including heterogeneous predator-prey systems, demonstrate that our method yields an accurate reconstruction of the underlying dynamics and is robust to noise.

💡 Research Summary

The paper tackles the ambitious problem of jointly learning three intertwined components of heterogeneous interacting particle systems (IPS): the underlying network topology, multiple interaction kernels, and the latent type assignments that dictate which kernel governs each pairwise interaction. Formally, the dynamics are described by a stochastic differential equation where the drift term aggregates contributions from Q distinct kernels Φ_q, each selected by an integer‑valued type matrix κ∈{1,…,Q}^{N×N}. The unknown parameters are the adjacency matrix a∈{0,1}^{N×N} (row‑normalized), the kernel coefficient matrix c∈ℝ^{K×Q} (expansion of each Φ_q in a known basis {ψ_k}), and the type matrix κ. Observations consist of M independent trajectories sampled at L time points, providing noisy estimates of particle positions and velocities.

Directly minimizing the empirical mean‑squared error over (a,c,κ) leads to a non‑convex mixed‑integer optimization problem that is computationally intractable. The authors circumvent this difficulty by re‑parameterizing the interaction structure of each agent i through a low‑rank matrix Z_i = Diag(p_i)·cᵀ, where p_i∈{0,1}^Q encodes the type of agent i. This representation reveals a rank‑Q structure across the stacked matrices {Z_i}, enabling a reduction to a low‑rank matrix sensing problem.

The proposed solution proceeds in three stages:

1. Matrix Sensing (Stage 1). Using the Alternating Least Squares (ALS) algorithm, the method estimates all Z_i simultaneously from the trajectory data. Each Z_i satisfies a linear relation B_i·Z_i ≈ \dot X_i, where B_i are sensing matrices constructed from observed positions and the basis functions. Under a Restricted Isometry Property (RIP) on the sensing operator, the authors prove that the ALS estimator achieves near‑optimal Frobenius‑norm error scaling as O(σ√((NQ+K)/M)), where σ denotes the noise level.

2. Clustering (Stage 2). The non‑zero rows of the estimated Z_i are normalized and clustered via K‑means. The geometric separation between rows belonging to different interaction types is quantified by angular distances. The paper establishes that if the minimal inter‑cluster angle exceeds a threshold that depends on the estimation error and the minimal row norm, K‑means exactly recovers the true type matrix κ, even when the number of types Q is unknown. An automatic model‑selection criterion based on intra‑ versus inter‑cluster angle ratios is provided.

3. Matrix Decomposition and Post‑Processing (Stage 3). With the cluster assignments fixed, the stacked Z_i are factorized to retrieve the adjacency matrix a and the kernel coefficient matrix c. This is done by solving a linear system under the row‑normalization constraint on a, a step that is computationally cheap. A final refinement pass of ALS is applied to improve robustness against measurement noise.

Theoretical contributions include (i) RIP‑based error bounds for the low‑rank sensing stage, (ii) exact recovery conditions for κ expressed via angular separation and minimal row norm, and (iii) a proof that the three‑stage pipeline yields consistent estimates of (a,c,κ) as the number of trajectories grows.

Empirical validation is performed on synthetic datasets, notably heterogeneous predator‑prey models with three distinct interaction kernels. Experiments vary the number of agents (N=50–200), kernel dimensionality (K=20–50), number of trajectories (M=20–200), and noise level (σ up to 0.1). Results demonstrate that the method recovers the network weights and kernel coefficients with Frobenius errors below 10⁻³ in noise‑free settings and below 10⁻² with moderate noise. Type assignment accuracy exceeds 98 % across all tested regimes, and a sharp phase transition is observed when the inter‑cluster angular separation falls below the theoretical threshold. The approach also works for complete graphs, where a dedicated variant skips the zero‑entry identification step.

In the context of related work, the paper distinguishes itself by providing a unified framework that simultaneously addresses graph inference, multi‑kernel learning, and latent type discovery—tasks that have traditionally been tackled separately and under restrictive assumptions (e.g., known types or dense graphs). By exploiting low‑rank structure and geometric clustering, the authors avoid the combinatorial explosion inherent in mixed‑integer formulations while still delivering strong statistical guarantees.

The authors conclude that their three‑stage algorithm offers a scalable, noise‑robust solution for learning heterogeneous IPSs and suggest future directions such as online learning, extensions to non‑Gaussian noise, application to real biological or ecological data, and theoretical extensions to hierarchical or time‑varying type matrices.