Reduced-Order Neural Operators: Learning Lagrangian Dynamics on Highly Sparse Graphs

Q1. Only explicit time integrator is considered. Please add motivation for this choice to paper.

A1. Indeed, the present work only considers explicit time integration. This is consistent with prior machine–learning–based time stepping methods. (See Learning to Simulate Complex Physics with Graph Networks). We agree that it is exciting to work on alternative time integration schemes, such as Runge-Kutta and implicit Euler, in the future.

Q2. Please clarify units of duration for Table 1.

A2. We have clarified these in the updated PDF. The time interval is 5e-3s per time-step.

Q3. Please add information about numerical scheme used in nclaw. Is it the same scheme used for training model?

A3. Yes, nclaw also uses an explicit time integration scheme. That said, the training data for the elastic examples (owl-shaped mesh from the LiCROM paper) uses an implicit integration scheme. Since our model works on both of these settings, we show that, our work is compatible with training data of any time integration scheme.

Q4. What about stability issues of numerical scheme used? How the time step is chosen?

A4. Consistent with other machine-learning-based methods (e.g., GNS, CROM), our approach remains stable even if we take a larger time step than the training data generated using classical numerical methods. However, as the reviewer pointed out, our approach also becomes unstable if too big of a time step is chosen. All our examples are generated with the same time step as the training data, whether they are generated using explicit (sand, nclaw dataset) or implicit scheme (elasticity, LiCROM dataset). We will further clarify in this paper and add a stability analysis (w.r.t. the time step size).

Q5. How does the developed approach compare with implicit numerical schemes in terms of stability and time-stepping requirements?

A5. As discussed in the previous answer, our approach can take the same time step size as the training data generated using implicit numerical schemes (elasticity, LiCROM dataset). It will be interesting for future work to explore the possibility of integrating implicit time integration with machine-learning-based force evaluation, which, to our knowledge, has not been done before. This will require differentiating the neural operator and graph neural network during runtime.

Major Revisions

We have revised the paper as suggested. However, we welcome any further feedback on any parts of the paper that require improvements in presentation.

1. The paper's novelty is either limited or not effectively communicated A. We have re-written the introduction section to better clarify our novelty and distinguish our method from prior approaches.

2. Sections 3, 4 - explanations of how they fit together, Section 5 - Additional explanation for experiments. A. We provided explanations at the start of sections 3 and 4, explaining how they relate to Figure 1. The introductory section further distinguishes how the parts explained in sections 3 and 4 differ from prior literature. We furthermore added additional details about the experiments in section 5.

3. To improve clarity, the authors should provide a more detailed rationale for each experiment, explaining the intentions behind the choices and the methodologies A. Each experiment has been updated with a paragraph explaining the experimental setup, evaluation metric, and the reasoning for the experiment.

Q1. Although the paper claims high fidelity across various initial conditions and resolutions, it does not address complex cases with extreme variability, such as materials undergoing phase transitions or highly turbulent flows. These more complex scenarios, which are generally more challenging to model with reduced-order methods, are not covered in the results presented.

A1. Indeed, our current approach can have issues with discontinuities and stress concentration (e.g., shear localization). We will add this to the discussion/limitation section of our paper. Increasing the sampling resolution can partially alleviate the issue, but tradeoffs between accuracy and the number of samples have to be made (See Figure 6). Importance sampling methods that increase sampling density in the singular geometry region can lead to more efficient results (See Contact-centric deformation learning, ACM SIGGRAPH 2022; Optimizing Cubature for Efficient Integration of Subspace Deformations, ACM SIGGRAPH 2008). The SPH-based Newtonian fluid datasets used in the paper exhibit chaotic behavior. We will release the dataset containing these trajectories, similar to what is presented in the GNS paper [Sanchez-Gonzales, 2020]. Modeling phase transitions is a consideration for future work which will solidify reduced-order modeling

Q2. What are the limitations of sparse sampling in terms of capturing fine-grained details? Can you provide information of the trade-offs between computational efficiency and the accuracy of capturing system details?

