Hierarchical Position Embedding of Graphs with Landmarks and Clustering for Link Prediction

We consider undirected graph $G=(V,E)$ where $V$ and $E\subseteq V\times V$ denote the set of vertices and edges, respectively. Let $N$ denote the number of nodes, or $N=|V|$. $\mathbf{A}\in\mathbb{R}^{N\times N}$ denotes the adjacency matrix of $G$. $d(v,u)$ denotes the geodesic (shortest-path) distance between $v$ and $u$. Node attributes are defined as $X=\{x_{1},...,x_{N}\}$ where $x_{i}\in\mathbb{R}^{n}$ denotes the feature vector of node $i$. We consider methods of embedding $X$ into the latent space as vectors $Z=\{z_{1},...,z_{N}\},z_{i}\in\mathbb{R}^{m}$. We study the node-pair-level task of estimating the link probability between $u$ and $v$ from $z_{u}$ and $z_{v}$.

The distances between nodes provide rich information on the graph structures. For example, a connected and undirected graph can be represented by a finite metric space with its vertex set and the inter-node distances, and the position of a node can be defined as its relative distances to other nodes. However, computing and storing shortest paths for all pairs of nodes incur high complexity.

Instead, we select a small number of representative nodes called landmarks. The landmarks are denoted by $\lambda_{1},\cdots,\lambda_{K}\in V$ where $K$ denotes the number of landmarks. For node $v$, we define $K$-tuple of distances to landmarks as distance vector (DV):

$D(v)$ can be used to represent the (approximate) position of node $v$. If the positions of $u$ and $v$ are to be used for link prediction, one can combine $D(u)$ and $D(v)$ (or the embeddings thereof) to estimate the probability of link formation.

The key question is, how much information $D(\cdot)$ has on the inter-node distances within the graph, so that positional information can be (approximately) recovered from $D(\cdot)$. For example, from triangle inequality, the distance between nodes $u$ and $v$ are bounded as

which states that the detour via landmarks, i.e., from $u$ to $\lambda_{i}$ to $v$, is longer than $d(u,v)$, but the shortest detour, or the minimum component of $D(u)+D(v)$, may provide a good estimate of $d(u,v)$. Indeed, it was empirically shown that the shortest detour via base nodes (landmarks) can be an excellent estimate of distances in Internet (Ng and Zhang, 2002).

We examine the achievable accuracy of distance estimates based on the detours via landmarks. Previous works on distance estimates using auxiliary nodes relied on heuristics (Potamias et al., 2009) or complex algorithms for arbitrary graphs (Linial et al., 1995; You et al., 2019). However, real-world graphs have structure, e.g., preferential attachment or certain distribution of node degrees. Presumably, one can exploit such structure to devise efficient algorithms. The key design questions are: how to select good landmarks, and how many of them? Drawing upon the theory of network science, we analyze a well-known class of random graphs, derive the detour distances via landmarks, and glean design insights from the analysis.

The framework by (Fronczak et al., 2004) provides a useful tool for analyzing path lengths for a wide range of classes of random graphs as follows. Consider a random graph where the probability of the existence of an edge for node $i$ and $j$ denoted by $q_{ij}$ is given by

where $h_{i}$ is tag information of node $i$, e.g., the connectivity or degree of the node, and $\beta$ is a parameter depending on the network model. From the continuum approximation (Albert and Barabási, 2002) for large graphs, tag information $h$ is regarded as a continuous random variable (RV) with distribution $\rho(\cdot)$. For some measurable function $f$, let $\langle f(h)\rangle$ denote the expectation of $f(h)$ of a node chosen at random, and $\langle f(h)\rangle_{Q}$ denote the expectation of $f(h)$ of landmarks chosen under some distribution $Q$.

The Erdős-Rényi (ER) model (Erdős et al., 1960) is a classical random graph in which every node pair is connected with a common probability. For ER graphs, model parameters $\beta$ and $h$ can be set as $\beta=\langle k\rangle N$ and $h\equiv\langle k\rangle$ respectively, where $\langle k\rangle$ denotes the mean degree of nodes (Boguná and Pastor-Satorras, 2003). From (3), we have $q_{ij}=\langle k\rangle/N$, i.e., the probability of edge formation between any node pair is constant. The node degree follows the Poisson distribution with mean $\langle k\rangle$ for large $N$.

