Node Embedding from Neural Hamiltonian Orbits in Graph Neural Networks
In the graph node embedding problem, embedding spaces can vary significantly
for different data types, leading to the need for different GNN model types. In
this paper, we model the embedding update of a node feature as a Hamiltonian
orbit over time. Since the Hamiltonian orbits generalize the exponential maps,
this approach allows us to learn the underlying manifold of the graph in
training, in contrast to most of the existing literature that assumes a fixed
graph embedding manifold with a closed exponential map solution. Our proposed
node embedding strategy can automatically learn, without extensive tuning, the
underlying geometry of any given graph dataset even if it has diverse
geometries. We test Hamiltonian functions of different forms and verify the
performance of our approach on two graph node embedding downstream tasks: node
classification and link prediction. Numerical experiments demonstrate that our
approach adapts better to different types of graph datasets than popular
state-of-the-art graph node embedding GNNs. The code is available at
\url{https://github.com/zknus/Hamiltonian-GNN}.
Qiyu Kang, Kai Zhao, Yang Song, Sijie Wang, Wee Peng Tay
# Node Embedding from Neural Hamiltonian Orbits in Graph Neural Networks

###### Abstract

In the graph node embedding problem, embedding spaces can vary significantly for different data types, leading to the need for different GNN model types. In this paper, we model the embedding update of a node feature as a Hamiltonian orbit over time. Since the Hamiltonian orbits generalize the exponential maps, this approach allows us to learn the underlying manifold of the graph in training, in contrast to most of the existing literature that assumes a fixed graph embedding manifold with a closed exponential map solution. Our proposed node embedding strategy can automatically learn, without extensive tuning, the underlying geometry of any given graph dataset even if it has diverse geometries. We test Hamiltonian functions of different forms and verify the performance of our approach on two graph node embedding downstream tasks: node classification and link prediction. Numerical experiments demonstrate that our approach adapts better to different types of graph datasets than popular state-of-the-art graph node embedding GNNs. The code is available at [https://github.com/zknus/Hamiltonian-GNN](https://github.com/zknus/Hamiltonian-GNN).

Machine Learning, Graph Neural Networks, Graph Neural Networks

## 1 Introduction

Graph neural networks (GNNs) (Yue et al., 2019; Ashoor et al., 2020; Kipf and Welling, 2017; Zhang et al., 2022; Wu et al., 2021) have shown remarkable inference performance on graph-structured data, including, but not limited to, social media networks, citation networks, and molecular graphs in chemistry. Most existing GNNs embed graph nodes in Euclidean spaces without further consideration of the dataset graph geometry. For some graph structures like the tree-like graphs (Liu et al., 2019), the Euclidean space may not be a proper choice for the node embedding. Recently, hyperbolic GNNs (Chami et al., 2019; Liu et al., 2019) propose to embed nodes into a hyperbolic space instead of the conventional Euclidean space. It has been shown that tree-like graphs can be inferred more accurately by hyperbolic GNNs.

Real-world graphs often have complex and varied structures, which can be more effectively represented by utilizing different geometric spaces. As shown in Figure 1, the Gromov \(\delta\)-hyperbolicity distributions1(Gromov, 1987) of various datasets exhibit a range of diverse values. This indicates that it is not optimal to embed each dataset with diverse geometries into a single globally homogeneous geometry. Works like (Zhu et al., 2020) have attempted to embed graph nodes in a mixture of the Euclidean and hyperbolic spaces, where the intrinsic graph local geometry is attained from the mixing weight. In several studies, including (Gu et al., 2018; Bachmann et al., 2020; Lou et al., 2020), researchers use (products of) constant curvature Riemannian spaces for graph node embedding where the spaces are assumed to be spherical, hyperbolic, or Euclidean. The work (Xiong et al., 2022) considers a special pseudo-Riemannian manifold named pseudo-hyperboloid that is of constant nonzero curvature and is diffeomorphic to the product manifolds of a unit sphere and the Euclidean space.

Footnote 1: If the \(\delta\)-hyperbolicity distribution is more concentrated at lower values, the more hyperbolic the graph dataset.

Embedding nodes in the aforementioned _restricted_ (pseudo-)Riemannian manifolds is achieved through the exponential map in closed forms, which is essentially a geodesic curve on the manifolds as the projected curve of the _cogeodesic orbits_ on the manifolds' cotangent bundles (Lee, 2013; Klingenberg, 2011). In our work, we propose to embed the nodes, via more general _Hamiltonian orbits_, into a general manifold, which generalizes the above graph node embedding works.

Manifolds have a diverse set of applications in physics, and their usage and development can be found interwoven throughout the literature. From the physics perspective, the cotangent bundles are the natural phase spaces in classical mechanics (De Leon and Rodrigues, 2011) where the physical system evolves according to the basic laws of physics modeled as differential equations on the phase spaces. In
this paper, we propose a new GNN paradigm based on Hamiltonian mechanics (Goldstein et al., 2001) with flexible Hamiltonian functions. Our objective is to design a new node embedding strategy that can automatically learn, without extensive tuning, the underlying geometry of any given graph dataset even if it has diverse geometries. We enable the node features to evolve on the manifold under the influence of neighbors. The learnable Hamiltonian function on the manifold guides the node embedding evolution to follow a learnable law analogous to basic physical laws.

**Main contributions.** Our main contributions are summarized as follows:

1. We consider the graph node embedding problem on an underlying manifold and enable node embedding through a learnable Hamiltonian orbit associated with the Hamiltonian scalar function on the manifold's cotangent bundle.
2. Our node embedding strategy can automatically learn, without extensive tuning, the underlying geometry of any given graph dataset even if it has diverse geometries. We empirically demonstrate its ability by testing on two graph node embedding downstream tasks: node classification and link prediction.
3. From empirical experiments, we observe that the over-smoothing problem of GNNs can be mitigated if the node features evolve through Hamiltonian orbits. By the conservative nature of the Hamiltonian equations, our model enables a stable training and inference process while updating the node features over time and layers.

## 2 Related Work

While our paper is related to Hamiltonian neural networks in the literature, we are the first, to our best knowledge, to model graph node embedding with Hamiltonian equations. We briefly review Hamiltonian neural networks, Riemannian manifold GNNs, and physics-inspired GNNs in Appendix A.

_Notations:_ We use the _Einstein summation convention_(Lee, 2013) for expressions with tensor indices. When using this convention, if an index variable appears twice in a term, once as a superscript and once as a subscript, it means a summation of the term over all possible values of the index variable. For example, \(a^{i}b_{i}\triangleq\sum_{i=1}^{d}a^{i}b_{i}\).

## 3 Motivations and Preliminaries

In this section, we briefly review the concepts of the geodesic curve on a (pseudo-)Riemannian manifold from the principle of stationary action in the form of Lagrange's equations. We then further generalize the geodesic curve to the Hamiltonian orbit associated with an energy function \(H\), which is a conserved quantity along the orbit. Our primary goal is to develop a more flexible and robust method for graph node embedding by leveraging the concepts of geodesic curves and Hamiltonian orbits on manifolds. We first summarize the motivation of our work as follows.

**Motivation I: from the exponential map to the Riemannian geodesic.** The geodesic curve gives rise to the exponential map that maps points from the tangent space to the manifold and has been utilized in (Chami et al., 2019; Bachmann et al., 2020; Xiong et al., 2022) to enable graph node embedding in some restricted (pseudo-)Riemannian manifolds where closed forms of the exponential map are obtainable. From this perspective, by using the geodesic curve, we generalize the graph node embedding to an arbitrary pseudo-Riemannian manifold with learnable local geometry \(g\) using Lagrange's equations.

**Motivation II: from geodesic to Hamiltonian orbit.** Despite the above conceptual generalization for node embedding using geodesic curves, the specific curve formulation involving minimization of curve length may result in a loss of generality for node feature evolution along the curve. We thus further generalize the geodesic curve to the Hamiltonian orbit associated with a more general energy function \(H\) that is conserved along the orbit. In Section 4, we propose graph node embedding without an explicit metric by using Hamiltonian orbits with learnable energy functions \(H\).

We kindly refer readers to Appendix B for a more comprehensive elucidation that may enhance their understanding after we introduce the related concepts in Section 3.2.

### Manifold and Riemannian Metric

**Manifold and local chart representation.** On a \(d\)-dimensional manifold \(M\), for each point on \(M\), there exists a triple \(\{q,U,V\}\), called a chart, such that \(U\) is an open

Figure 1: Gromov \(\delta\)-hyperbolicity distribution of datasets.
neighborhood of the point in \(M\), \(V\) is an open subset of \(\mathbb{R}^{d}\), and \(q:U\to V\) is a homeomorphism, which gives us a coordinate representation for a local area in \(M\).

**Tangent and cotangent vector spaces.** For any point \(q\) on \(M\) (we identify each point covered by a local chart on \(M\) by its representation \(q\)), we may assign two vector spaces named the tangent vector space \(T_{q}M\) and cotangent vector space \(T_{q}^{*}M\). The vectors from the tangent and cotangent spaces can be interpreted as representing the velocity and generalized momentum, respectively, of an object's movement in classical mechanics.

**Riemannian metric.** A Riemannian manifold is a manifold \(M\) equipped with a _Riemannian metric_\(g\), where we assign to any point \(q\in M\) and pair of vectors \(u,v\in T_{q}M\) an _inner product_\(\langle u,v\rangle_{g(q)}\). This assignment is assumed to be smooth with respect to the base point \(q\in M\). The length of a tangent vector \(u\in T_{q}M\) is then defined as

\[\|u\|_{g(q)}:=\langle u,u\rangle_{g(q)}^{1/2}. \tag{1}\]

* **Local coordinates representation:** In local coordinates with \(q=\left(q^{1},\ldots,q^{d}\right)^{\intercal}\in M,u=\left(u^{1},\ldots,u^{d} \right)^{\intercal}\in T_{q}M\) and \(v=\left(v^{1},\ldots,v^{d}\right)^{\intercal}\in T_{q}M\), the Riemannian metric \(g=g(q)\) is a real symmetric _positive definite_ matrix and the inner product above is given by \[\langle u,v\rangle_{g(q)}:=g_{ij}(q)u^{i}v^{j}\] (2)
* **Pseudo-Riemannian metric:** We may generalize the Riemannian metric to a metric tensor that only requires a _non-degenerate_ condition (Lee, 2018) instead of the stringent positive definiteness condition in the inner product. One example of a pseudo-Riemannian manifold is the Lorentzian manifold, which is important in applications of general relativity.

### Geodesic Curves and Exponential Maps

**Length and energy of a curve.** Let \(q:[a,b]\to M\) be a smooth curve.2 We use \(\dot{q}\) and \(\ddot{q}\) to denote the first and second order derivatives of \(q(t)\), respectively. We define the following:

Footnote 2: We abuse notations in denoting the chart coordinate map as \(q\) and the curve as \(q(t)\). It will be clear from the context which one is being referred to.

\(\bullet\) length of the curve:

\[\ell(q):=\int_{a}^{b}\|\dot{q}(t)\|_{g(q(t))}\,\mathrm{d}t. \tag{3}\]

\(\bullet\) energy of the curve:

\[E(q):=\frac{1}{2}\int_{a}^{b}\|\dot{q}(t)\|_{g(q(t))}^{2}\,\mathrm{d}t. \tag{4}\]

**Geodesic curves.** On a Riemannian manifold, geodesic curves are defined as curves that have a minimal length as given by (3) and with two fixed endpoints \(q(a)\) and \(q(b)\). However, computations based on minimizing the length to obtain the curves are difficult. It turns out that the minimizers of \(E(q)\) also minimize \(\ell(q)\)(Malham, 2016). Consequently, the geodesic curve formulation may be obtained by minimizing the energy of a smooth curve on \(M\).

**Principle of stationary action and Euler-Lagrange equation.** The Lagrangian function \(L(q(t),\dot{q}(t))\) minimizes the following functional (in physics, the functional is known as an **action**)

\[S(q)=\int_{a}^{b}L(q(t),\dot{q}(t))\mathrm{d}t \tag{5}\]

with two fixed endpoints at \(t=a\) and \(t=b\) only if the following _Euler-Lagrange equation_ is satisfied:

\[\frac{\partial L}{\partial q^{i}}(q(t),\dot{q}(t))-\frac{\mathrm{d}}{\mathrm{d }t}\frac{\partial L}{\partial\dot{q}^{i}}(q(t),\dot{q}(t))=0. \tag{6}\]

**Geodesic equation for geodesic curves.** The Euler-Lagrange equation derived from minimizing the energy (4) with local coordinates representation,

\[L=\frac{1}{2}\|\dot{q}(t)\|_{g(q(t))}^{2}=\frac{1}{2}g_{ik}(q)\dot{q}^{i}\dot{ q}^{k} \tag{7}\]

is expressed as the following ordinary differential equations called the _geodesic equation_:

\[\ddot{q}^{i}+\Gamma_{jk}^{i}\dot{q}^{j}\dot{q}^{k}=0, \tag{8}\]

for all \(i=1,\ldots,d\), where the Christoffel symbols \(\Gamma_{jk}^{i}:=\frac{1}{2}g^{ik}\left(\frac{\partial g_{ij}}{\partial q_{k}}+ \frac{\partial g_{kl}}{\partial q_{j}}-\frac{\partial g_{jk}}{\partial q_{t}}\right)\) and \([g^{ij}]\) denotes the inverse matrix of the matrix \([g_{ij}]\). The solutions to the geodesic equation (8) give us the geodesic curves.

**Exponential map.** Given the geodesic curves, at each point \(q\in M\), for velocity vector \(v\in T_{q}M\), the _exponential map_ is defined to obtain the point on \(M\) reached by the unique geodesic that passes through \(q\) with velocity \(v\) at time \(t=1\)(Lee, 2018). Formally, we have

\[\exp_{q}(v)=\gamma(1) \tag{9}\]

where \(\gamma(t)\) is the curve given by the geodesic equation (8) with initial conditions \(q(0)=q\) and \(\dot{q}(0)=v\).

With regards to **Motivation I**, we note that (Chami et al., 2019) considers graph node embedding over a homogeneous negative-curvature Riemannian manifold called hyperboloid manifold. In contrast, we generalize the embedding of nodes to an arbitrary pseudo-Riemannian manifold through the geodesic equation (8) with a learnable metric \(g\) that derives the local graph geometry from the nodes and their neighbors.
### From Geodesics to General Hamiltonian Orbits

The geodesic curves and the derived exponential map essentially come from (5) with \(L\) in (7) specified from the curve energy (4). However, the curves derived from this specific action may potentially sacrifice efficacy for the graph node embedding task since we do not know what a reasonable action formulation that guides the evolution of the node feature in this task is. Therefore, we follow the principle of stationary action but consider a learnable action that is more flexible than the length or energy of the curve. To better model the conserved quantity during the feature evolution, we reformulate the Lagrange equation to the Hamilton equation. This is our **Motivation II**.

**Hamiltonian function and equations.** The Hamiltonian orbit \((q(t),p(t))\) is given by the following _Hamiltonian equations_ with a _Hamiltonian function_\(H\):

\[\dot{q}^{i}=\frac{\partial H}{\partial p_{i}},\quad\dot{p}_{i}=-\frac{ \partial H}{\partial q^{i}}, \tag{10}\]

where \(q\) is the local chart coordinate on the manifold while \(p\) can be interpreted as a vector of generalized momenta in the cotangent vector space. In classical mechanics, the \(2d\)-dimensional pair \((q,p)\) is called _phase space_ coordinates that fully specify the state of a dynamic system with \(p\) guiding the movement direction and speed. Later, we consider the node feature evolution following the trajectory specified by the phase space coordinates.

**Hamiltonian function vs. Lagrangian function.** The Hamiltonian function can be taken as the Legendre transform of the Lagrangian function:

\[\begin{split} H(q,p)&=p_{i}\dot{q}^{i}-L(q,\dot{q} )\\ \text{with }\dot{q}&=\dot{q}(p)\text{ s.t. }p=\frac{ \partial L}{\partial\dot{q}}.\end{split} \tag{11}\]

If \(H\) is restricted to strictly convex functions, the Hamiltonian formalism is equivalent to a Lagrangian formalism (De Leon & Rodrigues, 2011).

**Geodesic equation reformulated as Hamiltonian function.** If \(H\) is set as

\[H(q,p)=\frac{1}{2}g^{ij}(q)p_{i}p_{j}, \tag{12}\]

where \([g^{ij}]\) denotes the inverse matrix of the matrix \([g_{ij}]\), we have the following Hamiltonian equations:

\[\dot{q}^{i}=g^{ij}p_{j},\quad\dot{p}_{i}=-\frac{1}{2}\partial_{i}g^{jk}p_{j}p_ {k}. \tag{13}\]

The Hamiltonian orbit \((q(t),p(t))\), as the solution of (13), gives us again the _geodesic curves_\(q(t)\) on the manifold \(M\) if we only look at the first \(d\)-dimensional coordinates.

**Theorem 1** (Conservation of energy (Da Silva & Da Salva, 2008)).: \(H(q(t),p(t))\) _is constant along the Hamiltonian orbit as solutions of (10)._

In physics, \(H\) typically represents the total energy of a system, and Theorem 1 indicates that the time evolution of the system follows the law of conservation of energy.

## 4 Proposed Framework

We consider an undirected graph \(\mathcal{G}=(\mathcal{V},\mathcal{E})\) consisting of a finite set \(\mathcal{V}\) of vertices, together with a subset \(\mathcal{E}\subset\mathcal{V}\times\mathcal{V}\) of edges. Our objective is to learn a mapping \(f\) that maps node features to embedding vectors:

\[f\left(\mathcal{V},\mathcal{E}\right)=Z\in\mathbb{R}^{|\mathcal{V}|\times d}.\]

The embedding \(Z\) is supposed to capture both semantic and topological information. Since the input node features in most datasets are sparse, a fully connected (FC) layer is first applied to compress the raw input node features. Let \({}^{n}q\) be the \(d\)-dimensional compressed node feature for node \(n\) after the FC layer.3 However, empirical experiments (see "MLP" results in Section 5.1) indicate that for the graph node embedding task, such simple raw compressing without any consideration of the graph topology does not produce a good embedding. Further graph neural network layers are thus required to update the node embedding.

Footnote 3: We put the node index to the left of the variable to distinguish it from the manifold dimension index.

We consider the node features \(\{^{n}q\}_{n\in\mathcal{V}}\) to be located on an embedding manifold \(M\) and take the node features as the chart coordinate representations for points on the manifold. In **Motivations I and II** in Section 3, we have provided the rationale for generalizing graph node embedding from the exponential map to the pseudo-Riemannian geodesic, and further to the Hamiltonian orbit. To enforce the graph node feature update on the manifold with well-adapted learnable local geometry, we make use of the concepts from Section 3.

### Model Architecture

**Node feature evolution along Hamiltonian orbits in a Hamiltonian layer.** As introduced in Section 3.3, the \(2d\)-dimensional _phase space_ coordinates \((q,p)\) fully specify a system's state. Consequently, for the node feature as a point on the manifold \(M\), we associate to each point \({}^{n}q\) a _learnable momentum vector_\({}^{n}p\). This vector guides the direction and speed of the feature update, allowing the node feature to evolve along Hamiltonian orbits on the manifold. More specifically, we set

\[{}^{n}p=Q_{\mathrm{net}}(^{n}q) \tag{14}\]

where \(Q_{\mathrm{net}}\)4 is instantiated by an FC layer. We consider a learnable Hamiltonian function \(H_{\mathrm{net}}\) that specifies the node feature evolution trajectory in the phase space via the
#### Hamiltonian equations

\[\dot{q}^{i}=\frac{\partial H_{\mathrm{net}}}{\partial p_{i}},\quad\dot{p}_{i}=- \frac{\partial H_{\mathrm{net}}}{\partial q^{i}} \tag{15}\]

with learnable Hamiltonian energy function

\[H_{\mathrm{net}}:(q,p)\mapsto\mathbb{R}. \tag{16}\]

The node features are updated along the Hamiltonian orbits, which are curves starting from each node \(({}^{n}q,{}^{n}p)\) at \(t=0\). In other words, they are the solution of (15) with the initial conditions \(({}^{n}q(0),{}^{n}p(0))=({}^{n}q,{}^{n}p)\) at \(t=0\). The solution of (15) on the phase space for each node \(n\in\mathcal{V}\) at time \(T\) is given by a differential equation solver (Chen et al., 2018), and denoted by \(({}^{n}q(T),{}^{n}p(T))\). The canonical projection \(\pi({}^{n}q(T),{}^{n}p(T))={}^{n}q(T)\) is taken to obtain the node feature embeddings on the manifold at time \(T\). The aforementioned operations are performed within one layer, and we call it the _Hamiltonian layer_.

**Neighborhood Aggregation.** After the node features update along the Hamiltonian orbits, we perform neighborhood aggregation on the features \(\{{}^{n}q^{(\ell)}(T)\}_{n\in\mathcal{V}}\), where \(\ell\) indicates the \(\ell\)-th layer. Let \(\mathcal{N}(n)=\{m:(n,m)\in\mathcal{E}\}\) denote the set of neighbors of node \(n\in\mathcal{V}\). We only perform a simple yet efficient aggregation (see Section 5) for node \(n\) as follows:

\[{}^{n}q^{(\ell+1)}={}^{n}q^{(\ell)}(T)+\frac{1}{|\mathcal{N}(n)|}\sum_{m\in \mathcal{N}(n)}{}^{m}q^{(\ell)}(T). \tag{17}\]

**Layer stacking for local geometry learning.** We stack up multiple Hamiltonian layers with neighborhood aggregation in between them. We first give an intuitive explanation for the case where \(H_{\mathrm{net}}\) is set as (12). A learnable metric \(g_{\mathrm{net}}\) for the manifold is involved (see Section 4.2.1 for more details) and the features are evolved following the geodesic curves with minimal length (see Section 3). Within each Hamiltonian layer, the metric \(g_{\mathrm{net}}\) that is instantiated by a smooth FC layer only depends on the local node position on a pseudo-Riemannian manifold that varies from point to point. Note that with layer stacking, these features contain information aggregated from their neighbors. The metric \(g_{\mathrm{net}}\), therefore, learns from the graph topology, and each node is embedded with a local geometry that depends on its neighbors. In contrast, (Chami et al., 2019; Bachmann et al., 2020; Xiong et al., 2022) consider graph node embedding using geodesic curves over some fixed manifold without adjustment of the local geometry.

At the beginning of Section 4, we have assumed the node features \(\{{}^{n}q\}_{n}\) to be located in local charts of a preliminary embedding manifold \(M\). The basic philosophy is that the embedding manifold evolves with a metric structure that adapts successively with neighborhood aggregation along multiple layers, whereas each node's features evolve to the most appropriate embedding on the manifold along the curves. For a general learnable \(H_{\mathrm{net}}\), the Hamiltonian orbit that starts from one node has aggregated information from its neighbors, which guides the learning of the curve that the node will be evolved along. Therefore, each node is embedded into a manifold with adaptation to the underlying geometry of any given graph dataset even if it has diverse geometries.

Figure 2: HamGNN architecture: in each layer, each node is assigned a learnable “momentum” vector \({}^{n}p\) (cf. (14)) at time \(t=0\), which initializes the evolution of the node feature. The node features evolve on a manifold following (15) to \(({}^{n}q(T),{}^{n}p(T))\) at the time \(t=T\). We only take \({}^{n}q(T)\) as the embedding and input it to the next layer. After \(L\) layers, we take \({}^{n}q^{(L)}(T)\) as the final node embedding.
**Conservation of \(H_{\rm net}\).** From Theorem 1, the feature updating through the orbit indicates that the \(H_{\rm net}\) is conserved along the curve.

**Model summary.** Our model is called _HamGNN_ as we use Hamiltonian orbits for node feature updating on the manifold. We summarize the HamGNN model architecture in Figure 2 and Algorithm 1. The forms of the Hamiltonian function \(H_{\rm net}\) are given in Section 4.2.

### Different Hamiltonian Orbits

We next propose different forms for \(H_{\rm net}\) from which the corresponding Hamiltonian orbit and its variations are obtained in our GNN model. The node features are updated along the Hamiltonian orbits, which are curves starting from each node \((^{n}q,^{n}p)\) at \(t=0\). Beginning with Motivation I from the paper, we design a learnable metric \(g_{\rm net}\) in (12) in Section 4.2.1. This approach relaxes the curve formulation constraint used in (Chami et al., 2019; Bachmann et al., 2020; Xiong et al., 2022) and enables learnable geodesic curves to guide feature evolution on arbitrary (pseudo-)Riemannian manifolds while learning local geometry from graph datasets. To further extend learnable geodesic curves to learnable Hamiltonian orbits on manifolds, as stated in Motivation II, we introduce a flexible \(H\) instantiated by an FC layer without constraints in Section 4.2.2. Subsequently, we present variations based on Section 4.2.2 in Sections 4.2.3 through 4.2.5 to examine their performance. In Section 4.2.3, we add constraints to ensure that \(H\) is a convex function, which allows us to equivalently test a more restricted Lagrangian formalism. In Section 4.2.4, we include less restricted Hamiltonian mechanics without strict constant energy constraints. Finally, in Section 4.2.5, we consider a more flexible representation using the symplectic 2-form in comparison to Section 4.2.2.

#### 4.2.1 Learnable Metric \(g_{\rm net}\)

In this subsection, we consider node embedding onto a pseudo-Riemannian and set \(H_{\rm net}\) as (12) where a learnable metric \(g_{\rm net}\) for the manifold is involved. Within each Hamiltonian layer, the metric \(g_{\rm net}\) instantiated by a smooth FC layer depends on the local node position on the pseudo-Riemannian manifold that varies from point to point. The output of \(g_{\rm net}\) at position \(q\) represents the inverse metric local representation \([g^{ij}]\). However, from (12), the space complexity is order \(d^{3}\) due to the partial derivative of \(g\)'s output being a \(d\times d\) matrix. We therefore only consider _diagonal metrics_ to mitigate the space complexity. More specifically, we now define

\[g_{\rm net}(q)={\rm diag}([\underbrace{-1,\ldots,-1}_{r},\underbrace{1,\ldots,1}_{s}]\odot h_{\rm net}(q)) \tag{18}\]

where \(h_{\rm net}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}\) consists of non-linear trainable layers and \(\odot\) denotes element-wise multiplication. To ensure non-degeneracy of the metric, the output of \(h_{\rm net}\) is set to be away from \(0\) with the final activation function of it being strictly positive. The vector \([-1,\ldots,-1,1,\ldots,1]\) controls the signature \((r,s)\)(Lee, 2018) of the metric \(g\) with \(r+s=d\), where \(r\) and \(s\) are the number of \(-1\)s and \(1\)s, respectively. The signature of the metric is set to be a hyperparameter. According to (13), we have

\[\dot{q}^{i}=g_{\rm net}^{ij}p_{j},\quad\dot{p}_{i}=-\frac{1}{2}\partial_{i}g_ {\rm net}^{jk}p_{j}p_{k}. \tag{19}\]

Intuitively, the node features evolve through the "shortest" curves on the manifold. The exponential map used in (Chami et al., 2019; Bachmann et al., 2020; Xiong et al., 2022) is essentially the geodesic curve on a hyperbolic manifold with an explicit formulation due to the manifold type restriction. We do not enforce any assumption here and let the model learn the embedding geometry.

#### 4.2.2 Learnable \(H_{\rm net}\)

Different from Section 4.2.1 where \(H\) is set as (12) with a pseudo-Riemannian metric, we choose a more flexible \(H\) instantiated by an FC layer and consider the Hamiltonian equations:

\[\dot{q}^{i}=\frac{\partial H_{\rm FC}}{\partial p_{i}},\quad\dot{p}_{i}=- \frac{\partial H_{\rm FC}}{\partial q^{i}}. \tag{20}\]

#### 4.2.3 Learnable Convex \(H_{\rm net}\)

As discussed in Section 3.3, if \(H_{\rm net}\) is restricted to strictly convex functions, the Hamiltonian formalism can degenerate to a _Lagrangian formalism_ through the Legendre transformation (11). We take the following restricted Hamiltonian equations

\[\dot{q}^{i}=\frac{\partial H_{\rm net}}{\partial p_{i}},\dot{p}_{i}=-\frac{ \partial H_{\rm net}}{\partial q^{i}},{\rm s.t.}\ H_{\rm net}\ \text{is convex}, \tag{21}\]

where a stationary action in (5) is achieved. To guarantee that \(H_{\rm net}\) is convex, we follow the work in (Amos et al., 2017) to set non-negative layer weights from the second layer in \(H_{\rm net}\), and all activation functions in \(H_{\rm net}\) to be convex and non-decreasing. This network design is shown to be able to approximate any convex functions in (Chen et al., 2018).

#### 4.2.4 Learnable \(H_{\rm net}\) with Relaxation

Different from Section 4.2.2, we enforce additional system biases along the curve as follows:

\[\dot{q}^{i}=\frac{\partial H_{\rm net}}{\partial p_{i}},\quad\dot{p}_{i}=- \frac{\partial H_{\rm net}}{\partial q^{i}}+f_{\rm net}(q). \tag{22}\]

Instead of keeping the energy during the feature update along the Hamiltonian orbit, we now also include an additional energy term during the node feature update.
#### 4.2.5 Learnable \(H_{\rm net}\) with a flexible symplectic form

Hamiltonian equations have a more flexible representation using the symplectic \(2\)-form (cf. Appendix G). The chart coordinate representation \((q,p)\) may not be the Darboux coordinate system for the symplectic \(2\)-form. Even if our learnable \(H_{\rm net}\) may be able to learn the energy representation under the chosen chart coordinate system, we consider a learnable symplectic \(2\)-form to act in concert with \(H_{\rm net}\). More specifically, following (Chen et al., 2021), we have the following 1-form

\[\theta^{1}_{\rm net}=f_{i,{\rm net}}{\rm d}q^{i},\]

where \(f_{i,{\rm net}}:M\rightarrow\mathbb{R}^{d}\) is the output's \(i\)-th component of the neural network parameterized function. The Hamiltonian orbit is then provided by (Chen et al., 2021) as follows:

\[\big{(}\dot{q}^{i},\dot{p}^{i}\big{)}=W^{-1}(q,p)\nabla H_{\rm net}(q,p), \tag{23}\]

where the skew-symmetric \(2d\times 2d\) matrix \(W\), whose elements are written in terms of \((\partial_{i}f_{j,{\rm net}}-\partial_{j}f_{i,{\rm net}})\), is given in (24) in the appendix due to space limitations.

## 5 Experiments

In this section, we implement the proposed HamGNNs with different settings as shown in Section 4.2 and Appendix C. We select **datasets with various geometries** including the three citation networks: Cora (McCallum et al., 2004), Citeseer (Sen et al., 2008), Pubmed (Namata et al., 2012); and two low hyperbolicity datasets (Chami et al., 2019): Disease and Airport (cf. Table 5). Furthermore, we create **new datasets with more complex geometry** by combining Disease and Cora/Citeseer so that the new datasets have a mixture of both hyperbolic and Euclidean local geometries.

Adhering to the experimental settings of (Chami et al., 2019; Gu et al., 2018; Bachmann et al., 2020; Lou et al., 2020; Xiong et al., 2022), we evaluate the effectiveness of the node embedding by performing _two downstream tasks: node classification and link prediction_ using the embeddings. Rather than beating all the existing GNNs on these two specific tasks, we want to demonstrate that our node embedding strategy is able to automatically learn, without extensive tuning, the underlying geometry of any given graph dataset even if its underlining geometry is very complex, e.g., a mixture of multiple geometries. Such examples can be found in Table 3. Due to space constraints, we refer the readers to Appendix E for the description of datasets and implementation details.

To fairly compare the performance of the proposed HamGNN, for the node classification tasks, we select several popular GNN models as the baseline. These include _Euclidean GNNs_: GCN (Kipf and Welling, 2017), GAT (Velickovic et al., 2018), SAGE (Hamilton et al., 2017), and SGC (Wu et al., 2019); _Hyperbolic GNNs_(Chami et al., 2019; Liu et al., 2019): HGNN, HGCN, HGAT and LGCN (Zhang et al., 2021); a _GNN that mixes Euclidean and hyperbolic embeddings_: GIL (Zhu et al., 2020); _(Pseudo-)Riemannian GNNs_: \(\kappa\)-GCN (Bachmann et al., 2020) and \(\mathcal{Q}\)-GCN (Xiong et al., 2022); as well as _Graph Neural Diffusions_: GRAND (Chamberlain et al., 2021) and GraphCON (Rusch et al., 2022). We also include the MLP baseline, which does not utilize the graph topology information. To further demonstrate the advantage of HamGNN, we also
include one vanilla ODE system, whose formulation is given in Appendix F.1. This vanilla ODE system neither includes the learnable "momentum" vector \(p\) nor adheres to the Hamiltonian orbits (26). For the link prediction task, we compare HamGNN to the _standard graph node embedding models_, including all the aforementioned baselines in the node classification task except the graph neural diffusion baselines. For the link prediction tasks, we report the best results from different versions of HamGNN: HamGNN (19) on the Disease dataset and (20) on the remaining datasets.

### Performance Results and Ablation Studies

**Node classification.** The node classification performance on the benchmark datasets using the baseline models and the proposed HamGNNs with different \(H_{\mathrm{net}}\)s is shown in Table 1. We observe that HamGNN adapts well to all datasets with various geometries. For the datasets such as Cora, Citeseer and Pubmed, which can be well embedded into Euclidean space, HamGNNs achieve comparable performance to the conventional Euclidean GNNs. For example, HamGNN (19) achieves the third-best performance on the Cora. For tree-like graph data Disease and Airport, HamGNNs outperform the other GNNs, including all the hyperbolic GNNs, which are tailored for hyperbolic datasets, and Riemannian GNNs which use a Cartesian product of spherical, hyperbolic, and Euclidean spaces. From the results in Table 1, we can see that Riemannian GNNs \(\kappa\)-GCN and \(\mathcal{Q}\)-GCN are biased towards embedding data into Euclidean space since they have the best performance on high hyperbolicity datasets but average performance on low hyperbolicity datasets.

Furthermore, the comparison between HamGNNs with the vanilla ODE GNN in (27) without any Hamiltonian mechanism implies the importance of the Hamiltonian layer.

**Link prediction.** In Table 2, we report the averaged ROC for the link prediction task. We observe HamGNN adapts well to all datasets and is the best performer on the Airport, Citeseer, and Cora.

**Comparison between different \(H_{\mathrm{net}}\).** We compare HamGNNs with different \(H_{\mathrm{net}}\)s as elaborated in Section 4.2 on the node classification task. In Section 3.3, we argue that the geodesic curves derived from the action of curve length may potentially sacrifice efficacy for the graph node embedding task since we do not know what a reasonable action formulation that guides the evolution of the node feature in this task is. We therefore also include a more flexible \(H_{\mathrm{net}}\) in (20) with other variations, e.g., (21) to (23). From Table 1, we observe that the original HamGNN (19) shows good adaptations of node embedding to various datasets. Simply using a FC layer for \(H\) in (20) has no obvious improvement. Further imposing convexity on \(H\) in (21) also has little positive impact on the performance. Including system bias in (22) achieves the best performance on Citeseer among all HamGNNs. The HamGNNs that achieve the best performance on Disease and Airport are (25) and (26). Overall, all HamGNN variants elaborated in Section 4.2 have well-adapted node embedding performance for various datasets even though some variants may perform slightly better. This observation indicates that good geometry adaptations may come from the HamGNN model architecture with the Hamiltonian orbit evolution. This is further verified by the node classification performance from vanilla ODE in (27), which does not have a well-adapted node embedding performance and is designed without the philosophy of Hamiltonian mechanics. Good experiment results in Section 5.2 using both (19) and (20) also provide further evidence to support our conclusion.

We next compare the HamGNNs using (20) and (23). The

\begin{table}
\begin{tabular}{c|c|c|c|c|c} \hline Method & Disease & Airport & Pubmed & Citeseer & Cora \\ \(\delta\)-hypercholcity & 0.0 & 1.0 & 3.5 & 4.5 & 11.0 \\ \hline MLP & 83.74 & 84.04 & 70.41 & 50.68 & 96.51 & 1.009 \\ HNN (Ganes et al., 2018) & 81.37 & 87.86 & 80.62 & 94.99 & 0.259 & 89.33 & 28.36 \\ GCN (Wolf & 2017) & 80.38 & 82.51 & 50.97 & 0.651 & 97.30 & 99.30 & 28.29 & 99.07 \\ GAT (Vehicular et al., 2018) & 82.03 & 82.04 & 82.51 & 90.63 & 83.01 & 99.30 & 63.34 & 99.34 \\ SAGE (Hamilton et al., 2017) & 80.02 & -0.43 & 91.04 & 80.89 & 95.61 & -0.263 & 37.05 & 89.24 & 94.00 \\ SOC (Wu et al., 2019) & 89.93 & -0.31 & 97.89 & 82.52 & 96.10 & 99.30 & 99.35 & 99.35 \\ \hline HGN (Liu et al., 2019) & 80.02 & -11.04 & -12.04 & -12.03 & -10.04 & -10.02 & -10.01 & -10.01 \\ HOCN (Chami et al., 2019) & 78.09 & -2.79 & -24.28 & -0.209 & 96.79 & 99.30 & 60.10 & 99.14 & 99.05 \\ HOGM (Zhang et al., 2012) & 76.32 & -3.41 & 94.61 & 99.64 & **80.49** & 99.30 & -0.259 & 99.36 \\ \hline GL (Zhu et al., 2020) & **89.74 & -0.21 & -12.49 & -12.04 & -12.04 & -10.01 & -10.01 & -10.05 \\ \hline \(\kappa\)-GCN (Bachmann et al., 2020) & **89.74 & -0.22 & -14.04 & -12.04 & -12.04 & -10.01 & -10.01 & -10.01 \\ Q-GCN (Xiong et al., 2022) & - & - & 96.30 & -0.22 & **96.56** & **10.37** & **97.11** & **93.30** \\ \hline HamGNN & **99.73 & -0.26 & **99.99 & -0.01** & -12.51 & 0.30 & **99.99** & **99.01** \\ \hline \end{tabular}
\end{table}
Table 1: Node classification accuracy(%). The best, second best, and third best results for each criterion are highlighted in **red**, **blue**, and **cyan**, respectively. “-” indicates the open source code and/or the result is unavailable.

\begin{table}
\begin{tabular}{c|c|c|c|c|c} \hline Method & Disease & Airport & Pubmed & Citeseer & Cora \\ \(\delta\)-hypercholcity & 0.0 & 1.0 & 3.5 & 4.5 & 11.0 \\ \hline MLP & 83.74 & 84.04 & 70.41 & 50.68 & 96.51 & 6.091 & 6.091 \\ HNN (Ganes et al., 2018) & 81.37 & 87.88 & 96.06 & -0.269 & 49.09 & -0.259 & 89.33 & 28.30 \\ GCN (Wolf & 2017) & 80.38 & -25.10 & 90.97 & 0.651 & 97.30 & 99.30 & 28.29 & 99.07 \\ GAT (Vehicular et al., 2018) & 82.03 & -25.10 & 90.63 & -0.31 & -0.633 & 99.30 & -0.593 & 99.34 \\ SAGE (Hamilton et al., 2017) & 80.02 & -0.43 & 91.04 & -0.089 & 95.61 & -0.263 & 99.37 & -0.389 \\ SOC (Wu et al., 2019) & 89.93 & -0.31 & 97.89 & -0.329 & -0.261 & 99.30 & -0.279 & 99.35 \\ \hline HGN (Liu et al., 2019) & 80.02 & -11.04 & -12.04 & -12.04 & -10.01 & -10.01 & -10.01 & -10.01 \\ HOCN (Chami et al., 2019) & 78.09 & -2.79 & -24.28 & -0.209 & 96.79 & -0.303 & -0.10.14 & 99.10 & -10.05 \\ HOGM (Zhang et al., 2012) & 76.32 & -3.41 & 94.61 & 99.64 & **80.49** & 99.30 & -0.259 & 99.36 \\ \hline GL (Zhu et al., 2020) & **89.74 & -0.22 & -14.04 & -12.04 & -12.04 & -10.01 & -10.01 & -10.01 \\ \hline \(\kappa\)-GCN (Bachmann et al., 2020) & **89.74 & -0.22 & -14.04 & -12.04 & -10.01 & -10.01 & -10.01 & -10.01 \\ Q-GCN (Xiong et al., 2022) & - & - & 96.30 & -0.22 & **96.56** & **10.37** & **97.01** & **93.30** \\ \hline HamGNN & **99.73 & -0.26 & **99.99 & -0.01** & -12.51 & 0.30 & **99.99** & **99.01** \\ \hline \end{tabular}
\end{table}
Table 2: Link prediction ROC(%). The best, second best, and third best results for each criterion are highlighted in **red**, **blue**, and **cyan**, respectively. “-” indicates the open source code and/or the result is not available.
difference between those two settings is that the symplectic form in (20) is set to be the special Poincare \(2\)-form while in (23), the symplectic form is a learnable one. We however observe that the two HamGNNs achieve similar performance and the more flexible symplectic form does not improve the model performance. This may be because of the fundamental Darboux theorem (Lee, 2013) in symplectic geometry, which states that we can always find a Darboux coordinate system to give any symplectic form the _Poincare \(2\)-form_. The feature compressing FC may have the network capacity to approximate the Darboux coordinate map, while the flexible learnable \(H_{\mathrm{net}}\) also has the network capacity to get the energy representation under the chosen chart coordinate system.

### Mixed Geometry Dataset

#### 5.2.1 Mixed Airport + Cora Dataset

To understand the node embedding capacity of HamGNN, we created a new mixed geometry graph dataset by combining the Airport (3188 nodes) and Cora (2708 nodes) datasets into a single large graph with 5896 nodes (cf. Table 5). The results are presented in the first row of Table 3.

This new dataset is composed of a mixture of hyperbolic and Euclidean geometries, as the Airport dataset is known to have a more hyperbolic geometry, and the Cora dataset is known to have a more Euclidean geometry, as shown in Figure 1. We choose these two datasets also because they have a similar number of nodes. To standardize the node feature dimension across datasets, we pad the node features with additional zeros. We do not create any _new_ edges in the new dataset so that the new graph contains two disconnected subgraphs: Airport and Cora. We use a 60%, 20%, 20% random split for training, validation, and test sets on this new dataset.

We test our HamGNN against several baselines including GCN, HGCN, GIL, and \(\mathcal{Q}\)-GCN. The results are shown in Table 3 from which we can observe that, in this more complex geometry setting, HamGNN performs the best.

#### 5.2.2 Mixed Airport + Citeseer Dataset

We also include experiments on another new mixed dataset created using Airport and Citeseer. The combination of these two datasets is the same as the previous one. This new dataset therefore has different geometry in each component from Airport or Citeseer. The results are reported in the second row of Table 3. We observe that our HamGNN still performs the best.

These two experiments show that our model can adapt well to the underlying complex geometry.

### Observation of Resilience to Over-Smoothing

As a side benefit of HamGNN, we observe from Table 4 that if more Hamiltonian layers are stacked, HamGNN is still able to distinguish nodes from difference classes, while the other GNNs suffer severe oversmoothing problem (Chen et al., 2020). This may be because when updating node features their energy is constrained to be alone in the Hamiltonian orbit. So, if two nodes are distinct in terms of their energy at the input of HamGNN, they will still be distinguishable by their energy at the output of HamGNN, no matter how many times the feature aggregation operation has been implemented.

**More Experiments.** We kindly refer readers to Appendix F where we include more empirical analysis and visualization.

## 6 Conclusion

In this paper, we have designed a new node embedding strategy from Hamiltonian orbits that can automatically learn, without extensive tuning, the underlying geometry of any given graph dataset even when multiple different geometries coexist. We have demonstrated empirically that our approach adapts better than popular state-of-the-art graph node embedding GNNs to various graph datasets on two graph node embedding downstream tasks.

## Acknowledgments

This research is supported by A*STAR under its RIE2020 Advanced Manufacturing and Engineering (AME) Industry Alignment Fund - Pre Positioning (IAF-PP) (Grant No. A19D6a0053) and the National Research Foundation, Singapore and Infocomm Media Development Authority under its Future Communications Research and Development Programme. The computational work for this article was partially performed on resources of the National Supercomputing Centre, Singapore ([https://www.nscc.sg](https://www.nscc.sg)).

\begin{table}
\begin{tabular}{c c c c c c} \hline \hline Dataset & Models & 3 layers & 5 layers & 10 layers & 20 layers \\ \hline \multirow{3}{*}{Cora} & GCN & 80.29\(\pm\)2.29 & 69.87\(\pm\)1.12 & 26.50\(\pm\)4.68 & 23.97\(\pm\)5.42 \\  & HGCN & 78.70\(\pm\)0.96 & 38.13\(\pm\)6.20 & 31.90\(\pm\)0.00 & 26.23\(\pm\)9.87 \\  & HamGNN (21) & **81.84\(\pm\)0.88** & **81.08\(\pm\)0.16** & **81.40\(\pm\)0.44** & **80.58\(\pm\)0.30** \\ \hline \hline \end{tabular}
\end{table}
Table 4: Node classification accuracy(%) when increasing the number of layers on the Cora dataset.

\begin{table}
\begin{tabular}{c|c|c|c|c|c|c} \hline datasets & GCN & HGCN & GIL & Q-GCN & HamGNN (19) & HamGNN (20) \\ \hline Air-Cora & 5.16\(\pm\)0.65 & 78.72\(\pm\)0.428 & 20.41\(\pm\)1.27 & 90.25\(\pm\)0.81 & **94.82\(\pm\)0.78** & 93.08\(\pm\)0.55 \\ Air-Cite & 0.97\(\pm\)0.347 & 4.65\(\pm\)0.407 & 77.3\(\pm\)1.578 & 76.3\(\pm\)0.28 & **90.03\(\pm\)0.29** & 85.94\(\pm\)0.75 \\ \hline \hline \end{tabular}
\end{table}
Table 3: Node classification on the Mixed Geometry Dataset. First row: mixed Airport+Cora dataset; Second row: mixed Airport+Citeseer dataset.