Port-Hamiltonian Architectural Bias for Long-Range Propagation in Deep Graph Networks

## Abstract

Abstract The dynamics of information diffusion within graphs is a critical open issue that heavily influences graph representation learning, especially when considering long-range propagation. This calls for principled approaches that control and regulate the degree of propagation and dissipation of information throughout the neural flow. Motivated by this, we introduce (port-)Hamiltonian Deep Graph Networks , a novel framework that models neural information flow in graphs by building on the laws of conservation of Hamiltonian dynamical systems. We reconcile under a single theoretical and practical framework both non-dissipative long-range propagation and non-conservative behaviors, introducing tools from mechanical systems to gauge the equilibrium between the two components. Our approach can be applied to general message-passing architectures, and it provides theoretical guarantees on information conservation in time. Empirical results prove the effectiveness of our port-Hamiltonian scheme in pushing simple graph convolutional architectures to state-of-the-art performance in long-range benchmarks.

## 1 Introduction

The conjoining of dynamical systems and deep learning has become a topic of great interest in recent years. In particular, neural differential equations (neural DEs) demonstrate that neural networks and differential equations are two sides of the same coin (Haber and Ruthotto, ; Chen et al., ; Chang et al., ). This connection has been pushed to the domain of graph learning (Bacciu et al., ; Wu et al., ), forging the field of differential-equations inspired Deep Graph Networks (DE-DGNs) (Poli et al., ; Chamberlain et al., ; Gravina et al., ; Han et al., ).

In this paper, we are interested in designing the information flow within a graph as a solution of a (port-)Hamiltonian system (Van der Schaft, ), which is a general formalism for physical systems that allows for both conservative and non-conservative dynamics, with the aim of allowing flexible long-range propagation in DGNs. Indeed, long-range propagation is an ongoing challenge that limits the power of the Message-Passing Neural Network (MPNNs) family (Gilmer et al., ), as their capacity to transmit information between nodes exponentially decreases as the distance increases (Alon and Yahav, ; Di Giovanni et al., ). This prevents DGNs from effectively solving real-world tasks, e.g., predicting anti-bacterial properties of peptide molecules (Dwivedi et al., ).
While recent literature proposes various approaches to mitigate this issue, such as graph rewiring (Gasteiger et al., ; Topping et al., ; Gutteridge et al., ) and graph transformers (Shi et al., ; Dwivedi and Bresson, ; Wu et al., ), here we aim to address this problem providing a theoretically grounded framework through the prism of Hamiltonian-inspired DE-DGNs. Therefore, we propose *(port-)Hamiltonian Deep Graph Network* (PH-DGN) a new message-passing scheme inspired by (port-)Hamiltonian dynamics. By design, our method introduces the flexibility to balance non-dissipative long-range propagation and non-conservative behaviors as required by the specific task at hand. Therefore, when using purely Hamiltonian dynamics, our method allows the preservation and propagation of long-range information by obeying the conservation laws. In contrast, when our method is used to its full extent, internal damping and additional forces can deviate from this purely conservative behavior, potentially increasing effectiveness in the downstream task. Leveraging the connection with Hamiltonian systems, we provide theoretical guarantees that information is conserved over time. Lastly, the general formulation of our approach can seamlessly incorporate any neighborhood aggregation function (i.e., DGN), thereby endowing these methods with the distinctive properties of our PH-DGN.

Our main contributions can be summarized as follows: (i) We introduce PH-DGN, a novel general DE-DGN inspired by (port-)Hamiltonian dynamics, which enables the balance and integration of non-dissipative long-range propagation and non-conservative behavior while seamlessly incorporating the most suitable aggregation function; (ii) We theoretically prove that, when pure Hamiltonian dynamic is employed, both the continuous and discretized versions of our framework allow for long-range propagation in the message passing flow, since node states retain their past; (iii) We introduce tools inspired by mechanical systems that deviate the Hamiltonian dynamic from its conservative behavior, thus facilitating a clear interpretation from the physics perspective; and (iv) We conduct extensive experiments to demonstrate the benefits of our method and the ability to stack thousands of layers. Our PH-DGN outperforms existing state-of-the-art methods on both synthetic and real-world tasks.

## 2 (Port-)Hamiltonian Deep Graph Network

We consider the problem of learning node embeddings for a graph $\mathcal{G}=(\mathcal{V},\mathcal{E})$, where $\mathcal{V}$ is a set of $n$ entities (the nodes) interacting through relations (i.e., edges) in $\mathcal{E}\subseteq\mathcal{V}\times\mathcal{V}$. Each node $u\in\mathcal{V}$ is associated to state ${\bf x}_{u}(t)\in\mathbb{R}^{d}$, that is the representation of the node at time $t$. The term ${\bf X}(t)\in{\mathbb{R}^{n\times d}}$ is the matrix of all node states’ in graph $\mathcal{G}$.

We introduce a new DE-DGN framework that designs the information flow within a graph as the solution of a (port-)Hamiltonian system (Van der Schaft, ).
Hamiltonian mechanics is a formalism for physical systems based on the Hamiltonian function $H(\mathbf{p},\mathbf{q},t)$, which represents the generalized energy of the system with position $\mathbf{q}$ and momentum $\mathbf{p}$. A classic example of a Hamiltonian system is that of a simple mass-spring pendulum, with mass $m$ attached to a spring with constant $k$ having position $q=x$ and momentum $p=m\dot{x}$. The Hamiltonian of the system is the total energy $H=K+P$ where $K$ is the kinetic component $K=\frac{p^{2}}{2m}=\frac{1}{2}m\dot{x}^{2}$ and $P$ the spring potential one $P=\frac{1}{2}kx^{2}$. Hamilton’s equations are then defined as:

$$ $\dot{p}=-\pdv{H}{q}=-kx,\quad\dot{q}=\pdv{H}{p}=\frac{p}{m},$ (1) $$

from which we recover the well-known mass-spring pendulum equation $m\ddot{x}=-kx$. The Hamiltonian formalism allows for an easy description of the dynamics of a system based on its energy and provides theoretical results that will allow us to guarantee relevant properties on our system,
such as energy preservation. It is widely used both in classical mechanics Arnold et al. () as well as in quantum mechanics Griffiths and Schroeter () due to its generality. Port-Hamiltonian systems let us introduce non-conservative phenomena in the system, such as internal dampening $D(q)p$ and external forcing $F(q,t)$, acting on the momentum equation as $\dot{p}=-\pdv{H}{q}-D(q)p+F(q,t)$.

Here, we show how the Hamiltonian formulation provides the backing to preserve and propagate long-range information between nodes, in adherence to the laws of conservation. The casting of the system in the more general port-Hamiltonian setting, then, introduces the possibility of
trading non-dissipation with non-conservative behaviors when needed by the task at hand.
Our approach is general, as it can be applied to any message-passing DGN, and frames in a theoretically sound way the integration of non-dissipative propagation and non-conservative behaviors.
In the following, we refer to a *Hamiltonian Deep Graph Network* (H-DGN) when the framework is instantiated to a purely conservative DGN, and *port-Hamiltonian Deep Graph Network* (PH-DGN) otherwise. Figure 1 shows our high-level architecture hinting at how the initial state of the system is propagated up to the terminal time $T$. While the state evolves preserving energy, internal dampening and additional forces can intervene to alter its conservative trajectory.

Figure: Figure 1: A high-level overview of the proposed (port-)Hamiltonian Deep Graph Network. It summarizes how the initial node state ${\bf x}_{u}(0)$ is propagated by means of energy preservation up until the terminal time $T$ (i.e., layer $L$), ${\bf x}_{u}(T)$. While the global system’s state $\mathbf{y}$ evolves preserving energy, external forces (i.e., dampening $D(\mathbf{y})$ and external control $F(\mathbf{y},t)$) can intervene to alter its conservative trajectory. The gray trajectories between the initial and final states represent the continuous evolution of the system. The discrete message passing step from layer $\ell$ to $\ell+1$, which is shown in middle of the figure, is given by the coupling of coordinates ${\bf q}$ and momenta ${\bf p}$ in terms of neighborhood aggregation $\Phi_{\mathcal{G}}$ and influence to adjacent neighbors $\Phi^{*}_{\mathcal{G}}$. Self-influence on both ${\bf q}$ and ${\bf p}$ from the previous step $\ell$ are omitted for simplicity.
Refer to caption: x1.png

In the following, we first derive our method from a purely Hamiltonian system and prove its conservative behavior theoretically. We then complete our method by incorporating non-conservative behaviors through port-Hamiltonian dynamics.

#### Hamiltonian-inspired message passing.

To inject the Hamiltonian dynamics into a DE-DGN, we start by considering the graph Hamiltonian system described by the following ODE

$$ $\odv{{\bf y}(t)}{t}=\mathcal{J}\nabla H_{\mathcal{G}}({\bf y}(t)),$ (2) $$

for time $t\in[0,T]$ and subject to an initial condition ${\bf y}(0)={\bf y}^{0}$. The term ${\bf y}(t)\in\mathbb{R}^{nd}$ is the vectorized view of ${\bf X}(t)$ that represents the global state of the graph at time $t$, with an even dimension $d$, following the notation of Hamiltonian systems (Hairer et al., ). $H_{\mathcal{G}}:\mathbb{R}^{nd}\to\mathbb{R}$ is a neural-parameterized Hamiltonian function capturing the energy of the system. The skew-symmetric matrix $\mathcal{J}=\begin{pmatrix}\bf 0&-{\bf I}_{nd/2}\\
{\bf I}_{nd/2}&\bf 0\end{pmatrix}$, with ${\bf I}_{nd/2}$ being the identity matrix of dimension $nd/2$, reflects a rotation of the gradient $\nabla H_{\mathcal{G}}$ and couples the position and momentum of the system.

Since we are dealing with a Hamiltonian system, the global state ${\bf y}(t)$ is composed by two components which are the momenta, ${\bf p}(t)=({\bf p}_{1}(t),\dots,{\bf p}_{n}(t))$, and the position, ${\bf q}(t)=({\bf q}_{1}(t),\dots,{\bf q}_{n}(t))$, of the system, thus ${\bf y}(t)=({\bf p}(t),{\bf q}(t))$. Therefore, from the node (local) perspective, each node state is expressed as ${\bf x}_{u}(t)=({\bf p}_{u}(t),{\bf q}_{u}(t))$.

Under this local node-wise perspective, Eq. 2 can be equivalently written as

$$ $\odv{{\bf x}_{u}(t)}{t}=\begin{pmatrix}\dot{{\bf p}}_{u}(t)\\ \dot{{\bf q}}_{u}(t)\end{pmatrix}=\begin{pmatrix}-\nabla_{{\bf q}_{u}}H_{ \mathcal{G}}({\bf p}(t),{\bf q}(t))\\ \phantom{-}\nabla_{{\bf p}_{u}}H_{\mathcal{G}}({\bf p}(t),{\bf q}(t))\end{ pmatrix},\,\,\forall u\in\mathcal{V}.$ (3) $$

With the aim of designing a Hamiltonian system based on message passing, we instantiate the Hamiltonian function $H_{\mathcal{G}}$ as

$$ $H_{\mathcal{G}}({\bf y}(t))=\sum_{u\in\mathcal{V}}\tilde{\sigma}({\bf W} \mathbf{x}_{u}(t)+\Phi_{\mathcal{G}}(\{{\bf x}_{v}(t)\}_{v\in\mathcal{N}_{u}}) +{\bf b})^{\top}\mathbf{1}_{d},$ (4) $$

where $\tilde{\sigma}(\cdot)$ is the anti-derivative of a monotonically non-decreasing activation function $\sigma$, $\mathcal{N}_{u}$ is the neighborhood of node $u$, and $\Phi_{\mathcal{G}}$ is a neighborhood aggregation permutation-invariant function. Terms ${\bf W}\in\mathbb{R}^{d\times d}$ and ${\bf b}\in\mathbb{R}^{d}$ are the weight matrix and the bias vector, respectively, containing the trainable parameters of the system; $\mathbf{1}_{d}$ denotes a vector of ones of length $d$.

By computing the gradient $\nabla_{{\bf x}_{u}}H_{\mathcal{G}}({\bf y}(t))$ we obtain an explicit version of Eq. 3, which can be rewritten
from the node-wise perspective
of the information flow as the sum of the self-node evolution influence and its neighbor’s evolution influence (referred to as $\Phi_{\mathcal{G}}^{*}$). More formally, for each node $u\in\mathcal{V}$

