Thank you for your insightful comments and constructive feedback. Below, we provide detailed responses to each point raised:

1. Quote Graph Laplacian pre-processing time Response: [Experiment: Precomputation Cost]:

While we provide a precomputation cost table for clarity, this data processing only occurs once and can be considered part of the dataset generation process.

The primary computational cost is in the eigen decomposition of the graph Laplacian, with a total precomputation time of about 3 minutes—negligible compared to training time. For the four datasets with static graph topologies, eigen decomposition is only needed at $t=0$ rather than at each time step.

For datasets with dynamic topologies, eigen decomposition is required at each time step but only needs to be done once, in batch, before training, and cached for reuse. This approach reduces on-the-fly calculation complexity, improving training efficiency. Importantly, inference does not require precomputed eigenvalues or eigenvectors, so inference speed remains unaffected.

2. Small Laplacian perturbations could lead to distinct eigenspaces. Could GFL hinder remeshing/generalization with varying graph connectivity? Response: Thank you for pointing this out. Our GFL remains robust under varying graph connectivity due to the Energy-Based Robustness: GFL operates on frequency-domain energy distributions ($E$) rather than relying directly on eigenvector precision. By focusing on mean energy across high- and low-energy components, it reduces sensitivity to small eigenvector perturbations caused by narrow eigenvalue gaps.

[Experiment: GFL Robustness Performance]:

3. Notable discrepancies exist between MGN results in your work and those in the original paper on the same datasets. Why? Response: We utilized the mgn_pytorch implementation to run MGN baseline, yielding the results shown in the table. There may be differences between our results and those in the original paper, possibly due to some variations between PyTorch and TensorFlow. Notice that we use RMSE(1E-2) to visualize and MGN use RMSE x$10^{-3}$, and FlagSimple's results seems substantially different because we used speed, whereas MGN used position. We note that other papers using MGN as a baseline[1][2][7] also show numbers that are quite substantially different.

Thank you for your insightful comments and constructive feedback. Below, we provide detailed responses to each point raised:

1. Reviewer: HPA introduces significant computational overhead; a detailed complexity analysis is missing. Response:Thank you for pointing this out. We have modified the paper and replaced the section 'Gradient Computation in HPA' with 'Complexity Analysis'. Complexity Analysis We calculate the FLOPs for the Scaled Hadamard-Product Attention (HPA) mechanism, with $N$ as the number of nodes, $s$ as the sequence length, and $d$ as the feature dimension. Constructing the key and value matrices $K$ and $V$ requires linear transformations across each sequence step, amounting to $2Nsd^2$ FLOPs. The calculation of attention weights $a$, including the element-wise product and softmax operation, requires $2Nsd$ FLOPs. The element-wise multiplication of $a$ and $V$, followed by a linear transformation to reshape the result back to dimension $d$, adds $Nsd + Nsd^2$ FLOPs. Summing these components, the total FLOPs for HPA is $3Nsd^2 + 3Nsd$. This complexity scales linearly with the number of nodes $N$, with a fixed sequence length $s$. Compared to a global self-attention mechanism, our decoder-only architecture avoids the $O(N^2)$ complexity, making it especially efficient and advantageous for large-scale simulations.

2.More attention-based baselines should be compared, such as [1], [2]. Response:Thank you for your suggestion. Currently, we use state-of-the-art message-passing GNNs(MGN, BSMS) and Transformer-based method(TIE) as our baselines. However, to enable a more comprehensive comparison, we have added Mesh-Transformer[1] and Graph MLP-MIXER[3] as additional baselines for testing. We additionally test HCMT[2] only on DeformingPlate dataset because its CMT and HMT modules were ultimately designed specifically for collision cases. Although MPT is slower, it outperforms all competitors in both short- and long-term predictions.

[Experiment: Additional Baslines Comparison]:

3. Figure 3a annotations are unclear; which line represents HPA? Response: Thank you for pointing this out. We apologize for the annotation mistake; the gray line represents Dot-Product Attention, while the red line represents our Hadamard-Product Attention. We have corrected this in the paper accordingly.

