Topology-aware Neural Flux Prediction Guided by Physics

Thank you for your constructive comments and questions. We address them in a Q&A format as follows.

Q1: How does PhyNFP perform on similar physical systems with different constraints?

A1: The PDEs adopted by PhyNFP (Saint-Venant (S-V) for hydrodynamics and Aw-Rascle (A-R) for traffic flow) are basic and general formulations applicable to diverse scenarios within river and traffic networks. For instance, the gravitational term in S-V equations, crucial for elevation variations, can be naturally attenuated for flat rivers by adjusting learned weights during training. This enables PhyNFP to adapt across varying physical constraints without structural changes. Moreover, our ablation study shows that PhyNFP is robust even without explicit PDE constraints, which confirms its generalizability by using the difference operators alone.

Q2: Can you give detailed dataset descriptions such as the number of training sequences and testing sequences?

A2: We have provided dataset descriptions in Section 4 (Datasets, page 5). In particular, for the river dataset derived from LamaH-CE, we used data from the period 2000–2017, where the data from the years of 2016 and 2017 were used as the test set, and the remaining years were used for training. We followed [1] to make this train/test split. In our revised manuscript and forthcoming code repository, we will include all dataset statistics and detailed train/test splits for all datasets used.

Q3: Can you provide detailed implementation of PhyNFP and baseline models?

A3: The implementation specifics of PhyNFP are presented in Section 3.3, including the construction of difference matrices $D_1$ and $D_2$ and the PDE integration mechanism. In our revised manuscript and subsequent code release, we will further provide complete architectural descriptions, hyperparameter settings, and detailed training settings. The baseline models used (GraphSAGE, GCN, GAT, GWN, MP-PDE Solver, MPNN, and GNO) are introduced in Section 4 (Competitors, page 6).

Q4: Comparison with Mesh-Based Benchmarks

A4: We clarify key differences between mesh-based methods and our graph-based approach:

1. Applicability to sparse, irregular networks. Mesh-based PDE methods require structured grids with regular spacing and clear geometric coordinates to ensure numerical accuracy. This does not match our datasets, which are sparse and irregular (such as river and traffic networks). In contrast, PhyNFP uses topology-based difference matrices that work well on irregular graphs and naturally capture directional flow.

2. Modeling topological structure. Mesh-based methods focus mostly on spatial resolution, but they do not model topological effects like upstream-downstream structure. Our method is designed to capture these topological patterns, which are important for real-world networks. This kind of topological sensitivity is hard to achieve using standard mesh-based solvers.

Q5: Efficiency / Scalability

A5: Although computational efficiency is not our core research contribution, we add an experiment below to show that the runtime of our method exhibits sub-linear scaling with larger graph size. As shown in the table above, when the number of nodes increases from 31 to 358 (more than 10×), the runtime per epoch only doubles, demonstrating that computational complexity grows significantly slower compared to the upscaling of graph size.

Table: Average Runtime per Epoch for Different Graph Structures

Reference:

[1] Nikolas Kirschstein and Yixuan Sun, The Merit of River Network Topology for Neural Flood Forecasting, ICML 2024.

Thank you for your constructive comments and questions. We address them in a Q&A format as follows.

Q1: How does PhyNFP scale to large networks and what are its computational bottlenecks?

Table: Average Runtime per Epoch for Different Graph Structures

A1: We add an experiment to show that the runtime of our method exhibits excellent sub-linear scaling with larger graph size. As shown in the table above, when the number of nodes increases from 31 to 358 (more than 10×), the runtime per epoch only doubles, demonstrating that computational complexity grows significantly slower than the graph size.

Q2: Can this framework be applied to link prediction tasks in dynamic graphs?

A2: Our framework is designed for flux prediction in physical systems, focusing on modeling directional flows governed by physical laws. Link prediction in dynamic graphs is a structurally different task, aiming to infer future or missing edges. Such problems are rarely encountered in our physical settings. We would like to defer this interesting setup by adapting our method for dynamic link prediction as a future work.

Q3: Is PhyNFP robust to incomplete or noisy data? How does it handle missing node attributes?

A3: Yes, our model is robust to both incomplete and noisy data. Our datasets include real-world measurement data, which naturally contain both noisy and missing entries. Our method can mitigate the effect of noise through its physics-guided inductive bias, which acts as implicit regularization to enhance robustness.