$$ $\odv{\mathbf{x}_{u}(t)}{t}=\mathcal{J}_{u}\Biggl{[}{\bf W}^{\top}\sigma({\bf W }\mathbf{x}_{u}(t)+\Phi_{\mathcal{G}}(\{{\bf x}_{v}(t)\}_{v\in\mathcal{N}_{u}} )+{\bf b})\\ +\underbrace{\sum_{v\in\mathcal{N}_{u}\cup\{u\}}{\pdv{\Phi_{\mathcal{G}}(\{{ \bf x}_{v}(t)\}_{v\in\mathcal{N}_{u}})}{\mathbf{x}_{u}{(t)}}^{\top}\sigma({\bf W }\mathbf{x}_{u}(t)+\Phi_{\mathcal{G}}(\{{\bf x}_{v}(t)\}_{v\in\mathcal{N}_{u}} )+{\bf b})}}_{\Phi^{*}_{\mathcal{G}}}\Biggr{]}.$ (5) $$

Here, $\mathcal{J}_{u}$ has the same structure as $\mathcal{J}$, but the identity blocks have dimension $d/2$ as it is applied to the single node $u$. We refer to the system in Eq. 5 as *Hamiltonian Deep Graph Network* (H-DGN) as it adheres solely to conservation laws.

Now, given an initial condition ${\bf x}_{u}(0)$ for a node $u$, and the other nodes in the graph, the ODE defined in Eq. 5 (i.e., H-DGN) is a continuous information processing system over a graph governed by conservation laws that computes the final node representation ${\bf x}_{u}(T)$. This is visually summarized in Figure 1 when dampening and external forcing are excluded.

Moreover, we observe that the general formulation of the neighborhood aggregation function $\Phi_{\mathcal{G}}(\{{\bf x}_{v}(t)\}_{v\in\mathcal{N}_{u}})$ allows implementing any function that aggregates nodes (and edges) information. Therefore, $\Phi_{\mathcal{G}}(\{{\bf x}_{v}(t)\}_{v\in\mathcal{N}_{u}})$ allows enhancing a standard DGNs with our Hamiltonian conservation. As a demonstration of this, in Section 3, we experiment with two neighborhood aggregation functions, which are the classical GCN aggregation (Kipf and Welling, ) and

$$ $\Phi_{\mathcal{G}}(\{{\bf x}_{v}(t)\}_{v\in\mathcal{N}_{u}})=\sum_{v\in \mathcal{N}_{u}}{\bf V}{\bf x}_{v}(t).$ (6) $$

Further details about H-DGN discretization are in Appendix A.3.

#### H-DGN allows long-range propagation.

We show that our H-DGN in Eq. 5 adheres to the laws of conservation, allowing long-range propagation in the message-passing flow.

As discussed in (Haber and Ruthotto, ; Gravina et al., ), non-dissipative (i.e., long-range) propagation is directly linked to the sensitivity of the solution of the ODE to its initial condition, thus to the stability of the system. Such sensitivity is controlled by the Jacobian’s eigenvalues of Eq. 5. Under the assumption that the Jacobian varies sufficiently slow over time and its eigenvalues are purely imaginary, then the initial condition is effectively propagated into the final node representation, making the system both stable and non-dissipative, thus allowing for long-range propagation.

###### Theorem 2.1 .

The Jacobian matrix of the system defined by the ODE in Eq. 5 possesses eigenvalues purely on the imaginary axis, i.e.,

$$ $\operatorname{Re}\left(\lambda_{i}\left(\pdv{}{{\bf x}_{u}}\mathcal{J}_{u} \nabla_{{\bf x}_{u}}H_{\mathcal{G}}({\bf y}(t))\right)\right)=0,\quad\forall i,$ $$

where $\lambda_{i}$ represents the $i$-th eigenvalue of the Jacobian.

We report the proof in Appendix B.1. Then, we take a further step and strengthen such result by proving that the nonlinear vector field defined by H-DGN is divergence-free, thus preserving information within the graph during the propagation process and helping to maintain informative node representations.
In other words, the H-DGN’s dynamics possess a non-dissipative behavior independently of both the assumption regarding the slow variation of the Jacobian and the position of the Jacobian eigenvalues on the complex plane.

###### Theorem 2.2 .

The autonomous Hamiltonian $H_{\mathcal{G}}$ of the system in Eq. 5 with learnable weights shared across time stays constant at the energy level specified by the initial value $H_{\mathcal{G}}({\bf y}(0))$, i.e.,

$$ $\odv{H_{\mathcal{G}}}{t}=0.$ (7) $$

H-DGN also possesses a divergence-free nonlinear vector field

$$ $\nabla\cdot\mathcal{J}_{u}\nabla_{\mathbf{x}_{u}}H_{\mathcal{G}}(\mathbf{y}(t) )=0,\quad t\in[0,T].$ (8) $$

See proof in Appendix B.2. This allows us to interpret the system dynamics as purely rotational, without energy loss, and demonstrates that H-DGN is governed by conservation laws.

We now provide a sensitivity analysis, following Chang et al. (); Gravina et al. (); Galimberti et al. (), to prove that H-DGN effectively allows for long-range information propagation. Specifically, we measure the sensitivity of a node state after an arbitrary time $T$ of the information propagation with respect to its previous state, $\|\partial{\bf x}_{u}(T)/\partial{\bf x}_{u}(T-t)\|$. In other words, we compute the backward sensitivity matrix (BSM). We now provide a theoretical bound of our H-DGN, with its proof in Appendix B.3.

###### Theorem 2.3 .

Considering the continuous system defined by Eq. 5, the backward sensitivity matrix (BSM) is bounded from below:

$$ $\left\|\pdv{\mathbf{x}_{u}(T)}{\mathbf{x}_{u}(T-t)}\right\|\geq 1,\quad\forall t \in[0,T].$ $$

The result of Theorem 2.3 indicates that the gradients in the backward pass do not vanish, enabling the effective propagation of previous node states through successive transformations to the final nodes’ representations. Therefore, H-DGN has a conservative message passing, where the final representation of each node retains its complete past.
We observe that Theorem 2.3 holds even during discretization when the Symplectic Euler method is employed (see Appendix A.3).

To give the full picture of the time dynamics of the gradients, we present a similar analysis and provide an upper bound of the BSM in Theorem A.1.

###### Theorem 2.1 .

###### Theorem 2.2 .

###### Theorem 2.3 .

#### Introducing dissipative components.

A purely conservative Hamiltonian inductive bias forces the node states to follow trajectories that maintain constant energy, potentially limiting the effectiveness of the DGN on downstream tasks by restricting the system’s ability to model all complex nonlinear dynamics.
To this end, we complete the formalization of our port-Hamiltonian framework by introducing tools from mechanical systems, such as friction and external control, to learn how much the dynamic should deviate from this purely conservative behavior. Therefore, we extend the dynamics in Eq. 5 to a *port*-Hamiltonian by including two new terms $D({\bf q})\in\mathbb{R}^{d/2\times d/2}$ and $F({\bf q},t)\in\mathbb{R}^{d/2}$, i.e.,

$$ $\odv{\mathbf{x}_{u}(t)}{t}=\left[\mathcal{J}_{u}-\begin{pmatrix}D({\bf q}(t))& \bf 0\\ \bf 0&\bf 0\end{pmatrix}\right]\nabla_{\mathbf{x}_{u}}H_{\mathcal{G}}({\bf y}( t))+\begin{pmatrix}F({\bf q}(t),t)\\ \bf 0\end{pmatrix},\quad\forall u\in\mathcal{V}.$ (9) $$

Depending on the definition of $D({\bf q}(t))$ we can implement different forces. Specifically, if $D({\bf q}(t))$ is positive semi-definite then it implements internal dampening, while a negative semi-definite implementation leads to internal acceleration. A mixture of dampening and acceleration is obtained otherwise. In the case of dampening, the energy is decreased along the flow of the system (Van der Schaft, ).
To further enhance the modeling capabilities, we integrate the learnable state- and time-dependent external force $F({\bf q}(t),t)$, which further drives node representation trajectories. Figure 1 visually summarizes how such tools can be plugged in our framework during node update.

Although $D({\bf q}(t))$ and $F({\bf q}(t),t)$ can be implemented as static (fixed) functions, in our experiments in Section 3 we employ neural networks to learn such terms. We provide additional details on the specific architectures in Appendix A.2.
In the following, we refer to a DGN following Eq. 9 as *port-Hamiltonian Deep Graph Network* (PH-DGN) to distinguish it from the purely conservative H-DGN in Eq. 5. We provide further details about the discretization of PH-DGN in Appendix A.3.

## 3 Experiments

We empirically verify both theoretical claims and practical benefits of our framework on popular graph benchmarks for long-range propagation. First (Section 3.1), we conduct a controlled synthetic test showing non-vanishing gradients even when thousands of layers are used. Afterward (Section 3.2), we run a graph transfer task inspired by Di Giovanni et al. () to assess the efficacy in preserving long-range information between nodes. Then, we assess our framework in popular benchmark tasks requiring the exchange of messages at large distances over the graph, including graph property prediction (Section 3.3) and the long-range graph benchmark (Dwivedi et al., ) (Section 3.4). We compare our performance to state-of-the-art methods, such as MPNN-based models, DE-DGNs (which represent a direct competitor to our method), higher-order DGNs, and graph transformers. Additional details on literature methods are in Appendix C.1. We investigate two neighborhood aggregation functions for our H-DGN and PH-DGN, which are the classical GCN aggregation and that proposed in Eq. 6. Our models are implemented in PyTorch (Paszke et al., ) and PyTorch-Geometric (Fey and Lenssen, ). We release openly the code implementing our methodology and reproducing our empirical analysis upon acceptance. Our experimental results were obtained using NVIDIA A100 GPUs and Intel Xeon Gold 5120 CPUs.

### 3.1 Numerical Simulations

Setup.
We empirically verify that our theoretical considerations on H-DGN hold true by an experiment requiring to propagate information within a Carbon-60 molecule graph without training on any specific task, i.e., we perform no gradient update step. While doing so, we measure the energy level captured in $H_{\mathcal{G}}({\bf y}(\ell\epsilon))$ in the forward pass and the sensitivity, $\|\partial{{\bf x}_{u}^{(L)}}/\partial{{\bf x}_{u}^{(\ell)}}\|$, from each intermediate layer $\ell=1,\dots,L$ in the backward pass.
We consider the 2-d position of the atom in the molecule as the input node features, fixed terminal propagation time $T=10$ with various integration step sizes $\epsilon\in\{0.1,0.01,0.001\}$ and $T=300$ with $\epsilon=0.3$. Note that the corresponding number of layers is computed as $L=T/\epsilon$, i.e., we use tens to thousands of layers. For the ease of the simulation, we use $\tanh$-nonlinearity, fixed learnable weights that are randomly initialized, and the aggregation function in Eq. 6.

Results.
In Figure 2(a), we show the energy difference $H_{\mathcal{G}}({\bf y}(\ell\epsilon))-H_{\mathcal{G}}({\bf y}(0))$ for different step sizes. For a fixed time $T$, a smaller step size $\epsilon$ is related to a higher number of stacked layers. We note that the energy difference oscillates around zero, and the smaller the step size the more accurately the energy is preserved. This supports our intuition of H-DGN being a discretization of a divergence-free continuous Hamiltonian dynamic, that allows for non-dissipative forward propagation, as stated in Theorem 2.1 and Theorem 2.2. Even for larger step sizes, energy is neither gained nor lost.

Regarding the backward pass, Figures 2(b), 2(c) assert that the lower bound $\|\partial{{\bf x}(L)}/\partial{{\bf x}(\ell)}\|\geq 1$ stated in Theorem 2.3 and its discrete version in Theorem A.2 leads to non-vanishing gradients. In particular, Figure 2(c) shows a logarithmic-linear increase of sensitivity with respect to the distance to the final layer, hinting at the exponential upper bound derived in Theorem A.1. This growing behavior can be controlled by regularizing the weight matrices, or by use of normalized aggregation functions, as in GCN (Kipf and Welling, ).

Figure: (a)
Refer to caption: extracted/5623546/figures/numerics/numerics_legend.png

### 3.2 Graph Transfer

Setup.
We build on the graph transfer tasks by Di Giovanni et al. () and consider the problem of propagating a label from a source node to a target node located at increasing distances $k$ in the graph. We use the same graph distributions as in the original work, i.e., line, ring, and crossed-ring graphs. To increase the difficulty of the task, we randomly initialize intermediate nodes with a feature uniformly sampled in $[0,0.5)$. Source and destination nodes are initialized with labels “1” and “0”, respectively. We considered problems at distances $k\in\{3,5,10,50\}$, thus requiring incrementally higher efficacy in propagating long-range information.
Given the conservative nature of the task, we focus on assessing the purely Hamiltonian H-DGN model.
More details on the task and the hyperparameters can be found in the Appendix C.2 and C.5.