A2. We have updated Figure 6 in the new Pdf to contrast accuracy, GPU usage, computation-time against sparsity and graph radius. We observe a decrease from ~7GB to <1GB on 0.031x sampled graphs, while achieving a rollout MSE of the order of 1e-4.

Q3. How does the model perform when interpolating or extrapolating to new materials or conditions not represented in the training data? It is mentioned that these cases are challenging but no empirical evidence is provided.

A3. We used a large number of trajectories and showed that the model can handle unseen trajectories with unseen conditions, generated by random initial velocity vector generated by random seeds. We also showed that the model can handle multi-material scenarios shown in Table 1. However, the model is unable to generalize to unseen materials because each material is encoded within the encoder and for unseen material, the encoding would be undefined or would be approximated to the nearest known material.

Q4. Can you provide a systematic study on the scaling performance with varying graph size for the same dynamical system?

A4. Yes, this has been provided in the updated figure 6. We observed lower GPU usage, faster computation time and comparable performance with higher levels of sparsity, upto 0.031x the full-order point cloud (78k particles). Increasing sparsity beyond 0.031x leads to degradation of performance.

Q5. Are there specific physics-based scenarios where this reduced-order approach might be less effective?

A5. When the system is not bounded and the number of particles in the system is not fixed, the reduced-order representations may fail to capture newer information that was previously not in the system. An example could be a fluid flow system with particles constantly entering and leaving the system of vortices.

Q1. The graph construction is computationally demanding.

A1. We show in Table 4 that using Graph layers allows us to achieve the best performance since our inputs are irregular and neural operators such as FNO, GNOT, DINO etc. are designed to work on regular grids. GINO does not account for inter-particle interactions, which is crucial for Lagrangian Dynamics. In the absence of which, the models fail to capture spatial information, leading to suboptimal performance. We created computation benchmarks in Figure 6 of the updated PDF which shows how graph size affects the GPU usage and computation time as a function of graph radius and sampling percentage.

Q2. Handling discontinuities/limitations with farthest point sampling

A2. Indeed, our current approach can have issues with discontinuities and stress concentration (e.g., shear localization). Increasing the sampling resolution can partially alleviate the issue, but tradeoffs between accuracy and the number of samples have to be made (See Figure 6 ). Importance sampling methods that increase sampling density in the singular geometry region can lead to more efficient results (See Contact-centric deformation learning, ACM SIGGRAPH 2022; Optimizing Cubature for Efficient Integration of Subspace Deformations, ACM SIGGRAPH 2008).

Q3. Complicated Pipeline, too many hyper-parameters.

A3. We agree with the reviewers that this model has several hyperparameters. However, we have shown in the config files that, once tuned, many of these hyperparameters can be considered universal hyperparameters, that work across physical systems, with standard deviation of the noise being one of the few features that changes across them. We require only little system-specific tuning. Also, its flexibility and modularity can also be seen as an advantage of being able to tailor it to more complicated systems (informed choices of subsampling, neighborhood selection, edge information, smoothness prior in the decoder).

Q4. Complexity of Encoder (with interaction operator) , Decoder

A4. We would like to clarify that the order reduction is done through farthest point sampling and not using learnable techniques. The encoder is used to capture the local spatial interactions. This has been clarified in the updated Figure 1. However, the use of Neural Field as the encoder would speed up the process by removing the message passing. As such we performed an ablation study with NF-NOT-NF setup on the elasticity dataset, where the order-reduction is learnable through the NF encoder. However, we observed that NF was not as efficient as the Interaction Operator in capturing spatial dependencies and interactions. We present the results below:

NF-NOT-NF: 1-step MSE: 1.26e-6, RolloutMSE: 1.35 IO-NOT-IO: 1-step MSE: 5.07e-10, RolloutMSE: 4e-4

Q5. Handling self-contact

