Learning Macroscopic Dynamics from Partial Microscopic Observations

Authors: mengyi  chen
qianxiao  li
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 16:14:59

Abstract: Macroscopic observables of a system are of keen interest in real applications such as the design of novel materials. Current methods rely on microscopic trajectory simulations, where the forces on all microscopic coordinates need to be computed or measured. However,  this can be computationally prohibitive for realistic systems.  In this paper, we propose a method to learn macroscopic dynamics requiring only force computations on a subset of the microscopic coordinates. Our method relies on a sparsity assumption: the force on each microscopic coordinate relies only on a small number of other coordinates. The main idea of our approach is to map the training procedure on the macroscopic coordinates back to the microscopic coordinates, on which partial force computations can be used as stochastic estimation to update model parameters. We provide a theoretical justification of this under suitable conditions. We demonstrate the accuracy, force computation efficiency, and robustness of our method on learning macroscopic closure models from a variety of microscopic systems, including those modeled by partial differential equations or molecular dynamics simulations.

Main paper: 
 introduction: 
 Macroscopic properties, including thestructural and dynamical properties, provide a way to describe and understand the collective behaviors of complex systems. In a wide range of real applications, researchers focus mainly on the macroscopic properties of a system, e.g., the viscosity and ionic diffusivity of liquid electrolytes for Li-ion batteries (Dajnowicz et al., 2022) . Macroscopic observables usually depend on the whole microscopic system, e.g., the calculation of mean squared displacement requires all microscopic coordinates during the simulation. With growing simulation and experimental data, data-driven learning of macroscopic properties from microscopic observations has become an active area of research (Zhang et al., 2018; Wang et al., 2019; Husic et al., 2020; Lee et al., 2020; Fu et al., 2023; Chen et al., 2024) .Accurate calculation of macroscopic properties requires large-scale microscopic simulation. However, accurate force computations on all microscopic coordinates for large systems are extremely expensive (Jia et al., 2020; Musaelian et al., 2023) . For example, in ab initio molecular simulations, accurate forces need to be calculated from density functional theory (DFT). The computational cost of DFT limits its application to relatively small systems, typically ranging from a few hundred atoms to several thousand atoms, depending on the level of accuracy and computation resources (Hafner et al., 2006; Luo et al., 2020) . This poses a dilemma: Accurate macroscopic properties are obtained from large-scale microscopic simulation, but the computation of forces on all the microscopic coordinates is extremely challenging.Figure 1 : Overview of our method. Left. Data generation workflow. For each configuration x, forces on a subset of all the microscopic coordinates are calculated by the microscopic force calculator. Right. Macroscopic dynamics identification. The macroscopic dynamics is mapped to the microscopic space first, then compared with the forces on a subset of the microscopic coordinates.To solve the dilemma, the corresponding question is: Can we still obtain accurate macroscopic observables even though only access to forces on a subset of the microscopic coordinates? In this work, we develop a method to learn the dynamics of the macroscopic observables directly, while only forces on a subset of the microscopic coordinates are needed. Efficient partial computation of microscopic forces relies on the sparsity assumption, where the computation cost of forces on a subset of microscopic coordinates does not scale with the microscopic system size. To learn the dynamics of the macroscopic observables, we first map the macroscopic dynamics back to the microscopic space, then compare it with the partial microscopic forces. Our key idea is summarized in Fig. 1 .Our main contributions are as follows:• We develop a novel method that can learn the macroscopic dynamics from partial computation of the microscopic forces. Our method can significantly reduce the computational cost for force computations.• We theoretically justify that forces on a subset of the microscopic coordinates can be used as stochastic estimation to update latent model parameters, even if the macroscopic observables depend on all the microscopic coordinates.• We empirically validate the accuracy, force computation efficiency, and robustness of our method through a variety of microscopic dynamics and latent model structures. 
related work: 
 Learning from Partial Observations Several works have sought to learn dynamics from partially observed state x utilizing machine learning (Ruelle & Takens, 1971; Sauer et al., 1991; Takens, 2006; Ayed et al., 2019; Ouala et al., 2020; Huang et al., 2020; Schlaginhaufen et al., 2021; Lu et al., 2022; Stepaniants et al., 2023) . For training, these methods reconstruct the unobserved state x first and model the dynamics of x " px, xq. Our work assumes full state x but partial forces f . Furthermore, we do not model the dynamics on state x directly, but rather on the latent space because the dimension of x would be extremely high for large systems.Reduced Order Models By modeling the dynamics in the latent space and then recovering the full states from them, reduced order models (ROMs) substitute expensive full order simulation with cheaper reduced order simulation (Schilders et al., 2008; Fresca et al., 2020; Lee & Carlberg, 2020; Hernandez et al., 2021; Fries et al., 2022) .Our method can be thought to fall in the range of closure modeling. Unlike ROMs, we aim to model the dynamics of some given macroscopic observables directly and are not interested in recovering the microscopic states from the latent states.Equation-free Framework The equation-free framework (EFF) has sought to simulate the macroscopic dynamics efficiently (Kevrekidis et al., 2003; Samaey et al., 2006; Liu et al., 2015) . The EFF is usually applied to partial differential equation (PDE) systems, and the macroscopic observables are chosen to be the solution of the PDE at the coarse spatial grid. In EFF, the macroscopic observables depend locally on the microscopic coordinates, allowing the macroscopic dynamics to be directly estimated from the microscopic simulations performed in small spatial domains. In contrast, the macroscopic observables may depend globally on the microscopic coordinates in our method, and the macroscopic dynamics may not be easily estimated from microscopic simulations performed in small spatial domains.Another difference is that our method explicitly learns the macroscopic dynamics, while EFF can bypass explicit derivation of macroscopic evolution law by coupling microscale and macroscale dynamics. During simulation, EFF still requires microscopic simulation to be performed in small spatial domains and for short times, but our method can enable fast macroscopic simulation without requiring any microscopic simulation However, for systems where the macroscopic evolution equations conceptually exist but are not available in closed form, EFF can efficiently handle such cases, but the learned dynamics in our method may involve approximation or statistical errors that are often challenging to estimate. 
