1 Introduction↩︎

Given high-dimensional microscopic observations of a physical system, accurately modeling its macroscopic dynamics is of broad interest across disciplines [1]. A central challenge is to construct a closed macroscopic model: the evolution of macroscopic quantities depends on microscopic degrees of freedom, and we therefore seek closure variables that encode the microscopic information. Recent data-driven approaches leverage neural networks to automatically discover low-dimensional descriptions that serve as closure variables. These methods have achieved notable success for many scientific applications, including learning reduced thermodynamic coordinates and predicting stretching length dynamics of polymers [2], discovering slow collective variables and reduced kinetic models for molecular conformational dynamics [3], [4], and learning macroscopic internal variables for history-dependent heterogeneous materials [5].

Most existing approaches learn closure variables using an autoencoder trained to minimize the pointwise reconstruction loss. The encoder and decoder are typically parameterized by the Multilayer Perceptron (MLP) [2], [6] or Convolutional Neural Network (CNN) [7], assuming that the microscopic degrees of freedom admit a fixed ordering to be represented as vectors or tensors (Fig. 1 (a)). However, many physical systems have no intrinsic ordering, such as a set of interacting particles. The desired closure variable must be permutation-invariant with respect to the input (Fig. 1 (b)). Applying existing methods, which require a canonical ordering on the input, can be prohibitively difficult in practice.

This gap arises because the autoencoders parameterized by MLP or CNN are not symmetric to permutation. The input and output are both indexed by an ordering, and the usual pointwise reconstruction loss fixes an index-by-index correspondence. Consequently, the learned latent vector encodes the input information with a specific ordering, and can change when the same microstate is presented under a different permutation. One may consider adopting models for sets (e.g., DeepSet [8] or Set Transformer [9]) as the encoder to produce a permutation-invariant latent variable. Yet, a fundamental difficulty remains on the decoding side: without a specified ordering, there is no natural mechanism to reconstruct the microstate from a single low-dimensional embedding using an index-wise loss. For instance, consider a microstate \(\{{\boldsymbol{x}}^i\}_{i=1}^N\) with \(N=3\) particles. If a decoder is trained with a pointwise MSE loss to output an ordered list \((\hat{{\boldsymbol{x}}}^1, \ldots, \hat{{\boldsymbol{x}}}^N)\), the loss will penalize \((\hat{{\boldsymbol{x}}}^1, \hat{{\boldsymbol{x}}}^2, \hat{{\boldsymbol{x}}}^3)\) and \((\hat{{\boldsymbol{x}}}^3, \hat{{\boldsymbol{x}}}^2, \hat{{\boldsymbol{x}}}^1)\) differently, even though they correspond to the same physical configuration. Crucially, the number of such physically equivalent permutations scales as \(N!\). Therefore, training often relies on explicit matching [10], [11] or permutation-invariant losses [12], [13]; both add substantial computational cost [14] and may introduce optimization instability [15], [16].

To overcome this difficulty, we propose an alternative reconstruction objective for closure modeling. Rather than reconstructing the microscopic degrees of freedom directly, we reconstruct their distribution. Concretely, we induce a mass distribution concentrated at the observed points, and train the decoder to reconstruct this density from the latent representation. This distributional reconstruction eliminates the need to align individual points or impose an ordering, thereby enabling learning of permutation-invariant latent variables. The learned latent variable summarizes the input distribution and serves as closure variables for modeling macroscopic dynamics. Our contributions can be summarized as follows.

• We propose a novel strategy for closure modeling that trains an autoencoder to reconstruct the distributional information instead of pointwise reconstruction.

• We demonstrate the effectiveness and robustness of the proposed approach across multiple microscopic settings, including learning the energy dynamics of interacting particle systems, the mixing dynamics of Lennard–Jones fluid and predicting the stretching dynamics of polymers from videos.

2 Related Work↩︎

Learning Reduced Dynamics Reduced-order models (ROMs) and closure modeling approximate high-dimensional dynamical systems using low-dimensional latent variables and learn latent dynamics, often with the goal of improving computational efficiency [17], [18]. ROMs typically aim to approximate the full-order state by evolving latent variables and decoding them back to the full-order space [19]–[21]. In contrast, closure modeling targets accurately modeling the dynamics of selected macroscopic observables and, in principle, does not have to reconstruct the full-order state [2], [22], [23]. Nevertheless, many data-driven closure approaches adopt a reconstruction objective to learn closure coordinates. The reconstruction loss offers a convenient objective that enforces the latent variables to encode microstate information and regularizes the latent variable [2], [24]–[26]. Accordingly, existing closure modeling methods often consist of two components: (i) an autoencoder for discovery of closure variables trained with a reconstruction loss and (ii) a reduced dynamical model for predicting the latent dynamics [2], [24]–[26]. Our method belongs to closure modeling, as we focus on modeling the dynamics of macroscopic quantities and leverage a reconstruction-based objective to learn closure variables. Prior work has mainly focused on dynamical systems whose state variables have a fixed indexing, such as discretized partial differential equation  (PDE) solutions on grids [21], [24], [27]. However, these methods are not directly applicable when the microstate has no canonical ordering. Our method addresses this challenge through a permutation-invariant architecture and a newly designed permutation-invariant training objective.

Learning Permutation-invariant Representations for Sets General machine learning models have been proposed to learn permutation-invariant representations of sets. Classic examples include pooling-based architectures such as DeepSet [8] and PointNet [28], as well as attention-based models such as the Set Transformer [9]. In our model, we adopt their architecture as the encoder. Unlike many previous works that use set encoders as standalone modules for supervised prediction [8], [9], [28]–[30], we instead learn permutation-invariant closure variables by pairing the set encoder with a decoder and optimizing an unsupervised reconstruction objective. Designing the proper decoder is non-trivial because the input has no ordering. A naive decoder that outputs individual elements raises the question of how to align reconstruction with the input. To do this, one has to explicitly handle element ordering, commonly via permutation-invariant set-matching losses (as commonly done in point-cloud reconstruction; see next paragraph). We avoid this complication by reconstructing the input distribution, thereby removing any dependence on the ordering.