We apply our analysis to ER graphs as follows. By plugging in the model parameters to (2), the average distance of the shortest detour via landmarks in ER graphs, denoted by $\bar{l}_{\textrm{ER}}$, is bounded above as

asymptotically in $N$. Meanwhile, the average inter-node distance in ER graphs denoted by $\bar{l}^{*}_{\textrm{ER}}$, is given by (Boguná and Pastor-Satorras, 2003; Fronczak et al., 2004)

By comparing (5) and (6), we observe that the detour via landmarks incurs the overhead of at most factor 2. This is because, the nodes in ER graphs appear homogeneous, and thus the path length to and from landmarks are on average similar to the inter-node distance. Thus, the average distance of a detour will be twice the direct distance. However, (5) implies that the minimum detour distance can be reduced by using multiple $(K(N)>1)$ landmarks. The reduction can be substantial, e.g., if $K(N)=N^{1-\varepsilon}$ for some $\varepsilon\in(0,1)$, we have, from (5),

For example, selecting $\sqrt{N}$ landmarks guarantees a 1.5-factor approximation of the shortest path distance. By making $\varepsilon$ close to 0, we get arbitrarily close to the shortest path distance.

The Barabási-Albert (BA) model (Barabási and Albert, 1999) generates random graphs with preferential attachment. BA graphs are characterized by continuous growth over time: initially there are $m$ nodes, and new nodes arrive to the network over time. Due to preferential attachment, the probability of a newly arriving node connecting to the existing node is proportional to its degree. The probability of an edge in BA graphs is shown to be (Boguná and Pastor-Satorras, 2003)

with $h_{i}=1/\sqrt{t_{i}}$ and $\beta=\frac{2}{m}$, where $t_{i}$ is the time of arrival of node $i$. Since the probability of a newly arriving node connecting to node $i$ is proportional to $h_{i}$, the degree of nodes with large $h_{i}$ is likely to be high. For large $N$, the distribution of $h$ is derived as (Fronczak et al., 2004)

By applying $\rho(\cdot)$ to (2), we bound the average path length with landmarks denoted by $\bar{l}_{\textrm{BA}}$ as (4). Importantly, unlike ER graphs, there exists a landmark selection strategy which achieves the asymptotically optimal distance, despite using a small number of landmarks relative to the network size.

The proof of Theorem 3 is in Appendix A.3. We claim the optimality of the strategy in Theorem 3 in the following sense. Since $K(N)$ is slowly increasing in $N$, e.g., $K(N)=O(\log N)$, a feasible strategy for Theorem 3 is to choose $O(\log N)$ landmarks from top–$O(\log^{2}N)$ degree nodes. For such $K(N)$, the numerator of the RHS of (8) is $\approx\log N$ for large $N$. Meanwhile, the average length of shortest paths in BA graphs is given by (Cohen and Havlin, 2003)

By comparing (8) and (9), we conclude that $\bar{l}_{\textrm{BA}}$ is asymptotically equal to $\bar{l}^{*}_{\textrm{BA}}$. This implies that the shortest detour via landmarks under our strategy has on average the same length as the shortest path in the asymptotic sense.

The key design questions we would like to address are: what kind of, and how many, landmarks should be selected?

In order to better capture local graph structure and evenly distribute landmarks over the graph, we propose to partition the graph into clusters such that the nodes within a cluster tend to be densely connected, i.e., close to one another. Next, we propose to pick the node with the highest degree within each cluster as the landmark, as suggested by the analysis. Each landmark represents the associated cluster, and the distance between nodes can be captured by distance vectors associated with the cluster landmarks. We empirically find that the combination of landmark selection and clustering yields improved results.

Graph clustering is a task of partitioning the vertex set into disjoint subsets called clusters, where the nodes within a cluster tend to be densely connected relative to those between clusters (Schaeffer, 2007). We will use FluidC graph clustering algorithm proposed in (Parés et al., 2017) as follows. Suppose we want to partition $G$ to $K$ clusters. FluidC initially selects $K$ random central nodes, assigns each node to clusters $C_{1},\cdots,C_{K}$, and iteratively assigns nodes to clusters according to the following rule. Node $u$ is assigned to cluster $P_{k^{*}(u)}$ such that