problem setup: 
 We consider a microscopic system consisting of n particles. Let the state of the microscopic system be x " px 1 , ¨¨¨, x n q P R N , x i P R m , N " mn, where x i P R m is some physical quantity associated with the i-th particle, such as the position and velocity. Assume the dynamics of the microscopic system can be characterized by an ordinary differential equation(ODE):dxptq dt " f pxptqq (1)where f pxq " pf 1 pxq, ¨¨¨, f n pxqq P R N . We will call x i the microscopic coordinate of the i-th particle and f i the force acting on the microscopic coordinate x i . In many real applications, we are interested in the dynamics of some macroscopic observables z ˚" φ ˚pxq. Here φ ˚is given beforehand and describes the functional dependence of z ˚on x. For example, z ˚can be chosen to be the instantaneous temperature or mean squared displacement in a Lennard-Jones system.The goal is to learn the dynamics of z ˚from microscopic simulation data. Existing methods that try to learn the macroscopic or latent dynamics require microscopic trajectories or forces on all the microscopic coordinates for training (Champion et al., 2019; Fries et al., 2022; Fu et al., 2023; Chen et al., 2024) . The problem is: When the microscopic system size N is very large such that the force computations on all the microscopic coordinates are impossible, these methods are no longer applicable. Instead, our method aims to learn from partial computation of microscopic forces.Consider we are given a microscopic force calculator S for computation of partial microscopic forces. Let the microscopic coordinate x be sampled from a distribution D. For each x, the microscopic force calculator S will first sample an n dimensional random variable Ipxq " pI 1 pxq, ¨¨¨, I n pxqq " P x P t0, 1u n according to a certain strategy. Next S will calculate the corresponding partial forces f Ipxq :" pf I1pxq , ¨¨¨, f Inpxq q. For notation simplicity sometimes we will simply write f Ipxq as f I . I i pxq indicate whether partial i is chosen to calculate the force or not. If particle i is chosen, I i pxq " 1, f Iipxq " f i , otherwise I i pxq " 0, f Iipxq " 0. We require the sampling strategy P x to satisfy:1. For each Ipxq " P x , exactly n ¨p items are equal to 1 and the rest are 0. 2. Each particle can be chosen with equal probability p, i.e. PpI i pxq " 1q " p, PpI i pxq " 0q " 1 ´p, i " 1, ¨¨¨, n.This means that the microscopic force calculator S can calculate forces on n ¨p microscopic coordinates, and 0 ă p ă 1 limits the computation capacity of the microscopic force calculator S.The above requirement is consistent with real applications since it is difficult to calculate all the microscopic forces due to computational cost. Furthermore, for efficient calculation of the partial microscopic forces, we will assume f satisfies the following sparsity assumption:Assumption: For a given error tolerance ϵ ą 0, there exists a constant M ! n, such that for any x " D and i P t1, ¨¨¨, nu, we can always find an index set Jpx i q Ă ti " 1, ¨¨¨, nu, |Jpx i q| ă M which satisfies:||f i px 1 , ¨¨¨, x n q ´f i ptx i u iPJpxiq q|| 2 ă ϵ (2)Intuitively, the assumption implies that the computational cost of force f i is independent of the microscopic system dimension N . Thus our microscopic force calculator S can compute partial forces in an efficient way. This assumption is prevalent in real-world applications. To better illustrate this, we give two examples here. The first example is about molecular dynamics. In molecular dynamics, each x i represents the position r i and velocities v i of the i-th atom, i.e., x i " pr i , v i q P R 6 . Then Eq. ( 1) becomes the Newton's law of motion:EQUATION)EQUATIONIt is common to limit the range of pairwise interactions to a cutoff distance (Allen et al., 2004; Zhou & Liu, 2022; Vollmayr-Lee, 2020) . To calculate the force on an atom, we only need to consider its interaction with other atoms that are within the cutoff. The second example is about systems modeled by partial differential equation (PDE). We consider a time-dependent PDE and apply finite difference scheme to discretize the spatial derivatives. Then, the resulting semi-discretized equation takes the form of Eq. (1), and each x i is the value at the i-th grid. f i only depends on those grids that are used for finite difference approximation of the spatial derivatives.Let the training data generated by the microscopic force calculator S be tx i , f Ipx i q u i"1,¨¨¨,K . The data generation procedure is provided in Algorithm 1. We will introduce in the next section how we can learn the macroscopic dynamics from the training data with partial forces. 