Point Clouds Reconstruction Various autoencoder architectures have been proposed for point cloud reconstruction [12], [13], [31]–[35]. These methods utilize permutation-invariant architectures as the encoder (as described above) and design the decoder that outputs an explicit set of points intended to match the input samples. To make the reconstruction invariant to how the points are indexed, they either introduce a combinatorial matching step [10], [11] or optimize permutation-invariant distances (e.g., Chamfer discrepancy, Earth Mover’s distance, and sliced Wasserstein distance) as reconstruction loss [12], [13]. More recently, [36] proposed a soft-matching reconstruction objective by representing both the input and output as point-centered Gaussian mixtures and measuring reconstruction accuracy via pairwise Gaussian overlaps. While our approach also learn permutation-invariant latent variables from a reconstruction objective, it differs in what is reconstructed. Rather than decoding an explicit point set and performing pointwise (hard or soft) alignment, we decode a density that describes the microstate. Concretely, we induce a mass distribution concentrated at the observed points and train the decoder to reconstruct this density from the latent variable. This formulation is inherently permutation-invariant and avoids the matching issue for unordered points.

3 Method↩︎

3.1 Problem Statement↩︎

Let \(X_t = \{{\boldsymbol{x}}^1_t, {\boldsymbol{x}}^2_t, \ldots, {\boldsymbol{x}}^n_t\} \in \mathcal{X}\) be the state of \(n\) unordered particles at time \(t\), where \({\boldsymbol{x}}^i_t \in \mathbb{R}^d\) is the state of particle with index \(i\). In many real-world systems, the number of particles \(n\) is very large, and the microscopic system is therefore high-dimensional. A deterministic function \(\bar \boldsymbol{\varphi}\) is given beforehand to extract the macroscopic quantities of interest from the microscopic state, \(\bar \boldsymbol{\varphi}: \mathcal{X} \rightarrow \bar{\mathcal{Z}}\) where \(\bar{\mathcal{Z}}\) is a space containing low-dimensional vectors. We restrict our discussion to intensive macroscopic observables such as the average energy of the interacting particle system, so that their magnitude does not scale with the system size \(n\). The goal of our method is to learn a permutation-invariant macroscopic dynamical model that can be applied to microscopic systems with varying numbers of particles.

To this end, an autoencoder is used to learn a low-dimensional representation that serves as the closure variable, together with a latent dynamical model for predicting macroscopic dynamics. The architecture is illustrated in Fig. 2. The encoder, \(\hat{\boldsymbol{\varphi}}: \mathcal{X} \rightarrow \hat{\mathcal{Z}}\), maps the microstate \(X\in \mathcal{X}\) to a low-dimensional latent vector \(\hat{{\boldsymbol{z}}}=\hat{\boldsymbol{\varphi}}\!\left(X\right)\). The latent \(\hat{{\boldsymbol{z}}}\) should be invariant to the ordering of \(X = \{{\boldsymbol{x}}^1, {\boldsymbol{x}}^2, \ldots, {\boldsymbol{x}}^n\}\), i.e., \[\label{eq:permutation-invariant} \hat{\boldsymbol{\varphi}}(\sigma X) = \hat{\boldsymbol{\varphi}}(X), \quad \forall\, \sigma \in S_n,\tag{1}\] where \(S_n\) denotes the symmetric group on \(n\) elements and \(\sigma X := \{{\boldsymbol{x}}^{\sigma(1)}, {\boldsymbol{x}}^{\sigma(2)}, \ldots, {\boldsymbol{x}}^{\sigma(n)}\}\) is the permuted particle set. In parallel, we extract the macroscopic quantity of interest \(\bar{{\boldsymbol{z}}}_{t}=\bar{\boldsymbol{\varphi}}\!\left(X_{t}\right)\), and form the closed macroscopic state \({\boldsymbol{z}}_{t}=[\bar{{\boldsymbol{z}}}_{t},\hat{{\boldsymbol{z}}}_{t}] \in {\mathcal{Z}}\). We denote the dimension of \(\hat{{\boldsymbol{z}}}\in \hat{\mathcal{Z}}\), \({\boldsymbol{z}}\in {\mathcal{Z}}\) by \(\hat{z}_\mathrm{dim}\), \({z}_\mathrm{dim}\). Next, we describe the closure modeling to learn \(\hat{{\boldsymbol{z}}}\) and the macroscopic dynamics modeling that predicts the evolution of \({\boldsymbol{z}}_t\).

3.2 Closure Modeling↩︎

As discussed in Sec. 1, a common approach is to train an autoencoder using a per-sample reconstruction loss and to use the learned latent representation as the closure variable. That is, they represent \(\{{\boldsymbol{x}}^1_t, {\boldsymbol{x}}^2_t, \ldots, {\boldsymbol{x}}^n_t\}\) by a vector \({\boldsymbol{y}}\in \mathbb{R}^{nd}\) according to an indexing and minimize an element-wise reconstruction error of the form: \[\textstyle \mathbb{E}_{{\boldsymbol{y}}}\left[ \sum_{i=1, 2, \dots, nd}\, l\Big( {\boldsymbol{y}}[i], \boldsymbol{\psi}(\hat{\boldsymbol{\varphi}}({\boldsymbol{y}}))[i] \Big) \right]~,\] where \([i]\) denotes the i-th entry in a vector and \(l\) is the loss such as squared error. However, when the input has no canonical ordering, the index-wise loss above is not permutation-invariant and can penalize reconstructions that are correct up to permutation. We resolve this by shifting the reconstruction objective from input particles themselves to a distributional representation of the input set.