Results.
Figure 3 reports the test mean-squared error (and std. dev.) of H-DGN compared to literature models. It appears that classical MPNNs do not effectively propagate information across long ranges as their performance decreases when $k$ increases. Differently, H-DGN achieves low errors even at higher distances, i.e., $k\geq 10$. The only competitor to our H-DGN is A-DGN, which is another non-dissipative method. Overall, H-DGN outperforms all the classical MPNNs baseline while having on average better performance than A-DGN, thus empirically supporting our claim of long-range capabilities while introducing a new architectural bias. Moreover, our results highlight how our framework can push simple graph convolutional architectures to
state-of-the-art performance when imbuing them with dynamics capable of long-range message exchange.

Figure: Figure 3: Information transfer performance on (a) Line, (b) Ring, and (c) Crossed-Ring graphs. Overall, baseline approaches are not able to transfer the information accurately as distance increase, while non-dissipative methods like A-DGN and our H-DGN achieve low errors.
Refer to caption: extracted/5623546/figures/GraphTransfer_logscale_merged.png

### 3.3 Graph Property Prediction

#### Setup.

We consider three graph property prediction tasks introduced in Corso et al. () under the experimental setting of Gravina et al. (). We investigate the performance of our (port-) Hamiltonian framework in predicting graph diameters, single source shortest paths (SSSP), and node eccentricity on synthetic graphs. Note that to effectively solve these tasks, it
is fundamental to propagate not only the information of direct neighborhoods, but also the information coming from nodes far away in the graph. Therefore, good performance in these tasks highlights long-range propagation capabilities.
In this experiment, we investigate the performance of our complete port-Hamiltonian framework, PH-DGN, and present the pure Hamiltonian H-DGN as an ablation study. More details on the task and the hyperparameters can be found in Appendix C.3 and C.5.

Results.
We present the results
on the graph property prediction tasks in Table 1, reporting $log_{10}(\text{MSE})$ as evaluation metric.
We observe that both our H-DGN and PH-DGN show a strong improvement with respect to baseline methods, achieving new state-of-the-art performance on all the tasks.
Indeed, our ablation reveals that the purely conservative H-DGN model has, on average, a $log_{10}(\text{MSE})$ that is 0.33 lower than the best baseline.
Such gap is pushed to 0.81 when the full port-Hamiltonian bias (i.e., PH-DGN) is employed, marking a significant decrease in the test loss. The largest gap is achieved by PH-DGN in the eccentricity task, where it improves the $log_{10}(\text{MSE})$ performance of the best baseline by 1.36.

Although our purely conservative H-DGN shows improved performance with respect to all baselines, it appears that relaxing such bias via PH-DGN is more beneficial overall,

**Table 1: Mean test $log_{10}(\text{MSE})$ and std average over 4 training seeds on the Graph Property Prediction. Our methods and DE-DGN baselines are implemented with weight sharing. The first, second, and third best scores are colored.**
| Model | Eccentricity | Diameter | SSSP |
| --- | --- | --- | --- |
| MPNNs |  |  |  |
| GCN | 0.8468±0.0028 | 0.7424±0.0466 | 0.9499${}_{\pm 9.2\cdot 10^{-5}}$ |
| GAT | 0.7909±0.0222 | 0.8221±0.0752 | 0.6951±0.1499 |
| GraphSAGE | 0.7863±0.0207 | 0.8645±0.0401 | 0.2863±0.1843 |
| GIN | 0.9504±0.0007 | 0.6131±0.0990 | -0.5408±0.4193 |
| GCNII | 0.7640±0.0355 | 0.5287±0.0570 | -1.1329±0.0135 |
| DE-DGNs |  |  |  |
| DGC | 0.8261±0.0032 | 0.6028±0.0050 | -0.1483±0.0231 |
| GraphCON | 0.6833±0.0074 | 0.0964±0.0620 | -1.3836±0.0092 |
| GRAND | 0.6602±0.1393 | 0.6715±0.0490 | -0.0942±0.3897 |
| A-DGN | 0.4296±0.1003 | -0.5188±0.1812 | -3.2417±0.0751 |
| Ours |  |  |  |
| H-DGN | -0.7248±0.1068 | -0.5473±0.1074 | -3.0467±0.1615 |
| PH-DGN | -0.9348±0.2097 | -0.5385±0.0187 | -4.2993±0.0721 |

leading to even greater improvements in long-range information propagation.
Our intuition is that such tasks do not require purely conservative behavior since nodes need to count distances while exchanging
more messages with other nodes, similar to standard algorithmic solutions such as Dijkstra (). Therefore, the energy may not be constant during the resolution of the task, hence benefiting from the non-purely conservative behavior of PH-DGN.

As for the graph transfer task, our results demonstrate that our PH-DGNs can effectively learn and exploit long-range information while pushing simple graph neural architectures to state-of-the-art performance when modeling dynamics capable of long-range propagation.

### 3.4 Long-Range Graph Benchmark

**Table 2: Results for Peptides-func and Peptides-struct averaged over 3 training seeds. The first, second, and third best scores are colored. Extended version of this table is provided in Appendix D.2. Note that all MPNN-based methods include structural and positional encoding.**
| Model | Peptides-func | Peptides-struct |
| --- | --- | --- |
| AP $\uparrow$ | MAE $\downarrow$ |  |
| MPNNs |  |  |
| GCN | 0.6860±0.0050 | 0.2460±0.0007 |
| GCNII | 0.5543±0.0078 | 0.3471±0.0010 |
| GINE | 0.6621±0.0067 | 0.2473±0.0017 |
| GatedGCN | 0.6765±0.0047 | 0.2477±0.0009 |
| Multi-hop DGNs |  |  |
| DIGL+MPNN+LapPE | 0.6830±0.0026 | 0.2616±0.0018 |
| MixHop-GCN+LapPE | 0.6843±0.0049 | 0.2614±0.0023 |
| DRew-GCN+LapPE | 0.7150±0.0044 | 0.2536±0.0015 |
| Transformers |  |  |
| Transformer+LapPE | 0.6326±0.0126 | 0.2529±0.0016 |
| SAN+LapPE | 0.6384±0.0121 | 0.2683±0.0043 |
| GraphGPS+LapPE | 0.6535±0.0041 | 0.2500±0.0005 |
| DE-DGNs |  |  |
| GRAND | 0.5789±0.0062 | 0.3418±0.0015 |
| GraphCON | 0.6022±0.0068 | 0.2778±0.0018 |
| A-DGN | 0.5975±0.0044 | 0.2874±0.0021 |
| Ours |  |  |
| H-DGN | 0.6961±0.0070 | 0.2581±0.0020 |
| PH-DGN | 0.7012±0.0045 | 0.2465±0.0020 |

#### Setup.

We assess the performance of our method on the real-world long-range graph benchmark (LRGB) from Dwivedi et al. (), focusing on the Peptide-func and Peptide-struct tasks. As in Section 3.3, we decouple our method into H-DGN and PH-DGN to provide an ablation study on the strictly conservative behavior. For our evaluation, we follow the experimental setting in Dwivedi et al. ().
Acknowledging the results from Tönshoff et al. (), we also report results with a 3-layer MLP readout. More details on the task and the hyperparameters can be found in Appendix C.3 and C.5.

Results.
We report results on the LRGB tasks in Table 2 (extended results are reported in Appendix D.2 due to space limits). Our results show that both H-DGN and PH-DGN outperform classical MPNNs, graph transformers, most of the multi-hop DGNs, and recent DE-DGNs (which represent a direct competitor to our method).
Overall, our (port-)Hamiltonian methods show great benefit in propagating long-range information without requiring additional strategies such as global position encoding, global attention mechanisms, or rewiring techniques that increase the overall complexity of the method. Consequently, our results reaffirm the effectiveness of our (port-)Hamiltonian framework in enabling efficient long-range propagation, even in simple DGNs characterized by purely local message exchanges.

## 4 Related works

DGN based on differential equations.
Recent advancements in the field of representation learning have introduced new architectures
that establish a connection between neural networks and dynamical systems.
Inspired by pioneering works on recurrent neural networks (Chen et al., ; Haber and Ruthotto, ; Chang et al., ; Galimberti et al., ), such connection has been pushed to the domains of DGNs (Han et al., ).
Indeed, works like GDE (Poli et al., ), GRAND (Chamberlain et al., ), PDE-GCN (Eliasof et al., ), DGC (Wang et al., ), GRAND++ (Thorpe et al., ) propose to interpret DGNs as discretization of ODEs and PDEs. The conjoint of dynamical systems and DGNs have found favorable consensus, as these new methods exploit the intrinsic properties of differential equations to extend the
characteristic of message passing
within DGNs. GRAND, GRAND++, and DGC bias the node representation trajectories to follow the heat diffusion process, thus performing a gradual smoothing of the initial node states. On the contrary, GraphCON (Rusch et al., ) used oscillatory properties to enable linear dynamics that preserve the Dirichlet energy encoded in the node features; PDE-GCN${{}_{\text{M}}}$ (Eliasof et al., ) uses an interpolation between anisotropic diffusion and conservative oscillatory properties; and HamGNN (Kang et al., ) leverages Hamiltonian dynamics to encode node input features, which are then fed into classical DGNs to enhance their conservative properties.
Similarly to HamGNN, HANG (Zhao et al., ) leverages Hamiltonian dynamics to improve robustness to adversarial attacks to the graph structure. More recently, A-DGN (Gravina et al., ) introduces an anti-symmetric mechanism that lead to node-wise non-dissipative dynamics.
Differently from HamGNN and HANG, our PH-DGN is the first port-Hamiltonian framework, which provides theoretical guarantees of both conservative and non-conservative behaviors while seamlessly incorporating the most suitable aggregation function for the task at hand, enabling long-range information propagation.

Long-range propagation on DGNs.
Effectively transferring information across distant nodes is still an open challenge in the graph representation learning community (Shi et al., ). Various strategies have been explored in recent years to address this challenge, such as graph rewiring. In this setting, methods like SDRF (Topping et al., ), GRAND (Chamberlain et al., ), BLEND (Chamberlain et al., ), and DRew (Gutteridge et al., ) (dynamically) alter the original edge set to densify the graph during preprocessing to facilitate node communication. Differently, Transformer-based methods (Shi et al., ; Dwivedi and Bresson, ; Ying et al., ; Wu et al., ) enable message passing between all node pairs. FLODE (Maskey et al., ) incorporates non-local dynamics by using a fractional power of the graph shift operator. Although these techniques are effective in addressing the problem of long-range communication, they can also increase the complexity of information propagation due to denser graph shift operators.

## 5 Conclusions

In this paper, we have presented (port-)Hamiltonian Deep Graph Network (PH-DGN), a general framework that gauges the equilibrium between non-dissipative long-range propagation and non-conservative behavior while seamlessly incorporating the most suitable neighborhood aggregation function. We theoretically prove that, when pure Hamiltonian dynamic is employed, both the continuous and discretized versions of our framework allow for long-range propagation in the message passing flow since node states retain their past. To demonstrate the benefits of including (port-) Hamiltonian dynamics in DE-DGNs, we conducted several experiments on synthetic and real-world benchmarks requiring long-range interaction. Our results show that our method outperforms state-of-the-art models and that the inclusion of data-driven forces that deviate from a purely conservative behavior is often key to maximize efficacy of the approach on tasks requiring long-range propagation.
Indeed, in practice, effective information propagation requires a balance between long-term memorization and propagation and the ability to selectively discard and forget information when necessary.
Looking ahead to future developments, our port-Hamiltonian dynamic can be extended to handle time-varying streams of graphs and can be evaluated with alternative discretization methods, e.g., adaptive multistep schemes (Rufai et al., ).

## Appendix A Additional Details of the (port-)Hamiltonian Framework

### A.1 Sensitivity Upper Bound

Although the sensitivity of a node state after a time $t$ with respect to its previous state can be bounded from below, allowing effective conservative message passing in H-DGN, we observe that it is possible to compute an upper bound on such a measure, which we provide in the following theorem. While the theorem shows that, theoretically, the sensitivity measure may grow (i.e., potentially causing gradient explosion), we emphasize that during our experiments we did not encounter such a problem.

###### Theorem A.1 .

Consider the continuous system defined by Eq. 5, if $\sigma$ is a non-linear function with bounded derivative, i.e. $\exists M>0,|\sigma^{\prime}(x)|\leq M$, and the neighborhood aggregation function is of the form $\Phi_{\mathcal{G}}=\sum_{v\in\mathcal{N}_{u}}\mathbf{V}{\bf x}_{v}$, the backward sensitivity matrix (BSM) is bounded from above:

$$ $\left\|\pdv{\mathbf{x}_{u}(T)}{\mathbf{x}_{u}(T-t)}\right\|\leq\sqrt{d}\,\text {exp}(QT),\quad\forall t\in[0,T],$ $$

where $Q=\sqrt{d}\,M\|\mathbf{W}\|_{2}^{2}+\sqrt{d}\,M\text{max}_{i\in[n]}|\mathcal{N
}_{i}|\|\mathbf{V}\|_{2}^{2}$.

We give the proof of this theorem in Appendix B.4.

###### Theorem A.1 .

### A.2 Architectural Details of Dissipative Components

As typically employed, we follow physics-informed approaches that learn how much dissipation and external control is necessary to model the observations [Desai et al., ]. In particular, we consider these (graph-) neural network architectures for the dampening term $D(\mathbf{q})$ and external force term $F(\mathbf{q},t)$, assuming for simplicity $\mathbf{q}_{u}\in\mathbb{R}^{\frac{d}{2}}$.

Dampening $D(\mathbf{q})$: it is a square $\frac{d}{2}\times\frac{d}{2}$ matrix block with only diagonal entries being non-zero and defined as:

- •
*param*: a learnable vector $\mathbf{w}\in\mathbb{R}^{\frac{d}{2}}$.
- •
*param+*: a learnable vector followed by a ReLU activation, i.e., $\text{ReLU}(\mathbf{w})\in\mathbb{R}^{\frac{d}{2}}$.
- •
*MLP4-ReLU*: a $4$-layer MLP with ReLU activation and all layers of dimension $\frac{d}{2}$.
- •
*DGN-ReLU*: a DGN node-wise aggregation layer from Eq. 6 with ReLU activation.

External forcing $F(\mathbf{q},t)$: it is a $\frac{d}{2}$ dimensional vector where each component is the force on a component of the system. Since it takes as input $\frac{d}{2}+1$ components it is defined as:

- •
*MLP4-Sin*: 3 linear layers of $\frac{d}{2}+1$ units with sin activation followed by a last layer with $\frac{d}{2}$ units.
- •
*DGN-tanh*: a single node-wise DGN aggregation from Eq. 6
followed by a tanh activation.

Note that dampening, i.e., energy loss, is only given when $D(\mathbf{q})$ represents a positive semi-definite matrix. Hence, we used ReLU-activation, except for *param*, which offers a flexible trade-off between dampening and acceleration learned by backpropagation.

### A.3 Discretization of (port-)Hamiltonian DGNs

As for standard DE-DGNs a numerical discretization method is needed to solve Eq. 5. However, as observed in Haber and Ruthotto [], Galimberti et al. [], not all standard techniques can be employed for solving Hamiltonian systems. Indeed, symplectic integration methods need to be used to preserve the conservative properties in the discretized system.

For the ease of simplicity, in the following we focus on the Symplectic Euler method, however, we observe that more complex methods such as Strömer-Verlet can be employed [Hairer et al., ].

The Symplectic Euler scheme, applied to our H-DGN in Eq. 5, updates the node representation at the $(\ell+1)$-th step as

$$ ${\bf x}_{u}^{(\ell+1)}=\begin{pmatrix}{\bf p}^{(\ell+1)}_{u}\[5.69046pt] {\bf q}^{(\ell+1)}_{u}\end{pmatrix}=\begin{pmatrix}{\bf p}^{(\ell)}_{u}\[5.69 046pt] {\bf q}^{(\ell)}_{u}\end{pmatrix}+\epsilon\mathcal{J}_{u}\begin{pmatrix}\nabla _{{\bf p}_{u}}H_{\mathcal{G}}({\bf p}^{(\ell+1)},{\bf q}^{(\ell)})\[5.69046pt ] \nabla_{{\bf q}_{u}}H_{\mathcal{G}}({\bf p}^{(\ell+1)},{\bf q}^{(\ell)})\end{ pmatrix},\quad\forall u\in\mathcal{V}.$ (10) $$

with $\epsilon$ the step size of the numerical discretization.
We note that Eq. 10 relies on both the current and future state of the system, hence marking an implicit scheme that would require solving a linear system of equations in each step. To obtain an explicit version of Eq. 10, we consider the neighborhood aggregation function in Eq. 6 and impose a structure assumption on ${\bf W}$ and ${\bf V}$, namely ${\bf W}=\begin{pmatrix}{\bf W}_{p}&\bf 0\\
\bf 0&{\bf W}_{q}\end{pmatrix}$ and ${\bf V}=\begin{pmatrix}{\bf V}_{p}&\bf 0\\
\bf 0&{\bf V}_{q}\end{pmatrix}$. We note that a comparable assumption can be made for other neighborhood aggregation functions, such as GCN aggregation.

Therefore, the gradients in Eq. 10 can be rewritten in the explicit form as

$$ $\displaystyle{\bf p}^{(\ell+1)}_{u}$ $\displaystyle={\bf p}^{(\ell)}_{u}-\epsilon\Biggl{[}{\bf W}_{q}^{\top}\sigma({ \bf W}_{q}{\bf q}^{(\ell)}_{u}+\Phi_{\mathcal{G}}(\{{\bf q}^{(\ell)}_{v}\}_{v \in\mathcal{N}_{u}})+{\bf b}_{q})$ $\displaystyle\hskip 85.35826pt+\sum_{v\in\mathcal{N}_{u}\setminus\{u\}}{{\bf V }_{q}^{\top}}\sigma({\bf W}_{q}{\bf q}^{(\ell)}_{v}+\Phi_{\mathcal{G}}(\{{\bf q }^{(\ell)}_{j}\}_{j\in\mathcal{N}_{v}})+{\bf b}_{q})\Biggr{]}$ (11) $\displaystyle{\bf q}^{(\ell+1)}_{u}$ $\displaystyle={\bf q}^{(\ell)}_{u}+\epsilon\Biggl{[}{\bf W}_{p}^{\top}\sigma({ \bf W}_{p}{\bf p}^{(\ell+1)}_{u}+\Phi_{\mathcal{G}}(\{{\bf p}^{(\ell+1)}_{v}\} _{v\in\mathcal{N}_{u}})+{\bf b}_{p})$ $\displaystyle\hskip 85.35826pt+\sum_{v\in\mathcal{N}_{u}\setminus\{u\}}{{\bf V }_{p}^{\top}}\sigma({\bf W}_{p}{\bf p}^{(\ell+1)}_{v}+\Phi_{\mathcal{G}}(\{{ \bf p}^{(\ell+1)}_{j}\}_{j\in\mathcal{N}_{v}})+{\bf b}_{p})\Biggr{]}.$ (12) $$

We observe that Eqs. 11 and 12 can be understood as coupling
two DGN layers. This discretization mechanism is visually summarized in the middle of Figure 1 where a message-passing step from layer $\ell$ to layer $\ell+1$ is performed.

In the case of PH-DGN in Eq. 9 the discretization employs the same step for $\mathbf{q}^{(\ell+1)}$ in Eq. 12 while Eq. 11 includes the dissipative components, thus it can be rewritten as

$$ $\mathbf{p}_{u}^{(\ell+1)}=\mathbf{p}_{u}^{(\ell)}+\epsilon\Biggl{[}-\nabla_{ \mathbf{q}_{u}}{{H}_{\mathcal{G}}}(\mathbf{p}^{(\ell)},\mathbf{q}^{(\ell+1)}) \\ -D_{u}(\mathbf{q}^{(\ell+1)})\nabla_{\mathbf{p}_{u}}{{H}_{\mathcal{G}}}( \mathbf{p}^{(\ell)},\mathbf{q}^{(\ell+1)})+F_{u}(\mathbf{q}^{(\ell+1)},t) \Biggr{]}.$ (13) $$

Lastly, it is important to acknowledge that properties observed in the continuous domain may not necessarily hold in the discrete setting due to the limitations of the discretization method. In the following theorem, we show that when the Symplectic Euler method is employed, then Theorem 2.3 holds.

###### Theorem A.2 .

Considering the discretized system in Eq. 10 obtained by Symplectic Euler discretization, the backward sensitivity matrix (BSM) is bounded from below:

$$ $\left\|\pdv{\mathbf{x}_{u}^{(L)}}{\mathbf{x}_{u}^{(L-\ell)}}\right\|\geq 1, \quad\forall\ell\in[0,L].$ $$

We provide the proof in Appendix B.5.
Again, this indicates that even the discretized version of H-DGN enables for effective propagation and conservative message passing.

###### Theorem A.2 .

## Appendix B Proofs of the Theoretical Statements

In this section, we provide the proofs of the theoretical statements in the main text and in appendix A. As for the rest of the paper, we will use the denominator notation, i.e., Jacobian matrices have output components on columns and input components on the rows.

### B.1 Proof of Theorem 2.1

###### Proof.

First, we note that

$$ $\pdv{}{{\bf x}_{u}}\mathcal{J}_{u}\nabla_{{\bf x}_{u}}H_{\mathcal{G}}(\mathbf{ y}(t))=\nabla^{2}_{{\bf x}_{u}}H_{\mathcal{G}}(\mathbf{y}(t))\mathcal{J}_{u}^{ \top},$ (14) $$

where $\nabla^{2}_{{\bf x}_{u}}H_{\mathcal{G}}$ is the symmetric Hessian matrix.
Hence, the Jacobian is shortly written as $\mathbf{AB}$, where $\mathbf{A}$ is symmetric and $\mathbf{B}$ is anti-symmetric. Consider an eigenpair of $\mathbf{AB}$, where the eigenvector is denoted by $\mathbf{v}$ and the eigenvalue by $\lambda\neq 0$. Then:

$$ $\displaystyle\mathbf{v}^{*}\mathbf{AB}$ $\displaystyle=\lambda\mathbf{v}^{*}$ $\displaystyle\mathbf{v}^{*}\mathbf{A}$ $\displaystyle=\lambda\mathbf{v}^{*}\mathbf{B}^{-1}$ $\displaystyle\mathbf{v}^{*}\mathbf{A}\mathbf{v}$ $\displaystyle=\lambda\left(\mathbf{v}^{*}\mathbf{B}^{-1}\mathbf{v}\right)$ $$

where $*$ represents the conjugate transpose. On the left-hand side, it is noticed that the $\left(\mathbf{v}^{*}\mathbf{A}\mathbf{v}\right)$ term is a real number. Recalling that $\mathbf{B}^{-1}$ remains anti-symmetric and for any real anti-symmetric matrix $\mathbf{C}$ it holds that $\mathbf{C}^{*}=\mathbf{C}^{\top}=-\mathbf{C}$, it follows that $\left(\mathbf{v}^{*}\mathbf{C}\mathbf{v}\right)^{*}=\mathbf{v}^{*}\mathbf{C}^{
*}\mathbf{v}=-\mathbf{v}^{*}\mathbf{C}\mathbf{v}$. Hence, the $\mathbf{v}^{*}\mathbf{B}^{-1}\mathbf{v}$ term on the right-hand side is an imaginary number. Thereby, $\lambda$ needs to be purely imaginary, and, as a result, all eigenvalues of $\mathbf{AB}$ are purely imaginary.
∎

###### Proof.

### B.2 Proof of Theorem 2.2

Our H-DGN has shared weights $\mathbf{W},\mathbf{V}$ across the layers of the DGN. This means that the Hamiltonian is autonomous and does not depend explicitly on time $H_{\mathcal{G}}(\mathbf{y}(t),t)=H_{\mathcal{G}}(\mathbf{y}(t))$ as we can see from Eq. 4. In such case, the energy is naturally conserved in the system it represents: this is a consequence of the Hamiltonian flow and dynamics, as we show in this theorem. For a more in-depth description and analysis of Hamiltonian dynamics in general we refer to Arnold et al. [].

###### Proof.

The time derivative of $H(\mathbf{y}(t))$ is given by means of the chain-rule:

$$ $\odv{H(\mathbf{y}(t))}{t}={\pdv{H(\mathbf{y}(t))}{\mathbf{y}(t)}}\cdot\odv{ \mathbf{y}(t)}{t}={\pdv{H(\mathbf{y}(t))}{\mathbf{y}(t)}}\cdot\mathcal{J}\pdv{ H(\mathbf{y}(t))}{\mathbf{y}(t)}=0,$ (15) $$

where the last equality holds since $\mathcal{J}$ is anti-symmetric. Having no change over time implies that $H(\mathbf{y}(t))=H(\mathbf{y}(0))=\text{const}$ for all $t$.