A5. The training data is generated by the material point method (MPM), which does not explicitly check collision and implicitly handles self-collision through a background grid, for both solids and fluids. We follow a similar approach as MPM as well as the GNS baseline by letting the neural network to implicitly these self contacts. This has been added to the discussion section. Better fine-grained self-contact sampling is an exciting future work direction (See Contact-centric deformation learning, ACM SIGGRAPH 2022). We will add it to the future work.

Q6. Figure 3: The figure is misleading. Why is the farthest point graph more dense than the full-order graph?

A6. We have updated Figure 3. to include N(odes), E(dges), R(adius) for each of the graphs. FPS based graph has more edges because it uses a larger radius than the Full-order system. The edges connect points that are farther away from each other. Delaunay graph on the other hand has a fixed number of edges, thus making it discretization-dependent. Having edges connecting points further away allows us to efficiently model long-range relationships, which due to the properties of graph layers, is retained if the density increases. However, importantly, if the density decreases, adding new edges by increasing the radius of the input during inference, allows the model to still infer the correct dynamics. This has been shown in Figure 6 (sparsity vs. accuracy). Furthermore, to justify that the speedup is not affected due to the addition of new edges, we show in Tables 9 and 10 that our setup outperforms other time-steppers in terms of speed even with more edges.

Thank you for your feedback.

Q1. In equation 4, how is the coefficient q derived in practice? Is it predicted by another neural network, or is it optimized directly via least-squares? If it is optimized online, then during the inference stage, it will require extra optimization for every reconstruction from reduced mesh to full mesh.

A1. Yes, it is optimized directly via a linear least squares. We follow equation 10 from https://arxiv.org/pdf/2310.15907. This boils down to solving a symmetric positive linear system using a single Cholesky factorization. Therefore there is no computationally expensive optimization involved (faster than a neural network evaluation), thereby introducing minimum overhead (See Table 9 upscale time, ~1e-4s). We have added this in the paper.

Q2. In table 4, several baselines are listed but there is no formal definition or a brief introduction of them, for example I could not find the definition of what is IPOT.

A2. We apologize for this oversight. These models have been defined in the updated Pdf, along with the relevant citations. GNN - Graph Neural Network, GAT - Graph Attention Network, GINO - Geometry Informed Neural Operator, GNOT - General Neural Operator Transformer, IPOT - Inducing Point Operator Transformer.

Q3. In table 3, what is the autoencoder? Is it a GNN or it's just a MLP? If it's a GNN then shouldn't it be permutation-equivariant?

A3. It is an MLP autoencoder, similar to the ones commonly used in standard reduced-order models. See [Model reduction of dynamical systems on nonlinear manifolds using deep convolutional autoencoders] by Lee et al. Our approach, on the other hand, guarantees permutation-equivariance and allows for evaluation at arbitrary locations in space. This has been clarified in the updated Pdf.

Q4. As shown in Figure 1, the reduction part does not contain any learnable parts but rather just rely on non-learnable sampling method, which can potentially result in information loss.

A4. Yes, information loss is possible. As we discussed in Figure 6 (left), there is a trade-off between accuracy and the number of points sampled. Indeed, it is an exciting future work to explore learnable methods to potentially reduce information loss. Furthermore, the “potential information loss” during training enhances the time-stepper model robustness and accelerates training. For inference, we can take any resolution—in particular, also omit subsampling.

Q5. When applying FNO as part of the latent dynamics model, how do you handle the empty voxels/regions?

A5. The FNO is applied to functions on an equidistant grid. The integral transform before and after the FNO allows us to change the discretization from a general point cloud, as done in GINO [Li et al., 2024]. In particular, function values on the latent grid originate from learned aggregations over neighborhoods with a fixed radius, i.e., having more neighbors if the particles are denser.

Q6. The paper compares with full-order model like GNS or other reduce-order model, how does the model compare to multi-grid model like Cao et al. [1]?

A6. The Interaction Operator presented in the paper can be converted to a multi-level graph model, which, as shown in "Multipole Graph Neural Operator for Parametric Partial Differential Equations, Li et al.", indeed behaves as a neural operator. It is an exciting future research direction to see how the improved time complexity of multi-grid operators can make learning Lagrangian dynamics faster and more efficient.