In our implementation (Fig. [fig:method95overview]), we parameterize the encoder \(\hat{\boldsymbol{\varphi}}\) using DeepSet to obtain a permutation-invariant representation \(\hat{{\boldsymbol{z}}}\) (Eq. 1 ). While we use DeepSet in this work, other permutation-invariant set encoders such as set transformer [9] can be used as well. For each \(X\in\mathcal{X}\), we induce a continuous density \(q_X(\mathbf{x})\) over the input space \(\mathcal{X}\) that corresponds to a mass distribution centered at the observation points, i.e., \[\label{eq:target95density95q} q_{X}(\mathbf{x}) = \frac{1}{|X|} \sum_{{\boldsymbol{x}}^i \in X} \delta_{\epsilon}(\mathbf{x}- {\boldsymbol{x}}^i),\tag{2}\] where \(\delta_\epsilon(\cdot)\) is an isotropic Gaussian kernel with variance \(\epsilon^2\). The \(q_{X}(\mathbf{x})\) serves as the reconstruction target. The parameter \(\epsilon\) acts as a bandwidth (smoothing scale) and therefore controls the resolution at which the set is represented: a smaller \(\epsilon\) yields sharper, more localized modes around individual particles, while a larger \(\epsilon\) produces a smoother density that discards more local details. More discussions on \(\epsilon\) are in Appendix 7.3 and the numerical results are in Appendix 8.2.

Given the learned \(\hat{{\boldsymbol{z}}}=\hat{\boldsymbol{\varphi}}(X)\) and the target density \(q_{X}(\mathbf{x})\), the decoder \(\boldsymbol{\psi}\) generates a conditional density \(p(\mathbf{x}| \hat{{\boldsymbol{z}}})\) via \[p_{\theta}(\mathbf{x}| \hat{{\boldsymbol{z}}}) = \boldsymbol{\psi}(\mathbf{x}; \hat{\boldsymbol{\varphi}}(X))~,\] where \({\theta}\) denotes all learnable parameters in \(\hat{\boldsymbol{\varphi}}\) and \(\boldsymbol{\psi}\). We implement \(p(\mathbf{x}| \hat{{\boldsymbol{z}}})\) with a conditional normalizing flow, though any conditional density model can be used [37], [38]. We refer readers to Appendix 7.2 for further implementation details. In this work, we use the KL divergence to measure the discrepancy between \(q_{X}(\mathbf{x})\) and \(p_{\theta}(\mathbf{x}| \hat{{\boldsymbol{z}}})\), although many alternative distributional losses could be used [39]: \[\label{eq:reconstruction95loss} \mathcal{L}_{\mathrm{rec}} = \mathbb{E}_{X}\left[ \mathrm{KL}\! \left( q_{X}(\mathbf{x})\;\|\;p_{\theta}(\mathbf{x}| \hat{{\boldsymbol{z}}})\right) \right]~.\tag{3}\] We estimate the KL divergence with Monte Carlo samples drawn from \(q_{X}(\cdot)\). Noting that \(q_{X}(\cdot)\) is a mixture of Gaussians sharing the same variance and weight, we can efficiently draw samples in parallel by firstly sampling the mixture component uniformly and then sampling from the local Gaussian.

This designed autoencoder brings several advantages for closure modeling. First, the learned \(\hat{{\boldsymbol{z}}}\) is permutation-invariant to input indexing. Second, the encoder admits variable-size inputs, enabling \(\hat{{\boldsymbol{z}}}\) to be computed for different particle counts. Third, the training objective avoids explicit set-to-set matching, which could be computationally expensive. Specifically, in our pipeline, the DeepSet encoder scales linearly with the input size, i.e., \(\mathcal{O}(n)\) with respect to the number of input particles \(n\), because the encoder processes each particle feature independently before aggregation. Given the latent variable generated by the encoder, the computational complexity of the density-based decoder is \(\mathcal{O}(1)\) with respect to the input size. This is because the reconstruction objective is evaluated using Monte Carlo samples from the target density, and the sampling procedure is independent of the number of particles under our construction. Therefore, unlike pointwise reconstruction losses that may require expensive pointwise matching, the computational cost of the proposed distribution-reconstruction objective depends on the number of Monte Carlo samples, the input feature dimension, and the cost of normalizing flow evaluations, rather than the number of input particles.

3.3 Macroscopic Dynamics Modeling↩︎

We model the dynamics of the augmented macroscopic state \({\boldsymbol{z}}_t=[\bar{{\boldsymbol{z}}}_t,\hat{{\boldsymbol{z}}}_t]\) by \[\label{eq:ode95or95sde} d{\boldsymbol{z}}_t = \boldsymbol{g}({\boldsymbol{z}}_t) \,dt + \boldsymbol{\Sigma}({\boldsymbol{z}}_t )\,d\boldsymbol{W}_t,\tag{4}\] where \(\boldsymbol{g}({\boldsymbol{z}}_t)\) is the drift term, \(\boldsymbol{\Sigma}({\boldsymbol{z}}_t)\) is the diffusion term, and \(\boldsymbol{W}_t\) denotes standard Brownian motion. Since the input consists of discrete frames, we learn a parameterized dynamics model by maximizing the likelihood of one-step transitions implied by an Euler–Maruyama discretization with step size \(\Delta t\): \[p_{\boldsymbol{g}, \boldsymbol{\Sigma}}({\boldsymbol{z}}_{t+1}| {\boldsymbol{z}}_{t}) = \mathcal{N}\!\Big({\boldsymbol{z}}_{t} + \boldsymbol{g}({\boldsymbol{z}}_{t})\,\Delta t,\;\; \Delta t\,\boldsymbol{\Sigma}({\boldsymbol{z}}_{t})\boldsymbol{\Sigma}({\boldsymbol{z}}_{t})^\top \Big)~,\] Where \(p_{\boldsymbol{g}, \boldsymbol{\Sigma}}\) indicates the likelihood depends on the parameters in \(\boldsymbol{g}\) and \(\boldsymbol{\Sigma}\). We therefore train \((\boldsymbol{g},\boldsymbol{\Sigma})\) by minimizing the negative conditional log-likelihood \[\label{eq:dyn95nll} \mathcal{L}_{\mathrm{dyn}} = \mathbb{E}_{{\boldsymbol{z}}_t, {\boldsymbol{z}}_{t+1}} [ - \log p_{\boldsymbol{g}, \boldsymbol{\Sigma}}\!\left({\boldsymbol{z}}_{t+1}| {\boldsymbol{z}}_{t}\right)]~.\tag{5}\] In the deterministic case, we set \(\boldsymbol{\Sigma}({\boldsymbol{z}})\equiv 0\) and minimizes a standard one-step squared-error loss.