For missing node attributes, we adopt the preprocessing step used in~[1], where nodes with missing features are excluded from the network.

Q4: What does the "undirected" in Figure 1 mean and how it was implemented?

A4: In the “undirected” setting we remove edge directions by symmetrizing the adjacency matrix (i.e., using $A_{\text{undirected}} = A + A^\top$), so each node aggregates messages from its upstream and downstream neighbors. In this setting, the model does not differentiate between upstream and downstream connections. In contrast, the forward setting aggregates only from upstream neighbors, and the reverse setting flips all edge directions, so aggregation only considers those originally downstream nodes.

[1] Nikolas Kirschstein and Yixuan Sun, The Merit of River Network Topology for Neural Flood Forecasting, ICML 2024.

Thank you, and we address your questions in a Q&A fashion as follows

Q1: Why the reverse problem will amplify numerical errors? Is it a concern for the forward problem?

A1: Solving inverse problems is ill-posed, as many different upstream conditions can lead to the same downstream fluxes. When inferring upstream states from downstream inputs, the mapping becomes one-to-many and numerically unstable, small noise at downstream nodes can cause large errors upstream. In our inverse experiment (Q9), we observe that GCN suppresses such noise entirely, showing no upstream variation. This supports our hypothesis that GCN cannot distinguish forward from reverse settings and lacks sensitivity to edge direction. Further details are provided in our response to Q9.

This is not a concern in the forward setting because the PDE describing physical structure can enforce a stable and one-to-one mapping from upstream to downstream nodes.

Q2: Generalize to broader systems cannot be discretized?

A2: Our discretization does not rely on any specific PDE; rather, it is based on the topological structure describing the system. In fact, we can construct multiple difference operators [1] from the graph through its spatial or functional adjacency, and our difference matrix derived from local spatial variations and directional flow is one possible instantiation.

Q3: How were Eqs.(12) and (13) derived from Eqs.(8) and (11)?

A3: We will supplement a step-by-step derivation, as presented in (https://anonymous.4open.science/r/PhyNFP-D88F/Q3.pdf). The main idea is to use independent MLPs that take raw node inputs to learn representations of physical quantities. Difference matrices $D_1$ and $D_2$ are applied to guide the representation learning, so to align with the original PDE. The intuitions behind $D_1$ and $D_2$ are in our response to Q4.

Q4: Are D1 and D2 just introduced to enable Eqs.(12) and (13) derivation?

A4: $D_1$ and $D_2$ are not just for derivation, they capture key physical quantities. $D_1$ models horizontal gradients on the graph, reflecting convection and spatial variation. $D_2$ captures vertical gradients from elevation, tied to gravity-driven flow. These matrices allow our framework to generalize to other PDE systems involving similar terms.

Q5: The text following Eq. (11) does not with the equation.

A5: We will fix this typo, as also suggested by Reviewer RhQb.

Q6: Why is $\Delta t$ (and $\hat{g}$) set as learnable scalar but not fixed?

A6: The amount of $\Delta t$ determines either a message-passing layer will update the node embeddings by reusing the output from previous layer (small $\Delta t$) or it will incorporate the information for updating the physical status of PDE (large $\Delta t$). For example, a large $\Delta t$ in Eq.(12) means that the embedding update will mostly rely on using the convection and gravity terms derived from the Eq.(6). By learning it from data, we make $\Delta t$ adaptable to real conditions.

We validate through two experiments: 1) Larger $\Delta t$ improves $DS$, as shown in Table 1 in the anonymous link ( https://anonymous.4open.science/r/PhyNFP-D88F/README.md), although inverse problems typically prefer smaller steps [2]. With fixed layers, small $\Delta t$ limits each update. A learnable $\Delta t$ adapts to network depth without this concern. 2) Figure 2 in the link plots the variation of $\Delta t$ w.r.t. the number of epochs, where it converges to small values in the reverse setting.

Q7: Handle input features with 24 time steps?

A7: Concatenation. Please refer to our response to Q1 of Reviewer RhQb.

Q8: Misunderstanding in DS?

A8: Indeed a high reverse MSE can produce high DS performance. However, our model enjoys the highest DS = +.0105 by attaining the lowest reverse MSE = .0906. The positive and high DS value indicates a better distinguishability of edge directions.

Q9: 1) Why ResGCN is not in Table 1? 2) How perturbations reflect in reverse setting?

A9: Answer to 1) is in Q3 of Reviewer RhQb.