where $\mathcal{N}(u)$ denotes the neighbors of $u$. This assignment rule prefers community candidates of small sizes (denominator), which results in well-balanced cluster sizes. The rule also prefers communities already containing many neighbors of $u$ (numerator), which makes clusters locally dense. As a result, FluidC generates densely-connected cohesive clusters of relatively even sizes. In addition, FluidC has complexity $O(|E|)$, and thus is highly scalable (Parés et al., 2017). Importantly, the number of clusters $K$ can be specified beforehand, which is crucial because $K$ is a hyperparameter in our method.

For each cluster, we select the node with the highest degree as the landmark, where $\lambda_{k}$ denotes the landmark of cluster $C_{k}$. For node $v$, we compute the distances to landmarks $\lambda_{1},\cdots,\lambda_{K}$ to yield distance vector $D(v)$. Since $G$ may be disconnected, we extend the definition of distance $d(v,\lambda)$ such that, in case there exists no path from $v$ to $\lambda$, $d(v,\lambda)$ is set to $\bar{d}+1$ where $\bar{d}$ is the maximum diameter among the connected components. As stated earlier, we set number of clusters $K=\eta\log N$ where integer $\eta$ is a hyperparameter.

A supporting argument for incorporating clustering into analysis can be made for BA graphs. Theorem 3 states that, choosing $K(N)=O(\log N)$ landmarks from top-$(\log N)\cdot K(N)=O(\log^{2}N)$ nodes in degrees is optimal. We show that, $\log N$ landmarks chosen after FluidC clustering are indeed within top-$(\log N)^{2}$ degree centrality through simulation. Table 1 compares the rank of degrees of landmarks selected after FluidC clustering versus the rank of top-$(\log N)^{2}$ degree nodes in BA graphs. We observe that the landmarks selected after clustering have sufficiently large degrees for optimality, i.e., all of their degrees are within top-$(\log N)^{2}$.

The result can be explained as follows. Power-law graphs are known to have modules which are densely connected subgraphs and are sparsely connected to each other (Goh et al., 2006). A proper clustering algorithm is expected to detect modules. Moreover, each cluster contains a relatively large number of nodes, because the number of clusters is relatively small $(O(\log N))$. Thus, each cluster is likely to contain nodes with sufficiently high degrees which are, according to Theorem 3, good candidates for landmarks.

One can expect similar effects from FluidC clustering for real-world graphs. Given the small number $(O(\log N))$ of clusters as input, the algorithm will detect locally dense clusters resulting in landmarks which have high degrees and tend to be evenly distributed over the graph. Thus, landmark selection after clustering can be a good heuristic motivated by the theory.

For the nodes within the same cluster, we augment the nodes’ embeddings with information identifying that they are members of the same community, which we call membership. The membership is extracted from landmarks, and is based on the relative positional information among landmarks. Thus, not only the nodes within the same cluster have the same membership, but also the nodes of neighboring clusters have “similar” membership. The membership of a node is encoded into a membership vector (MV) which is combined with DV in computing the node embeddings.

We consider the graph consisting only of landmark nodes, and use graph Laplacian (Belkin and Niyogi, 2003) to encode their relative positions as follows. The landmark graph uses the edge weight between landmarks $u$ and $v$ which is set to $e_{uv}=\exp(-d(u,v)^{2}/T)$ where $T$ denotes normalizing parameter of heat kernel. Let $\hat{A}\in\mathbb{R}^{K\times K}$ denote the weighted adjacency matrix. The normalized graph Laplacian is given by $L=I-\Delta^{-\frac{1}{2}}\hat{{A}}\Delta^{-\frac{1}{2}}$ where degree matrix $\Delta$ is given by $\Delta_{ii}=\sum_{j}\hat{A}_{ij}$. The eigenvectors of $L$ are used as MVs, i.e., for node $v\in C_{k}$, the MV of node $v$, denoted by $M(v)$, is the eigenvector of $L$ associated with landmark $\lambda_{k}$. MV provides positional information in addition to DV, using the graph of upper hierarchy, i.e., landmarks. We use a random flipping of the signs of eigenvectors to resolve ambiguity (Dwivedi et al., 2020). A problem with Laplacian encoding is the time complexity: one needs the eigendecomposition whose complexity is cubic in network size. Thus, previous approaches had limitations, i.e., used only a subset of eigenspectrum (Dwivedi et al., 2020; Kreuzer et al., 2021). However, in our method, the landmark graph has small size, i.e., $K=O(\log N)$. Thus, we can exploit the full eigenspectrum, enabling more accurate representation of landmark graphs.