Learning the dynamics over \({\boldsymbol{z}}\) naturally yields the macroscopic dynamics of \(\bar{{\boldsymbol{z}}}\), which is our quantity of interest. In this work, we parameterize the drift \(\boldsymbol{g}\) and diffusion \(\boldsymbol{\Sigma}\) if it is present with MLPs, but any dynamical model can be used as well, such as the OnsagerNet [40]. Once trained, the encoder is used only to construct the initial state \({\boldsymbol{z}}_0\). Then the learned dynamics model generates the trajectory autoregressively, i.e., it predicts future macroscopic states \({\boldsymbol{z}}_1, {\boldsymbol{z}}_2, \ldots, {\boldsymbol{z}}_t\) given only an initial state \({\boldsymbol{z}}_0\).

Overall Objective We train the full framework with the total loss \[\mathcal{L} = \mathcal{L}_{\mathrm{rec}} + \lambda_{\mathrm{dyn}}\mathcal{L}_{\mathrm{dyn}},\] where \(\lambda_{\mathrm{dyn}}\) balances distribution reconstruction and dynamical consistency. Following existing approaches [24], [41], we use a two-stage training procedure to learn the latent states and their corresponding latent dynamical models separately. The pseudocode of training and inference procedures is provided in Appendix 6.1.

4 Experiments↩︎

We compare our method against prior closure modeling approaches to assess its effectiveness. As noted, most existing works assume ordered microstates and employ an autoencoder with MLPs for the encoder and decoder, trained by minimizing a pointwise reconstruction loss [6], [24], [42]. In contrast, we consider permutation-invariant closure modeling from unordered microstates. To enable a fair comparison, we adopt the common strategy of sampling random input permutations during training [43]. During training, we randomly permute the order of the input anchors within each mini-batch before feeding them to the flattened MLP autoencoder. Combining this permutation-augmentation strategy with a standard MLP autoencoder trained using pointwise MSE reconstruction loss gives our first baseline, denoted as (AE-Aug). Besides, we also replace the MLP encoder with a symmetric architecture to learn permutation-invariant closure variables. Concretely, we use a DeepSet encoder for fair comparison paired with an MLP decoder, and train the autoencoder using a pointwise MSE reconstruction loss. We denote this baseline as (AE-InvE). Moreover, since point clouds reconstruction (see Sec. 2) can also generate permutation-invariant latent variables, we include an additional baseline aligned with that literature for a more comprehensive comparison. Concretely, we use the same DeepSet encoder and MLP decoder but train the autoencoder with the Chamfer distance as the reconstruction loss [12]. We denote this baseline as (AE-InvE-CD). Finally, because closure modeling does not necessarily require a reconstruction objective, we also consider the strategy of directly training an encoder jointly with a latent dynamics model. In practice, however, this optimization can be ill-posed: the latent representation may collapse, with the encoder producing constant or otherwise non-informative vectors [44], [45]. We nevertheless implement this approach as a baseline (InvE), where a DeepSet encoder generates the permutation-invariant latent variables, which are then concatenated with the macroscopic features of interest and passed to the dynamical model.

We evaluate our method and baselines on macroscopic variable prediction. Concretely, given an initial macroscopic state \({\boldsymbol{z}}_0\) generated by \(\bar{\boldsymbol{\varphi}}\) and the trained \(\hat{\boldsymbol{\varphi}}\) on \(X_0\), we apply the trained latent dynamics model (Eq. 4 ) to recursively predict future states \({\boldsymbol{z}}_1, {\boldsymbol{z}}_2, ..., {\boldsymbol{z}}_t\). We then compare the predicted and the ground-truth \(\bar{{\boldsymbol{z}}}_1, \bar{{\boldsymbol{z}}}_2, ..., \bar{{\boldsymbol{z}}}_t\). For deterministic dynamics, we report the mean relative error (MRE) [42]; for stochastic dynamics, we use the maximum mean discrepancy (MMD) [46], [47]. Details are provided in Appendix 6.2. We ensure that all models have the same latent-variable dimension for \(\hat{{\boldsymbol{z}}}\) per experiment to enable a fair comparison.

4.1 Visualizing Distributional Reconstruction↩︎

We first analyze the reconstruction of our distribution-aware autoencoder as we vary \(\epsilon\) in the target distribution \(q_X(\mathbf{x})\) and the encoder dimension \(\hat{z}_\mathrm{dim}\). To this end, we construct a 1D synthetic dataset comprising multiple particle sets. For each choice of \(\epsilon\), we induce the target distribution and minimize Eq. 3 to learn \(\hat{{\boldsymbol{z}}}\) under different values of \(\hat{z}_\mathrm{dim}\) Specifically, we generate \(3300\) sets, using \(3000\) for training and \(300\) for validation. Each set consists of \(200\) particles sampled from a Gaussian mixture distribution, with mixture means and variances varying across sets. For each \(\epsilon\), we train our model with different \(\hat{z}_\mathrm{dim}\) to approximate the induced target distribution, and the results are reported in Fig. 3.

Results indicate that the encoder bottleneck \(\hat{z}_\mathrm{dim}\) strongly affects reconstruction quality. When \(\hat{z}_\mathrm{dim}\) is small (e.g., \(\hat{z}_\mathrm{dim}=1\) in this case), the model frequently underfits the target: the learned density is overly smooth and fails to match the multimodal structure. Increasing \(\hat{z}_\mathrm{dim}\) improves reconstruction, enabling the learned distribution closely reconstruct the target density. We also observe a clear effect that a smaller \(\epsilon\) yields a more local and sharper target density (top row). Larger \(\hat{z}_\mathrm{dim}\) is required to reconstruct these details. As \(\epsilon\) grows (middle and bottom rows), the target becomes smoother, and correspondingly, the reconstruction becomes easier. A moderate latent dimension \(\hat{z}_\mathrm{dim}\) can achieve near-target fits. Overall, the figure demonstrates how \(\epsilon\) controls the smoothness of the target distribution, and a small \(\hat{z}_\mathrm{dim}\) is enough for a smooth target density.