Since the Hessian $\nabla^{2}H_{\mathcal{G}}(\mathbf{y}(t))$ is symmetric, it follows directly

$$ $\displaystyle\nabla\cdot\mathcal{J}_{u}\nabla_{{\bf x}_{v}}H_{\mathcal{G}}( \mathbf{y}(t))$ $\displaystyle=\sum_{i=1}^{d}{-\pdv[style-var=multipl]{H_{\mathcal{G}}(\mathbf{ y}(t))}{q_{v}^{i},p_{v}^{i}}+\pdv[style-var=multipl]{H_{\mathcal{G}}(\mathbf{y }(t))}{p_{v}^{i},q_{v}^{i}}}=0$ $$

∎

###### Proof.

### B.3 Proof of Theorem 2.3

In order to prove the lower bound on the BSM, we need a technical lemma that describes the time evolution of the BSM itself, which extends the result from Galimberti et al. [].

###### Lemma B.1 .

Given the system dynamics of the ODE in Eq. 2 governing the H-DGN, we have that

$$ $\odv{}{t}\pdv{\mathbf{y}(T)}{\mathbf{y}(T-t)}=\mathcal{J}\left.\pdv{H}{\mathbf {y}}\right|_{\mathbf{y}(T-t)}\pdv{\mathbf{y}(T)}{\mathbf{y}(T-t)}$ (16) $$

as in [Galimberti et al., ]. The same applies, with a slightly different formula, for each node $u$, that is the BSM satisfies

$$ $\odv{}{t}\pdv{\mathbf{x}_{u}(T)}{\mathbf{x}_{u}(T-t)}=\left.\pdv{\mathbf{y}}{ \mathbf{x}_{u}}\right\rvert_{(T-t)}\left.\pdv{f_{u}}{\mathbf{y}}\right\rvert_{ \mathbf{y}(T-t)}\pdv{\mathbf{x}_{u}(T)}{\mathbf{x}_{u}(T-t)}=F_{u}\pdv{\mathbf {x}_{u}(T)}{\mathbf{x}_{u}(T-t)}$ (17) $$

where $f_{u}$ is the restriction of $f=\mathcal{J}\pdv{H}{\mathbf{y}}$ to the components corresponding to $\mathbf{x}_{u}$, that is the dynamics of node $u$, which can be written as

$$ $f_{u}=L_{u}f=L_{u}\mathcal{J}\pdv{H}{\mathbf{y}}$ (18) $$

where $L_{u}$ is the readout matrix, of the form

$$ $L_{u}=\begin{bmatrix}0_{\frac{d}{2}\times\frac{d}{2}(u-1)}&I_{\frac{d}{2} \times\frac{d}{2}}&0_{\frac{d}{2}\times\frac{d}{2}(n-u)}&0_{\frac{d}{2}\times \frac{d}{2}(u-1)}&0_{\frac{d}{2}\times\frac{d}{2}}&0_{\frac{d}{2}\times\frac{d }{2}(n-u)}\\ 0_{\frac{d}{2}\times\frac{d}{2}(u-1)}&0_{\frac{d}{2}\times\frac{d}{2}}&0_{ \frac{d}{2}\times\frac{d}{2}(n-u)}&0_{\frac{d}{2}\times\frac{d}{2}(u-1)}&I_{ \frac{d}{2}\times\frac{d}{2}}&0_{\frac{d}{2}\times\frac{d}{2}(n-u)}\end{bmatrix}$ (19) $$

which is a projection on the coordinates of a single node $u$.
Notice as well that, in denominator notation

$$ $\pdv{\mathbf{y}}{\mathbf{x}_{u}}=L_{u}$ (20) $$

We only show Eq. 17, as it is related to graph networks and is actually a harder version of Eq. 16, with the latter being already proven in Galimberti et al. []. We also show the last part of the proof in a general sense, without using the specific matrices of the Hamiltonian used.

###### Proof.

Following Galimberti et al. [], the solution to the ODE $\odv{\mathbf{x}_{u}}{t}=f_{u}(\mathbf{y}(t))$ can be written in integral form as

$$ $\mathbf{x}_{u}(T)=\mathbf{x}_{u}(T-t)+\int_{T-t}^{T}f_{u}(\mathbf{y}(\tau)) \textnormal{d}\tau=\mathbf{x}_{u}(T-t)+\int_{0}^{t}f_{u}(\mathbf{y}(T-t+s)) \textnormal{d}s$ (21) $$

Differentiating by the solution at a previous time $\mathbf{x}_{u}(T-t)$ we obtain

$$ $\begin{split}\pdv{\mathbf{x}_{u}(T)}{\mathbf{x}_{u}(T-t)}&=I_{u}+\pdv{\int_{0} ^{\top}f_{u}(\mathbf{y}(T-t+s))\textnormal{d}s}{\mathbf{x}_{u}(T-t)}=\\ &=I_{u}+\int_{0}^{t}{\pdv{f_{u}(\mathbf{y}(T-t+s))}{\mathbf{x}_{u}(T-t)}}\\ &=I_{u}+\int_{0}^{t}{\pdv{\mathbf{y}(T-t+s)}{\mathbf{x}_{u}(T-t)}\left.\pdv{f_ {u}}{\mathbf{y}}\right\rvert_{\mathbf{y}(T-t+s)}\textnormal{d}s}\end{split}$ (22) $$

where in the second equality we brought the derivative term under the integral sign and in the third we used the chain rule of the derivative (recall we are using denominator notation). Considering a slight perturbation in time $\delta$, we consider $\pdv{\mathbf{x}_{u}(T)}{\mathbf{x}_{u}(T-t-\delta)}$ as this will be used to calculate the time derivative of the BSM. Using again the chain rule for the derivative and the formula above with $T-t-\delta$ instead of $T-t$, we have that

$$ $\begin{split}\pdv{\mathbf{x}_{u}(T)}{\mathbf{x}_{u}(T-t-\delta)}&=\pdv{\mathbf {x}_{u}(T-t)}{\mathbf{x}_{u}(T-t-\delta)}\pdv{\mathbf{x}_{u}(T)}{\mathbf{x}_{u }(T-t)}\\ &=\left(I_{u}+\int_{0}^{\delta}{\pdv{\mathbf{y}(T-t-\delta+s)}{\mathbf{x}_{u}( T-t-\delta)}\left.\pdv{f_{u}}{\mathbf{y}}\right\rvert_{\mathbf{y}(T-t-\delta+s )}\textnormal{d}s}\right)\pdv{\mathbf{x}_{u}(T)}{\mathbf{x}_{u}(T-t)}\end{split}$ (23) $$

This way, we have expressed $\pdv{\mathbf{x}_{u}(T)}{\mathbf{x}_{u}(T-t-\delta)}$ in terms of $\pdv{\mathbf{x}_{u}(T)}{\mathbf{x}_{u}(T-t)}$. To calculate our objective, we want to differentiate with respect to $\delta$. We first calculate the difference:

$$ $\pdv{\mathbf{x}_{u}(T)}{\mathbf{x}_{u}(T-t-\delta)}-\pdv{\mathbf{x}_{u}(T)}{ \mathbf{x}_{u}(T-t)}=\left(\int_{0}^{\delta}{\pdv{\mathbf{y}(T-t-\delta+s)}{ \mathbf{x}_{u}(T-t-\delta)}\left.\pdv{f_{u}}{\mathbf{y}}\right\rvert_{\mathbf{ y}(T-t-\delta+s)}\textnormal{d}s}\right)\pdv{\mathbf{x}_{u}(T)}{\mathbf{x}_{u} (T-t)}$ (24) $$

We can now divide by $\delta$ and take the limit $\delta\rightarrow 0$

$$ $\begin{split}&\lim_{\delta\rightarrow 0}\frac{1}{\delta}\left(\pdv{\mathbf{x}_ {u}(T)}{\mathbf{x}_{u}(T-t-\delta)}-\pdv{\mathbf{x}_{u}(T)}{\mathbf{x}_{u}(T-t )}\right)\\ &=\lim_{\delta\rightarrow 0}\left(\frac{1}{\delta}\int_{0}^{\delta}{\pdv{ \mathbf{y}(T-t-\delta+s)}{\mathbf{x}_{u}(T-t-\delta)}\left.\pdv{f_{u}}{\mathbf {y}}\right\rvert_{\mathbf{y}(T-t-\delta+s)}\textnormal{d}s}\right)\pdv{\mathbf {x}_{u}(T)}{\mathbf{x}_{u}(T-t)}\\ &=\left.\pdv{y}{\mathbf{x}_{u}}\right\rvert_{(T-t)}\left.\pdv{f_{u}}{y}\right \rvert_{\mathbf{y}(T-t)}\pdv{\mathbf{x}_{u}(T)}{\mathbf{x}_{u}(T-t)}\end{split}$ (25) $$

Where in the final equality we used the fundamental theorem of calculus. Finally

$$ $\odv{}{t}\pdv{\mathbf{x}_{u}(T)}{\mathbf{x}_{u}(T-t)}=\left.\pdv{\mathbf{y}}{ \mathbf{x}_{u}}\right\rvert_{(T-t)}\left.\pdv{f_{u}}{\mathbf{y}}\right\rvert_{ \mathbf{y}(T-t)}\pdv{\mathbf{x}_{u}(T)}{\mathbf{x}_{u}(T-t)}=F_{u}\pdv{\mathbf {x}_{u}(T)}{\mathbf{x}_{u}(T-t)}$ (26) $$

giving us the final result.
∎

We are now ready to prove Theorem 2.3. First, we calculate that

$$ $\pdv{f_{u}}{\mathbf{y}}=\pdv{}{y}\left(L_{u}\mathcal{J}\pdv{H}{\mathbf{y}} \right)=\pdv[order={2}]{H}{\mathbf{y}}\mathcal{J}^{\top}L_{u}^{\top}=S\mathcal {J}^{\top}L_{u}^{\top}$ (27) $$

so that $F_{u}=L_{u}S\mathcal{J}^{\top}L_{u}^{\top}$. This will be helpful in the following matrix calculations.

###### Proof.

For brevity, we call $\left[\pdv{\mathbf{x}_{u}(T)}{\mathbf{x}_{u}(T-t)}\right]=\Psi_{u}(T,T-t)$, which will be indicated simply as $\Psi_{u}$. When $t=0$, $\Psi_{u}$ is just the Jacobian of the identity map $\Psi_{u}(T,T)=I_{u}$ and the result $\Psi_{u}^{\top}\mathcal{J}_{u}\Psi_{u}=\mathcal{J}_{u}$ is true for $t=0$. Calculating the time derivative on $\Psi_{u}^{\top}\mathcal{J}_{u}\Psi_{u}$ we have that

$$ $\begin{split}\odv{}{t}[\Psi_{u}^{\top}\mathcal{J}_{u}\Psi_{u}]&=\dot{\Psi}_{u} ^{\top}\mathcal{J}_{u}\Psi_{u}+\Psi_{u}^{\top}\mathcal{J}_{u}\dot{\Psi}_{u}\\ &=(F_{u}\Psi_{u})^{\top}\mathcal{J}_{u}\Psi_{u}+\Psi_{u}^{\top}\mathcal{J}_{u} F_{u}\Psi_{u}=\\ &=\Psi_{u}^{\top}L_{u}\mathcal{J}S^{\top}L_{u}^{\top}\mathcal{J}_{u}\Psi_{u}+ \Psi_{u}^{\top}\mathcal{J}_{u}L_{u}S\mathcal{J}^{\top}L_{u}^{\top}\Psi\\ &=\Psi_{u}^{\top}\left(L_{u}\mathcal{J}SL_{u}^{\top}\mathcal{J}_{u}+\mathcal{J }_{u}L_{u}S\mathcal{J}^{\top}L_{u}^{\top}\right)\Psi_{u}\end{split}$ (28) $$

where in the second equality we used the result from lemma B.1.
We just need to show that the term in parentheses is zero so that the time derivative is zero. Using the relations $\mathcal{J}^{\top}L_{u}^{\top}=-L_{u}^{\top}\mathcal{J}_{u}$ and $J_{u}L_{u}=L_{u}\mathcal{J}$ we easily see that, finally

$$ $\odv{}{t}\left(\left[\pdv{\mathbf{x}_{u}(T)}{\mathbf{x}_{u}(T-t)}\right]^{\top }\mathcal{J}_{u}\left[\pdv{\mathbf{x}_{u}(T)}{\mathbf{x}_{u}(T-t)}\right] \right)=\Psi_{u}^{\top}\left(L_{u}\mathcal{J}SL_{u}^{\top}\mathcal{J}_{u}+L_{u }\mathcal{J}S(-L_{u}^{\top}\mathcal{J}_{u})\right)\Psi_{u}=0$ (29) $$