method: 
 Existing works for macroscopic dynamics identification consist of two parts: dimension reduction and macroscopic dynamics identification (Fu et al., 2023; Chen et al., 2024) . We will follow these two parts. We start with most standard parts of closure modeling with an autoencoder, next, we turn to the main difficulty of macroscopic dynamics identification from partial forces.We will use an autoencoder to find the closure ẑ " φpxq to z ˚" φ ˚pxq such that z " pz ˚, ẑq forms a closed system. Here we define z as forming a closed system if its dynamics 9 z depends only on z, not any external variables. Note that in z ˚" φ ˚pxq, φ ˚is determined beforehand and contains no trainable parameters. This ensures z ˚represents the desired macroscopic observables and remains unchanged during the training of the autoencoder. Denote the encoder by φ " pφ ˚, φq and the decoder by ψ, we will minimize the following reconstruction loss:L rec " 1 K ř K i"1 ∥x i ´ψ ˝φpx i q∥ 2 2 (5)We also want φ 1 pxqφ 1 pxq T to be well-conditioned (see Section 4.2), then we impose constraints on the condition number of φ 1 pxqφ 1 pxq T :L cond " 1 K ř K i"1 ∥κpφ 1 px i qφ 1 px i q T q ´1∥ 2 2 (6)By enforcing φ 1 px i qφ 1 px i q T to be well-conditioned, we are also enforcing φ 1 px i q P R dˆN , d ! N to have full row rank, which will be used later. The overall loss to train the autoencoder is :EQUATIONAlgorithm 1 Data generation.Require: D:configuration distribution S: microscopic force calculator K: training data size P: partial index sampling strategy 1: for i " 1 to K do 2:x i " D 3:Ipx i q " P x i 4: calculate f Ipx i q px i q using S return tx i , f Ipx i q px i qu i"1,¨¨¨,K Calculate L x,p in Eq. ( 12) with tx i , f x i px i qu iPJ 4: θ t`1 Ð opt(θ t , ∇L x,p ) return model parameter λ cond is a hyperparameter to adjust the ratio of L cond and is chosen to be quite small in our experiments, e.g., 10 ´5 or 10 ´6. The aim of training the decoder ψ is to help the discovery of the closure variables. The decoder ψ will not be used for further macroscopic dynamics identification. To facilitate comparison between models tainted with all and partial forces, we will train the autoencoder first and freeze it for macroscopic dynamics identification.We now address the difficulty of learning from data with partial microscopic forces. Substitute z " φpxq into equation Eq. ( 1) and make use of chain rule, we get the dynamics of z:dz dt " φ 1 pxqf pxq, zp0q " φpx 0 q (8)here we use φ 1 pxq to denote the Jacobian of ∇ x φpxq for notation simplicity. If the dynamics of z is closed, the right-hand side of Eq. (8) will only depend on z, and we parametrize it using a neural network g θ pzq « φ 1 pxqf pxq. Since we are only interested in macroscopic dynamics identification, the loss would be naturally defined on the macroscopic coordinates:EQUATIONEq. ( 9) is used commonly in existing work (Champion et al., 2019; Fries et al., 2022; Bakarji et al., 2022; Park et al., 2024) .The main difficulty of L z is that it includes the matrix-vector product φ 1 pxqf pxq. Note that the i-th entry of φ 1 pxqf pxq can be written as ř n j"1 φ 1 ij pxqf j pxq, and it is difficult to find an unbiased estimation of the i-th entry using a subset of tf j pxqu j"1,¨¨¨,N . Thus the accurate calculation of L z requires the forces tf j pxqu j"1,¨¨¨,N on all the microscopic coordinates.The main idea of our method is to map the loss on the macroscopic coordinates back to the microscopic coordinates:L x pθq " 1 K ř K i"1 ||f px i q ´pφ 1 px i qq : g θ pz i q|| 2 (10)where pφ 1 px i qq : P R N ˆd is the Moore-Penrose inverse. Since φ 1 px i q is of full row rank, pφ 1 px i qq : is in fact the right inverse of φ 1 px i q, i.e., φ 1 px i qpφ 1 px i qq : is an identity matrix. Below we will show our main theoretical result: Theorem 1. Assume for any x " D, the eigenvalues of φ 1 pxqφ 1 pxq T are lower bounded by b 1 and upper bounded by b 2 , 0 ă b 1 ď b 2 . Then:EQUATIONhere C does not depend on θ hence does not affect the optimization.The proof relies on singular value decomposition of φ 1 pxq and we provide the full proof in Appendix A.1. Theorem 1 states that by minimizing L x pθq, we are actually narrowing the range of L z pθq. Hence we want b 1 and b 2 to be as close as possible, this is the reason why we constrain the condition number of φ 1 pxqφ 1 pxq T in Eq. ( 6). In the very special case where b 1 " b 2 , minimizing L x pθq is just equivalent to minimizing L z pθq.Note that in loss L x , the term ||f pxq ´pφ 1 pxqq : g θ pzq|| can be rewritten as ř n j"1 ||f j pxq ṕφ 1 pxqq : j g θ pzq||, and 1 p ř jPIpxq ||f j pxq ´pφ 1 pxqq : j g θ pzq|| can be regarded as its unbiased stochastic estimation. Then, we can introduce our loss defined for partial forces:L x,p pθq " 1 pK ř K i"1 ∥f Ipx i q px i q ´pφ 1 px i qq : Ipx i q g θ pz i q∥ 2 2 (12)Here a constant 1{p is multiplied to L x,p to guarantee:EQUATIONWe provide the full proof of Eq. ( 13) in Appendix A.2. By training with the model with L x,p , we can use data with partial forces as stochastic estimation to update model parameters. The full training procedure is provided in Algorithm 2.Note that in Algorithm 1 during the data generation, Ipx i q is also sampled from its distribution. Thus, L x,p is deterministic once the samples tx i , f Ipx i q px i qu i"1,¨¨¨,K are generated. In the limit, the estimation is unbiased: Theorem 2 (informal). Let Lx pθq " EL x pθq, θ ˚P arg min θ Lx pθq, θ K,p P arg min θ L x,p pθq, then under certain conditions: Lx pθ K,p q ´L x pθ ˚q a.s.ÝÑ 0The proof utilized Rademacher complexity and a crucial assumption used in the proof is the uniform boundedness of L x,p . The formal version of Theorem 2 and the complete proof is provided in Appendix A.3. Theorem 2 theoretically justifies the expected risk Lx pθ K,p q at the optimal parameter found by L x,p , converges to the optimal expected risk Lx pθ ˚q as K goes to infinity. 
experiments: 
 In this section, we experimentally validate the accuracy, force computation efficiency, and robustness through a variety of microscopic dynamics.We first consider a one-dimensional spatiotemporal Predator-Prey system, mainly to validate the correctness and force computations efficiency of our method.Predator-Prey System The simplified form of the Predator-Prey system with diffusion (Murray, 2003) isEQUATIONwhere u, v denote the dimensionless populations of the prey and predator respectively, a, b, D are three parameters. The complex dynamics of Predator-Prey interaction, including the pursuit of the predator and the evasion of the prey in ecosystems, can be described by Eq. ( 15).We discretize the spatial domain of Eq. ( 15) into 50 uniform grids with x i " pi ´1 2 q∆x, ∆x " 0.02, 1 ď i ď 50. Let uptq " pupx 1 , tq, ¨¨¨, upx 50 , tqq, vptq " pvpx 1 , tq, ¨¨¨, vpx 50 , tqq, then puptq, vptqq P R 100 are treated as the microscopic states. After approximating the spatial derivatives in Eq. ( 15) with the finite difference method, we consider the semi-discrete equation which is an N " 100 dimensional ODE as our microscopic evolution law (see Appendix B.1).We choose the macroscopic observable of interest to be z ˚" pū, vq, the spatial average of the predator and the prey's population:ū " 1 50EQUATIONWe find another 2 closure variables using the autoencoder, then the total dimension of the latent space z is 4. We choose D to be the trajectory distribution of the state x. For the data generation with partial forces, given x we randomly choose forces on 100 ¨p microscopic coordinate. For example, if p " 1{5, then for each configuration, forces on 20 coordinates are calculated for training. Results Fig. 2 shows the results on the test dataset which consists of 100 trajectories (for more detials, see Appendix Table 2 ). For models trained with partial microscopic forces, we report the equivalent number of training data with full forces throughout the paper. For example, for a model trained with L x,p pp " 1{5q on 3 ˆ10 3 data, we will report the number of training data to be 3 ˆ10 3 ˆ0.2 " 600. The test error is defined to be the mean relative error of the macroscopic observables between the ground truth and the predicted trajectories as in Appendix Eq. ( 43).First, we observe that the mean relative error of all the models is around 10 ´4 when the number of training data is large enough. This tells us that training with L x,p is correct and accurate, which is consistent with Theorem 2. The predicted trajectories fit quite well with the ground truth trajectories (see Appendix Fig. 6 and Fig. 7 ).We can conclude from Fig. 2 that, under the same number of training data, L x,p with smaller p p1{4, 1{5q performs better. Similarly, to achieve the same performance, L x,p with smaller p requires less training data. We set the error tolerance to be e tol " 3 ˆ10 ´3 and investigate how much training data is required to reach the error tolerance. In Fig. 2 the x-coordinate of the intersection point between the black dashed line and the other curves indicates the minimum data size required If we arrange each model according to their minimal required training data, then L x,p pp " 1{5q « L x,p pp " 1{4q ă L x,p pp " 1{2q ă L x,p pp " 3{4q ă L x . Model trained with L x,p pp " 1{4, 1{5q requires less data to reach e tol , or equivalently, less force computations. This validates the force computation efficiency of our method. One explanation could be that there are many redundant information in the forces acting on all the microscopic coordinates. By using partial microscopic forces, L x,p can explore more configurations x given the same size of training data, thus can make use of more useful information. Another observation from Fig. 2 is that as the training data size increases, the gap between models trained with different p narrows down. This is because as more data are provided, these data can contain almost all the information of the Predator-Prey system, thus more information will not lead to significant improvement.Having validated the correctness and force computation efficiency of our method, we are ready to apply our method to a variety of latent structures. We will tackle the Lennard-Jones system in this subsection, and validate the robustness of our method to different latent model structures. Three latent model structures are considered: MLP, OnsagerNet (Yu et al., 2021) , GFINNs (Zhang et al., 2022) . Both the OnsagerNet and GFINNs endow the latent dynamical model with certain thermodynamic structure to ensure stability and interpretability. The specific implementations of these two models are slightly different.The Lennard-Jones system is widely used in molecular simulation to study phase transition, crystallization and macroscopic properties of a system (Hansen & Verlet, 1969; Bengtzelius, 1986; Lin et al., 2003; Luo et al., 2004) . The Lennard-Jones potential describes the interaction between two atoms i and j through the potential of the following form:EQUATIONTable 1 : Summary of the results on each system. Results of the Predator-Prey and Lennard-Jones (small) system are taken from Section 5.1, Section 5.2. For each system, L z and L x,p are trained with the same size of data. All the results in this experiment will be shown in reduced Lennard-Jones units. We consider a threedimensional Lennard-Jones fluid with N atoms " 800 atoms of the same type. The microscopic state consists of the positions and velocities of the 800 atoms, thus the microscopic dimension is N " 4800. We simulate the Lennard-Jones system under NVE ensemble using LAMMPS (Thompson et al., 2022) .We choose the instantaneous temperature (T ) as our macroscopic observables:T " 2 3pN atoms ´1q ˆNatoms ÿ i"1 m i v 2 i 2 (18)here v i is the velocity of the i-th atom and m i " 1. We find another 31 closure variables using the autoencoder, then the latent dimension is 32. In our experiment, we also adopt the trajectory distribution of microscopic statex for D. For data generation with partial forces, we choose p " 1{16. For each x we randomly choose 50 atoms for force computations.Results All the models are trained with the same size of data. Fig. 3 shows the test error on 10 test trajectories. The test errors of using L x,p pp " 1{16q are relatively small (" 10 ´3), which validates the accuracy of our model on the Lennard-Jones system. It is easy to observe from Fig. 3 that for all the latent model structures, models trained with L x,p can always outperform those trained with L z . This validates L x,p is robust over different latent model structures.We have already validated the accuracy and force computation efficiency of L x,p on the Predatory-Prey system and the Lennard-Jones system, but their microscopic dimension is still not big enough.In this subsection, we focus on the robustness of our method to different systems including those with much larger microscopic dimension. We will consider two large systems: the Allen-Cahn system and a larger Lennard-Jones system with 51200 atoms.The Allen-Cahn equation is widely used to model the phase transition process in binary mixtures (Allen & Cahn, 1979; Del Pino et al., 2008; Shen & Yang, 2010; Kim et al., 2021; Yang et al., 2023) . We consider the 2-dimensional Allen-Cahn equation with zero Neumann boundary condition on a bounded domain:EQUATIONwhere vp´1 ď v ď 1q denotes the difference of the concentration of the two phases. F pvq is usually chosen to be the double potential taking the form of F pvq " 1 4 pv 2 ´1q 2 . The Allen-Cahn equation is the L 2 gradient flow of the free energy functional Epvq P R in Eq. ( 20) (Bartels, 2015) .EQUATIONThe free energy functional z ˚" Epvq is a macroscopic observable of wide interest, hence we choose Epvq as our target macroscopic observable.The spatial domain is discretized into 200 ˆ200 grids, then the dimension of the microscopic system is N " 40000. We find another 31 closure variables using the autoencoder hence the total dimension of the macroscopic system is 32, which is much smaller compared to the dimension of the microscopic system. We consider D to be the trajectory distribution of x. We choose p " 1{25, each time the forces on 1600 grids are calculated for training L x,p .Lennard-Jones System (large) To further demonstrate the capacity of our method, we scale up the Lennard-Jones system in Section 5.2 to encompass 51200 atoms, then N " 307200. The size of the simulation box is increased from 10 ˆ10 ˆ10 to 40 ˆ40 ˆ40 to keep the density unchanged.Results For a summary of the experiments and the results, we refer the reader to Table 1 . Note that for the Lennard-Jones system (large), we still use the forces on 50 atoms for training, thus p " 1{1024. For the training of L x,p pp " 1{1024q, only 5000 configurations with partial forces are used due to memory limit, which is equivalent to 5000{1024 « 5 training data with forces on all the atoms. Obviously, training data with size 5 is way too small, hence the results of L z are not reported for this system.From the results shown in Table 1 , one can observe that for a variety of problems, including those modeled by partial differential equations or molecular dynamics simulations, L x,p can always outperform L z . This shows the robustness of our method to a variety of systems. Moreover, the success of our method on the Lennard-Jones system (large) demonstrates the ability and efficiency of our method when scaled to very large systems. We also compare the number of force computations that are required for models trained with L z and L x,p to reach test error e tol " 3 ˆ10 ´3. Fig. 4 shows the results on the Lennard-Jones system with different sizes. Lennard-Jones system with 800, 2700, 6400, 21600 atoms are considered and the density is 0.8 for all the systems. The number of force computations here refers to the total number of forces on atoms that are used. For example, if L x,p uses 100 configurations to train, and for each configuration, forces on 50 atoms are calculated, then the number of force computations would be 100 ˆ50 " 5000.From Fig. 4 we can observe that as the system size increases, the number of force computations required by L z continues increasing. In the experiment, we find that the number of training data does not change a lot. The increase in the number of force computations is mainly due to the system size increases, then for each configuration, more force computations are required. However, the number of force computations even decreases a bit. One possible explanation is that, as the system size increases, the dynamics become less fluctuating and are easy to learn. 
conclusion: 
 We present a framework for modeling the dynamics of macroscopic observables from partial computation of microscopic forces. We theoretically and experimentally demonstrate the accuracy, force computation efficiency, and robustness of our method through different problems. Finally, we apply our method to a very large Lennard-Jones system which contains 51200 atoms.While our method can learn the macroscopic dynamics from partial computation of microscopic forces, it relies on the sparsity assumption. For systems that do not satisfy the sparsity assumption, the calculation of partial computation of microscopic forces is not efficient, thus it is not beneficial to learn from partial forces computation. For example, in the McKean-Vlasov system, the force on each microscopic coordinate depends on the collective behavior of all the other coordinates (Méléard, 1996) .Another limitation is the structure of the autoencoder. For particle systems such as the Lennard-Jones system, an ideal encoder should be permutation-invariant. Currently, we use MLP for the encoder, which can be improved. Additionally, our method assumes the microscopic state to be sampled from a distribution D. We choose D to be trajectory distribution in the experiments, but in reality trajectory distribution of large systems may be impossible to obtain. Active learning is commonly applied to efficiently select microscopic configurations for training (Ang et al., 2021; Zhang et al., 2019; Farache et al., 2022; Kulichenko et al., 2023; Duschatko et al., 2024) . It is of interest to combine active learning and our proposed method to overcome the difficulty of the choice of D. Furthermore, our method assumes the microscopic dynamics to be deterministic, another future direction could be generalizing our method to stochastic systems.(21)In our method, we map the training procedure from the macroscopic coordinates to microscopic coordinates and use partial computation of the microscopic forces to train. We treat L z pθq as baseline and give detailed theoretical analysis of the possible error introduced by using L x,p pθq. The error can be controlled by two parts: (i) Convert loss from z space to x space, i.e., the error introduced by using L x pθq (ii) Use partial labels, i.e., the error between L x,p pθq and L x pθq.We will analyze the first part and the second part of the error accordingly (Appendix A). More experimental details are provided in Appendix B.Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?Answer: [Yes] Justification: Details of the experiment setting is provided in Section 5.Guidelines:• The answer NA means that the paper does not include experiments.• The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them. • The full details can be provided either with the code, in appendix, or as supplemental material. 
experiment statistical significance: 
 Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?Answer: [Yes] Justification: We report the error bar for all the experiments.Guidelines:• The answer NA means that the paper does not include experiments.• The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper. • The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions). • The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors).• It should be clear whether the error bar is the standard deviation or the standard error of the mean. • It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified. • For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates). • If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text. 