Thank you for your thoughtful review and valuable feedback. We appreciate the time you have taken to evaluate our work. We would like to ask if there is anything else we can do to improve the rating of our submission. Your insights would be greatly appreciated as we continue to refine our approach.

We conducted robustness tests on GFL, further addressing your concern in Weakness 1 regarding the robustness of GFL: Our GFL remains robust under varying graph connectivity due to the Energy-Based Robustness: GFL operates on frequency-domain energy distributions ($E$) rather than relying directly on eigenvector precision. By focusing on mean energy across high- and low-energy components, it reduces sensitivity to small eigenvector perturbations caused by narrow eigenvalue gaps.

[Experiment: GFL Robustness Performance on Airfoil]:

Additionally, in our discussions with Reviewer ZnqH, we further optimized the computational efficiency of our method, significantly reducing the increase in time and memory usage compared to MGN. This means that the wall clock time and memory footprint of our model are no longer limiting factors you mentioned in Weakness 2. As shown in the table below, the optimized version of HMIT achieves substantial improvements in both training and inference time, while also reducing memory usage:

[Experiment: Memory Consumption and Computational Speed]:

Despite the promising advancements offered by the Historical Message-Passing Integration Transformer (HMIT) in simulating physical systems, its initial implementation faced notable limitations in computational speed and memory consumption. To address these challenges, we incorporated two key optimizations. First, KV Cache eliminates redundant computations by caching key and value matrices, reducing the attention complexity from $O(3Nsd^2 + 3Nsd)$ to $O(Nsd^2 + 3Nsd)$. Second, Dynamic Weighted Value Selection dynamically selects the first $m$ rows of the weighted value matrix $w \in \mathbb{R}^{s \times d}$, where $m$ corresponds to the current message-passing step, further enhancing computational efficiency. As shown in the above table, these optimizations significantly reduce training and inference times while decreasing memory requirements.

Thank you for your insightful comments and constructive feedback. Below, we provide detailed responses to each point raised:

1. The computational cost for precomputing Laplacian eigenvectors Response:

[Experiment: Precomputation Cost]:

While we provide a precomputation cost table for clarity, this data processing only occurs once and can be considered part of the dataset generation process.

The primary computational cost is in the eigen decomposition of the graph Laplacian, with a total precomputation time of about 3 minutes—negligible compared to training time. For the four datasets with static graph topologies, eigen decomposition is only needed at $t=0$ rather than at each time step.

For datasets with dynamic topologies, eigen decomposition is required at each time step but only needs to be done once, in batch, before training, and cached for reuse. This approach reduces on-the-fly calculation complexity, improving training efficiency. Importantly, inference does not require precomputed eigenvalues or eigenvectors, so inference speed remains unaffected.

2. The precomputation of Laplacian eigenvectors is not feasible for dynamic graphs with frequent topology changes Response: Thank you for pointing this out. We clarify that our GFL can be applied to cases with changing graph topology, as long as GFT is performed during preprocessing for each time step. Relevant sections of the paper have been revised accordingly.

3. The ablation study in Table 2 is limited to the CylinderFlow dataset. Response:

[Experiment: ablation in FlagSimple]

In the ablation study focusing on the FlagSimple dataset, the integration of both HPA and GFL demonstrated the highest improvement in reducing error rates across all RMSE measures.

4. In Figure 3(a), the graph label is confusing, especially for demonstrating Hadamard-Product Attention effectiveness. Response: Thank you for pointing this out. We apologize for the annotation mistake; the gray line represents Dot-Product Attention, while the red line represents our Hadamard-Product Attention. We have corrected this in the paper accordingly.

5. What is meant by "weight of the edge" in Line 131? Response: We have modified the paper’s wording to make it clearer: The adjacency matrix of $G$ is denoted by $A$, where $A_{ij} = 1$ if there is an edge between vertices $i$ and $j$, and $A_{ij} = 0$ otherwise.