which means that $\left[\pdv{\mathbf{x}_{u}(T)}{\mathbf{x}_{u}(T-t)}\right]^{\top}\mathcal{J}_{u
}\left[\pdv{\mathbf{x}_{u}(T)}{\mathbf{x}_{u}(T-t)}\right]$ is constant and equal to $\mathcal{J}_{u}$ for all $t$, that is our thesis.
Now, the bound on the gradient follows by considering any sub-multiplicative norm $\|\cdot\|$:

$$ $\left\|\mathcal{J}_{u}\right\|=\left\|\left[\pdv{\mathbf{x}_{u}(T)}{\mathbf{x} _{u}(T-t)}\right]^{\top}\mathcal{J}_{k}\left[\pdv{\mathbf{x}_{u}(T)}{\mathbf{x }_{u}(T-t)}\right]\right\|\leq\left\|\pdv{\mathbf{x}_{u}(T)}{\mathbf{x}_{u}(T- t)}\right\|^{2}\left\|\mathcal{J}_{u}\right\|$ $$

and simplifying by $\left\|\mathcal{J}_{u}\right\|=1$.
∎

###### Lemma B.1 .

###### Proof.

###### Proof.

### B.4 Proof of Theorem A.1

To prove the upper bound, we use the following technical lemma:

###### Lemma B.2 ( Galimberti et al. [ 2023 ] ) .

Consider a matrix $\mathbf{A}\in\mathbb{R}^{n\times n}$ with columns $\mathbf{a}_{i}\in\mathbb{R}^{n}$, i.e., $\mathbf{A}=\left[\begin{array}[]{llll}\mathbf{a}_{1}&\mathbf{a}_{2}&\cdots&
\mathbf{a}_{n}\end{array}\right]$, and assume that $\left\|\mathbf{a}_{i}\right\|_{2}\leq$ $\gamma^{+}$ for all $i=1,\ldots,n$. Then, $\|\mathbf{A}\|_{2}\leq\gamma^{+}\sqrt{n}$.

This lemma gives a bound on the spectral norm of a matrix when its columns are uniformly bounded in norm.
Therefore, our proof strategy for Theorem A.1 lies in bounding each column of the BSM matrix.

###### Proof.

Consider the ODE in Eq. 26 from Lemma B.1 and split $\pdv{{\bf x}_{u}(T)}{{\bf x}_{u}(T-t)}$ into columns $\pdv{{\bf x}_{u}(T)}{{\bf x}_{u}(T-t)}=\left[\begin{array}[]{llll}\mathbf{z}_{
1}(\mathbf{t})&\mathbf{z}_{2}(t)&\ldots&\mathbf{z}_{d}(t)\end{array}\right]$. Then, Eq. 26 is equivalent to

$$ $\dot{\mathbf{z}}_{i}(t)=\mathbf{A}_{u}(T-t)\mathbf{z}_{i}(t),\quad t\in[0,T],i =1,2\ldots,d,$ (30) $$

subject to $\mathbf{z}_{i}(0)=e_{i}$, where $e_{i}$ is the unit vector with a single nonzero entry in position $i$. The solution of the linear system of ODEs in Eq. 30 is given by the integral equation

$$ $\mathbf{z}_{i}(t)=\mathbf{z}_{i}(0)+\int_{0}^{t}\mathbf{A}_{u}(T-s)\mathbf{z}_ {i}(s)ds,\quad t\in[0,T].$ (31) $$

By assuming that $\|\mathbf{A}_{u}(\tau)\|_{2}\leq Q$ for all $\tau\in[0,T]$, and applying the triangular inequality in Eq. 31, it is obtained that:

$$ $\left\|\mathbf{z}_{i}(t)\right\|_{2}\leq\left\|\mathbf{z}_{i}(0)\right\|_{2}+Q \int_{0}^{t}\left\|\mathbf{z}_{i}(s)\right\|_{2}ds=1+Q\int_{0}^{t}\left\| \mathbf{z}_{i}(s)\right\|_{2}ds,$ $$

where the last equality follows from $\left\|\mathbf{z}_{i}(0)\right\|_{2}=\left\|e_{i}\right\|_{2}=1$ for all $i=1,2,\ldots,d$. Then, applying the Gronwall inequality, it holds for all $t\in[0,T]$

$$ $\left\|\mathbf{z}_{i}(t)\right\|_{2}\leq\exp(QT).$ (32) $$

By applying Lemma B.2, the general bound follows.

Lastly, we characterize $Q$ by bounding the norm $\|\mathbf{A}_{u}(\tau)\|_{2}\;\forall\tau\in[0,T]$.
From Lemma B.1 $\mathbf{A}_{v}$ can be expressed as $\mathbf{A}_{u}=\mathbf{L}_{u}\mathbf{S}\mathcal{J}^{\top}\mathbf{L}_{u}^{\top}$, which is equivalently $\mathbf{A}_{u}=\nabla^{2}_{{\bf x}_{u}}{H}_{\mathcal{G}}(\mathbf{y})\mathcal{J
}_{u}^{\top}$, since $\mathcal{J}^{\top}\mathbf{L}_{u}^{\top}=\mathbf{L}_{v}^{\top}\mathcal{J}_{u}^{\top}$.
The Hessian $\nabla^{2}_{{\bf x}_{u}}{H}_{\mathcal{G}}(\mathbf{y})$ is of the form:

$$ $\displaystyle\nabla^{2}_{{\bf x}_{u}}H_{\mathcal{G}}(\mathbf{y})=\mathbf{W}^{ \top}\operatorname{diag}(\sigma^{\prime}(\mathbf{W}{\bf x}_{u}+\Phi_{u}+b)) \mathbf{W}+\sum_{v\in\mathcal{N}_{u}}\mathbf{V}^{\top}\operatorname{diag}( \sigma^{\prime}(\mathbf{W}{\bf x}_{v}+\Phi_{v}+b))\mathbf{V}.$ $$

After noting that $\|\nabla^{2}_{{\bf x}_{u}}{H}_{\mathcal{G}}(\mathbf{y})\mathcal{J}_{u}^{\top}
\|_{2}\leq\|\nabla^{2}_{{\bf x}_{u}}{H}_{\mathcal{G}}(\mathbf{y})\|_{2}\|
\mathcal{J}_{u}^{\top}\|_{2}$, the only varying part is the Hessian $\nabla^{2}_{{\bf x}_{u}}{H}_{\mathcal{G}}(\mathbf{y})$ since $\|\mathcal{J}_{u}^{\top}\|_{2}=1$.
By Lemma B.2 $\operatorname{diag}(\sigma^{\prime}(x))\leq\sqrt{d}\,M$ and noting that $\|\mathbf{X}^{\top}\|=\|\mathbf{X}\|$ for any square matrix $\mathbf{X}$, then

$$ $\left\|\nabla^{2}_{{\bf x}_{u}}H_{\mathcal{G}}(\mathbf{y})\right\|_{2}\leq \sqrt{d}\,M\,\|\mathbf{W}\|_{2}^{2}+\sqrt{d}\,M\,\operatorname{max}_{i\in[n]}| \mathcal{N}_{i}|\;\|\mathbf{V}\|_{2}^{2}\eqqcolon Q.$ $$

This also justifies our previous assumption that $\|\mathbf{A}_{u}(\tau)\|_{2}$ is bounded.
∎

###### Lemma B.2 ( Galimberti et al. [ 2023 ] ) .

###### Proof.

### B.5 Proof of Theorem A.2

###### Proof.

In the discrete case, since the semi-implicit Euler integration scheme is a symplectic method, it holds that:

$$ $\left[\pdv{{\bf x}_{u}^{(\ell)}}{{\bf x}_{u}^{(\ell-1)}}\right]^{\top}\mathcal {J}_{u}\left[\pdv{{\bf x}_{u}^{(\ell)}}{{\bf x}_{u}^{(\ell-1)}}\right]= \mathcal{J}_{u}$ (33) $$

Further, by using the chain rule and applying Eq. 33 iteratively we get:

$$ $\left[\pdv{{\bf x}_{u}^{(L)}}{{\bf x}_{u}^{(L-\ell)}}\right]^{\top}\mathcal{J} _{u}\left[\pdv{{\bf x}_{u}^{(L)}}{{\bf x}_{u}^{(L-\ell)}}\right]=\left[\prod_{ i=L-\ell}^{L-1}\pdv{{\bf x}_{u}^{(i+1)}}{{\bf x}_{u}^{(i)}}\right]^{\top} \mathcal{J}_{u}\left[\prod_{i=L-\ell}^{L-1}\pdv{{\bf x}_{u}^{(i+1)}}{{\bf x}_{ u}^{(i)}}\right]=\mathcal{J}_{u}$ $$

Hence, the BSM is symplectic at arbitrary depth and we can conclude the proof with:

$$ $\left\|\mathcal{J}_{u}\right\|=\left\|\left[\pdv{{\bf x}_{u}^{(L)}}{{\bf x}_{u }^{(L-\ell)}}\right]^{\top}\mathcal{J}_{u}\left[\pdv{{\bf x}_{u}^{(L)}}{{\bf x }_{u}^{(L-\ell)}}\right]\right\|\leq\left\|\pdv{{\bf x}_{u}^{(L)}}{{\bf x}_{u} ^{(L-\ell)}}\right\|^{2}\left\|\mathcal{J}_{u}\right\|.$ $$

∎

###### Proof.

## Appendix C Experimental Details

### C.1 Employed Baselines

In our experiments, the performance of our (port-)Hamiltonian DGNs is compared with various state-of-the-art DGN baselines from the literature. Specifically, we consider:

- •
classical MPNN-based methods, i.e., GCN [Kipf and Welling, ], GraphSAGE [Hamilton et al., ], GAT [Veličković et al., ], GatedGCN [Bresson and Laurent, ], GIN [Xu et al., ], GINE [Hu et al., ], and GCNII [Chen et al., ];
- •
DE-DGNs, i.e., DGC [Wang et al., ], GRAND [Chamberlain et al., ], GraphCON [Rusch et al., ], and A-DGN [Gravina et al., ];
- •
Graph Transformers, i.e., Transformer [Vaswani et al., , Dwivedi and Bresson, ], SAN [Kreuzer et al., ], and GraphGPS [Rampášek et al., ];
- •
Higher-Order DGNs, i.e., DIGL [Gasteiger et al., ], MixHop [Abu-El-Haija et al., ], and DRew [Gutteridge et al., ].

Note that since DE-DGNs belong to the same family as our proposed method, they are a direct competitor to our (port-)Hamiltonian DGNs.

### C.2 Graph Transfer Task

The graph transfer task consists of the problem of propagating a label from a source node to a target node located at increased distance $k$ in the graph. We develop this task based on the work of [Di Giovanni et al., ]. We employ the same graph distributions as in the original work, i.e., line, ring, and crossed ring graphs (see Figure 4 for a
visual exemplification of the three types of graphs when the distance between the source and target nodes is $5$). However, to make the task more challenging, we initialize node input features with random values uniformly sampled in the range $[0,0.5)$. In each graph, the source node is assigned with label “1”, and a target node is assigned with label “0”. We use an input dimension of $1$, and with a source-target distance equal to $3,5,10,$ and $50$. We generate $1000$ graphs for training, $100$ for validation, and $100$ for testing.

Figure: Figure 4: Three topologies for Graph Transfer. Left) Line. Center) Ring. Right) Crossed-Ring. The distance between source and target nodes is equal to 5. Nodes marked with S are source nodes, while the nodes with a T are target nodes.
Refer to caption: extracted/5623546/figures/rings.png

We design each model using three main components. First is the encoder, which maps the node input features into a latent hidden space. Second is the graph convolution component (i.e., H-DGN or other baselines). Third is the readout, which maps the convolution’s output into the output space. The encoder and readout have the same architecture across all models in the experiments.

We perform hyperparameter tuning via grid search, optimizing the Mean Squared Error (MSE) computed on the node features of the whole graph. We train the models using Adam optimizer for a maximum of 2000 epochs and early stopping with patience of 100 epochs on the validation loss. For each model configuration, we perform 4 training runs with different weight initialization and report the average and standard deviation of the results. We report in Table 3 the grid of hyperparameter employed by each model in our experiment.

### C.3 Graph Property Prediction