4.2 Energy Evolution of Interacting Particle Systems↩︎

We study the energy evolution in particle systems. For this, we adopt the 2D interacting particle system from [48], where each particle interacts with all others through a pairwise step-force law \(F(r)=\tanh(a(1-r))+b .\) We set \(a=4\) and \(b=0.1\), and simulate trajectories using a forward Euler integrator with time step \(\Delta t=0.002\). The initial particle positions are sampled from a two-component Gaussian mixture model, with component means and isotropic covariances randomly drawn per trajectory. The 2D particle positions serve as microstates \(X\), and the \(\bar{\boldsymbol{\varphi}}\) computes the normalized pairwise interaction energy, which is the total pairwise interaction energy divided by the number of particle pairs (see Appendix 6.3 for the energy computation). Fig. 4 (a) visualizes a examplary trajectory: particles evolve from the initial configuration to a spatial pattern. Correspondingly, the normalized interaction energy decreases over time and converges to a steady value.

We generate \(10\)k simulated trajectories and split them into training/validation sets with an \(8/2\) split. Each trajectory contains 300 particles recorded at 301 time steps. Since this macroscopic evolution is deterministic, we model \({\boldsymbol{z}}_t\) using an ODE (i.e., setting \(\boldsymbol{\Sigma}({\boldsymbol{z}})=0\) in Eq. 4 ). We evaluate macroscopic prediction in three test regimes: First, we generate \(1\)k trajectories with the same number of particles as training, and sample initial positions also from a two-component GMM with independently resampled component means and isotropic covariances. This is for in-distribution evaluation, and we denote it as in-dst. Secondly, we generate \(1\)k trajectories with the same number of particles as training, but sample initial positions from a three-component GMM. This evaluates whether the learned macroscopic dynamics is generalizable to different initial patterns, and we denote it as diff-init. Third, we generate \(1\)k trajectories initialized from the same two-component GMM pattern but with more particles (\(400\)), denoted as diff-N. It evaluates whether the learned models are generalizable to different input sizes.

We first test the macroscopic variable prediction by our model in Fig. 4 (b). Specifically, we select three trajectories from the in-dst test set, provide only their initial macroscopic state \({\boldsymbol{z}}_0\), and recursively predict the future states. The prediction closely matches the ground truth throughout the rollout, capturing both the rapid initial decay and the convergence to a steady value. Furthermore, we test whether the learned macroscopic variable is permutation invariant. We take one in-dst test trajectory, apply three random permutations of the particle ordering, and compare the macroscopic prediction from our model and AE-Aug in Fig. 4 (c). AE-Aug produces noticeably different energy predictions under different permutations, even with permutation augmentation during training. In contrast, our method is exactly permutation-invariant: the three prediction curves overlap, appearing as a single curve. We then compare our method to baselines quantitatively in Table 1. On in-dst, our method achieves the lowest error, substantially improving over the plain AE-Aug baseline. On the particle-count shift (diff-N), our method again attains the best performance, suggesting that the learned macroscopic dynamics transfers across different system sizes. Under the initialization-pattern shift test (diff-init), our method remains competitive and improves over AE-Aug and AE-InvE, while AE-InvE-CD attains the lowest error. We also examine the influence of \(\hat{{\boldsymbol{z}}}\) dimension and the smoothing scale \(\epsilon\), and the sensitivity test is provided in Appendix 8.2.

4.3 Mixing of Binary Particles↩︎

We next study a stochastic multi-species system where the macroscopic state quantifies the mixing of two types of particles. Specifically, we consider a two-dimensional binary particle system mixing under a pairwise Lennard-Jones potential, a standard model in atomistic binary fluids [49], [50]. We simulate \(529\) particles in a square domain of side length \(l=32\) with reflective boundaries in both \(x\) and \(y\). At initialization, particles are randomly placed, and a vertical boundary is placed at a random location to partition the domain into two regions: particles on one side are labeled type A and those on the other side are labeled type B (Fig. 5 (a) left). We run the NVT simulation at temperature \(1.0\) for \(10^5\) steps and record particle positions every \(200\) steps as microstates \(X_t\). Fig. 5 (a) shows a representative trajectory, where the initially separated particles gradually interpenetrate until reaching a well-mixed configuration. We compute the mixing ratio as the macroscopic variable \(\bar{{\boldsymbol{z}}}\), i.e., \(\bar{{\boldsymbol{z}}} = (R_{AB}, R_{BA})\), where \(R_{AB}\) is the probability that a neighbor of a type-A particle is type-B (and vice versa for \(R_{BA}\)), computed with a fixed cutoff radius. The mixing ratio is equivalent to the short-range order measure commonly used in chemistry [51]. Intuitively, \(\bar{{\boldsymbol{z}}}\) increases as mixing progresses and converges when the two species become uniformly mixed. We generate \(10\)k trajectories with sampling the initial separation boundary uniformly in \([10,22]\), and use the \(8/2\) split for training and validation. Similar to the previous experiment, we evaluate the rollout prediction in three regimes. First, we generate 200 test trajectories with the separation boundary sampled in the same range \([10,22]\), denoted as in-dst. Then, we generate test data diff-N with a reduced domain height while keeping density fixed, resulting in \(437\) particles per simulation. The separation boundary is also sampled in \([10,22]\). Finally, we simulate 200 trajectories with the separation boundary sampled outside the training range (uniformly in \([8,10] \cup [22, 24]\)), denoted as diff-init.

To accommodate two particle types, our method applies a shared encoder and a shared decoder to each species separately, and then aggregates these two embeddings to form the closure variable \(\hat{{\boldsymbol{z}}}\) for the full system. For a fair comparison, we configure each baseline to use particle-type information whenever its architecture permits; otherwise, we train it on the full particle set without distinguishing types (see Appendix 6.4 for details). Because of the thermal fluctuation, we model the dynamics of \({\boldsymbol{z}}_t\) with an SDE (i.e., learning both \(\boldsymbol{g}\) and \(\boldsymbol{\Sigma}\) in Eq. 4 ).