experiments compute resources: 
 Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?Answer: [Yes] Justification: We mention in Appendix C.Guidelines:• The answer NA means that the paper does not include experiments.• The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage. • The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute. • The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper). 
code of ethics: 
 Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https://neurips.cc/public/EthicsGuidelines?Answer: [Yes] Justification: Our research conform with the NeurIPS Code of Ethics.Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 
broader impacts: 
 Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?Answer: [NA] Justification: no societal impacts Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.• The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: no such risks Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?Answer: [Yes] Justification: Existing assets are properly cited.Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. 
 
 Appendix: A.1 Proof of Theorem 1Theorem. Assume for any x " D, the eigenvalues of φ 1 pxqφ 1 pxq T are lower bounded by b 1 and upper bounded by b 2 , 0 ă b 1 ď b 2 . Then:EQUATIONhere C does not depend on θ hence does not affect the optimization.Proof. We can write φ 1 pxq in the following form by leveraging singular value decomposition:EQUATIONwhere φ 1 pxq P R dˆN , Σ P R dˆN is a rectangular diagonal matrix, U P R dˆd and V P R N ˆN are two orthogonal matrices, d ! N . U, Σ, V actually depends on x but for simplicity we omit the dependence in notation. During the training of the autoencoder, we enforce φ 1 pxq to have full row rank, then the diagonal items are all nonzero, i.e., Σ i,i " λ i pxq ‰ 0, i " 1, ¨¨¨, d.We define Σ : P R N ˆd to be the rectangular diagonal matrix such that the diagonal items are pΣ : q i,i " λ ´1 i pxq ‰ 0, i " 1, ¨¨¨, d and all the remaining items are zero. Actually Σ : is the Moore-Penrose inverse of Σ.Then pφ 1 pxqq : can be calculated by:pφ 1 pxqq : " VΣ : U T (24)If we denote the i-th column of V by v i , the i-th column of U by u i , then V " pv 1 , ¨¨¨, v N q, U " pu 1 , ¨¨¨, u d q, and we can rewrite L z pθq and L x pθq:L z pθq " 1 K ř px,yqPD K ∥φ 1 pxqy ´gθ pzq∥ 2 2 " 1 K ř px,yqPD K ∥UΣV T y ´gθ pzq∥ 2 2 " 1 K ř px,yqPD K ∥UΣV T y ´UU T g θ pzq∥ 2 2 " 1 K ř px,yqPD K pΣV T y ´UT g θ pzqq T U T UpΣV T y ´UT g θ pzqq " 1 K ř px,yqPD K ∥ΣV T y ´UT g θ pzq∥ 2 2 " 1 K ř d i"1 ř px,yqPD K ∥λ i pxqv T i y ´uT i g θ pzq∥ 2 2 (25) L x pθq " 1 K ř px,yqPD K ∥y ´pφ 1 pxqq : g θ pzq∥ 2 2 " 1 K ř px,yqPD K ∥y ´VΣ : U T g θ pzq∥ 2 2 " 1 K ř px,yqPD K ∥VV T y ´VΣ : U T g θ pzq∥ 2 2 q " 1 K ř px,yqPD K pV T y ´Σ: U T g θ pzqq T V T VpV T y ´Σ: U T g θ pzq " 1 K ř px,yqPD K ∥V T y ´Σ: U T g θ pzq∥ 2 2 " 1 K ř d i"1 ř px,yqPD K ∥v T i y ´λ´1 i pxqu T i g θ pzq∥ 2 2 `1 K ř N i"d`1 ř px,yqPD K ∥v T i y∥ 2 2 " 1 K ř d i"1 ř px,yqPD K λ ´2 i pxq∥λ i pxqv T i y ´uT i g θ pzq∥ 2 2 `1 K ř N i"d`1 ř px,yqPD K ∥v T i y∥ 2 2 (26) We define: Lx pθq " 1 K ř d i"1 ř px,yqPD K λ ´2 i pxq∥λ i pxqv T i y ´uT i g θ pzq∥ 2 2 C " ´1 K ř N i"d`1 ř px,yqPD K ∥v T i y∥ 2 2 (27)then L x pθq " Lx pθq ´C, C does not depend on θ and:EQUATION)Comparing L z pθq and Lx pθq, we observe that the only difference between L z pθq and Lx pθq is that for every term ∥λ i pxqv T i f pxq ´uT i g θ pzq∥ 2 2 , there is a constant λ ´2 i pxq multiplied to it. Hence Lx pθq is a weighted version of L z pθq. Note that if the eigenvalues of φ 1 pxq are λ i pxq, i " 1, ¨¨¨, d, then the eigenvalues of φ 1 pxqφ 1 pxq T are λ 2 i pxq, i " 1, ¨¨¨, d. Since the eigenvalues of φ 1 pxqφ 1 pxq T are lower bounded by b 1 ą 0 and upper bounded by b 2 , i.e., @x P X , 0 ă b 1 ď λ 2 i pxq ď b 2 , then:EQUATIONor equivalently, b 1 pL x pθq `Cq ď L z pθq ď b 2 pL x pθq `Cq (30)By minimizing Lx pθq, we are actually narrowing the region of L z pθq. Hence we want b 1 and b 2 to be as close as possible. In the extreme case where b 1 " b 2 , minimizing Lx pθq is just equivalent to minimizing L z pθq. Another observation is that Lx pθq is a weighted version of L x pθq, and if there exists i such that λ i pxq 2 is too small compared to the others, the weighted sum will be dominated by it in Eq. ( 27).The above insights guide us to constrain the condition number of φ 1 pxqφ 1 pxq T during the training of autoencoder, i.e., we require φ 1 pxqφ 1 pxq T to be well-conditioned through Eq. ( 6).A.2 Proof of Eq. ( 13)EQUATIONNote that θ ˚P arg min θ Lx pθq P arg min θ Lx,p pθq, then 0 ď Lx pθ K,p q ´L x pθ ˚q " Lx,p pθ K,p q ´Lx,p pθ K,p q loooooooooooooomoooooooooooooon a.s.`Lx,p pθ K,p q ´Lx,p pθ ˚q looooooooooooomooooooooooooon ď0 `Lx,p pθ ˚q ´L x,p pθ ˚q loooooooooooomoooooooooooon a.s.(39) thus Lx,p pθ K,p q ´L x,p pθ ˚q a.s.ÝÑ 0. The Neumann boundary condition is commonly used in mathematical models of ecosystems, restricting any movement of the species out of the boundary.We approximate the spatial derivatives in Eq. ( 15) with finite difference method:EQUATION)Let h u pu, vq " up1 ´u ´vq, h v pu, vq " avpu ´bq, and h u pu, vq, h v pu, vq denote the element-wise application of h u , h v to each upx i , tq, vpx i , tq,EQUATIONhere A P R 50ˆ50 is a matrix defined according to Eq. ( 42). Hence the predator-prey system after spatial discretization can be written in the form of Eq. (1).In our experiment, we choose a " 3, b " 0.4, λ " 0. The training parameter set T train of pairs pµ, σq are sampled uniformly from r0, 0.2s ˆr0.4, 0.6s. For testing, T test is also sampled uniformly from r0, 0.2s ˆr0.4, 0.6s , but with a different random seed from T train . The mean relative error is defined as :EQUATIONhere we use zp¨; µ, σq to denote the dependency of the solution on the initial condition. The microscopic equation is solved with a uniform time step ∆t " 0.01 from t " 0 to t " 30 using the Euler method. We subsampled every tenth snapshot for training. During testing, the microscopic evolution equation is solved with the same ∆t " 0.01 using Runge-Kutta 4-order RK4 solver. Then we encode the microscopic trajectories to obtain the ground truth latent trajectories. The predicted latent trajectories are obtained by encoding the initial microscopic state first, then solved using RK4 solver with ∆t " 0.1, ∆ " 0.5 on r0, 30s. Fig. 6 and Fig. 7 show the true and predicted trajectories. Figure 7 : Latent trajectories with initial condition µ " 0.02, σ " 0.52 and ∆t " 0.5 in the Predator-Prey system.In our experiment, we consider the initial condition of a torus (Kim et al., 2021) :vpx, y, 0q " ´1 `tanhp r 1 ´dpx, yq ? 2ϵ q ´tanhp r 2 ´dpx, yq ? 2ϵ q (44)here dpx, yq " a px ´0.5q 2 `py ´0.5q 2 , r 1 P r0.3, 0.4s is the circumscribed circle radius and r 2 P r0.1, 0.15s is the inscribed circle radius. The initial condition is visualized in 5 (b).The free energy in Eq. ( 20) also tends to decrease with time, following the energy dissipation law in Eq. ( 45). Then the minimization of the free energy drives the evolution of the system towards equilibrium.BEpvq Bt " ´żµ ∥B t v∥ 2 2 dx dy (45)Data Generation For both training and testing, the microscopic evolution law is solved using RK4 method with ∆t " 1{N " 2.5 ˆ10 ´5 from t " 0 to t " minpt f , 1q.Here t f is the time when the Allen-Cahn system reaches the equilibrium. We subsample every hundredth snapshots for training. We choose ϵ in Eq. ( 19) and Eq. ( 44) to be 10ˆ200 2 ? 2 tanh ´1p0.9q as in (Kim et al., 2021) . For testing, 50 parameter points pr 1 , r 2 q are chosen uniformly from r0.3, 0.4s ˆr0.1, 0.15s, and we report the mean relative error. The test parameter set T test contains 50 parameter points pr 1 , r 2 q are chosen uniformly from r0.3, 0.4s ˆr0.1, 0.15s, but are sampled with a different random seed.We consider Lennard-Jones systems containing different atoms in this paper: N atoms " 800, 2700, 6400, 21600, 51200. We use periodic boundary conditions and fix the density to be 0.8, then the corresponding box side lengths are 10, 15, 20, 30, 40. We simulate the Lennard-Jones system under the NVE ensemble using the LAMMPS (Thompson et al., 2022) . In our experiment, ϵ ij " σ ij " 1, @i, j , and r cut " 2.5. The integration step is 0.001 and each trajectory is integrated for 250 steps. We sample the initial temperature randomly from r0.5, 1.5s. The initial velocities are then sampled from the Maxwell-Boltzmann distribution. For each system, the initial configuration has the same atom positions and velocity direction. For testing, the initial temperatures are also randomly sampled from r0.5, 1.5s but with a different random seed to the training data.All the experiments are run on a single NVIDIA GeForce RTX 3090 GPU. For all the experiments, we use the multilayer perceptron (MLP) for both the encoder and the decoder. The autoencoders are trained with L AE in Eq. ( 7). The condition number is the maximal eigenvalue λ max divided by the minimal eigenvalue λ min of φ 1 pxqφ 1 pxq T : κ `φ1 pxqφ 1 pxq T ˘" |λ max pφ 1 pxqφ 1 pxq T q| |λ min pφ 1 pxqφ 1 pxq T q| ě 1Since φ 1 pxq P R dˆN , φ 1 pxqφ 1 pxq T P R dˆd , and d is small, the condition number κpφ 1 pxqφ 1 pxq T q can be calculated efficiently. In our experiments, we calculate λ max and λ min with torch.linalg.svd. To better compare L x and L x,p , once finish the training of the autoencoders, we freeze them and use the encoder for macroscopic dynamics identification. For the macroscopic dynamics identification, MLP and GFINNs are used for the latent model in Section 5.2. For the rest experiments, we adopt the structure of OnsagerNet to enhance the stability for latent dynamics prediction for g θ (Yu et al., 2021) .D.1 Loss Curve of L z , L x , L x,pTo give the readers a better idea of the behaviors of the loss L z , L x , L x,p pp " 1{4q trained on the same number of training data. Fig. 8 shows the training and test loss curve of different training metrics. Note that the training metrics of these models are different, but they are tested with the same metric Eq. ( 43).To evaluate the influence of the hyperparameter λ cond on the performance of the loss L x,p , we conduct experiments with different values of λ cond and show the test error in Table 3 .From Table 3 we can observe when λ cond increases from 0 to 10 ´6, the test error gradually decrease. When λ cond further increases from 10 ´6 to 10 ´2, the test error gradually decrease. Among all the λ cond that we tried, the test error has the minimal value when λ cond " 10 ´6.Theoretically, if λ cond is too low, since L AE " L rec `λcond L cond , there may not have enough constraint on L cond and the condition number of φ 1 pxqφ 1 pxq T may be very large. By Theorem 1, • Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted.38th Conference on Neural Information Processing Systems (NeurIPS 2024). 