The tasks consist of predicting two node-level (i.e., single source shortest path and eccentricity) and one graph-level (i.e., graph diameter) properties on synthetic graphs generated by different graph distributions. In our experiments, we follow the data generation and experimental procedure outlined in Gravina et al. []. Therefore, graphs contains between 25 and 35 nodes and are randomly generated from multiple distributions. Nodes are randomly initialized with a value uniformly sampled in the range $[0,1)$, while target values represent single source shortest path, eccentricity, and graph diameter depending on the selected task. 5120 graphs are used as the training set, 640 as the validation set, and the rest as the test set.

As in the original work, each model is designed as three components, i.e., encoder, graph convolution component (H-DGN, PH-DGN, or baselines), and readout.
We perform hyperparameter tuning via grid search, optimizing the Mean Square Error (MSE), training the models using Adam optimizer for a maximum of 1500 epochs, and early stopping with patience of 100 epochs on the validation error. For each model configuration, we perform 4 training runs with different weight initialization and report the average of the results. We report in Appendix C.5 the grid of hyperparameters employed by each model in our experiment.

### C.4 Long-Range Graph Benchmark

We consider the Peptide-func and Peptide-struct tasks from the Long-Range Graph Benchmark [Dwivedi et al., ]. Both tasks contain graphs corresponding to peptide molecules for a total of 15,535 graphs, 2.3 million nodes and 4.7 million edges. Each graph has an average of 150 nodes, 307 edges, and an average diameter of 57. Peptide-func is a multi-label graph classification task with the aim of predicting the peptide function. Peptide-struct is a multi-dimensional graph regression task, whose objective is to predict structural properties of the peptides.

We follow the experimental setting of Dwivedi et al. [], thus we consider a stratified splitting of the data, with 70% of graphs as the training set, 15% as the validation set, and 15% as the test set. We perform hyperparameter tuning via grid search, optimizing the Average Precision (AP) in the Peptide-func task and the Mean Absolute Error (MAE) in the Peptide-struct task, training the models using AdamW optimizer for a maximum of 300 epochs. For each model configuration, we perform 3 training runs with different weight initialization and report the average of the results. Moreover, we stay within 500K parameter budget as proposed by Dwivedi et al. []. In our experiments, we also consider the setting of Tönshoff et al. [] that propose to employ a 3-layer MLP as readout. We report in Appendix C.5 the grid of hyperparameters employed by each model in our experiment.

### C.5 Explored Hyperparameter Space

In Table 3 we report the grid of hyperparameters employed in our experiments by each method.

**Table 3: The grid of hyperparameters employed during model selection for the graph transfer tasks (*Transfer*), graph property prediction tasks (*GraphProp*), and Long Range Graph Benchmark (*LRGB*). We refer to Appendix A.2 for more details about dampening and external force implementations.**
| Hyperparameters | Values |  |  |
| --- | --- | --- | --- |
| *Transfer* | *GraphProp* | *LRGB* |  |
| Optimizer | Adam | Adam | AdamW |
| Learning rate | 0.001 | 0.003 | 0.001,  0.0005 |
| Weight decay | 0 | $10^{-6}$ | 0 |
| embedding dim | 64 | 10, 20, 30 | 195,  300 |
| N. layers ($L$) | 3,  5,  10,  50, | 1,  5,  10,  20,  30 | 32,  64 |
|  | 100,  150 |  |  |
| Termination time ($T$) | $L\epsilon$ | 0.1,  1,  2,  3,  4 | 3,  5,  6 |
| $\epsilon$ | 0.5,  0.2,  0.1, | $T/L$ | $T/L$ |
|  | 0.05,  0.01,  $10^{-4}$ |  |  |
| $\Phi_{\bf p}$ | Eq. 6,  GCN | Eq. 6,  GCN | Eq. 6,  GCN |
| $\Phi_{\bf q}$ | Eq. 6,  GCN | Eq. 6,  GCN | GCN |
| Readout input | ${\bf p}$, ${\bf q}$, ${\bf p}\|{\bf q}$ | ${\bf p}$, ${\bf q}$, ${\bf p}\|{\bf q}$ | ${\bf p}$, ${\bf q}$, ${\bf p}\|{\bf q}$ |
| $\sigma$ | tanh | tanh | tanh |
| Dampening | – | *param*, *param+*, | *param* |
|  |  | *MLP4-ReLU*, *DGN-ReLU* |  |
| External Force | – | *MLP4-Sin*, *DGN-tanh* | *DGN-tanh* |
| N. readout layers | 1 | 1 | 1,  3 |

## Appendix D Additional Experimental Results

### D.1 Graph Property Prediction Runtimes

In Table 4 we report runtimes of both our H-DGN and PH-DGN as well as baseline methods on the graph property prediction task in Section 3.3. To obtain a fair comparison, we configured all models with $20$ layers and an embedding dimension of $30$.
We averaged the time (in seconds) over 4 random weight initializations for a complete epoch, which includes the forward and backward passes on the training data, as well as the forward pass on the validation and test data.

Table 4 shows that both H-DGN and PH-DGN have improved or comparable runtimes compared to MPNNs. Notably, H-DGN is on average 5.92 seconds faster than GAT and 5.19 seconds faster than GCN. Compared to DE-DGN baselines, our methods show longer execution times, which are inherently caused by the sequential computation of both ${\bf p}$ and ${\bf q}$ explicited in Section A.3 and non-conservative components (i.e., *param* and *MLP4-ReLU* in Table 4).

**Table 4: Average time per epoch (measured in seconds) and std, averaged over 4 random weight initializations. Each time is obtained by employing $20$ layers and an embedding dimension equal to $30$. Our PH-DGN employs *param* and *MLP4-ReLU* as dampening and external force, respectively. The evaluation was carried out on an AMD EPYC 7543 CPU @ 2.80GHz. First, second, and third best results.**
| Model | Diameter | SSSP | Eccentricity |
| --- | --- | --- | --- |
| MPNNs |  |  |  |
| GCN | 32.45±2.54 | 17.44±3.85 | 11.78±2.43 |
| GAT | 20.20±5.18 | 26.41±8.34 | 17.28±1.92 |
| GraphSAGE | 13.12±2.99 | 13.12±2.99 | 8.20±0.75 |
| GIN | 6.63±0.28 | 21.16±2.33 | 14.22±3.17 |
| GCNII | 13.13±6.85 | 14.96±7.17 | 15.70±3.92 |
| DE-DGNs |  |  |  |
| DGC | 8.97±9.07 | 12.54±1.62 | 7.21±11.10 |
| GRAND | 133.84±42.57 | 109.15±27.49 | 202.46±85.01 |
| GraphCON | 9.26±0.47 | 7.76±0.05 | 7.80±0.05 |
| A-DGN | 8.42±2.71 | 7.86±2.11 | 13.18±9.07 |
| Ours |  |  |  |
| H-DGN | 15.49±0.05 | 15.28±0.02 | 15.34±0.04 |
| PH-DGN | 17.18±0.04 | 17.12±0.07 | 17.13±0.06 |

### D.2 Complete LRGB Results

In Table 5 we report the complete results for the LRGB tasks including more multi-hop DGNs and ablating on the scores obtained with the original setting from Dwivedi et al. [] and the one proposed by Tönshoff et al. [], which leverage 3-layers MLP as a decoder to map the DGN output into the final prediction. Since we are unable to compare the validation scores of the baselines directly, we cannot determine the best method between [Dwivedi et al., ] and Tönshoff et al. []. Therefore, in Table 5, we present all the results.

**Table 5: Results for Peptides-func and Peptides-struct averaged over 3 training seeds. Baseline results are taken from Dwivedi et al. [] and Gutteridge et al. []. Re-evaluated methods employ the 3-layer MLP readout proposed in Tönshoff et al. []. Note that all MPNN-based methods include structural and positional encoding. The first, second, and third best scores are colored.**
| Model | Peptides-func | Peptides-struct |
| --- | --- | --- |
| AP $\uparrow$ | MAE $\downarrow$ |  |
| MPNNs |  |  |
| GCN | 0.5930±0.0023 | 0.3496±0.0013 |
| GINE | 0.5498±0.0079 | 0.3547±0.0045 |
| GCNII | 0.5543±0.0078 | 0.3471±0.0010 |
| GatedGCN | 0.5864±0.0077 | 0.3420±0.0013 |
| Multi-hop DGNs |  |  |
| DIGL+MPNN | 0.6469±0.0019 | 0.3173±0.0007 |
| DIGL+MPNN+LapPE | 0.6830±0.0026 | 0.2616±0.0018 |
| MixHop-GCN | 0.6592±0.0036 | 0.2921±0.0023 |
| MixHop-GCN+LapPE | 0.6843±0.0049 | 0.2614±0.0023 |
| DRew-GCN | 0.6996±0.0076 | 0.2781±0.0028 |
| DRew-GCN+LapPE | 0.7150±0.0044 | 0.2536±0.0015 |
| DRew-GIN | 0.6940±0.0074 | 0.2799±0.0016 |
| DRew-GIN+LapPE | 0.7126±0.0045 | 0.2606±0.0014 |
| DRew-GatedGCN | 0.6733±0.0094 | 0.2699±0.0018 |
| DRew-GatedGCN+LapPE | 0.6977±0.0026 | 0.2539±0.0007 |
| Transformers |  |  |
| Transformer+LapPE | 0.6326±0.0126 | 0.2529±0.0016 |
| SAN+LapPE | 0.6384±0.0121 | 0.2683±0.0043 |
| GraphGPS+LapPE | 0.6535±0.0041 | 0.2500±0.0005 |
| DE-DGNs |  |  |
| GRAND | 0.5789±0.0062 | 0.3418±0.0015 |
| GraphCON | 0.6022±0.0068 | 0.2778±0.0018 |
| A-DGN | 0.5975±0.0044 | 0.2874±0.0021 |
| Re-evaluated |  |  |
| GCN | 0.6860±0.0050 | 0.2460±0.0007 |
| GINE | 0.6621±0.0067 | 0.2473±0.0017 |
| GatedGCN | 0.6765±0.0047 | 0.2477±0.0009 |
| DRew-GCN+LapPE | 0.6945±0.0021 | 0.2517±0.0011 |
| GraphGPS+LapPE | 0.6534±0.0091 | 0.2509±0.0014 |
| Ours |  |  |
| H-DGN | 0.6961±0.0070 | 0.2581±0.0020 |
| PH-DGN | 0.7012±0.0045 | 0.2465±0.0020 |

## Appendix E Limitations and Broader Impact

This paper aims to contribute to the field of Machine Learning, specifically focusing on advancing differential-equation inspired Deep Graph Networks (DE-DGNs). The research presented herein has the potential to positively impact the ongoing exploration and applications of DGNs.

Although our work offers a principled framework designed through the lens of port-Hamiltonian systems, we believe it is important to note that our presented inductive bias relies on the assumption of long-range interaction to solve the downstream task, and hence, we focused on data settings that follow this regime. Moreover, the numerical integration of all DE-DGNs intrinsically introduces a computation-accuracy trade-off, which is commonly explored through hyperparameter selection. As our experiments show, our message-passing architectures offer an approach to capture global interactions, but if one aims to obtain state-of-the-art results regardless of computational constraints and theoretical guarantees, other approaches may become superior by relaxing the given
graph topology by means of Transformers or rewiring techniques. The latter, in particular, do not contrast with the methodology we present here, and might be used synergistically with our approach in future developments. Lastly, we note that the upper bound discussed in Appendix A.1 shows that the sensitivity of node states with respect to previous states may increase, potentially leading to gradient growth (although gradient explosion was not observed in our experiments), thus requiring additional care by the user, e.g., using normalized neighborhood aggregation functions.

## Peer Reviews

confidence: 3
contribution: 2
initial_rating: 6
presentation: 3
questions: 1. While the paper shows good empirical performance, it's unclear what types of problems would theoretically benefit most from a port-Hamiltonian approach versus standard message passing. Could the authors provide analysis or insights about which properties of the underlying data generation process would suggest using their method?

2. The authors demonstrate that adding dissipative components often improves performance, but how does the balance between conservative (Hamiltonian) and dissipative parts relate to the nature of the task? It would be valuable to see an analysis of when pure conservation might be preferable to including dissipation, and how practitioners should choose this balance for new problems.

3. Given that the method is inspired by physical Hamiltonian systems, it's surprising that there are no experiments on problems with actual Hamiltonian dynamics (e.g., molecular dynamics simulations). Such experiments could help validate whether the method's conservation properties provide advantages for physically meaningful conservation laws, beyond just improving general information flow.
soundness: 4
strengths: Overall, the idea of using Hamilton's dynamics for GNN is attractive, though not entirely new. 
Nevertheless, the paper is solid and is more general than existing approaches. 
The improved efficiency over previous Hamiltonian GNN approaches and the practical benefits for long-range propagation make it a useful contribution to the field.
The experimental results are also interesting. 