The quantitative comparisons between our methods and baselines are reported in Table 2. Our method achieves the lowest MMD on all test settings. Comparing diff-init to in-dst, MMD increases for all methods because the macrostates are unseen during training, but our method still performs best. The MMD is lower for diff-N than for in-dst for all methods because the system of diff-N is a subset of the system used in training. Moreover, we evaluate macroscopic predictions in three representative settings by placing the initial separation boundary on the left, middle, or right of the domain at \(x \in {12,16,20}\). For each setting, we generate 100 trajectories with randomly sampled initial particle positions. Starting from the initial macrostates \({\boldsymbol{z}}_0\), we compare the predicted \({\boldsymbol{z}}_t\) to the ground-truth. We compare our method to baseline InvE, which is the reconstruction-free method. The predictions are as shown in Fig. 5 (b). Across all three cases, the means of our predictions are close to the ground-truth means for both \(R_{AB}\) and \(R_{BA}\), whereas InvE underestimates the mixing ratios and tends to plateau prematurely. These results support the observation that closure modeling does not necessarily require reconstruction in principle, but directly learning latent dynamics can be challenging because the representation may collapse [44], [45]. Careful initialization of learnable parameters or additional stabilization techniques may alleviate this issue, and we leave it for future work.

4.4 Polymer Extension↩︎

Besides particle systems, many scientific applications provide data in image format [52], [53]. We therefore test whether our method can learn closure variables for macroscopic dynamics modeling from videos of microstate trajectories, even though it is not specifically designed for image inputs. To this end, we consider the polymer dynamics dataset from [2], which simulates the dynamics of polymer chains using the Brownian dynamics approach in a planar elongational flow. We convert each polymer frame (3D bead coordinates) into a 2D image by rendering a Gaussian blob at each bead’s (x, y) coordinate, with the blob width determined by the z coordinate. Additional details on the original dataset and our image construction procedure are provided in Appendix 6.5.

To apply our method, we treat each non-white pixel (i.e., pixels occupied by the polymer) as a “particle”, using its image coordinates together with its grayscale intensity as features (\({\boldsymbol{x}}^i \in \mathbb{R}^3\)). The collection of these non-white pixels defines the microstate \(X\), and we compute the stretching length \(\bar{{\boldsymbol{z}}}\) from \(X\). Note that the number of non-white pixels varies across time steps, and our model can deal with variable-sized microstates.

We evaluate the learned macroscopic dynamics on the three test scenarios provided in [2], referred to as Test Fast, Test Medium, and Test Slow. Fig. 6 (a) visualizes the microstates of the first trajectory at different time steps, indicating the different extension speeds. Fig. 6 (b) compares our prediction to the ground-truth. It shows that our model captures the stretching-length dynamics for Test Fast and Test Medium, but overestimates the stretching rate in Test Slow. A likely reason is that the pattern of the initial configuration in Test Slow appears visually similar to that in Test Fast, making it challenging for the learned closure variable to distinguish these two. We also test the state-of-the-art image models, including CNN [54] and Vision Transformers [55], to learn the latent variables by reconstructing input images. They exhibit similar results to our method, as they also succeed on Test Fast and Test Medium but fail on Test Slow. More details are provided in Appendix 6.5.

5 Conclusion and Discussion↩︎

This work revisits a classical problem of learning macroscopic dynamics from high-dimensional microscopic observations. The novel insight delivered by it is that, instead of minimizing pointwise reconstruction loss, we can train an autoencoder to reconstruct the distributional information of the input. Following this idea, we design the decoder trained to recover a smooth density representation of macrostates, making the training objective inherently permutation-invariant and allowing inputs of varying size.

A few previous works also use the decoder as the density on modeling point sets [56]–[58]. In these approaches, the observed points are treated as samples from an underlying distribution, and the decoder is trained by maximizing the (conditional) likelihood of the point set \(\{{\boldsymbol{x}}_1, \ldots, {\boldsymbol{x}}_n\}\), i.e., by fitting \(p_\theta(\mathbf{x}|\hat{{\boldsymbol{z}}})\) to the empirical samples so that the model can generate new samples. In contrast, our goal is not microstate generation, but the learning of latent closure variables. We induce, for each input set \(X\), a target density \(q_X\) and train the decoder by minimizing \(\mathrm{KL}(q_X\|p_\theta)\). By varying \(\epsilon\) in \(q_X\), we can tune the resolution at which the decoder approximates.

A practical limitation is that the approach can struggle for stiff systems, where small microscopic perturbations produce large macroscopic changes. We believe this is not a limitation specific to our method, but a challenge common to all closure modeling methods. For example, consider heating a solid: throughout the process, the atoms remain close to equilibrium configurations, yet the temperature can vary substantially. The method can fail to learn a meaningful closure variable because the underlying microstates are nearly indistinguishable. In such regimes, learning a stable mapping from microstate distributions to useful closure variables may require much denser coverage of the microstate space. As a result, the current method is best suited to systems where macroscopic evolution is accompanied by substantial changes in microscopic configurations, so that the induced target distributions vary sufficiently to train the autoencoder effectively.

Several extensions are promising, including exploring other set encoders such as attention-based Set Transformers, investigating alternative density models such as score-based generative models, and improving the training strategy to couple closure modeling with dynamical prediction more tightly. We expect to apply the method to solve more real-world problems.

Acknowledgements↩︎