6. How is the weighted value matrix 𝑤 reshaped back to 𝑑 dimensions, and is an MLP used for this? Response: The weighted value matrix $w \in \mathbb{R}^{d \times s}$ is reshaped to $d$ dimensions via a learned linear layer $W_{\text{linear}} \in \mathbb{R}^{s \times 1}$, applying a transformation $w W_{\text{linear}}$. This is a simple linear layer without non-linear activation, not an MLP.

7. During each step of message passing, are the networks $f_1, f_2, f_3, f_4, f_5$ shared or independent? Response: The five networks are independent of each other.

• $f_1$ and $f_2$ take node features (e.g., velocity) and edge features (e.g., length) as inputs, respectively, and output 128-dimensional vectors.
• $f_3$ is an independent MLP for each MP layer. For example, if $\text{MP} = 15$, then $f_3$ generates 15 MLPs, each with an input size of 128, an output size of 128.
• $f_4$ uses the same weights across different MP layers and implements Multihead Hadamard-Product Attention.
• $f_5$ decodes the 128-dimensional node latent features back to the real space.

Thank you for your answers.

Regarding W1. Can you analyze the computational complexity when using Hadamard-Product Attention compared to regular message-passing and Dot-Product Attention?

Regarding Table 1. Why is there a huge discrepancy in the results of MGN and BSMS compared to the results in their original paper? One may even find some results of baselines in their original paper are better than the results of the proposed MPT (eg. RMSE-all of MGN on Airfoil and BSMS on Plate).

I also recommend that the authors upload the revised paper to address some concerns in the rebuttal.

Dear reviewer, this is a friendly reminder that the public discussion period for Submission 2831 is ending soon. If you have any questions or concerns about our rebuttal, we’d be happy to address them. Thank you for your time and feedback!

Thank you for your insightful comments and constructive feedback. Below, we provide detailed responses to each point raised:

1. Wall clock time and memory comparisons are missing. A detailed time and memory analysis per component would aid future research. Response: We have noted this as a limitation in the paper and also discussed it in section F of the appendix. Additionally, we have included a graph Laplacian precomputation cost table and replace "Gradient Computation in HPA" section with "Complexity Analysis" in the paper. [Experiment: Precomputation Cost]:

While we provide a precomputation cost table for clarity, this data processing only occurs once and can be considered part of the dataset generation process.

The primary computational cost is in the eigen decomposition of the graph Laplacian, with a total precomputation time of about 3 minutes—negligible compared to training time. For the four datasets with static graph topologies, eigen decomposition is only needed at $t=0$ rather than at each time step.

For datasets with dynamic topologies, eigen decomposition is required at each time step but only needs to be done once, in batch, before training, and cached for reuse. This approach reduces on-the-fly calculation complexity, improving training efficiency. Importantly, inference does not require precomputed eigenvalues or eigenvectors, so inference speed remains unaffected.

Complexity Analysis We calculate the FLOPs for the HPA mechanism, with $N$ as the number of nodes, $s$ as the sequence length, and $d$ as the feature dimension. Constructing the key and value matrices $K$ and $V$ requires linear transformations across each sequence step, amounting to $2Nsd^2$ FLOPs. The calculation of attention weights $a$, including the element-wise product and softmax operation, requires $2Nsd$ FLOPs. The element-wise multiplication of $a$ and $V$, followed by a linear transformation to reshape the result back to dimension $d$, adds $Nsd + Nsd^2$ FLOPs. Summing these components, the total FLOPs for HPA is $3Nsd^2 + 3Nsd$. This complexity scales linearly with the number of nodes $N$, with a fixed sequence length $s$. Compared to a global self-attention mechanism, our decoder-only architecture avoids the $O(N^2)$ complexity, making it especially efficient and advantageous for large-scale simulations.

2. The derivation of GFL, including the expression of $\alpha$ and the segment rate $s_r$, is unclear, making it difficult to extend. Response:

• The concept of the segment rate $s_r$ was introduced to divide the frequency domain energy into high-energy and low-energy components, ensuring that the model focuses on the most relevant frequency components during training.
• The adjustment factor $\alpha$, derived from the ratio of mean energy in high- and low-energy partitions, adaptively balances their contributions to the loss.

• Without the adjustment factor $\lambda$, $\alpha$ would excessively boost low-energy components to match high-energy levels, amplifying noise and reducing performance. Introducing $\lambda$ allows adaptively controlled scaling in low-energy regions.

3. Section 1, line 47-48: Should the phrase about replacing MLP with Hadamard-Product Attention instead reference replacing Scaled Dot Product Attention? Response: Thank you for your feedback.

• In lines 47-48, "MLP" refers to using MLPs in each Message Passing step to model node-edge interactions, similar to MGN and GNS, where concatenated node and edge features are input and the next step’s node features are output.
• In our architecture, however, the FFN applies non-linearity to HPA-modulated node information, with nodes as both input and output.

Dear reviewer, this is a friendly reminder that the public discussion period for Submission 2831 is ending soon. If you have any questions or concerns about our rebuttal, we’d be happy to address them. Thank you for your time and feedback!

Thank you for your thoughtful review and valuable feedback. We appreciate the time you have taken to evaluate our work. We would like to ask if there is anything else we can do to improve the rating of our submission. Your insights would be greatly appreciated as we continue to refine our approach.

In our discussions with Reviewer ZnqH, we further optimized the computational efficiency of our method, significantly reducing the increase in time and memory usage compared to MGN. This means that the wall clock time and memory footprint of our model are no longer limiting factors. As shown in the table below, the optimized version of HMIT achieves substantial improvements in both training and inference time, while also reducing memory usage:

Despite the promising advancements offered by the Historical Message-Passing Integration Transformer (HMIT) in simulating physical systems, its initial implementation faced notable limitations in computational speed and memory consumption. To address these challenges, we incorporated two key optimizations. First, KV Cache eliminates redundant computations by caching key and value matrices, reducing the attention complexity from $O(3Nsd^2 + 3Nsd)$ to $O(Nsd^2 + 3Nsd)$. Second, Dynamic Weighted Value Selection dynamically selects the first $m$ rows of the weighted value matrix $w \in \mathbb{R}^{s \times d}$, where $m$ corresponds to the current message-passing step, further enhancing computational efficiency. As shown in the above table, these optimizations significantly reduce training and inference times while decreasing memory requirements.

Thank you for your insightful comments and constructive feedback. Below, we provide detailed responses to each point raised:

1. Additional landmark works Response: Thank you for your suggestion. We have included The Universal Physics Transformers by Alkin et al. and Poseidon: Efficient Foundation Models for PDEs by Herde et al. as additional citations, to acknowledge recent foundational work in this area. As our model is a GNN-based architecture, we have avoided direct comparisons with non-GNN-based methods. However, we have added three additional GNN-based baselines in response to the question below.

2. The paper's evaluation and literature integration are lacking. Response: [Experiment: Additional Dataset Test]:

We additionally conducted an experiment on the Dam Problem from the 2023 CFDBench [2] benchmark, which models the rapid release of water from a column collapse and represents complex free-surface flows with varying velocities. Our method outperforms previous methods (MGN, BSMS, and TIE) across all three RMSE metrics—RMSE-1, RMSE-50, and RMSE-all—by 33%-57%.

[Experiment: Additional Baslines Comparison]:

Currently, we use state-of-the-art message-passing GNNs(MGN, BSMS) and Transformer-based method(TIE) as our baselines. However, to enable a more comprehensive comparison, we have added Mesh-Transformer[1] and Graph MLP-MIXER[3] as additional baselines for testing. We additionally test [2] only on DeformingPlate dataset because its CMT and HMT modules were ultimately designed specifically for collision cases.