To answer 2), we add local perturbation in reverse setting (and the Traffic Network) as shown in Fig 1 (and Fig 3) in the anonymous link.

Q10: Performance downgrade in traffic network?

A10: Mainly due to cycles in the traffic network. Using Johnson's algorithm, we find 80 cycles in it, but the river network has zero. These cycles allow messages to propagate in both directions, making it more challenging for models to distinguish forward and reverse flows.

Q11: Details on the model architecture and training settings?

A11: We use 19 layers based on the graph diameter in our dataset. In experiment we observe that our implementation is robust across different layer numbers. We will publish our code with all hyperparameter settings for reproducibility.

[1] Grady, et al. Discrete calculus: Applied analysis on graphs for computational science[M]. London: Springer, 2010.

[2] Bertero, et al. The stability of inverse problems[M]. Springer Berlin Heidelberg, 1980.

Thank you for your constructive comments and questions. We address them in a Q&A format as follows.

Q1: How was the flux prediction task modeled? Is it autoregressive for the past 24 time steps?

A1: The flux prediction task is formulated as a supervised node regression problem rather than an autoregressive one. For each node, a fixed-length temporal window of 24 past time steps $W=24$ is concatenated into one feature vector, following the prior research~[1]. This vector is used to predict the future flux $y$ at time $t+6$. The base learner is graph neural network that integrates discretized difference matrices and physical priors in its message passing process.

Q2: Notation inconsistency in Eq.(11) and their elaborations.

A2: This is a typo, and we will make revision in our camera-ready as follows:

$$\rho_i^{t+1} = \rho_i^t - \alpha \left( \rho_i^t (\hat{D} u^t)_i + u_i^t (\hat{D} \rho^t)_i \right),$$

where $(\hat{D}u^t)_i$ represents velocity differences, and $(\hat{D}\rho^t)_i$ encodes density-driven effects. These terms approximate the spatial derivatives of velocity and density, respectively, using a difference operator $\hat{D}$. The scalar factor $\alpha$ is defined as $\Delta t / \Delta x$.

Q3: What is the main finding of Fig.2? Why is ResGCN good/bad given its reaction to the "change" (what is this change?) in v-1? Why is ResGCN not included as one of the GNN competitors?

A3: The "change" refers to a manually injected spike in the input flux at node $v_1$, simulating a sudden perturbation in the river network. Fig.2 examines how this perturbation propagates downstream through the model predictions. The full analysis of the main finding in Fig.2 is provided in RQ4 of the main text. We will further clarify in the final version.

The result of ResGCN in Figure 2 is bad compared to our method because it fails to capture the correct flow dynamics. Specifically, when the perturbation is applied to a node, it should cause an increase in its downstream flux, while GCN shows almost no change in downstream nodes. Although ResGCN propagates this perturbation signal to some extent, it produces incorrect trends; for example, predicting a decrease at node $v_2$ where an increase is expected.

We select GCN due to its wider application over its residual variant. In fact, ResGCN exhibits similar performance to standard GCN with shallow layers. Even after careful tuning (e.g., adding more layers), ResGCN only performs on a par with GWN by having MSE (F) = .1114, MSE (R) = .1139, and RDS = -76.2%. We shall include those results in our camera ready.

Q4: How robust is PhyNFP to data incompleteness and how are the missing entries handled?

A4: Our datasets include real-world measurement data, which naturally contain both noisy and missing entries. Our method can mitigate the effect of noise through its physics-guided inductive bias, which acts as implicit regularization to enhance robustness.

For missing node attributes, we adopt the preprocessing step used in~[1], where nodes with missing features are excluded from the network.

Q5: What are the settings of forward vs reverse vs undirected?

A5: In our work, these terms describe different ways of using edge directions during graph construction and message passing. In the forward setting, we use the original directed graph where edges follow the real physical flow (i.e., from upstream to downstream). This reflects the correct direction of information propagation.

In the reverse setting, all edge directions are flipped (i.e., the adjacency matrix is transposed), so information flows in the opposite direction (i.e., from downstream to upstream), which violates the physical rules.

In the undirected setting, edge directions are ignored. Each node exchanges messages with both upstream and downstream neighbors, which is common in standard GNNs.

Reference: [1] Nikolas Kirschstein and Yixuan Sun, The Merit of River Network Topology for Neural Flood Forecasting, ICML 2024.