Decision: The reviewers find this work studying an important problem and encourage the authors to incorporate the suggestions and developments during the review process.

Reviews: review: paper_summary: paper_summary: The authors describe a method to efficiently obtain aggregate information about forces acting on all particles in a system, in an effort to compute dynamics of macroscopic (aggregated) quantities. The key idea described in the paper is to sub-sample the particles to a small set, and only compute the forces on the small set (instead of computing all forces). The forces are then used to train a vector field of a latent space variable (through the chain rule). The latent space is obtained by training an auto-encoder on the microscopic states, independent of the force computations. The authors demonstrate the efficiency on multiple numerical systems, and include theoretical insights into the loss functions on micro- and macro (latent) states.
strength_weakness: strength_weakness: 1) The methods proposed in the paper are (a) sub-sampling particle indices to reduce number of force computations and (b) using a standard auto-encoder network to obtain a latent space. The former idea - especially when combined with the chain rule for the macroscopic dynamics - is noteworthy. Still, the general advancement of the state of the art is not large enough to warrant acceptance, in my opinion. In particular, the choice of auto-encoder is not suitable for the given problem (see question 6). 

2) Some parts of the literature are not cited.

[A] l.65: reduced order modeling has not started in 2020 with deep learning methods. Appropriate literature from the decades before should be cited. A review can be found here:

[A1] Schilders, Wilhelmus H.A., Joost Rommes, and Henk A. van der Vorst, eds. Model Order Reduction: Theory, Research Aspects and Applications. Mathematics in Industry. Springer Berlin Heidelberg, 2008. https://doi.org/10.1007/978-3-540-78841-6.

[B] l.55: similarly, "learning from partial observations" is a decades (if not centuries) old problem. Key mathematical contributions are from Ruelle and Takens over 40 years ago, refined by Yorke:

[B1] Takens, Floris. “Detecting Strange Attractors in Turbulence.” Lecture Notes in Mathematics, 1981, 366–81. https://doi.org/10.1007/bfb0091924.
[B2] Ruelle, David, and Floris Takens. “On the Nature of Turbulence.” Commun. Math. Phys. 20, no. 3 (September 1971): 167–92. https://doi.org/10.1007/bf01646553.
[B3] Sauer, Tim, James A. Yorke, and Martin Casdagli. “Embedology.” Journal of Statistical Physics 65, no. 3 (1991): 579–616. https://doi.org/10.1007/BF01053745.

[C] In general, the idea of "sparse sampling" of the microscopic dynamics has also been published before. This must not just be cited but critically compared to, especially [C2].

[C1] Samaey, Giovanni, Ioannis G. Kevrekidis, and Dirk Roose. “Patch Dynamics with Buffers for Homogenization Problems.” Journal of Computational Physics 213, no. 1 (March 2006): 264–87. https://doi.org/10.1016/j.jcp.2005.08.010.
[C2] Liu, Ping, Giovanni Samaey, C. William Gear, and Ioannis G. Kevrekidis. “On the Acceleration of Spatially Distributed Agent-Based Computations: A Patch Dynamics Scheme.” Applied Numerical Mathematics 92 (June 2015): 54–69. https://doi.org/10.1016/j.apnum.2014.12.007.

[D] The idea of creating a latent space with an auto-encoder is of course also not new. The authors cite a few papers in l.143, but should do so in section 4.1.


3) l141: it is misleading to write "$g_\theta = \phi'(x)f(x)$". g depends on z, not x (correct?). It is not clear if the authors mean "$g_\theta(z) = \phi'(\phi^{-1}(z))f(\phi^{-1}(z))$, i.e., as a combination of $\phi'$ (encoder Jacobian), $\phi^{-1}$ (decoder) and $f$. At this point in the paper, it could also be that training of g happens directly on dynamics of z and avoid $\phi$. The latter is probably the case, looking at equation 9, but it is not clear from the text.


Minor:
1) l74: "where" instead of "here", and it is $f(x)\in\mathbb{R}^N$, not $f\in\mathbb{R}^N$ (the former is a vector, the latter a function).
2) l157: "we constrict the condition number..." (I think "constrain" is meant?) is not ideal wording; because the loss in eq. 6 does not constrain the condition number at all, it just penalizes it.
questions: questions: 1) What could be done with quantities that do not rely on individual molecule evaluations (e.g. force), but on integral quantities (e.g. "density")?
2) I assume the authors require all particles for the force computation because even if a particle is not chosen for force computaiton, it still needs to be present because other particles require its position for "their" force computation? This would also answer question 1.
3) What is the core difference between reference C2 and the approach by the authors?
4) l159: "we get rid of the matrix-vector product...": why is this beneficial? The new loss also has a matrix vector product with phi' and g.
5) Equation 2 is not precise enough (and contains two typos: the parenthesis and no period at the end). What is the "almost equivalent" here? What is assumed, precisely? How accurate does the approximate force need to be?
6) The autoencoder described in section 4.1 seems to have a very simple architecture. How is it possible to train it for such high-dimensional microscopic input states, if the encoder is not permutation invariant? For example, if the model is trained using a fixed index for each molecule, and then the indices are shuffled at random, the input to the autoencoder completely changes while the conformation (independent of the index) does not change at all. Since the autoencoder is not permutation invariant, it would need to learn all conformations including permutations, which is an incredibly high-dimensional space (not possible to sample in reasonable time).
limitations: limitations: The authors do not discuss the limitations of the auto-encoder they use to map to a latent space, in particular regarding permutation invariance of the molecular indices (see my question).
Societal concerns are not addressed at all; e.g. that speeding up force computations can speed up insights into materials that are harmful (not just beneficial).
paper_summary: paper_summary: The authors aims to use ML to compute macroscopic dynamics of a system from partially observed microscopic dynamics. The paper defines the macroscopic dynamics to be the lower dimensional latent space of an autoencoder that encodes the microscopic system. Given the autoencoder they train a macroscopic dynamics model using the partially observed forces of the microscopic dynamics that are projected to the macroscopic space via the encoder. Experiments are carried out on toy systems and lennard jones potential systems that are scaled to a large number of particles.
strength_weakness: strength_weakness: 1.  Experiment system choices: The authors carry out experiments on small toy systems and large lennard jones systems of the same particles. It seems to me that the problem should become significantly harder with the "homogeneity" of the modeled system decreasing. It further seems to me that this would be the case in real scientific applications where we potentially have different atom types and multiple forces beyond lennard jones potential forces. Is this assessment correct? It would be great if you could put into perspective how close the experiments you carry out and the systems you choose are to actual scientific applications of interest. Thanks!
2. Presentation: you briefly mention closure modeling in the related work. Then in the methods section you say you follow closure modeling and define quantities such that a closed system arises. At this point I ask myself, what does closed mean? I hoped it becomes clear later in the paper but find there to be no sufficient explanation. Could you please let me know how you train the autoencoder in such a manner that the latent space captures specific desired macroscopic dynamics?
Minor:
1. Presentation: Motivation: What are some compelling application examples in which learning from partially observed microscopic dynamics is valuable?
questions: questions: 1. What if the particles to observe are not chosen uniformly at random? Is that also a common scenario where we do not have a uniform subsampling but our observations are biased by the specific subregion that we observe.
limitations: limitations: The paper prominently discusses two limitations (sparsity assumption and sampled distribution of microscopic dynamics) which are significant and put the paper into perspective of the field of work. I do not see further meaningful limitations that should be discussed apart from the possible limitation in terms of realistic evaluation which I mention in the weaknesses section.
paper_summary: paper_summary: This paper introduces an efficient framework for learning macroscopic dynamics of complex systems. The training data consists of microscopic configurations with only partial dynamics, i.e., time derivatives of only a small subset of the microscopic variables. The paper shows how this information nevertheless allows unbiased estimation of dynamics at the level of low-dimensional macroscopic dynamics. Experiments show the data-efficiency of the scheme compared to learning on full microscopic dynamics.
strength_weakness: strength_weakness: * Although the theoretical results are nice and the experimental results empirically support them, I am still skeptical of the efficacy of $L_x$ as a proxy loss for $L_z$. To bolster the claims, I strongly recommend that authors present some empirical analysis on the spectrum of the Jacobian $\psi’$, or even better, plot loss curves over the course of a training run showing the behavior of $L_z$ versus $L_x$. This will give the reader a much better intuition for how these different loss functions behave.
* Without an understanding of how significant the final test errors are, it is not particularly insightful to merely show test errors after an arbitrarily selected amount of training data. Instead it would be nice to see scaling plots like in Figure 2 for the other systems in the analysis.
* Table 1 and Figure 2 are redundant and I would recommend eliminating the former.
* Theorem 2 could be further clarified. What is $\arg \min L_{x,p}$ since $L_{x,p}$ is also a random variable?
* Is it difficult or tricky to tune the weight of the conditioning number loss term? What happens if the weight on this term is not high enough? What if it is too high? Ablation analyses of these factors would improve the paper.
questions: questions: See above
limitations: limitations: Yes