This research is supported by the National Research Foundation, Singapore under its AI Singapore Programme (AISG Award No: AISG3-RP-2022-028). The computational work for this article was partially performed on resources of the National Supercomputing Centre, Singapore (https://www.nscc.sg). We thank Zhuoyuan Li and Aiqing Zhu for helpful discussions, and Fushuai Wang and Prof. Beatrice W. Soh for kindly sharing the script to generate polymer images from simulated polymer chains.

Impact Statement↩︎

This work addresses a gap in closure modeling to learn permutation-invariant macroscopic dynamics from unordered microscopic observations. This can benefit multiscale modeling in particle-based domains, improve the scientific surrogate models and accelerate exploratory simulation and analysis. The main risk is that the prediction can become inaccurate under distribution shift or new boundary conditions that never appear in the training data. We recommend a careful out-of-distribution check when it is used for making decisions.

6 Experiments Details↩︎

6.1 Pseudocode of Train and Inference↩︎

6.2 Evaluation Metrics↩︎

For the ODE dynamics, the evaluation metric is the mean relative error (MRE) which is defined as: \[\text{MRE}(\mathcal{T}_{\text{test}}) = \frac{1}{|\mathcal{T}_{\text{test}}|} \sum_{ \mathcal{T}_{\text{test}} } \left( \frac{ \sum_t \| \bar{{\boldsymbol{z}}}^{\text{true}}_t - \bar{{\boldsymbol{z}}}^{\text{pred}}_t \|_2^2 }{ \sum_t \| \bar{{\boldsymbol{z}}}^{\text{true}}_t \|_2^2 } \right)~,\] where \(\mathcal{T}_{\text{test}}\) represents the trajectories in the test data and \(|\mathcal{T}_{\text{test}}|\) is the number of test trajectories. \(\{\bar{{\boldsymbol{z}}}^{\text{pred}}_{t}\}_{t=1}^{T}\) and \(\{\bar{{\boldsymbol{z}}}^{\text{true}}_{t}\}_{t=1}^{T}\) denote the predicted and ground-truth states for a testing trajectory.

For the SDE macrodynamics, we also predict the whole trajectory for each initial state autoregressively in the testing data and evaluate the discrepancy between the predicted and true trajectory distributions using the time-average multi-kernel maximum mean discrepancy (MMD). Again, let \(\{\bar{{\boldsymbol{z}}}^{\text{pred}}_{t}\}_{t=1}^{T}\) and \(\{\bar{{\boldsymbol{z}}}^{\text{true}}_{t}\}_{t=1}^{T}\) denote the predicted and true states, and let \(\bar{{\boldsymbol{z}}}^{\text{pred},(i)}_{t}\) and \(\bar{{\boldsymbol{z}}}^{\text{true},(i)}_{t}\) be the \(i\)-th trajectory in \(\mathcal{T}_{\text{test}}\) (suppose \(|\mathcal{T}_{\text{test}}|=N\)) at time \(t\). The per-timestep multi-kernel MMD is computed as \[\begin{align} \widehat{\mathrm{MMD}}^{2}_{t} &= \frac{1}{|\Gamma|}\sum_{\sigma\in\Gamma} \Bigg[ \frac{1}{N^{2}}\sum_{i=1}^{N}\sum_{j=1}^{N} k_{\sigma}\!\left(\bar{{\boldsymbol{z}}}^{\text{pred},(i)}_{t},\bar{{\boldsymbol{z}}}^{\text{pred},(j)}_{t}\right) + \frac{1}{N^{2}}\sum_{i=1}^{N}\sum_{j=1}^{N} k_{\sigma}\!\left(\bar{{\boldsymbol{z}}}^{\text{true},(i)}_{t},\bar{{\boldsymbol{z}}}^{\text{true},(j)}_{t}\right)\\ &- \frac{2}{N^{2}}\sum_{i=1}^{N}\sum_{j=1}^{N} k_{\sigma}\!\left(\bar{{\boldsymbol{z}}}^{\text{pred},(i)}_{t},\bar{{\boldsymbol{z}}}^{\text{true},(j)}_{t}\right) \Bigg], \end{align}\] where the RBF kernel is \[k_{\sigma}({\boldsymbol{a}},{\boldsymbol{b}})=\exp\!\left(-\frac{\|{\boldsymbol{a}}-{\boldsymbol{b}}\|_2^{2}}{2\sigma^{2}}\right), \qquad \Gamma=\{0.01,\,0.05,\,0.25,\,1.25\}.\] Finally, averaging over timesteps and test trajectories yields \[\mathrm{MMD}(\mathcal{T}_{\text{test}}) = \frac{1}{|\mathcal{T}_{\text{test}}|}\sum_{\mathcal{T}_{\text{test}}} \left( \frac{1}{T}\sum_{t=1}^{T}\widehat{\mathrm{MMD}}^{2}_{t} \right).\]

6.3 Interaction Particle System↩︎

The total pairwise energy of the interaction particle studied in Sec. 4.2 is computed as \[E_\text{total}=\sum_{i<j}\Big[\frac{1}{a}\log\cosh\big(a(1-r_{ij})\big)+b(1-r_{ij})\Big],\] where \(r_{ij}\) is the distance between particle \(i\) and particle \(j\). We set \(a=4\) and \(b=0.1\). We divide this total energy by \(n \times (n-1) / 2\) to get the normalized pairwise interaction energy.

6.4 Binary Particle Mixing↩︎

Because the numbers of type-A and type-B particles vary across trajectories, baselines with an MLP decoder cannot reconstruct each species separately. Therefore, we train AE-Aug, AE-InvE, and AE-InvE-CD to reconstruct the states of all particles jointly, without conditioning on type. For baselines (AE-InvE, AE-InvE-CD and InvE) that can incorporate type information, we use a shared encoder applied to each type and aggregate the resulting embeddings into a single closure variable, matching the setup of our method for a fair comparison. The encoder of baseline AE-Aug just takes all particles as the input to avoid the varying numbers of type-A and type-B particles.

6.5 Polymer Image Dataset↩︎

The original polymer dataset is from [2]. In detail, they simulate the dynamics of polymer chains using Brownian dynamics approach in a planar elongational flow. Each polymer chain consisted of 300 beads moving in (x,y,z) dimensions. Every trajectory contains 1001 time steps. They provide 610 trajectories for training, 110 trajectories for validation and three testing cases (“test fast”, “test medium” and “test slow”). In each testing case, they simulate 100 trajectories of the same initial configuration of the 300 beads.

Given the polymer simulation, we construct the video dataset in the following way. Each polymer frame (3D bead coordinates) is projected into a 2D image by placing a Gaussian blob at each bead’s \((x, y)\) position on a fixed 500×100 grid. The Gaussian width depends on the bead’s \(z\)-offset from the frame’s mean \(z\) (larger \(|z|\), wider blob), then all blobs are summed, normalized to the frame’s max intensity, and quantized to 8-bit grayscale. Fig.9 visualizes an arbitrary trajectory in the video dataset.

We test the state-of-the-art image models, CNN and ViT, for learning the closure variables. We train CNN and ViT autoencoders on image reconstruction with MSE loss and use the latent vector as the closure variables. The macroscopic dynamics learned with their closure variables are shown in Fig. 10.

7 Model Implementation Details↩︎

7.1 Encoder↩︎

The encoder implements \(\hat{\boldsymbol{\varphi}}: \mathcal{X} \rightarrow \hat{{\mathcal{Z}}}\) (Sec. 3.2). We choose DeepSet [8] but any permutation-invariant architecture can also be the choice. In detail, the elementwise function \(\phi\) is applied on the input feature of every particle and outputs a latent representation. We use mean pooling to aggregates the latent representations to the set-level representation. Finally, the set-level neural network \(\rho\) transfers the aggregated set-level representation to the final \(\hat{\boldsymbol{z}}\). We implement the elementwise function \(\phi\) and the set-level function \(\rho\) by MLPs.

7.2 Decoder↩︎

The decoder \(\boldsymbol{\psi}\) parameterizes a conditional density \(p(\mathbf{x}| \hat{{\boldsymbol{z}}})\), i.e., \(p(\mathbf{x}| \hat{{\boldsymbol{z}}}) = \boldsymbol{\varphi}(\mathbf{x}; \hat{{\boldsymbol{z}}})\). We instantiate \(p(\mathbf{x}| \hat{{\boldsymbol{z}}})\) with a stack of conditional autoregressive rational quadratic spline flow layers  [59] since it is powerful to approximate complex densities and allows us to compute the likelihood conveniently. But any conditional density model is applicable. In detail, the conditional normalizing flow defines an invertible mapping \(\mathbf{x}= f(\mathbf{u};\hat{{\boldsymbol{z}}})\) from a simple base random variable \(\mathbf{u}\) (e.g., \(\mathbf{u}\sim \mathcal{N}(\mathbf{0},\mathbf{I})\)) to the particle feature \(\mathbf{x}\). By the change-of-variables formula, the conditional log-likelihood \(\log p(\mathbf{x}| \hat{{\boldsymbol{z}}})\) can be computed exactly, allowing for training by maximizing likelihood.

7.3 The Induced Density↩︎

The \(q_{X}(\mathbf{x})\) constructed in Sec. 3.2 can be viewed as a kernel density estimate of the empirical measure associated with \(X\). As \(\epsilon\) goes close to 0 (with \(|X|\) fixed), \(q_X\) concentrates and converges in the weak sense to the empirical measure \(\frac{1}{|X|}\sum_j \delta_{{\boldsymbol{x}}_j}\) (i.e., a sum of point masses) [60]. When the decoder family \(p(\cdot| \hat{\boldsymbol{z}})\) is restricted to smooth densities (e.g., a normalizing flow), we keep \(\epsilon>0\) and interpret \(p(\cdot| \hat{\boldsymbol{z}})\approx q_X\) as distributional recovery at resolution \(\epsilon\), rather than exact pointwise reconstruction of every particle in the limit \(\epsilon\to 0\).

8 Supplementary Numerical Results↩︎

8.1 Compatibility with Encoders Beyond DeepSet↩︎

To demonstrate that the proposed framework is not restricted to DeepSet, we replaced the DeepSet encoder with a Set Transformer [9] in the interacting particle system experiment while keeping all other settings unchanged. The results are reported in Table 3.

The Set Transformer achieves lower rollout prediction errors in the in-dst and diff-N settings, whereas DeepSet is better under the diff-init distribution shift setting. These results suggest that more expressive permutation-equivariant encoders can improve performance when the testing distribution is close to the training, but may also learn more specialized latent representations that are less robust under distribution shift. Overall, the experiment empirically demonstrates that the proposed framework is compatible with different encoders and that the main contribution is orthogonal to the encoder choice.

8.2 Sensitivity to the Latent Dimension \(\hat{z}_\mathrm{dim}\) and Smoothing Scale \(\epsilon\)↩︎

In general, closure modeling aims to compress the high-dimensional microscopic information into a low-dimensional latent variable, which can be used as the state vector for modeling the system’s macroscopic dynamics. In our method, the latent dimension \(\hat{z}_\mathrm{dim}\) and smoothing scale \(\epsilon\) both affect the information compression, i.e., how much microscopic information is encoded into the closure variable. Specifically, \(\hat{z}_\mathrm{dim}\) controls the information compression capacity, while \(\epsilon\) controls the information compression resolution. A larger \(\epsilon\) smooths out fine-scale details, whereas a smaller \(\epsilon\) preserves detailed information.

To study their influence, we perform a sensitivity analysis on the in-dst setting of the interacting particle system. Table 4 shows that the rollout prediction error decreases as the latent dimension increases and becomes stable once \(z^{\mathrm{dim}} \geq 5\), suggesting that the latent representation has sufficient capacity to encode the microscopic information needed for modeling macroscopic dynamics beyond this point. The large error with small \(\hat{z}_\mathrm{dim}\) (e.g., \(\hat{z}_\mathrm{dim}= 1\)) is under expectation because the too small \(\hat{z}_\mathrm{dim}\) poses an overly restrictive bottleneck to encode microscopic information.

Table 5 shows that the method remains relatively stable across a broad range of smoothing scales. Interestingly, even extremely small values of \(\epsilon\) (e.g., \(10^{-8}\)) still produce useful closure variables for macroscopic dynamics prediction. The reason is as follows. In the limit \(\epsilon \rightarrow 0\), the reconstruction objective (Eq. 3 ) effectively reduces to maximizing likelihood on the observed point set. In this regime, the learned density no longer accurately reconstructs the smoothed target density, but still preserves sufficient underlying distributional structure of the particles, which remains informative for predicting macroscopic dynamics.

References↩︎