{'fraction_non_alphanumeric': 0.037437672559442865, 'fraction_numerical': 0.01532945976016813, 'mean_word_length': 4.542423204293708, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Weisfeiler and Leman go Hyperbolic: Learning Distance Preserving Node Representations

In recent years, graph neural networks (GNNs) have emerged as a promising tool for solving machine learning problems on graphs. Most GNNs are members of the family of message passing neural networks (MPNNs). There is a close connection between these models and the Weisfeiler-Leman (WL) test of isomorphism, an algorithm that can successfully test isomorphism for a broad class of graphs. Recently, much research has focused on measuring the expressive power of GNNs. For instance, it has been shown that standard MPNNs are at most as powerful as WL in terms of distinguishing non-isomorphic graphs. However, these studies have largely ignored the distances between the representations of nodes/graphs which are of paramount importance for learning tasks. In this paper, we define a distance function between nodes which is based on the hierarchy produced by the WL algorithm, and propose a model that learns representations which preserve those distances between nodes. Since the emerging hierarchy corresponds to a tree, to learn these representations, we capitalize on recent advances in the field of hyperbolic neural networks. We empirically evaluate the proposed model on standard node and graph classification datasets where it achieves competitive performance with state-of-the-art models.


Introduction
Over the past few years, graph neural networks (GNNs) have been applied with great success to machine learning problems on graphs in various application domains. For instance, in chemistry, much attention has been devoted to deep learning systems for drug screening or design where molecules are represented as graphs [Kearnes et al., 2016]. Likewise, in biology, an issue of high interest is the prediction of the functions of proteins modeled as graphs [You et al., 2021] So far, the field of GNNs has been largely dominated by message passing architectures. Indeed, most of these models share the same basic idea, and can be reformulated into a single common framework, so-called message passing neural networks (MPNNs) [Gilmer et al., 2017]. Specifically, these models follow a message passing procedure, where each node updates its feature vector by aggregating the feature vectors of its neighbors. This iterative scheme runs for a number of steps, and then to compute a feature vector for the entire graph, MPNNs typically employ some permutation invariant readout function such as summing or averaging the feature vectors of all the nodes of the graph.
Recently, much research has focused on measuring the expressive power and limitations of MPNNs. One line of research focused on establishing a connection between MPNNs and the Weisfeiler-Leman (WL) test of graph isomorphism [Weisfeiler and Leman, 1968], a powerful heuristic that can successfully test isomorphism for a broad class of graphs [Babai and Kucera, 1979]. It was shown that standard MPNNs do not have more power in terms of distinguishing non-isomorphic graphs than the WL algorithm [Morris et al., 2019;Xu et al., 2019]. Based on those findings, in the past years, considerable effort was devoted to the development of models that are more powerful than the WL algorithm [Morris et al., 2019;Maron et al., 2019;Chen et al., 2019;Morris et al., 2020b]. Most of the above studies investigate the power of GNNs in terms of distinguishing non-isomorphic graphs. However, in graph classification/regression problems, we are not that much interested in testing whether two (sub)graphs are isomorphic to each other, and it has been observed that stronger GNNs (in the aforementioned sense) do not necessarily outperform weaker GNNs [Dwivedi et al., 2020]. On the other hand, we would like to learn representations which could preserve the similarities or distances of subgraphs/graphs. Deep neural networks are known for being very sensitive to their input. For instance, it has been reported that GNNs are extremely vulnerable to adversarial attacks on the graph structure [Dai et al., 2018]. Such adversarial attacks can mislead a GNN and significantly decrease its performance.
To increase the robustness of standard neural network layers, previous studies proposed to constrain their local Lipschitz constants, i. e., to make their output preserve distances as much as possible for some given distance function [Virmaux and Scaman, 2018]. Furthermore, functions that vary at a slower rate are usually considered simpler and generalize better [Gouk et al., 2021]. Thus, in many settings, it is of paramount importance to learn representations that preserve distances of nodes.
Present Work. To achieve that, in this paper, we propose a new MPNN that learns node representations that respect the distances between nodes, as those are defined by the WL algorithm. The WL algorithm iteratively updates a given node's color (or label) by aggregating the colors of its neighbors. The color classes produced in a given iteration of the algorithm are at least as fine as those produced in the previous iteration. Thus, the sequence of node colors gives rise to a family of nested subsets, which can naturally be represented by a hierarchy (i. e., a tree). It is well-known that trees can be embedded with arbitrarily low distortion into the hyperbolic space, while Euclidean space cannot achieve such a low distortion [Sarkar, 2011]. Therefore, Euclidean MPNNs cannot encode accurately the information contained in the WL hierarchy. Thus, in this paper, we propose a MPNN which delivers the "best of both worlds" from Euclidean space and hyperbolic space. Our proposed model, so-called Weisfeiler-Leman Hyperbolic Network (WLHN), takes into account the distance of the input nodes according to the hierarchy induced by WL, but also according to the representations that emerge from the neighborhood aggregation procedure. To embed the WL tree hierarchy into the hyperbolic space, we capitalize on recent advances in hyperbolic representations using the Poincaré ball model [Ganea et al., 2018b;Chami et al., 2021]. Specifically, we propose an embedding algorithm (DiffHypCon) which generalizes Sarkar's algorithm [Sarkar, 2011] such that the model is differentiable and end-to-end trainable. Experiments on standard benchmark datasets demonstrate that the proposed model either outperforms or achieves performance comparable with that of state-of-the-art models in node and graph classification tasks.

Related Work
GNNs. The first GNN models were proposed several years ago [Sperduti and Starita, 1997;Scarselli et al., 2009;Micheli, 2009]. However, it was not until recently that this family of models attracted significant attention, mainly due to the advent of deep learning [Bruna et al., 2014;Li et al., 2015;Defferrard et al., 2016;Kipf and Welling, 2017;Hamilton et al., 2017;Zhang et al., 2018]. Broadly speaking, these models can be categorized into spectral and spatial approaches depending on which domain the convolutions (neighborhood aggregations) are applied to. Interestingly, the majority of these models follow the same principle and can be reformulated into a single common framework.
These models are known as Message Passing Neural Networks (MPNNs) [Gilmer et al., 2017]. These models employ an iterative message passing procedure, where each node combines the representations of its neighbors with its own representation to compute a new representation. This procedure is typically followed by a readout (or pooling) phase where a feature vector for the entire graph is produced using some permutation invariant function. Several works have proposed extensions and improvements to the message passing procedure of MPNNs. For instance, some works have proposed more expressive or learnable aggregation functions [Murphy et al., 2019;Seo et al., 2019;Dasoulas et al., 2021], schemes that incorporate different local structures or high-order neighborhoods [Abu-El-Haija et al., 2019;Jin et al., 2020;Nikolentzos et al., 2020], approaches that utilize node positional information [You et al., 2019], while others have focused on efficiency [Gallicchio and Micheli, 2020]. Fewer works have focused on the readout phase and have proposed more sophisticated pooling functions Gao and Ji, 2019]. Note that there also exist GNNs that are not MPNNs [Nikolentzos et al., 2022].
Expressive power of GNNs. Recently, a line of research has started exploring the expressive power of GNNs. Several of those studies have investigated how GNNs are related to the WL test of isomorphism and its higher-order variants. For instance, it was shown that standard GNNs are at most as powerful as the WL algorithm in terms of distinguishing non-isomorphic graphs [Morris et al., 2019;Xu et al., 2019]. Other studies proposed families of GNNs whose message passing scheme is equivalent to high-order variants of the WL algorithm, and can thus distinguish more pairs of non-isomorphic graphs than standard MPNNs [Morris et al., 2019[Morris et al., , 2020b. Maron et al. [2019] introduced a class of GNNs which are at least as powerful as the folklore variant of the k-WL graph isomorphism test in terms of distinguishing non-isomorphic graphs, while Chen et al. [2019] also proposed a GNN that is more powerful than the 2-WL algorithm. Barceló et al. [2020] characterized the expressive power of GNNs in terms of classical logical languages based on a connection between first-order logic and the WL algorithm. Sato et al. [2021] showed that random features make standard MPNNs more expressive. Some works propose neighborhood aggregation schemes that take into account all possible node permutations and produce universal graph representations [Murphy et al., 2019;Dasoulas et al., 2020]. However, due to the dramatically high complexity, approximation schemes are necessary. For a comprehensive overview of the expressive power of GNNs, the interested reader is referred to the survey by Sato [2020]. The above studies investigate the power of GNNs in terms of distinguishing non-isomorphic graphs or in terms of how well they can approximate combinatorial problems. On the other hand, the proposed model learns representations that capture the distance of nodes which is largely ignored by the above models.
Hyperbolic Embedding Algorithms and Neural Networks. In the past years, there has been a growing interest in algorithms that can learn hyperbolic embeddings. Such algorithms can be very effective in embedding hierarchical graphs. Nickel and Kiela [2017] proposed a new algorithm for learning hierarchical representations of symbolic data by embedding them into the hyperbolic space using the Poincaré-ball model. Later, they proposed a new optimization approach based on the Lorentz model of hyperbolic space for learning, while they found that learning embeddings in the Lorentz model is more efficient than in the Poincaré-ball model [Nickel and Kiela, 2018]. Ganea et al. [2018a] proposed an algorithm for embedding directed acyclic graphs in the Poincaré-ball model, while Balazevic et al. [2019] proposed an algorithm for embedding hierarchical multi-relational data also in the Poincaré-ball model. Chami et al. [2020] introduced a class of hyperbolic embedding models that can capture both hierarchical and logical patterns. Some other studies have extended deep learning methods to the hyperbolic space by deriving hyperbolic versions of common neural network layers such as recurrent layers and feed-forward layers [Ganea et al., 2018b;Shimizu et al., 2021]. The works closest to ours are the ones proposed by Chami et al. [2019] and by Liu et al. [2019], where the authors generalize Euclidean MPNNs to the hyperbolic space. However, both these models assume that there exists a hierarchical structure in the input data, and they aim to capture such latent hierarchies. On the other hand, in our setting, we have access to an explicit hierarchy produced by the WL algorithm and our model learns representations that respect that hierarchy.

Preliminaries
Notation. Let N denote the set of natural numbers, i. e., {1, 2, . . .}. Then, [n] = {1, . . . , n} ⊂ N for n ≥ 1. Let also { {} } denote a multiset, i. e., a generalized concept of a set that allows multiple instances for its elements. Let G = (V, E) be an undirected graph, where V is the vertex set and E is the edge set. We will denote by n the number of vertices and by m the number of edges, i. e., n = |V | and m = |E|. Let N (v) denote the the neighborhood of vertex v, i. e., the set WL Algorithm. The WL algorithm (also known as color refinement) has been revisited a lot in the past years since it found its way into several machine learning applications [Morris et al., 2021]. The algorithm iteratively computes a coloring of the set of nodes of a graph. Formally, given a graph G = (V, E), a coloring of the set of nodes V is a mapping c : V → N. In other words, a coloring c assigns a number (or color) to every node of the graph. The WL algorithm runs for a number of iterations and is associated with a sequence of colorings c 0 , c 1 , . . . , c T . If the initial coloring c 0 of V is not specified, we either assume a monochromatic coloring, i. e., all vertices have the same color or we define a coloring that maps nodes to their degrees, i. e., c 0 (v) = deg(v). The refinement of a coloring c i is a new coloring c i+1 defined as follows. For every vertex v, collect and lexicographically sort the multiset of colors of its neighbors { {c i (u)|u ∈ N (v)} }. Let M i (v) denote the sequence of colors that emerges from the above operation. Then, define (c i (v), M i (v)) and assign a new color c i+1 (v) to node v by employing a one-to-one mapping from tuples such as (c i (v), M i (v)) to new colors. The new color c i+1 (v) of v depends on the colors of its neighbors and the previous color of v. Thus, in iteration i + 1, vertices v and u receive different colors, if they already had different colors in iteration i, or if the multisets of colors of their neighbors in iteration i are different. Observe also that each coloring c 0 , c 1 , . . . , c T partitions V into color classes, i. e., sets of nodes with the same color. Since the implication c i (v) = c i (u) ⇒ c i+1 (v) = c i+1 (u) holds for all v, u ∈ V and any i ∈ N, the color classes produced in iteration i + 1 are at least as fine as those produced in iteration i. Hence, the sequence of colorings c 0 , c 1 , . . . , c T gives rise to a family of nested subsets, which can naturally be represented by a hierarchy. Figure 1 illustrates the color refinement procedure and the hierarchy that emerges from that procedure.
Poincaré Ball Model. A hyperbolic space is a non-Euclidean space with constant negative curvature. Specifically, the Poincaré ball model of hyperbolic space corresponds to the Riemannian manifold (D n , g D x ) where D n = {x ∈ R n : ||x|| < 1}, i. e., an open unit ball. A Riemannian manifold is a real and smooth manifold equipped with an inner product which is called a Riemannian metric. The above manifold is equipped with the metric g D x = λ 2 x g E where λ x = 2 1−||x|| 2 is known as the conformal factor and g E is the Euclidean metric tensor defined as g E = diag([1, 1, . . . , 1]). The distance between two points x, y ∈ D n is computed as d D (x, y) = (1−||x|| 2 )(1−||y|| 2 ) . For any point x ∈ D n , let T x D n denote the associated tangent space, which is a subset of the Euclidean space. We use the exponential map exp x : T x D n → D n and the logarithmic map log x : D n → T x D n to map points from the tangent space to the hyperbolic space and from the hyperbolic space to the tangent space, respectively. For the tangent vector v = 0, we where ⊕ is the Mobius addition for any x, y ∈ D n (hyperbolic analogous to vector addition in Euclidean space): x ⊕ y = (1+2 x,y +||y|| 2 )x+(1−||x|| 2 )y 1+2 x,y +||x|| 2 ||y|| 2 . Similar to previous works, we choose the origin as the point on the manifold in whose tangent space we operate.  Figure 1: An illustration of (a) a graph G with uniform initial colors c 0 and refined colors c i for i ∈ [4], and of (b) its corresponding WL tree hierarchy H G . The nodes of G are the leaves of the hierarchy.

Embedding WL hierarchies
As already discussed, prior GNNs largely ignore the distance between nodes. It is not thus clear whether two nodes that have structurally dissimilar neighborhoods will obtain dissimilar representations. Likewise, there are no guarantees that two nodes whose neighborhoods are structurally very similar will be mapped close to each other by the intermediate layers of a model. However, capturing such distances between nodes is of paramount importance for machine learning applications since in node classification, similar nodes usually belong to the same class, while in node regression, similar nodes are associated with similar target values. A natural question then is: how is the distance between two nodes defined? Unfortunately, there is no clear answer to the above question. Several different similarity and distance functions were proposed over the past decades for comparing graphs, subgraphs or nodes (i. e., subgraphs centered at nodes), however, most of those functions are hard to compute. For instance, the maximum common subgraph problem is known to be NP-hard [Garey and Johnson, 1979]. Since non-polynomial time computable functions are of no practical use, we focus on a distance function of nodes which can be derived from the hierarchy generated by the WL algorithm. This distance function has already been introduced in previous studies [Kriege et al., 2019;Togninalli et al., 2019]. Formally, let H = (V, E) be a rooted tree representing the hierarchy produced by the WL algorithm.
Producing representations that preserve the distances between the nodes of the hierarchy is then equivalent to embedding the hierarchy into some space. Unfortunately, it is not possible to embed trees into any Euclidean space with arbitrarily low distortion [Linial et al., 1995]. This highlights a limitation of Euclidean GNNs which might fail to accurately encode into the learned node representations the information contained in the hierarchy produced by the WL algorithm. On the other hand, it has been shown that trees can be embedded into the Poincaré disk D 2 with arbitrarily low distortion [Sarkar, 2011]. This motivates the use of models that embed nodes into the hyperbolic space and which can naturally encode the information contained into hierarchical structures. Note, however, that most existing approaches for learning hyperbolic embeddings assume that there exists a latent hierarchy in the input data and their objective is to learn that hierarchy. In our setting, the hierarchy is known in advance, and the objective is to directly embed the hierarchy into some space.
Sarkar's construction and limitations. There exists a construction proposed by Sarkar [2011] which we can leverage to embed the nodes of the hierarchy into the Poincaré disk D 2 and preserve the WL distance with arbitrarily low distortion. However, Sarkar's construction (more details are given in the Appendix) is purely combinatorial and involves no learning (thus cannot be part of an end-to-end learning model). Furthermore, the number of bits of precision used to represent components of the embedded nodes scales linearly with the maximum path length [Sala et al., 2018]. Thus, Sarkar's construction might experience numerical instabilities in the case of trees that contain long paths. Sala et al. [2018] generalized Sarkar's construction from the Poincare disk D 2 to the Poincare ball D d to deal with such problems, but still the new construction is also combinatorial. The above algorithms aim to maintain an equal distance between the children of a given node. In other words, these algorithms equally space the children of a node around a circle. However, in our setting, the hierarchy that is produced by the WL algorithm is not entirely unordered. For instance, given three nodes v, u, w, if deg(v) < deg(u) < deg(w), then the three nodes belong to different color classes and v is more similar to u than to w. Such relationships between the different color classes can be captured by a MPNN model which iteratively updates the representations of the nodes.
Proposed embedding approach. To take into account both the hierarchy produced by the WL algorithm and the aforementioned relationships between the nodes, and to also allow the embedding algorithm to be integrated into the pipeline of a learning model, we generalize Sarkar's algorithm. The proposed construction which is shown in Algo-rithm 1 uses a node placement step which is based on the representations generated by a MPNN (line 4). Specifically, we first normalize those representations such that they lie on the surface of the unit hypersphere (line 6). Then, the emerging representations are scaled by τ and the children are placed in the corresponding positions in the hypersphere (lines 10 and 12). In contrast to previous works [Sarkar, 2011;Sala et al., 2018], we do not space out children of a given node, but we allow them to be placed close to each other in the hyperbolic space. Unfortunately, while angles are preserved between hyperbolic and Euclidean space, distances are not preserved. Thus, the distances derived from the Euclidean representations produced by the MPNN are not preserved. However, as we show next, relationships between distances are preserved.
Proposition 4.1. Let v 1 , v 2 , v 3 denote the children of some node in the WL hierarchy, and let v 1 , v 2 , v 3 ∈ R n denote some transformation of their Euclidean representations such Based on the above result, and since isometric reflections across geodesics preserve hyperbolic distances, the final hyperbolic representations of the children of a given node preserve their corresponding distance relations from the Euclidean space. It is also necessary to separate the children of a given node from the node's parent reflected representation. We can obtain an angle at least equal to θ = π 2 between the parent of the node and each child as follows: we can ensure that the components of the Euclidean representations of the children contain non-negative values by applying the ReLU function. We can then set the representation of the parent equal to [−1/ ∈ R d and rotate the hypersphere so that this vector is located at node's parent actual representation (line 8).

Model Architecture
Neighborhood aggregation. We assume that the initial feature vectors of the nodes are from a countable universe. For finite graphs, node representations at deeper layers of a MPNN are also from a countable universe [Xu et al., 2019]. Thus, each node representation can be mapped to a color. Then, the node representations that emerge at the different iterations of the message passing procedure correspond to different colors and from those colors we can construct a hierarchy. In this paper, we use the neighborhood aggregation scheme of the GIN model [Xu et al., 2019] which is defined as follows: For each depth t ∈ [T ] of the hierarchy, each node belongs to one of the color classes of that depth. Two nodes v, u belong to the same color class of depth t if h Algorithm 1 Proposed Differentiable Hyperbolic Construction (DiffHypCon) 1: Input: Euclidean representations of current iteration H (t) , hyperbolic representations of previous two iterations representations of children of a } 5: { store hyperbolic representations in corresponding rows of Z (t) } 14: end for 15: Output: Embedded vectors Z (t) in D n holds. Note that the GIN model can be as powerful as the WL algorithm in terms of distinguishing non-isomorphic graphs, and thus the emerging hierarchy could potentially be identical to the one produced by the WL algorithm (i. e., in case proper weights are found).
Note that we do not construct a single hierarchy consisting of the color classes of all the nodes of all graphs. Instead, a different hierarchy is constructed and embedded for each batch of samples. This does not pose any problem since the proposed construction does not space out children and hence, it does not make any assumptions about the number of children of a given node. Let H note the set that contains the representations of all the nodes of the graphs contained in a batch B after the t-th iteration of the neighborhood aggregation procedure where V B is the set of nodes of the graphs contained in B. These representations form the color classes of the hierarchy produced by the WL algorithm after its t-th iteration. Specifically, there is a one-to-one mapping between those vectors and the nodes of the hierarchy whose path distance from the root node is equal to t (or t + 1 in the case of non-monochromatic initial coloring). Furthermore, |H (t) | ≤ |H (t+1) | holds, and each representation has a parent. For instance, the parent of h . Therefore, we can use the construction presented above to embed those representations into the hyperbolic 1: Input: Adjacency matrix A and matrix of node features X of graph G, number of iterations T , scaling factor τ 2: Z (−1) ← 0 { root of the tree } 3: { initial hyperbolic representations } 5: for i = 1 to T do 6: is the hyperbolic representation of node v in the t-th iteration of the neighborhood aggregation procedure. Thus, each neighborhood aggregation operation is followed by an embedding phase where the emerging node representations are mapped to the Poincaré ball using the proposed construction. The proposed method is illustrated in Algorithm 2.
Readout. To perform graph classification, we map the output of the last neighborhood aggregation layer of the model to the tangent space using the logarithmic map and aggregate the node representations using the sum operator.
Then, the graph representation h G is fed to further layers (e. g., Euclidean multinomial logistic regression) to produce the output. Alternatively, the rest of the computations can also be realized in the hyperbolic space, i. e., we can directly classify samples on the hyperboloid manifold using the hyperbolic multinomial logistic loss [Ganea et al., 2018b]. In this case, we need to map the tangent vector back to the hyperbolic space using the exponential map: z G = exp 0 h G . Thus, this model first maps node representations to the tangent space, it performs the node aggregation in the tangent space, and then maps the graph representation back to the hyperbolic space. In preliminary experiments, we observed that this approach performs similarly to Euclidean classification.
Note that if the generated node representations accurately capture the distances of nodes, then the emerging graph representations also capture the distances of graphs.
Proposition 4.2. Let G 1 = (V 1 , E 1 ), G 2 = (V 2 , E 2 ) denote two graphs. Without loss of generality, we assume that the two graphs have the same number of nodes, i. e., |V 1 | = |V 2 | = n. Let v i , u j , i, j ∈ [n] denote the representations of the nodes of the two graphs. Let also h G1 and h G2 denote the vector representations of G 1 and G 2 , respectively, which emerge by applying the sum operator to the representations of the nodes of the two graphs. Then, we have that: where T ∈ {0, 1} n×n , T 1 = 1 and 1 T = 1.
The above Proposition implies that given two sets of vectors, the difference of the output of the sum aggregator for the two sets can be bounded via the stable matching between the two sets. Thus, if two graphs consist of similar nodes (in terms of their representations), the graph representations will also be similar to each other. Model depth. It has been reported that MPNNs usually do not benefit from more than few neighborhood aggregation layers. This is attributed to different phenomena such as over-smoothing . We should note that in contrast to Euclidean MPNNs, a large number of neighborhood aggregation layers does not have a negative impact on the performance of the proposed model since the learned node representations respect the structure of the WL hierarchy, and thus distances between nodes are preserved no matter how large the depth of the hierarchy is. Furthermore, no residual connections are required since the final node representations encode the entire structure of the WL hierarchy, i. e., the whole history of previous node representations.
Computational complexity. The proposed model uses the neighborhood aggregation mechanism of GIN to update the representations of the nodes and then uses the DiffHypCon algorithm to embed the nodes into the hyperbolic space. DiffHypCon places points into the unit hypersphere in linear time. Therefore, the computational complexity of the model is comparable with that of other standard MPNNs and is in the order of O(T md 2 ) where d denotes the dimension of the node features. We empirically measured the running time of the proposed model on real-world datasets and we report the results in the Appendix.

Synthetic Datasets
Qualitative results. To empirically verify that the proposed model can encode the hierarchy H G created by the WL algorithm, we applied it to the graph shown in Figure 1. Specifically, we randomly initialized the parameters of the model and we performed a feedforward pass (no training was performed) to generate the graph representations that lie in the hyperbolic space (τ was set equal to 1). Then, given the hyperbolic representations of the nodes of the hierarchy, we computed the distance between the different nodes. The emerging distances were then compared against the WL distance defined in section 4. We computed the Pearson correlation coefficient which was found to be equal to 0.98. We also retrieved the Euclidean representations of the nodes  (produced by the GIN model). We computed the Euclidean distance between those representations and we also compared the emerging distances against the WL distance. In this case, the correlation turned out to be much smaller and equal to 0.35. Hence, the emerging hyperbolic representations lead to distances that are much more correlated with the distances derived from the WL hierarchy compared to those that emerge from the Euclidean representations. We also produced a heatmap (shown in Figure 2) that illustrates the distances between the nodes of the hierarchy which were computed based on the nodes' hyperbolic representations.
We can see that the obtained distances are similar to the corresponding WL distances, while the structure is preserved since the nodes closest to a given node are its direct neighbors in the hierarchy. For example, the blue and the red nodes are neighbors in the hierarchy H G while the distance of their hyperbolic representations is small.
Quantitative results. We also constructed two node regression datasets where the task is to predict the density of the nodes' ego-networks and their effective size, respectively. The first dataset contains 10 Barabási-Albert graphs. The number of nodes of each graph is set to 1, 000 and the number of edges to attach from a new node to existing nodes is set to 5 or 10 (i. e., m = 5 or m = 10). We compute the target of each node as follows: we extract the nodes' ego-networks of radius 2 and compute their densi-ties (the density of a graph consisting of n nodes and m edges is equal to m/n). The second dataset contains 10 Erdös-Rényi graphs. The number of nodes is set to 1, 000 and the probability for edge creation is set to either 0.008 or 0.01 (i. e., p = 0.008 or p = 0.01). The target of each node is set equal to its effective size which for a node v is defined as e where t is the number of edges between neighbors of v. We split the nodes of each graph into training/validation/test sets with a ratio of 60%/20%/20%. We trained the models by minimizing the mean squared error loss and we report in Table 1 the obtained results. We can see that the proposed model outperforms the two baselines (GCN [Kipf and Welling, 2017] and GIN [Xu et al., 2019]) under all setting. The two predicted properties (i. e., density and effective size) are related to the structure of the neighborhood of each node, and thus the experimental results verify our claim that WLHN can better capture the structural distance of nodes.

Real-World Datasets
Datasets. We evaluated the proposed model on five node classification datasets: Cornell, Texas, Wisconsin, Squirrel, and Actor. The first three datasets are extracted from the WebKB dataset and the rest of the datasets from Wikipedia, and have been employed in previous studies [Pei et al., 2020;Rozemberczki et al., 2021]. We also evaluated the proposed model on standard graph classification datasets derived from the TUD repository [Morris et al., 2020a]. For datasets that contain graphs whose nodes are annotated with discrete labels, we map those labels to one-hot vectors. We also evaluated the proposed model on two graph classification datasets from the Open Graph Benchmark (OGB) [Hu et al., 2020], a collection of large-scale datasets. Specifically, we used the following two molecular property prediction datasets: ogbg-molhiv and ogbg-molpcba.
Experimental setup. For the node classification datasets, we randomly split each dataset into training/validation/test sets with a ratio of 60%/20%/20%, while the whole process was repeated 10 times. For the standard graph classification datasets, we perform 10-fold cross-validation, where within each fold a model's hyperparameters are selected based on a 90%/10% split of the training set. We use the evaluation framework from Errica et al. [2020], thus employing the   same split for each evaluated method. For the two datasets from OGB, we use the standard splits associated with those datasets. A detailed description of the hyperparameter selection approach is given in the Appendix.
Results. Table 2 presents the performance of the different models on the five node classification datasets. We observe that WLHN is the best performing method since it outperforms the baselines on 4 out of the 5 datasets. On some datasets, the proposed model improves significantly over the baseline models. For instance, on the Cornell, Wisconsin and Actor datasets, the proposed model yields respective absolute improvements of 3.78%, 2.74% and 4.20% in accuracy over the best competitor. Table 3 illustrates the average accuracies and corresponding standard deviations for the TUD datasets. We observe that the proposed model outperforms the baselines on 6 out of the 10 datasets. On MUTAG, PROTEINS and ENZYMES, the proposed model provides significant performance gains compared to the baseline models. The proposed model embeds the WL tree hierarchy produced by the GIN model into the hyperbolic space to produce expressive node representations. On the other hand, our two main competitors, GIN and HGCN, have some serious limitations. GIN operates in the Euclidean space where hierarchical structures cannot be accurately captured. HGCN embeds nodes in the hyperbolic space but it does not have explicit access to the WL hierarchy, thus leading to less informed representations. Finally, Table 4 illustrates the performance of the proposed model on the two OGB datasets. On ogbg-molhiv, the proposed model outperforms all the baselines, most of them by wide margins. On the other hand, on ogbg-molpcba, while the proposed model is still the best-performing method, some of the baselines achieve an average precision close to that of WLHN. Overall, the results indicate that the proposed model exhibits competitive performance in both node and graph classification tasks.
Additional Experiments. We further investigate the robustness of WLHN under structural and feature noise. We observe that noise does not have a very large impact on the model's performance, while WLHN also significantly outperforms the baselines. We also experimentally demonstrate that WLHN does not suffer from over-smoothing even if the number of neighborhood aggregation layers is large, since the proposed hyperbolic embedding algorithm preserves the distances of the nodes. We present the results in the Appendix.

Conclusion
In this paper, we proposed a new model which preserves distances of nodes/graphs, a problem that has not been investigated thoroughly so far. Our model learns node representations in the hyperbolic space which respect the hierarchy generated by MPNNs. We evaluated the proposed model on synthetic and real-world datasets. Our results demonstrate that the proposed model can indeed encode meaningful distances in the learned representations, while it achieves high levels of performance in node and graph classification tasks. The Appendix is organized as follows. In Sections A and B, we prove Proposition 4.1 and Proposition 4.2, respectively. In Section C, we give more details about Sarkar's construction. In Sections D and E, we present a set of additional experiments we performed (measuring running time, robustness to noise and to the problem of oversmoothing) and describe the hyperparameters we used in our experiments, respectively. In Section F, we provide statistics and descriptions of the datasets we used in our experiments.

B Proof of Proposition 4.2
Let G 1 = (V 1 , E 1 ) and G 2 = (V 2 , E 2 ) denote two graphs. Without loss of generality, we assume that the two graphs have the same number of nodes, i. e., |V 1 | = |V 2 | = n. Let H 1 and H 2 denote the set of vector representations of the nodes of G 1 and G 2 , respectively. Then, H 1 and H 2 have the same cardinality, thus |H 1 | = |H 2 | holds. Let also H 1 = {v 1 , . . . , v n } and H 2 = {u 1 , . . . , u n }. Then, we have: which concludes the proof. Note that the right part of the last equality above is the solution of the following bipartite matching problem. We construct a complete bipartite graph where the first and second partition contain nodes of G 1 and G 2 , respectively, while the weight of an edge is some positive value inversely proportional to the Euclidean distance of the vector representations of the two endpoints. Then, matrix T is the solution of the maximum bipartite matching on the above graph, thus it can be seen as a function that maps nodes of G 1 to matched nodes of G 2 .

C Sarkar's Construction
There exists a construction proposed by Sarkar which can embed the nodes of a tree into the Poincaré disk D 2 and preserve the distances between nodes with arbitrarily low distortion [Sarkar, 2011]. Sarkar's construction, shown in 3, is purely combinatorial. The inputs to the Algorithm are a scaling factor τ , a node a (of degree deg(a)) from the tree and its parent node b. These nodes have already been embedded into D 2 , and the objective is to to also embed the deg(a) − 1 children of a into that space. Specifically, to place the children c 1 , c 2 , . . . , c deg(a)−1 of node a into D 2 , the algorithm performs the following steps: (1) it reflects the embeddings of nodes a and b (f (a) and f (b), respectively) across a geodesic such that f (a) is mapped onto the origin (i. e., 0) and f (b) is mapped onto some point z; (2) it places the children of node a to vectors y 1 , y 2 , . . . , y deg(a)−1 equally spaced around a circle with radius exp(τ )−1 /exp(τ)+1 (which is a circle of radius τ in the hyperbolic metric), and maximally separated from node b's reflected embedding z; and (3) it reflects all of the points back across the geodesic, thus the origin is mapped back onto f (a), while z is mapped onto f (b). To embed the entire tree, we follow a recursive approach: the root is placed at the origin and its children in a circle around it, and then the children of all nodes are embedded in the space until no more nodes are left to be placed. This complexity of the entire algorithm in linear in the number of nodes of the tree.
One limitation of Sarkar's construction is that the number of bits of precision used to represent components of the embedded nodes scales linearly with the maximum path length [Sala et al., 2018]. Thus, Sarkar's construction might experience numerical instabilities in the case of trees that contain long paths. Sala et al. [2018] generalized Sarkar's construction from the Poincare disk D 2 to the Poincare ball D d to deal with such problems. The new construction can produce embeddings of higher quality in the case of bushy trees, i. e., trees whose maximum degree is large.

D.1 Running Time
In Table 5, we present the average running time per epoch of the proposed model and two baselines (GIN [Xu et al., 2019] and HGCN [Chami et al., 2019]) on the 10 graph classification datasets (experiments performed on an Nvidia Titan Xp GPU). As discussed in Section 4, the time complexity of the model is in the same order of magnitude as those of standard MPNNs. We observe that the running time of WLHN is comparable with those of GIN and HGCN and thus not prohibitive for real-world problems.

D.2 Robustness Experiments
We have performed experiments investigating the performance of our model after applying different types of perturbations on the input. Specifically, we use the Cornell, Texas, Wisconsin datasets and we apply Gaussian noise with µ = 2 and σ = 0.2 to the node features of the training graphs. We report in the Table 6 the mean classification accuracy and the standard deviation when 10% or 20% of the nodes are affected by the noise, across 10 different splits. We observe that noise   has not a very large impact on the performance of the proposed model, while WLHN also significantly outperforms the baselines on all 3 datasets.
We further investigate the performance of our model after applying structural noise. Specifically, we sample randomly nodes and we add k edges between nodes with different labels. We report the results in Table 7 for k ≈ 0.1|E| and k ≈ 0.2|E|, where |E| is the number of the edges in the graph, across 10 different random splits. We compare the performance of WLHN against that HGCN (Hyperboloid). We observe that structural noise also does not have a large impact on the performance of the proposed model, and thus that WLHN is more robust to noise than HGCN (Hyperboloid). These results partially validate our intuition that the proposed model is more robust to noise than existing models that do not capture the nodes' structural distance in the generated representations.

D.3 Oversmoothing Experiments
From the theoretical point of view, it is known from the proposed construction that the distance between two nodes in the hyperbolic space cannot decrease as more neighborhood aggregation steps are performed. Thus, over-smoothing cannot occur, i.e., the node representations of structurally different nodes (according to WL) will never converge to indistinguishable vectors. We also experimentally demonstrate that our model can have increased depth without suffering from over-smoothing.
We performed experiments on the MUTAG and PROTEINS datasets and report the results in Figure 3. We observe that even with a large number of layers, there is no decrease in the performance of our model.

E Hyperparameters
For a fair comparison, we used the same hyperparmaters setup between our model and HGCN [Chami et al., 2019]. For all datasets, the hidden dimension size of the GNN layers was chosen from {32, 64, 128}. The number of neighborhood aggregation layers was chosen from {1, . . . , 5} and we used batch normalization. The emerging graph representations were fed to a multi-layer perceptron consisting of two hidden layers with a hidden-dimension size of 128 and 64, respectively. We used the ReLU as activation function. We used dropout with ratio chosen from {p = 0.0, 0.5} between the hidden layers. We trained the model by minimizing the cross-entropy loss. We choose between two optimizers, Adam and RiemannianAdam, with learning rate in {10 −3 , 10 −2 }. We trained the networks for 300 epochs. The experiments were run on a Intel(R) Xeon(R) CPU E5-1607 v2 @ 3.00GHz processor with 128GB ram and a Nvidia TITAN Xp GPU.

F Datasets
We evaluated the proposed model on 10 publicly available graph classification datasets including 5 bio/chemo-informatics datasets: MUTAG, D&D, NCI1, PROTEINS and ENZYMES, as well as 5 social interaction datasets: IMDB-BINARY, IMDB-MULTI, REDDIT-BINARY, REDDIT-MULTI-5K and COLLAB [Morris et al., 2020a]. A summary of the 10 datasets is given in Table 8  is to predict whether or not each chemical compound has mutagenic effect on the Gram-negative bacterium Salmonella typhimurium [Debnath et al., 1991]. ENZYMES contains 600 protein tertiary structures represented as graphs obtained from the BRENDA enzyme database. Each enzyme is a member of one of the Enzyme Commission top level enzyme classes (EC classes) and the task is to correctly assign the enzymes to their classes [Borgwardt et al., 2005]. NCI1 contains more than four thousand chemical compounds screened for activity against non-small cell lung cancer and ovarian cancer cell lines [Wale et al., 2008]. PROTEINS contains proteins represented as graphs where vertices are secondary structure elements and there is an edge between two vertices if they are neighbors in the amino-acid sequence or in 3D space. The task is to classify proteins into enzymes and non-enzymes [Borgwardt et al., 2005]. D&D contains over a thousand protein structures. Each protein is a graph whose nodes correspond to amino acids and a pair of amino acids are connected by an edge if they are less than 6 Ångstroms apart. The task is to predict if a protein is an enzyme or not [Dobson and Doig, 2003]. IMDB-BINARY and IMDB-MULTI were created from IMDb, an online database of information related to movies and television programs. The graphs contained in the two datasets correspond to movie collaborations. The vertices of each graph represent actors/actresses and two vertices are connected by an edge if the corresponding actors/actresses appear in the same movie. Each graph is the ego-network of an actor/actress, and the task is to predict which genre an ego-network belongs to [Yanardag and Vishwanathan, 2015]. REDDIT-BINARY and REDDIT-MULTI-5K contain graphs that model the social interactions between users of Reddit. Each graph represents an online discussion thread. Specifically, each vertex corresponds to a user, and two users are connected by an edge if one of them responded to at least one of the other's comments. The task is to classify graphs into either communities or subreddits [Yanardag and Vishwanathan, 2015]. COLLAB is a scientific collaboration dataset that consists of the ego-networks of several researchers from three subfields of Physics (High Energy Physics, Condensed Matter Physics and Astro Physics). The task is to determine the subfield of Physics to which the ego-network of each researcher belongs [Yanardag and Vishwanathan, 2015].
We also evaluated the proposed kernel on two datasets from the Open Graph Benchmark (OGB) [Hu et al., 2020], a collection of large-scale and diverse benchmark datasets for machine learning on graphs. They are both molecular property prediction datasets that are adopted from the MoleculeNet . ogbg-molhiv consists of 41, 127 molecules and corresponds to a binary classification dataset where the task is to predict whether a molecule inhibits HIV virus replication or not. The average number of vertices per graph is equal to 25.5, while the average number of edges is equal to 27.5. The dataset is split into training/validation/test sets with a ratio of 80/10/10. For ogbg-molpcba dataset, the class balance is very skewed (only 1.4% of data is positive) and the dataset contains 128 classification tasks. Therefore, we use the Average Precision (AP) averaged over the tasks as the evaluation metric. It contains 437, 929 molecule graphs with average number of vertices per graph equals to 26.0, average number of edges equals to 28.1. For both datasets, the molecules in the training, validation and test sets are divided using a scaffold splitting procedure that splits the molecules based on their two-dimensional structural frameworks. The scaffold splitting attempts to separate structurally different molecules into different subsets.
Finally, we evaluated the proposed model on five node classification datasets: Cornell, Texas, Wisconsin, Squirrel, and Actor. The first three datasets are extracted from the WebKB dataset and the rest of the datasets from Wikipedia, and have been employed in previous studies [Pei et al., 2020;Rozemberczki et al., 2021]. A summary of the three datasets is given in Table 9.