Comments: title: A neat trick indeed! Why reject because permutation sensitive architecture in toy experiment?
comment: Thank you very much for your valuable suggestions that helped us to improve our paper.  

We will revise our descriptions and terminology according to your suggestion to make our paper easier to understand: 

-  For the 'more comments part', given the same number of data points, the larger the $p$ is, the more force computation we will have, and theoretically the error will be smaller. Here the error is controlled through the force computation.Small condition number of $\varphi^{\prime}(x) \varphi^{\prime}(x)^T$ controls the error through the encoder.  We will try to explain this in more detail in the paper.
-  We agree that calling $f$ as a vector field is more accurate. We will show the reason why we call $f$ the 'force' in the paper. 
As shown in Eq 3 and Eq 4, the Newton's law of motion in molecular system is :
$$\dot{r_i} = v_i$$
$$\dot{v_i} = \frac{1}{m_i}F_i$$
here $f_i$ is the position, $v_i$ is the velocity, $m_i$ is the mass and $F_i$ is the force on the $i-th $ atom. In this case the microscopic state $x = (r_i, v_i)$, and its dynamics $f = (v_i, \frac{1}{m_i} F_i)$. $\frac{1}{m_i} F_i$ is part of the $f$ in molecular dynamics, then we call $f$ 'force' in other cases such as the PDE case even though $f$ is not a 'force'. Thank you for pointing out that calling $f$ 'force' is difficult to understand. For notation simplicity, how about adding an explanation for the reason why we call $f$ the 'force' after Eq 1 in the paper as the following?  

*Note that the $f(x(t))$ is indeed a vector field. In molecular dynamics, ... (The above-mentioned molecular dynamics example)...for notation simplicity, we will call $f$ the force in the rest of the paper.*

- You are correct that the 'black-box solver' only computes the 'forces', rather than the trajectory. We will modify the terminology to 'microscopic force calculator' or 'vector field calculator'. Our method does not need the trajectories of the microscopic system, but we assume we know a distribution of the microscopic state $x$. If the distribution is unknown, there are still several possible ways for our method to work, such as 1. using active learning to add the most informative configuration and 2. using cheaper coarse-grain or Markov chain Monte Carlo (MCMC) methods to sample an approximate distribution. 

- We will change Figure 1 to make it easier to understand. In the Data generation part of Figure 1, we will remove all the black arrows on the left of the 'black-box solver', and only preserve the pink arrows on the right of the 'black-box solver'.
Thank you for pushing me on this. I had a very hard time understanding the core idea of your method, also because there are several issues with the terminology used in the manuscript. I strongly urge you to revise most of the descriptions and terminology in the paper based on my feedback. I will raise my score to 5, but not higher, because I cannot read the revised version and the submitted version has too many issues with terminology and missing explanations for a higher score. I will outline what I understand now, regarding the setting, the main idea, and give hints on the terminology below.

## Main setting
The proposed method is designed to approximate a vector field ${g}$ on macroscopic variables - including some of interest, and some to guarantee a closed model. The former are given, the latter are the output of a trained encoder network. All macroscopic variables $z$ considered are functions of microscopic variables $x$, approximated by an encoder network $\phi$ s.t. $\phi(x)=z$.

## Main idea
Many types of macroscopic variables can be *estimated* from a small subset $I$ of the microscopic variables $x$ (e.g., the mean of 1000 numbers can be estimated by randomly sub-sampling only 100 of them and dividing by 100 instead of 1000). This is at the heart of many coarse-graining approaches in MD. The core idea of the authors - as I understand it now - is that randomly choosing microscopic variables and computing their individual effect on the macroscopic variables (with reweighting based on the sample size) can also be done between the vector fields $f$ and $g$. Instead of computing all entries of $f(x)$, the authors only compute entries $f_{I(x)}(x)$, but then reweight them so that it can be used to estimate $g$. I understand this in terms of coarse-graining: the fraction $p\in(0,1)$ of sampled particles are "copied" 1/p times (by reweighting), so that the total influence is still as before. Using identical copies of a subset, instead of all particles, introduces an error. That error is reduced based on the "averaging effects" introduced by computing the macroscopic variables $z$. As some of these macroscopic variables are *learned*, the map $\phi$ can even be adjusted so that the error is reduced.

## More comments
I would have liked to connect the error reduction through averaging with the requirement on "small condition number" of $\phi'(x)\phi'(x)^T$, maybe the authors can do this.
The fact that *no microscopic simulation* is required at all makes this method extremely appealing for MD applications. As soon as reasonable configurations of particles can be obtained, and the vector field (e.g. including forces) can be evaluated on a small sample, it is possible to approximate the macroscopic states (the encoder) and their dynamics.

## Terminology
1) The function $f$ in Equation 1 is not a force, but this term is used throughout the paper. It is a *vector field*, and if $x$ is a position, then it can be called a velocity. Only if $\ddot{x}=f(x)$ would $f$ be a force, but then all derivations lateron for $\dot{z}$ etc. would need to be adjusted. It is clear why in most settings used in the paper, the complication to compute the vector field is in the computation of the forces, but the notation is misleading and must be cleared up.

2) The fact that you call it a "solver" in many places, while what you actually mean is a subroutine that computes forces - and does not "solve" anything - was the main reason I did not understand the appeal of your method. In my opinion, a solver "solves" something, usually a differential equation, so that the output is a solution (e.g. a trajectory). In the paper, the "black-box solvers" in figure 1 and algorithm 1 are actually not solving the microscopic dynamics, they are just computing the "forces" (or, rather, "evaluate the vector field"). This is the main reason I now think the method is great (as I discuss above), but it was extremely difficult for me to understand because of the mismatch of terminology (solver vs. vector field evaluation).

3) Similar to 2, Figure 1 also did not help to mitigate my main concerns. On the left there are many black arrows indicating the vector field (or "force") on each particle. But your method does not even need any of this in the "data set", correct? The figure suggests that the "black-box solver" must compute all the black arrows, and then only a subset of them is used in training. This is how an actual solver for the microscopic dynamics would do it, too - because the microscopic state must be advanced over time, and an actual solver can only do that after computing all forces for all particles (maybe subsampling which neighbors are used in the process). For your method, none of that is relevant, which I now understand, because you do not even need any trajectories on the microscopic scale for training the macroscopic vector field.
Thank you for the suggestion! We will revise Section 4.1 to state this idea more explicitly and make it easier to understand.
Thank you! This explanation fully addresses my fundamental concern of "how do you make sure that the learned dynamics correspond to the desired dynamics" and allows me to fully assess the paper. I am not sure whether I am the only reader who needs this more explicitly stated, and I leave it up to your assessment whether or not it would be beneficial to change the text. 