We propose cluster-group encoding as an additional model optimization as follows. Several neighboring clusters are further grouped into a macro-cluster. The embeddings of nodes in a macro-cluster use an encoder specific to that macro-cluster. The motivation is that, using a separate encoder per local region may facilitate capturing attributes specific to local structures or latent features of the community, which is important for link prediction.

For cluster $C_{k}$, let $\Phi_{i(k)}$ denote the macro-cluster which contains $C_{k}$ where $i(k)$ denote the index of the macro-cluster. The node embedding $z_{v}$ for node $v$ is computed as

where $x_{v}$ is the input node feature, $M(v)$ and $D(v)$ are MV and DV, $\phi$ is membership/distance encoder, and $\oplus$ denotes concatenation. $f_{i(\cdot)}$ denotes the encoder specific to macro-cluster $\Phi_{i(\cdot)}$. MLPs are used for $\phi(\cdot)$ and $\Phi_{i(\cdot)}$. The number of macro-cluster, denoted by $R$, is set such that $K$ is a multiple of $R$. Each macro-cluster contains $K/R$ clusters. Cluster-group encoding requires $R$ encoders, one per macro-cluster. To limit the model size, we empirically set $R=\min(15,\lfloor K/\eta\rfloor)$. For ease of implementation, $G$ is first partitioned into $R$ macro-clusters, then each macro-cluster is partitioned into $K/R$ clusters. This is simpler than, say, partitioning $G$ into $K$ (smaller) clusters and grouping them back to macro-clusters.

Next, $z_{v}$ is input to GNN whose $l$-th layer output $h^{l}_{v}$ is given by

where $h^{0}_{v}=z_{v}$. HPLC is a plug-in method and thus can be combined with different types of GNNs. The related study is provided in Sec. 7.2. Finally, we concatenate node-pair embeddings and compute scores by $y_{v,u}=\sigma(h^{L}_{v}\oplus h^{L}_{u})$ where $y_{v,u}$ denotes the link prediction score between node $v$ and $u$, and $\sigma(\cdot)$ is an MLP.

The time complexity of HPLC is incurred mainly from computing DVs and MVs. We first consider the complexity of computing $D(v)$ for all $v\in V$. There are $\eta\log N$ landmarks, and for each landmark, computing the distances from all $v\in V$ to the landmark requires $O(|E|+N\log N)$ using Fibonacci heap. Thus, the overall complexity is $O(|E|\log N+N(\log N)^{2})$. Next, we consider the complexity of computing $M(v)$. We compute Laplacian eigenvectors of landmark graph, where its time complexity is $O(K^{3})=O((\log N)^{3})$. Importantly, the computation of $D(v)$ and $M(v)$ can be done during the preprocessing stage.

The space complexity of HPLC is mainly from the GNN models for computing the node embeddings. Additional space complexity of HPLC is from the membership/distance encoder $\phi(\cdot)$ which is $O\left((F+R)H_{\textrm{in}}\right)$, and cluster-group encoding which is $O\left(H_{\textrm{in}}H_{\textrm{out}}R\right)$ respectively, where $F$ denotes the node feature dimension, $R=\min(15,\lfloor K/\eta\rfloor)$ denotes the number of macro-clusters, $H_{\textrm{in}}$ denotes the hidden dimension of $\phi(\cdot)$, and $H_{\textrm{out}}$ denotes the hidden dimension of cluster-group encoders.

Overall, the time and space complexity of HPLC is low, and our experiments shows that HPLC handles large or dense graphs well. We also measured the actual resource usage, and the related results are provided in Appendix D.