Q7. Is it necessary to do the sampling of full-order mesh at every timestep, and why not stay in the latent space with a fixed set of particles sampled in the beginning?

A7. During training, to increase robustness, the inputs are partitioned into slices of 6 time-steps, 5 for input and 1 for output. All 6 frames have the same set of particles. However, other slices may have a different set of particles. During inference, the model predicts a fixed set of particles sampled at the initial time-step over the entire temporal sequence. An additional feature supported by the model, is that the second (decoder) integral transform layer can infer particles that were not in the input graph. This theoretically allows us to have varying set of particles (between input and output) at each time-step. This feature was, however not tested because we saw no application of it in the scenarios shown in the paper.

Thank you for your feedback

Q1. Is there any reason why you didn't use existing 3D datasets, e.g., from Sanches-Gonzalez et al. (2020)? The only subset that is not there is a 3D multi-material one. Looking at the baseline results in your Table 4, GNS and the proposed model perform similarly. Comparability with the baseline numbers from the GNS paper would have been very useful.

A1. The NCLAW datasets were generated with complex geometries, allowing us to test discretization invariance and geometry invariance. 3D GNS datasets offer low flexibility for testing these invariance features since their 3D datasets do not have complex geometries.

Q2. You say on lines 351-352 that random sampling has a similar performance to FPS, but on L. 329, you say that you do use FPS after all. Having some experience with FPS, I know that it is a sequential process over the point cloud, with the number of iterations being the number of subsampled points. And this sequential nature becomes a bottleneck when working with large point clouds as the presented 3D ones. Why didn't you just use random sampling?

A2. We show in Table 13 that the two sampling strategies have comparable performance. While our point-clouds were small enough (<100k points), we used FPS for even distribution of points, on larger point-clouds, random sampling may be used without loss of accuracy.

Q3. What do you mean by POD and an autoencoder in Table 3? Do you explain this somewhere? I don't think any of these two are representative baselines. In addition, Table 3 has multiple identical lines, all of which could be summarized in one sentence. Please consider removing this table altogether.

A3. For reduced-order methods (ROM), POD and autoencoders are standard baselines. See [Model reduction of dynamical systems on nonlinear manifolds using deep convolutional autoencoders] by Lee et al. We have added more explantions and background citations for these two methods. We apologize for the identical lines. We have removed the table from the paper.

Q4. Which dataset do you use for Table 3?

A4. We have updated the paragraph to include the dataset used for computing the metrics. We used The Owl Dataset (Elasticity) for all the experiments shown in (now removed) Table 3.

Q5. Table 3: Did you consider discussing "ROM baselines" (L.431) and "Discretization-agnostic ROM" (L.514) next to each other, so that Table 3 is close to both?

A5. Thank you for the recommendation. We have updated it in the revised pdf.

Q6. You cite the DINo paper (Yin et al., 2023) and a few similar ones, but I didn't see a discussion of how your approach improves on them. To my knowledge, you could have used the DINo encoder and decoder on the velocity field (and Euler-integrate once), and you would only have had to add an operator transformer in between. Am I missing something?

A6. DINo works on grid based inputs. Therefore, it is necessary to include integral transform operation to convert from unstructured pointcloud to a structured grid, before passing it to DINo. We performed experiments with DINo as the processor on the elasticity dataset, in place of NOT, and observed similar performance as NOT, showing that it is a valid time-stepper.

IO-DINo-IO: 1-step MSE: 6.25e-10, Rollout MSE: 7e-4

Q7.

The following changes have been made in the updated pdf.

L. 310, "Ma et al. (2023)" -> "(Ma et al., 2023)" L. 412: please add some space between the captions of subplots (a) and (b). Currently, it is unclear that there are two subplots. Table 13: please reduce the font size of this table L. 1048: "our model is at least 5x faster than [GNNs]" doesn't match with the numbers in Tables 9 and 10. It is rather a 2-4x speedup. Please reformulate.</b>