I am changing my overall score from 6 to 7 and my confidence has increased through the clarifications of the authors.
Thanks for your insightful question. Let us address the question more clearly. 

In l 124, we mention '*we will use an autoencoder to find the closure $\hat{z} = \hat{\varphi}(x)$ to $z^{\ast}$ such that $z = (z^{\ast}, \hat{z})$ forms a closed system*'. 
Here $z^{\ast} = \varphi^{\ast}(x)$ is the macroscopic observables that we are interested in, and the dynamics of $z^{\ast}$ is the DESIRED macroscopic dynamics. 
The function $\varphi^{\ast}$ **is determined beforehand** and contains no trainable parameters. The closure $\hat{z} = \hat{\varphi}(x)$ is learned by the autoencoder. During the training of the autoencoder, $\varphi^{\ast}$ is fixed all the time, and only the parameters of $\hat{\varphi}$ are updated. 

For example, in the Lennard-Jones experiment, we choose the macroscopic observable $z^{\ast}$ as the Temperature $T$.  $z^{\ast} = \varphi^{\ast}(x)$ is defined through Eq 18:
$$
T = \frac{2}{3(N_{atoms}-1)} \times \sum_{i=1}^{N_{atoms}} \frac{m_iv_i^2}{2}
$$
In l 250, we mention we find another 31 closure variables using the autoencoder. 
We **directly concatenate** the macroscopic observable $z^{\ast}$ and the closure $\hat{z}$ to get the latent variable $z$:
$$z = (z^{\ast}, \hat{z})$$
Then the first dimension of $z$ is the desired macroscopic observable. 

The method mentioned above is common in closure modeling. For example,  the authors in [1] want to learn the stretching dynamics of polymer. Then they fix $z^{\ast}$ to be the length of the polymer and find another two-dimensional closure variable $\hat{z}$.  As for the interpretability, $z^{\ast}$ is interpretable since it is exactly the macroscopic observables we want. Usually, the closure variables $\hat{z}$ are not easy to interpret, since they are learned by the neural network. For simple systems like a nonlinear pendulum system, it may be possible to interpret $\hat{z}$ learned by the autoencoder [2, 3]. 

We hope that this explanation clarifies your question.


---
References:

[1] Chen, X., Soh, B.W., Ooi, ZE. et al. Constructing custom thermodynamics using deep learning. Nat Comput Sci 4, 66–85 (2024). https://doi.org/10.1038/s43588-023-00581-5


[2] Champion, K., Lusch, B., Kutz, J. N., Brunton, S. L. (2019). Data-driven discovery of coordinates and governing equations. Proceedings of the National Academy of Sciences, 116(45), 22445-22451.

[3] Evangelou, N., Giovanis, D. G., Kevrekidis, G. A., Pavliotis, G. A., & Kevrekidis, I. G. (2023). Machine Learning for the identification of phase-transitions in interacting agent-based systems. arXiv preprint arXiv:2310.19039.
Thank you for the reply. We wish to clarify certain aspects of our paper.
First, our paper is centered on machine learning for the following reasons.

- The key contribution of our paper is the novel training loss $L_{x,p}$. 
    By training with $L_{x,p}$, we map the training procedure on the macroscopic coordinates back to the microscopic coordinates, on which partial force computations can be used as stochastic estimation to update model parameters. From this perspective, our method is directly related to ML,
    as it resolves a key bottleneck in learning reduced macroscopic (closure) models for large systems.
- Our method enables learning macroscopic dynamics from partial microscopic observations using ML.
    In our paper, we validate our method in several experiments, not only including microscopic dynamics described by MD,
    but also PDE based microscopic models.
    Hence, our method applies not only to MD, but rather it addresses some common challenges in *AI for science*,
    namely learning macroscopic dynamics from data.

Next, we would like to clarify the distinction between our method and the sub-sampling techniques employed in the MD solvers you mentioned:
- The sub-sampling method employed in MD solvers sub-samples a subset of the interactions between particles. For example, the random-batch molecular dynamics (RBMD) in [1] sub-samples the non-bonded interaction for each integration step. The force on all the atoms is still calculated. However, during the calculation of the force on each atom, only a subset of the interaction with other particles is considered. In our paper, when we refer to 'partial computation of microscopic forces', we only calculate the forces on a subset of particles.
    In contrast, our method does not at any instant compute the forces on all particles, since this is not
    necessary to obtain a (stochastic) approximation of the effective force on the macroscopic variables.

- The sub-sampling method is applied to microscopic MD simulation. The computational cost scales with the system size $N$. Our method instead learns the macroscopic dynamics and simulates the dynamics in macroscale. 

--- 
References:

[1] Gao, W., Zhao, T., Guo, Y., Liang, J., Liu, H., Luo, M., ... Xu, Z. (2024). RBMD: A molecular dynamics package enabling to simulate 10 million all-atom particles in a single graphics processing unit. arXiv preprint arXiv:2407.09315.
As I see it, the main idea of the paper is to utilize a small sample of all microscopic forces $f$ from a particle system to learn the vector field of macroscopic dynamics $g$ on learned latent variables $z$. The latter part is not performed with the right tools (GNNs could work, indeed). I agree that that is not directly related to the main idea, and changing to another architecture could make the examples more general. However, the main idea is not related to ML, but rather to molecular dynamics (MD). In fact, any reasonably efficient MD solver would probably already subsample forces to compute the microscopic dynamics, especially in the cases that the authors mention (e.g. forces that decay quickly with distance). I still fail to see how the main idea advances the state of the art in ML. For a conference on MD, this contribution may be more suitable.
I agree with d4vx that microscopic space loss is interesting and a neat trick (if I gathered that opinion correctly)!

It may be a weakness of the paper (I think it is) that it only evaluates on rather toy tasks (one is quite a large scale but, as far as I understand, should still be considered toy because there is no heterogeneity).\
However, in the large scale toy task, the permutation equivariance seems very much beside the point - generalization along this axis is not what we are concerned in evaluating. There is no reason to believe that a GNN wouldn't work if that generalization is required. \
The paper's experiments are not meant to show that their implementation is optimal or an "advancement to networks" and addresses the concrete task (toy task) as well as possible, but rather provides evaluation and a sanity check of a neat new method. I think that can be considered interesting ML.
Sorry for the late response - I reread the paper and will be faster to respond from now on.
(and thanks for the careful responses and explanations) 

I think your answers that you put under these headings best relate to my main concern:

1. "Q: Could you please let me know how you train the autoencoder in such a manner that the latent space captures specific desired macroscopic dynamics?"
2. "Q: What are some compelling application examples in which learning from partially observed microscopic dynamics is valuable?"

I think your answer to the first question misses my concern a little bit - I should have explained my question better. I understand that you train an autoencoder and a model operating on its latent space to capture the dynamics of the latent space. 

The concern is: why would the latent space capture the DESIRED macroscopic dynamics. If it is closed, it will capture all dynamics - sure. Also, it will be lower dimensional and capture SOME macroscopic dynamics. But why would the latent dimensions necessarily correspond to any of the macroscopic observables that "we" as the practitioners were interested in?

For instance, considering your example that you gave for question 2 where we care about the viscosity of the system: what guarantees us that any of the latent variables correspond to viscosity? It seems to me that I am missing the entire point of the paper and that the latent representations remain uninterpretable and we cannot ensure that they correspond to certain macroscopic observables that we were interested in.
I thank the authors for their answers and clarification. The discussion of the literature is adequate. Still, my concern about major advancement of the field of networks (and ML in general) stands - enforced by the clarification that the auto-encoder indeed cannot deal with permutations of the individual atoms. I will keep the review score at 3.
I appreciate the authors' detailed response to my concerns and will raise the score to 7. However as I am not as familiar with existing literature on reduced order models as other reviewers, I will decrease my confidence to 3.