Datasets. Experiments were conducted on 7 datasets widely used for evaluating link prediction. For experiments on small graphs, we used PubMed, Cora, Citeseer, and Facebook. For experiments on dense or large graphs, we chose DDI, COLLAB, and CITATION2 provided by OGB (Hu et al., 2020). Detailed statistics and evaluation metrics of the datasets are provided in Table 9 in Appendix E.
Baseline models. We compared HPLC with Adamic Adar (AA) (Adamic and Adar, 2003), Matrix Factorization (MF) (Koren et al., 2009), Node2Vec (Grover and Leskovec, 2016), GCN (Kipf and Welling, 2017), GraphSAGE (Hamilton et al., 2017), GAT (Veličković et al., 2018), P-GNN (You et al., 2019), NBF-net (Zhu et al., 2021), and plug-in type approaches like JKNet (Xu et al., 2018), SEAL (Zhang and Chen, 2018), GCN+DE (Li et al., 2020), GCN+LPE (Dwivedi et al., 2020), GCN+LRGA (Puny et al., 2021), Graph Transformer+LPE (Dwivedi and Bresson, 2021), and PEG-DW+ (Wang et al., 2022). All methods except for AA and GAE are computed by the same decoder, which is a 2-layer MLP. For a fair comparison, we use GCN in most plug-in type approaches: SEAL, GCN+DE, GAE, JKNet, GCN+LRGA, GCN+LPE, and HPLC.
Evaluation metrics. Link prediction was evaluated based on the ranking performance of positive edges in the test data over negative ones. For COLLAB and DDI, we ranked all positive and negative edges in the test data, and computed the ratio of positive edges which are ranked in top-$k$. We did not utilize validation edges for computing node embeddings when we predicted test edges on COLLAB. In CITATION2, we computed all positive and negative edges, and calculated the reverse of the mean rank of positive edges. Due to high complexity when evaluating SEAL, we only trained 2% of training set edges and evaluated 1% of validation and test set edges respectively, as recommended in the official GitHub of SEAL. The metrics for DDI, COLLAB and CITATION2 are chosen according to the official OGB settings  (Hu et al., 2020). For Cora, Citeseer, PubMed, and Facebook, Area Under ROC Curve (AUC) is used as the metric similar to prior works (Kipf and Welling, 2017; Wang et al., 2022). If applicable, we calculated the average and harmonic mean (HM) of the measurements. HM penalizes the model for very low scores, thus is a useful indicator of robustness.
Hyperparameters. We used GCN as our base GNN encoder. MLP is used in decoders, except GAE. None of the baseline methods used edge weights. During training, negative edges are randomly sampled at a ratio of 1:1 with positive edges. The details of hyperparameters are provided in Appendix F.

Experimental results are summarized in Table 2. HPLC outperformed the baselines on most datasets. HPLC achieved large performance gains over GAE combined with GCN on all datasets, which are 88.9% on DDI, 27.0% on COLLAB, 12.7% on Citeseer, 7.1% on Cora, and 1.4% on CITATION2. HPLC showed superior performance over SEAL, achieving gains of 167% on DDI, and 12.0% on Citeseer. We compare HPLC with other distance-based methods. Compared to GCN+DE which encodes distances from a target node set whose representations are to be learned, or to P-GNN which uses distances to random anchor sets, HPLC achieved higher performance gains by a large margin. HPLC outperformed other positional encoding methods, e.g., GCN+LPE, Graph Transformer+LPE, and PEG-DW+. The results show that the proposed landmark-based representation can be effective for estimating the positional information of nodes.

SEAL and NBF-net performed poorly on DDI which is a highly dense graph. Since the nodes of DDI have a large number of neighbors, the enclosing subgraphs are both very dense and large, and the models struggle with learning the representations of local structures or paths between nodes. By contrast, HPLC achieved the best performance on DDI, demonstrating its effectiveness on densely connected graphs.

Finally, we computed the average and harmonic means of performance measurements except for the methods with OOM problems. Although the averages are taken over heterogeneous metrics and thus do not represent specific performance metrics, they are presented for comparison purposes. In summary, HPLC achieved the best average and harmonic mean of performance measurements, demonstrating both its effectiveness and robustness.

Table 3 shows the performances in the ablation analysis for the model components. The components denoted by“DV”, “CE” and “MV” columns in Table 3 indicate the usage of distance vector, cluster-group encoding and membership vector, respectively. The results show that all of the hierarchical components, i.e., “DV”, “CE” and “MV”, indeed contribute to the performance improvement.

HPLC can be combined with different GNN encoders. We experimented the combination with three types of widely-used GNN encoders. Table 4 shows that, HPLC enhances the performance of various types of GNNs.