1. Technical soundness:
- Clear theoretical analysis of conservation properties
- Explicit connection between Hamiltonian dynamics and message passing
- Thorough experimental validation

2. Practical value:
- More efficient than previous Hamiltonian GNN approaches
- Good performance on long-range tasks
- Can incorporate different message passing schemes

3. Clarity:
- Well-structured presentation
- Good balance of theory and empirics
- Clear comparisons to prior work
summary: This paper introduces port-Hamiltonian Deep Graph Networks (PH-DGN), a new framework for graph neural networks that addresses the challenge of long-range information propagation. The approach embeds message passing within a port-Hamiltonian dynamical system framework, where node states are split into position (q) and momentum (p) components coupled through skew-symmetric matrices. The Hamiltonian function incorporates graph structure through neighborhood aggregation functions, and the port-Hamiltonian extension introduces learnable dissipative components (internal dampening and external forces) that allow the network to balance conservative information preservation with task-dependent information modification. The authors provide theoretical analysis of the framework's properties, including bounds on sensitivity and gradient behavior, and demonstrate empirically that their method outperforms existing approaches on several tasks requiring long-range information propagation, including synthetic graph property prediction tasks and real-world molecular property prediction. The framework can incorporate different message passing schemes and provides improved efficiency compared to previous Hamiltonian-based graph neural networks, with experimental results showing that while the purely conservative version performs well, including the dissipative components often leads to better task performance.
weaknesses: 1. Novelty is incremental:
- Builds on existing ideas (Hamiltonian GNNs, message passing)
- Main contribution is combining these effectively rather than fundamentally new concepts

2. Technical questions:
- The derivation of the discretization scheme could use more justification
- Some assumptions about the structure of $W$ and $V$ matrices for explicit updates feel restrictive
- Could better explain why port-Hamiltonian framework is more appropriate than simpler alternatives

3. Empirical:
- Some ablation studies could be stronger (e.g., analyzing impact of different dissipative terms)
- Could better justify hyperparameter choices

---

confidence: 2
contribution: 2
initial_rating: 5
presentation: 3
questions: **Questions**

- l.161: What is the definition of anti-derivative of an activation function $\sigma$?
- l.223, l.228: Which norm is used in $\| \partial \mathbf{b}_u(T) / \partial \mathbf{b}_u(T-t) \|$? 
- l.259: *Therefore, [...]. , it holds the capability of PH-GDN to perform long-range propagation effectively*: the authors compare the upper bounds of sensitivity for the existing model (Theorem 2.3) and the proposed model (Theorem 2.4), then argue that since the latter is larger than the former, the proposed model is more effective in performing long-range propagation. However, it is difficult to claim so because there is the possibility that Theorem 2.4 only shows a looser bound than Theorem 2.3, which does not exclude the possibility that the upper bound does not properly reflect the sensitivity of the PH-GDN. It should be shown theoretically or experimentally that this upper bound adequately reflects the sensitivity of the PH-GDN.
- l.1129: I want to clarify the setup of the Graph Transfer Task: If I understand correctly, the feature vectors of all nodes are 1-dimensional and randomly sampled from $\mathrm{Unif}([0, 0.5))$. The target value is 0 for the source node and 1 for the target node. Assuming that this problem setup is correct, I need help understanding why this task is solvable because the model cannot distinguish source or target nodes from other nodes from feature vectors.


**Minor Comments**

- l.196: *non-dissipative (i.e., long-range)*: I think this is a slightly misleading expression. My understanding is that this paper uses non-dissipative in the sense of energy-perserving. Although non-dissipative implies long-range propagation (Theorem 2.3), they are not equivalent. In fact, this paper argues that non-dissipative PH-DGN performs better than conservative PH-DGN in predicting the LRGB dataset in some numerical experiments.
- l.436: The use of position encoding is only mentioned in the caption of Table 2 and should be explicitly stated in the text, specifically in the setup of Section 3.4.
- l.436: The correspondence between Table 2 and Table 5 needs to be clarified. For example, the method titled *MPNNs* in Table 2 is titled *re-evaluated* in Table 5. However, GCNII is not listed as *re-evaluated*, but as *MPNNs*. Furthermore, it is difficult to tell from the captions of Table 5 whether the authors test each listed method by themselves or is a citation of existing research. The reference should be indicated for each method in Table 5, for example, by adding a column such as *reference* to Table 5.
- l.443: *Overall, our port-Hamiltonian framework [...] shows great benefit [...] without requiring additional strategies such as global position encoding, global attention mechanism, or rewiring techniques [...].*: I suggest clarifying how the usefulness of each strategy is shown by comparing it with existing methods. More specifically:
  - The superiority of the proposed method over global position encoding is justified by the comparison with MPNN-based models using position encoding.    　　
  - The superiority over the global attention mechanism is by comparison with the Transformer-based method.
  - The superiority of rewiring methods by comparison with Drew.
- l.1153: The reference to Adam (Kingma & Ba, 2015) should be cited.
- l.1343: n. layers -> the number of layers
- l.1353: Table 6 claims that PH-DGN achieves both Hamiltonian conservation and Learnable driving forces. However, I think this is misleading because to achieve Hamiltonian conservation, the dissipative component must be removed, resulting in $\mathrm{PH-DGN}_{\mathrm{C}}$. It is not possible to achieve both characteristics at the same time in one architecture. Instead, two variants that achieve only one of the two are proposed.
soundness: 2
strengths: - The background knowledge of the port-Hamiltonian system is explained carefully, increasing accessibility for readers unfamiliar with this topic.
- For conservative PH-DGN, the ability of long-range interaction is theoretically shown by giving the lower bounds of sensitivity.
- The validity of the proposed methods for synthesis datasets is properly demonstrated.
summary: This paper proposes PH-DGN, a new GNN based on the port-Hamiltonian system to develop GNN that can solve graph learning tasks that require long-range dependencies. Two variants of PH-DGN are proposed: a conservative PH-DGN based on the Hamiltonian system and a PH-DGN based on the port-Hamiltonian system by introducing learnable dissipative terms. 
The theoretical analyses show that the conservative PH-DGN is stable and energy-preserving as a dynamical system, and derive a lower bound for the sensitivity, implying the possibility of long-range interaction.
Numerical experiments show that the conservative PH-DGN is energy-preserving without gradient vanishing using synthesis datasets (Section 3.1) and that the long-range interaction can be achieved in graph learning tasks that require long-range interactions (Sections 3.2, 3.3). Also, the usefulness of two variants of PH-DGN is evaluated on long-range graph benchmarks of real datasets.
weaknesses: - As the title indicates, the theme of this paper is the relationship between GNN and long-range propagation based on the port-Hamiltonian system. However, no theoretical guarantees on long-range propagation are given for general PH-DGNs with dissipative components.
- The tables of experimental results could be clearer (in particular Tables 2 and 5).
- For $\mathrm{PH-GDN}_{\mathrm{C}}$, which has theoretical guarantees, the prediction performance on the real dataset (long-range graph benchmark) does not outperform existing methods, and hence its practical usefulness is limited.

---

confidence: 3
contribution: 3
initial_rating: 8
presentation: 3
questions: - How would classical GCN aggregation interact with Theorem 2.4? Can the bound be easily extended for that case?
- In Section 3.1, you mention that the growing behavior can be controlled by regularizing the weight matrices or using normalized aggregation functions. Did you try this? How are the empirical results?
- Have you examined the interpretability of a trained PH-DGN? In particular, do the learned non-conservative forces make sense for the associated problem?
soundness: 4
strengths: - There is very clear motivation to this work, and it builds nicely upon other references.
- The proposed port-Hamiltonian approach is an original and clever way to allow for non-conservative dynamics in a graph network while still maintaining long-range message passing.
- The theoretical results for the conservative case are strong, and the motivation and interpretation for these are presented nicely.
- The experimental setup is superb; the care taken to ensure ease of replication is applauded. Model details are presented very clearly and choices are explained well for each step of the setup.
- A strong suite of models are compared against, with many different competitors and different approaches. The consistently favorable results provide a great deal of strength to claims of the proposed method's performance.
- The appendices are comprehensive for both proofs and experimental setup, and made clear many of the questions I had on an initial read.
summary: This work provides a novel methodology for the incorporation of port-Hamiltonian dynamics in graph representation learning. The central model, called a port-Hamiltonian Deep Graph Network (PH-DGN), is introduced first in a purely conservative setting using only the Hamiltonian and no non-conservative terms. Several theorems are developed to show that this conservative case leads to long-range data propagation between graph nodes, where the graph dynamics exhibit no energy loss and gradients do not vanish as the backward sensitivity matrix is bounded below. Dissipitive forces are then added to the full port-Hamiltonian model, which are two additional networks that may be trained to capture non-conservative forces. Several experiments follow, including a showcase of energy conservation and sensitivity to empirically verify theoretical work, and a graph transfer problem and graph property prediction problem to compare performance on benchmark tasks against other graph models in settings which require long-range information propagation.
weaknesses: - The majority of the theoretical results are developed for the conservative case. This makes sense in the context, as conservative long-range message passing is stated as a goal, but I would also be quite interested to see what could be proven for the fully general port-Hamiltonian case.
- In the explanation of Theorem 2.3, the statement that "the final representation of each node retains its complete past" seems somewhat strong. While I understand that the BSM result shows the influence of the entire past history on the current state, this statement as written seems to imply something stronger, and perhaps could be made more clear.
- The dissipitive force terms are added in some experiments to great success, but the explanations of their performance are more intuitive and are not supported by hard data in the paper. There may be a great opportunity here for visualization to support the intuitive claims.

There are two very minor typos:
- In the first paragraph of Section 2, "node states' in graph" has an unnecessary apostrophe.
- In Appendix D.3, "on a grid of n. layers" has an unnecessary period.

Meta-review: **(a) Scientific Claims and Findings:**
The paper introduces a novel framework called port-Hamiltonian Deep Graph Networks (pH-DGNs). This framework models neural information flow in graphs by leveraging principles from Hamiltonian dynamical systems, aiming to address challenges in long-range information propagation within graph representation learning. By incorporating both non-dissipative and non-conservative behaviors, the approach seeks to balance information conservation and dissipation, providing theoretical guarantees on information preservation over time. Empirical results demonstrate that pH-DGNs enhance the performance of simple graph convolutional architectures, achieving state-of-the-art results on benchmarks requiring long-range propagation.

**(b) Strengths:**
* Innovative Framework: The introduction of port-Hamiltonian systems into graph neural networks offers a fresh perspective on managing information flow, potentially addressing limitations in existing architectures concerning long-range dependencies.
* Theoretical Foundations: The framework is grounded in well-established principles from Hamiltonian dynamics, providing a solid theoretical basis for the proposed approach.
* Empirical Performance: The proposed method demonstrates superior performance on benchmarks that require long-range information propagation, indicating its practical effectiveness.
* Applicability: The approach can be integrated into general message-passing architectures, suggesting broad applicability across various graph-based learning tasks. 

**(c) Weaknesses:**
* Complexity: The incorporation of port-Hamiltonian systems may introduce additional complexity into the model, potentially impacting computational efficiency and implementation.
* Scope of Evaluation: While the empirical results are promising, the evaluation is primarily focused on benchmarks requiring long-range propagation. Assessing the framework's performance across a wider range of tasks and datasets would provide a more comprehensive understanding of its capabilities. The empirical improvements are limited in some cases.
* Practical Implementation Details: The paper could benefit from a more detailed discussion on the practical aspects of implementing the proposed framework, including computational requirements and potential challenges in real-world applications.

**(d) Reasons for Acceptance:**
After a thorough evaluation of the paper, I recommend acceptance based on the following considerations:
1. Novel Contribution: The paper presents a unique integration of port-Hamiltonian systems into graph neural networks, offering a new approach to addressing challenges in long-range information propagation.
2. Theoretical Rigor: The proposed framework is underpinned by solid theoretical foundations from Hamiltonian dynamics, enhancing the credibility and potential impact of the work.
3. Empirical Validation: The method demonstrates state-of-the-art performance on relevant benchmarks, providing empirical evidence of its effectivenes. Comparisons for the Minesweeper task from Platonov et al. (2023) were added during the rebuttal.
4. Broader Impact: The approach's applicability to general message-passing architectures suggests it could influence a wide range of graph-based learning tasks, contributing to advancements in the field.
While there are areas for improvement, such as expanding the scope of evaluation and providing more practical implementation details, the paper's strengths and contributions to the field warrant its acceptance.

Decision: Accept (Poster)