Link Prediction Based on Graph Neural Networks
Muhan Zhang Yixin Chen



User Association for Offloading in Heterogeneous Network Based on Matern Cluster Process
    Yuxuan Xie1, Xuefei Zhang1, Qimei Cui1 and Yanyan Lu1
1Key Laboratory of Universal Wireless Communication, Ministry of Education

Beijing University of Posts and Telecommunications, Beijing, 100876, China

Email: 1{cuiqimei, zhangxuefei}@bupt.edu.cn

    
==============================================================================================================================================================================================================================================================






Link prediction is a key problem for network-structured data. Link prediction heuristics use some score functions, such as common neighbors and Katz index, to measure the likelihood of links. They have obtained wide practical uses due to their simplicity, interpretability, and for some of them, scalability. However, every heuristic has a strong assumption on when two nodes are likely to link, which limits their effectiveness on networks where these assumptions fail. In this regard, a more reasonable way should be learning a suitable heuristic from a given network instead of using predefined ones. By extracting a local subgraph around each target link, we aim to learn a function mapping the subgraph patterns to link existence, thus automatically learning a “heuristic” that suits the current network. In this paper, we study this heuristic learning paradigm for link prediction. First, we develop a novel γ-decaying heuristic theory. The theory unifies a wide range of heuristics in a single framework, and proves that all these heuristics can be well approximated from local subgraphs. Our results show that local subgraphs reserve rich information related to link existence. Second, based on the γ-decaying theory, we propose a new method to learn heuristics from local subgraphs using a graph neural network (GNN). Its experimental results show unprecedented performance, working consistently well on a wide range of problems.








§ INTRODUCTION

Link prediction is to predict whether two nodes in a network are likely to have a link <cit.>. Given the ubiquitous existence of networks, it has many applications such as friend recommendation <cit.>, movie recommendation <cit.>, knowledge graph completion <cit.>, and metabolic network reconstruction <cit.>.

One class of simple yet effective approaches for link prediction is called heuristic methods. Heuristic methods compute some heuristic node similarity scores as the likelihood of links <cit.>. Existing heuristics can be categorized based on the maximum hop of neighbors needed to calculate the score.
For example, common neighbors (CN) and preferential attachment (PA) <cit.> are first-order heuristics, since they only involve the one-hop neighbors of two target nodes. Adamic-Adar (AA) and resource allocation (RA) <cit.> are second-order heuristics, as they are calculated from up to two-hop neighborhood of the target nodes. We define h-order heuristics to be those heuristics which require knowing up to h-hop neighborhood of the target nodes. There are also some high-order heuristics which require knowing the entire network. Examples include Katz, rooted PageRank (PR) <cit.>, and SimRank (SR) <cit.>. Table <ref> in Appendix <ref> summarizes eight popular heuristics.



Although working well in practice, heuristic methods have strong assumptions on when links may exist. For example, the common neighbor heuristic assumes that two nodes are more likely to connect if they have many common neighbors. This assumption may be correct in social networks, but is shown to fail in protein-protein interaction (PPI) networks – two proteins sharing many common neighbors are actually less likely to interact <cit.>.

In fact, the heuristics belong to a more generic class, namely graph structure features. Graph structure features are those features located inside the observed node and edge structures of the network, which can be calculated directly from the graph. Since heuristics can be viewed as  predefined graph structure features, a natural idea is to automatically learn such features from the network. <cit.> first studied this problem. They extract local enclosing subgraphs around links as the training data, and use a fully-connected neural network to learn which enclosing subgraphs correspond to link existence. Their method called Weisfeiler-Lehman Neural Machine (WLNM) has achieved state-of-the-art link prediction performance. The enclosing subgraph for a node pair (x,y) is the subgraph induced from the network by the union of x and y's neighbors up to h hops. Figure <ref> illustrates the 1-hop enclosing subgraphs for (A,B) and (C,D). These enclosing subgraphs are very informative for link prediction – all first-order heuristics such as common neighbors can be directly calculated from the 1-hop enclosing subgraphs. 



However, it is shown that high-order heuristics such as rooted PageRank and Katz often have much better performance than first and second-order ones <cit.>. To effectively learn good high-order features, it seems that we need a very large hop number h so that the enclosing subgraph becomes the entire network. This results in unaffordable time and memory consumption for most practical networks. But do we really need such a large h to learn high-order heuristics? 

Fortunately, as our first contribution, we show that we do not necessarily need a very large h to learn high-order graph structure features. We dive into the inherent mechanisms of link prediction heuristics, and find that most high-order heuristics can be unified by a γ-decaying theory. We prove that, under mild conditions, any γ-decaying heuristic can be effectively approximated from an h-hop enclosing subgraph, where the approximation error decreases at least exponentially with h. This means that we can safely use even a small h to learn good high-order features. It also implies that the “effective order” of these high-order heuristics is not that high.





Based on our theoretical results, we propose a novel link prediction framework, SEAL, to learn general graph structure features from local enclosing subgraphs. SEAL fixes multiple drawbacks of WLNM. First, a graph neural network (GNN) <cit.> is used to replace the fully-connected neural network in WLNM, which enables better graph feature learning ability. Second, SEAL permits learning from not only subgraph structures, but also latent and explicit node features, thus absorbing multiple types of information. We empirically verified its much improved performance.









Our contributions are summarized as follows. 1) We present a new theory for learning link prediction heuristics, justifying learning from local subgraphs instead of entire networks. 2) We propose SEAL, a novel link prediction framework based on GNN (illustrated in Figure <ref>). SEAL outperforms all heuristic methods, latent feature methods, and recent network embedding methods by large margins. SEAL also outperforms the previous state-of-the-art method, WLNM.







































































































§ PRELIMINARIES



Notations  
Let G = (V, E) be an undirected graph, where V is the set of vertices and E ⊆ V × V is the set of observed links. Its adjacency matrix is A, where A_i,j = 1 if (i,j) ∈ E and A_i,j = 0 otherwise. 

For any nodes x, y ∈ V, let Γ(x) be the 1-hop neighbors of x, and d(x,y) be the shortest path distance between x and y. A walk w = ⟨ v_0, ⋯, v_k ⟩ is a sequence of nodes with (v_i, v_i+1) ∈ E. We use |⟨ v_0, ⋯, v_k ⟩| to denote the length of the walk w, which is k here.

Latent features and explicit features  
Besides graph structure features, latent features and explicit features are also studied for link prediction. Latent feature methods <cit.> factorize some matrix representations of the network to learn a low-dimensional latent representation/embedding for each node. Examples include matrix factorization <cit.> and stochastic block model <cit.> etc. Recently, a number of network embedding techniques have been proposed, such as DeepWalk <cit.>, LINE <cit.> and node2vec <cit.>, which are also latent feature methods since they implicitly factorize some matrices too <cit.>. 

Explicit features are often available in the form of node attributes, describing all kinds of side information about individual nodes. It is shown that combining graph structure features with latent features and explicit features can improve the performance <cit.>.



Graph neural networks   Graph neural network (GNN) is a new type of neural network for learning over graphs <cit.>). Here, we only briefly introduce the components of a GNN since this paper is not about GNN innovations but is a novel application of GNN. 
A GNN usually consists of 1) graph convolution layers which extract local substructure features for individual nodes, and 2) a graph aggregation layer which aggregates node-level features into a graph-level feature vector. Many graph convolution layers can be unified into a message passing framework <cit.>. 






Supervised heuristic learning   
There are some previous attempts to learn supervised heuristics for link prediction. The closest work to ours is the Weisfeiler-Lehman Neural Machine (WLNM) <cit.>, which also learns from local subgraphs. However, WLNM has several drawbacks. Firstly, WLNM trains a fully-connected neural network on the subgraphs' adjacency matrices. Since fully-connected neural networks only accept fixed-size tensors as input, WLNM requires truncating different subgraphs to the same size, which may lose much structural information. Secondly, due to the limitation of adjacency matrix representations, WLNM cannot learn from latent or explicit features. Thirdly, theoretical justifications are also missing. We include more discussion on WLNM in Appendix <ref>.
Another related line of research is to train a supervised learning model on different heuristics' combination. For example, the path ranking algorithm <cit.> trains logistic regression on different path types' probabilities to predict relations in knowledge graphs. <cit.> propose to incorporate heuristic features into tensor factorization models. However, these models still rely on predefined heuristics – they cannot learn general graph structure features.








§ A THEORY FOR UNIFYING LINK PREDICTION HEURISTICS
 
In this section, we aim to understand deeper the mechanisms behind various link prediction heuristics, and thus motivating the idea of learning heuristics from local subgraphs. Due to the large number of graph learning techniques, note that we are not concerned with the generalization error of a particular method, but focus on the information reserved in the subgraphs for calculating existing heuristics. 
 (Enclosing subgraph)
For a graph G = (V, E), given two nodes x, y ∈ V, the h-hop enclosing subgraph for (x,y) is the subgraph G^h_x,y induced from G by the set of nodes { i  |  d(i,x) ≤ h  or  d(i,y)≤ h  }.
























The enclosing subgraph describes the “h-hop surrounding environment" of (x,y). Since G^h_x,y contains all h-hop neighbors of x and y, we naturally have the following theorem. 


Any h-order heuristic for (x, y) can be accurately calculated from G^h_x,y.


For example, a 2-hop enclosing subgraph will contain all the information needed to calculate any first and second-order heuristics. 
However, although first and second-order heuristics are well covered by local enclosing subgraphs, an extremely large h seems to be still needed for learning high-order heuristics. Surprisingly, our following analysis shows that learning high-order heuristics is also feasible with a small h. We support this first by defining the γ-decaying heuristic. We will show that under certain conditions, a γ-decaying heuristic can be very well approximated from the h-hop enclosing subgraph. Moreover, we will show that almost all well-known high-order heuristics can be unified into this γ-decaying heuristic framework.




 (γ-decaying heuristic)
A γ-decaying heuristic for (x,y) has the following form:

    ℋ(x,y) = η∑_l=1^∞γ^l f(x,y,l),

where γ is a decaying factor between 0 and 1, η is a positive constant or a positive function of γ that is upper bounded by a constant, f is a nonnegative function of x,y,l under the the given network.


Next, we will show that under certain conditions, a γ-decaying heuristic can be approximated from an h-hop enclosing subgraph, and the approximation error decreases at least exponentially with h.


Given a γ-decaying heuristic ℋ(x,y) = η∑_l=1^∞γ^l f(x,y,l), if f(x,y,l) satisfies:




  * (property 1) f(x,y,l) ≤λ^l where λ < 1/γ; and



  * (property 2) f(x,y,l) is calculable from G^h_x,y for l=1,2,⋯, g(h), where g(h)=ah + b with a,b ∈ℕ and a>0, 



then ℋ(x,y) can be approximated from G^h_x,y and the approximation error decreases at least exponentially with h.



We can approximate such a γ-decaying heuristic by summing over its first g(h) terms.

    ℋ(x,y) := η∑_l=1^g(h)γ^l f(x,y,l).

The approximation error can be bounded as follows.

    |ℋ(x,y) - ℋ(x,y) | = η∑_l=g(h)+1^∞γ^l f(x,y,l) ≤η∑_l=ah+b+1^∞γ^l λ^l = η (γλ)^ah+b+1 (1-γλ)^-1.


In practice, a small γλ and a large a lead to a faster decreasing speed. Next we will prove that three popular high-order heuristics: Katz, rooted PageRank and SimRank, are all γ-decaying heuristics which satisfy the properties in Theorem <ref>. First, we need the following lemma.

Any walk between x and y with length l ≤ 2h+1 is included in G^h_x,y.



Given any walk w = ⟨ x, v_1, ⋯, v_l-1, y ⟩ with length l, we will show that every node v_i is included in G^h_x,y. Consider any v_i. Assume d(v_i,x) ≥ h+1 and d(v_i,y) ≥ h+1. Then, 2h+1 ≥ l = |⟨ x, v_1,⋯,v_i⟩| + |⟨ v_i,⋯,v_l-1,y⟩| ≥ d(v_i,x) + d(v_i,y) ≥ 2h+2, a contradiction. Thus, d(v_i,x) ≤ h or d(v_i,y) ≤ h. By the definition of G^h_x,y, v_i must be included in G^h_x,y. 


Next we will analyze Katz, rooted PageRank and SimRank one by one.




 §.§ Katz index


The Katz index <cit.> for (x, y) is defined as

    Katz_x,y = ∑_l = 1^∞β^l  |walks^⟨ l⟩(x,y)| = ∑_l=1^∞β^l [A^l]_x,y,

where walks^⟨ l⟩(x,y) is the set of length-l walks between x and y, and A^l is the l^th power of the adjacency matrix of the network. Katz index sums over the collection of all walks between x and y where a walk of length l is damped by β^l (0 < β < 1), giving more weight to shorter walks.

Katz index is directly defined in the form of a γ-decaying heuristic with η=1,γ=β, and f(x,y,l) = |walks^⟨ l⟩(x,y)|. According to Lemma <ref>, |walks^⟨ l⟩(x,y)| is calculable from G^h_x,y for l≤ 2h+1, thus property 2 in Theorem <ref> is satisfied. Now we show when property 1 is satisfied.


For any nodes i,j, [A^l]_i,j is bounded by d^l, where d is the maximum node degree of the network.



We prove it by induction. When l=1, A_i,j≤ d for any (i,j). Thus the base case is correct. Now, assume by induction  that [A^l]_i,j≤ d^l for any (i,j), we have



    [A^l+1]_i,j = ∑_k=1^|V| [A^l]_i,kA_k,j≤ d^l∑_k=1^|V| A_k,j≤ d^l d 
    = d^l+1.




Taking λ = d, we can see that whenever d < 1/β, the Katz index will satisfy property 1 in Theorem <ref>. In practice, the damping factor β is often set to very small values like 5E-4 <cit.>, which implies that Katz can be very well approximated from the h-hop enclosing subgraph.








 §.§ PageRank


The rooted PageRank for node x calculates the stationary distribution of a random walker starting at x, who iteratively moves to a random neighbor of its current position with probability α or returns to x with probability 1-α. Let π_x denote the stationary distribution vector. Let [π_x]_i denote the probability that the random walker is at node i under the stationary distribution. 

Let P be the transition matrix with P_i,j = 1/|Γ(v_j)| if (i,j) ∈ E and P_i,j = 0 otherwise. Let 𝐞_x be a vector with the x^th element being 1 and others being 0. The stationary distribution satisfies


    π_x    = α Pπ_x + (1-α)𝐞_x.


When used for link prediction, the score for (x,y) is given by [π_x]_y (or [π_x]_y + [π_y]_x for symmetry). To show that rooted PageRank is a γ-decaying heuristic, we introduce the inverse P-distance theory <cit.>, which states that [π_x]_y can be equivalently written as follows:

    [π_x]_y = (1-α)∑_w: x ⇝ y P[w]  α^len(w),

where the summation is taken over all walks w starting at x and ending at y (possibly touching x and y multiple times). For a walk w = ⟨ v_0, v_1, ⋯, v_k ⟩, len(w) := |⟨ v_0, v_1, ⋯, v_k ⟩| is the length of the walk. The term P[w] is defined as ∏_i=0^k-11/|Γ(v_i)|, which can be interpreted as the probability of traveling w. Now we have the following theorem.


The rooted PageRank heuristic is a γ-decaying heuristic which satisfies the properties in Theorem <ref>.


We first write [π_x]_y in the following form.

    [π_x]_y = (1-α)∑_l=1^∞∑_w: x ⇝ y 
    len(w) = l  P[w]  α^l.

Defining f(x,y,l) := ∑_w: x ⇝ y 
len(w) = l  P[w] leads to the form of a γ-decaying heuristic. Note that f(x,y,l) is the probability that a random walker starting at x stops at y with exactly l steps, which satisfies ∑_z∈ Vf(x,z,l) = 1. Thus, f(x,y,l) ≤ 1 < 1/α (property 1). According to Lemma <ref>, f(x,y,l) is also calculable from G^h_x,y for l≤ 2h+1 (property 2).





 §.§ SimRank


The SimRank score <cit.> is motivated by the intuition that two nodes are similar if their neighbors are also similar. It is defined in the following recursive way: if x=y, then s(x,y) := 1; otherwise,

    s(x,y) := γ∑_a∈Γ(x)∑_b∈Γ(y) s(a,b)/|Γ(x)| · |Γ(y)|

where γ is a constant between 0 and 1. According to <cit.>, SimRank has an equivalent definition:

    s(x,y) = ∑_w: (x,y) ⊸ (z,z)P[w]  γ^len(w),

where w: (x,y) ⊸ (z,z) denotes all simultaneous walks such that one walk starts at x, the other walk starts at y, and they first meet at any vertex z. For a simultaneous walk w = ⟨ (v_0, u_0), ⋯, (v_k, u_k) ⟩, len(w) = k is the length of the walk. The term P[w] is similarly defined as ∏_i=0^k-11/|Γ(v_i)||Γ(u_i)|, describing the probability of this walk. Now we have the following theorem.


SimRank is a γ-decaying heuristic which satisfies the properties in Theorem <ref>.


We write s(x,y) as follows.

    s(x,y) = ∑_l=1^∞∑_w: (x,y) ⊸ (z,z)
    len(w) = l P[w]  γ^l,

Defining f(x,y,l) := ∑_w: (x,y) ⊸ (z,z)
len(w) = l P[w] reveals that SimRank is a γ-decaying heuristic. Note that f(x,y,l) ≤ 1 < 1/γ. It is easy to see that f(x,y,l) is also calculable from G^h_x,y for l≤ h.


Discussion  
There exist several other high-order heuristics based on path counting or random walk <cit.> which can be as well incorporated into the γ-decaying heuristic framework. We omit the analysis here. 
Our results reveal that most high-order heuristics inherently share the same γ-decaying heuristic form, and thus can be effectively approximated from an h-hop enclosing subgraph with exponentially smaller approximation error. We believe the ubiquity of γ-decaying heuristics is not by accident – it implies that a successful link prediction heuristic is better to put exponentially smaller weight on structures far away from the target, as remote parts of the network intuitively make little contribution to link existence.
Our results build the foundation for learning heuristics from local subgraphs, as they imply that local enclosing subgraphs already contain enough information to learn good graph structure features for link prediction which is much desired considering learning from the entire network is often infeasible. To summarize, from the small enclosing subgraphs extracted around links, we are able to accurately calculate first and second-order heuristics, and approximate a wide range of high-order heuristics with small errors. Therefore, given adequate feature learning ability of the model used, learning from such enclosing subgraphs is expected to achieve performance at least as good as a wide range of heuristics. There is some related work which empirically verifies that local methods can often estimate PageRank and SimRank well <cit.>. Another related theoretical work <cit.> establishes a condition of h to achieve some fixed approximation error for ordinary PageRank.



§ SEAL: AN IMPLEMETATION OF THE THEORY USING GNN










In this section, we describe our SEAL framework for link prediction. SEAL does not restrict the learned features to be in some particular forms such as γ-decaying heuristics, but instead learns general graph structure features for link prediction.
It contains three steps: 1) enclosing subgraph extraction, 2) node information matrix construction, and 3) GNN learning.
Given a network, we aim to learn automatically a “heuristic” that best explains the link formations. Motivated by the theoretical results, this function takes local enclosing subgraphs around links as input, and output how likely the links exist. To learn such a function, we train a graph neural network (GNN) over the enclosing subgraphs. Thus, the first step in SEAL is to extract enclosing subgraphs for a set of sampled positive links (observed) and a set of sampled negative links (unobserved) to construct the training data. 

A GNN typically takes (A,X) as input, where A (with slight abuse of notation) is the adjacency matrix of the input enclosing subgraph, X is the node information matrix each row of which corresponds to a node's feature vector. The second step in SEAL is to construct the node information matrix X for each enclosing subgraph. This step is crucial for training a successful GNN link prediction model. In the following, we discuss this key step. The node information matrix X in SEAL has three components: structural node labels, node embeddings and node attributes.





















 §.§ Node labeling

The first component in X is each node's structural label. A node labeling is function f_l: V →ℕ which assigns an integer label f_l(i) to every node i in the enclosing subgraph. The purpose is to use different labels to mark nodes' different roles in an enclosing subgraph:
1) The center nodes x and y are the target nodes between which the link is located. 2) Nodes with different relative positions to the center have different structural importance to the link. 
A proper node labeling should mark such differences. If we do not mark such differences, GNNs will not be able to tell where are the target nodes between which a link existence should be predicted, and lose structural information. 

Our node labeling method is derived from the following criteria: 1) The two target nodes x and y always have the distinctive label “1”. 2) Nodes i and j have the same label if d(i,x) = d(j,x) and d(i,y) = d(j,y). The second criterion is because, intuitively, a node i's topological position within an enclosing subgraph can be described by its radius with respect to the two center nodes, namely (d(i,x), d(i,y)). Thus, we let nodes on the same orbit have the same label, so that the node labels can reflect nodes' relative positions and structural importance within subgraphs. 


Based on the above criteria, we propose a Double-Radius Node Labeling (DRNL) as follows. First, assign label 1 to x and y. Then, for any node i with (d(i,x), d(i,y)) = (1,1), assign label f_l(i)=2. Nodes with radius (1,2) or (2,1) get label 3. Nodes with radius (1,3) or (3,1) get 4. Nodes with (2,2) get 5. Nodes with (1,4) or (4,1) get 6. Nodes with (2,3) or (3,2) get 7. So on and so forth. In other words, we iteratively assign larger labels to nodes with a larger radius w.r.t. both center nodes, where the label f_l(i) and the double-radius (d(i,x),d(i,y)) satisfy

1) if d(i,x) + d(i,y) ≠ d(j,x) + d(j,y), then d(i,x) + d(i,y) < d(j,x) + d(j,y) ⇔  f_l(i) < f_l(j);

2) if d(i,x) + d(i,y) = d(j,x) + d(j,y), then d(i,x)d(i,y) < d(j,x)d(j,y) ⇔  f_l(i) < f_l(j). 
 
One advantage of DRNL is that it has a perfect hashing function

    f_l(i) = 1 + min(d_x, d_y) + (d / 2)[(d / 2) + (d % 2) - 1],

where d_x := d(i,x), d_y := d(i,y), d := d_x + d_y, (d / 2) and (d % 2) are the integer quotient and remainder of d divided by 2, respectively. This perfect hashing allows fast closed-form computations.

For nodes with d(i,x)=∞ or d(i,y)=∞, we give them a null label 0. Note that DRNL is not the only possible way of node labeling, but we empirically verified its better performance than no labeling and other naive labelings. We discuss more about node labeling in Appendix <ref>. After getting the labels, we use their one-hot encoding vectors to construct X.






 §.§ Incorporating latent and explicit features

Other than the structural node labels, the node information matrix X also provides an opportunity to include latent and explicit features. By concatenating each node's embedding/attribute vector to its corresponding row in X, we can make SEAL simultaneously learn from all three types of features.

Generating the node embeddings for SEAL is nontrivial. Suppose we are given the observed network G=(V,E), a set of sampled positive training links E_p ⊆ E, and a set of sampled negative training links E_n with E_n ∩ E = ∅. If we directly generate node embeddings on G, the node embeddings will record the link existence information of the training links (since E_p ⊆ E). We observed that GNNs can quickly find out such link existence information and optimize by only fitting this part of information. This results in bad generalization performance in our experiments.
Our trick is to temporally add E_n into E, and generate the embeddings on G' = (V, E ∪ E_n). This way, the positive and negative training links will have the same link existence information recorded in the embeddings, so that GNN cannot classify links by only fitting this part of information. 
We empirically verified the much improved performance of this trick to SEAL. We name this trick negative injection. 

We name our proposed framework SEAL (learning from Subgraphs, Embeddings and Attributes for Link prediction), emphasizing its ability to jointly learn from three types of features.







§ EXPERIMENTAL RESULTS

We conduct extensive experiments to evaluate SEAL. Our results show that SEAL is a superb and robust framework for link prediction, achieving unprecedentedly strong performance on various networks. We use AUC and average precision (AP) as evaluation metrics. We run all experiments for 10 times and report the average AUC results and standard deviations. We leave the the AP and time results in Appendix <ref>.
SEAL is flexible with what GNN or node embeddings to use. Thus, we choose a recent architecture DGCNN <cit.> as the default GNN, and node2vec <cit.> as the default embeddings. 
The code and data are available at https://github.com/muhanzhang/SEAL.

Datasets   The eight datasets used are: USAir, NS, PB, Yeast, C.ele, Power, Router, and E.coli (please see Appendix <ref> for details). We randomly remove 10% existing links from each dataset as positive testing data. Following a standard manner of learning-based link prediction, we randomly sample the same number of nonexistent links (unconnected node pairs) as negative testing data. We use the remaining 90% existing links as well as the same number of additionally sampled nonexistent links to construct the training data.










Comparison to heuristic methods   We first compare SEAL with methods that only use graph structure features.  We include eight popular heuristics (shown in Appendix <ref>, Table <ref>): common neighbors (CN), Jaccard, preferential attachment (PA), Adamic-Adar (AA), resource allocation (RA), Katz, PageRank (PR), and SimRank (SR). We additionally include Ensemble (ENS) which trains a logistic regression classifier on the eight heuristic scores. We also include two heuristic learning methods: Weisfeiler-Lehman graph kernel (WLK) <cit.> and WLNM <cit.>, which also learn from (truncated) enclosing subgraphs. We omit path ranking methods <cit.> as well as other recent methods which are specifically designed for knowledge graphs or recommender systems <cit.>. As all the baselines only use graph structure features, we restrict SEAL to not include any latent or explicit features. In SEAL, the hop number h is an important hyperparameter. Here, we select h only from {1,2}, since on one hand we empirically verified that the performance typically does not increase after h≥ 3, which validates our theoretical results that the most useful information is within local structures. On the other hand, even h=3 sometimes results in very large subgraphs if a hub node is included. This raises the idea of sampling nodes in subgraphs, which we leave to future work.
The selection principle is very simple: If the second-order heuristic AA outperforms the first-order heuristic CN on 10% validation data, then we choose h=2; otherwise we choose h=1. For datasets PB and E.coli, we consistently use h=1 to fit into the memory. We include more details about the baselines and hyperparameters in Appendix <ref>.

F>p2.5em
G>p2.8em
H>p5.2em












Table <ref> shows the results. Firstly, we observe that methods which learn from enclosing subgraphs (WLK, WLNM and SEAL) generally perform much better than predefined heuristics. This indicates that the learned “heuristics” are better at capturing the network properties than manually designed ones. Among learning-based methods, SEAL has the best performance, demonstrating GNN's superior graph feature learning ability over graph kernels and fully-connected neural networks. From the results on Power and Router, we can see that although existing heuristics perform similarly to random guess, learning-based methods still maintain high performance. This suggests that we can even discover new “heuristics” for networks where no existing heuristics work.



L0.61

Comparison with latent feature methods (AUC).
0.602!

Data       MF        SBM       N2V       LINE      SPC       VGAE      SEAL  
 USAir      94.08±0.80     94.85±1.14     91.44±1.78     81.47±10.71     74.22±3.11     89.28±1.99     97.09±0.70 
 
NS         74.55±4.34     92.30±2.26     91.52±1.28     80.63±1.90      89.94±2.39     94.04±1.64     97.71±0.93 

PB         94.30±0.53     93.90±0.42     85.79±0.78     76.95±2.76      83.96±0.86     90.70±0.53     95.01±0.34 

Yeast      90.28±0.69     91.41±0.60     93.67±0.46     87.45±3.33      93.25±0.40     93.88±0.21     97.20±0.64 

C.ele      85.90±1.74     86.48±2.60     84.11±1.27     69.21±3.14      51.90±2.57     81.80±2.18     89.54±2.04 

Power      50.63±1.10     66.57±2.05     76.22±0.92     55.63±1.47      91.78±0.61     71.20±1.65     84.18±1.82 

Router     78.03±1.63     85.65±1.93     65.46±0.86     67.15±2.10      68.79±2.42     61.51±1.22     95.68±1.22 

E.coli     93.76±0.56     93.82±0.41     90.82±1.49     82.38±2.19      94.92±0.32     90.81±0.63     97.22±0.28 
 





Comparison to latent feature methods   Next we compare SEAL with six state-of-the-art latent feature methods: matrix factorization (MF), stochastic block model (SBM) <cit.>, node2vec (N2V) <cit.>, LINE <cit.>, spectral clustering (SPC), and variational graph auto-encoder (VGAE) <cit.>. Among them, VGAE uses a GNN too. Please note the difference between VGAE and SEAL: VGAE uses a node-level GNN to learn node embeddings that best reconstruct the network, while SEAL uses a graph-level GNN to classify enclosing subgraphs. Therefore, VGAE still belongs to latent feature methods. For SEAL, we additionally include the 128-dimensional node2vec embeddings in the node information matrix X. Since the datasets do not have node attributes, explicit features are not included. 

Table <ref> shows the results. As we can see, SEAL shows significant improvement over latent feature methods. One reason is that SEAL learns from both graph structures and latent features simultaneously, thus augmenting those methods that only use latent features.

We observe that SEAL with node2vec embeddings outperforms pure node2vec by large margins. This implies that network embeddings alone may not be able to capture the most useful link prediction information located in the local structures. 
It is also interesting that compared to SEAL without node2vec embeddings (Table <ref>), joint learning does not always improve the performance.
More experiments and discussion are included in Appendix <ref>.























§ CONCLUSIONS

Learning link prediction heuristics automatically is a new field. In this paper, we presented theoretical justifications for learning from local enclosing subgraphs. In particular, we proposed a γ-decaying theory to unify a wide range of high-order heuristics and prove their approximability from local subgraphs. Motivated by the theory, we proposed a novel link prediction framework, SEAL, to simultaneously learn from local enclosing subgraphs, embeddings and attributes based on graph neural networks. Experimentally we showed that SEAL achieved unprecedentedly strong performance by comparing to various heuristics, latent feature methods, and network embedding algorithms. We hope SEAL can not only inspire link prediction research, but also open up new directions for other relational machine learning problems such as knowledge graph completion and recommender systems.








  §.§.§ Acknowledgments

The work is supported in part by the III-1526012 and SCH-1622678 grants from the National Science Foundation and grant 1R21HS024581 from the National Institute of Health. 
















§ MORE ABOUT THE THREE TYPES OF FEATURES FOR LINK PREDICTION

In this section, we discuss more about the difference among the three types commonly used features for link prediction: graph structure features, latent features, and explicit features. 





Graph structure features locate inside the observed node and edge structures of the network, which can be directly observed and computed. Link prediction heuristics belong to graph structure features. We show eight popular heuristics in Table <ref>. In addition to link prediction heuristics, node centrality scores (degree, closeness, betweenness, PageRank, eigenvector, hubs etc.), graphlets, network motifs etc. all belong to graph structure features. Although effective in many domains, these predefined graph structure features are handcrafted – they only capture a small set of structure patterns, lacking the ability to express general structure patterns underlying different networks. Considering deep neural networks' success in feature learning, a natural question to ask is whether we can automatically learn such features, no longer relying on predefined ones.

Graph structure features are inductive, meaning that these features are not associated with a particular node or network. For example, the common neighbor heuristic between any pair of nodes x and y is consistently calculated by counting the number of their common one-hop neighbors, invariant to where x and y are located. Thus, graph structure features are transferrable to new nodes and new networks. This is in contrast to latent features, which are often transductive – the changing of network structure will require a complete retraining to get the latent features again.


L[1]>
m#1


Latent features are latent properties or representations of nodes, often obtained by factorizing a specific matrix derived from a network, such as the adjacency matrix or the Laplacian matrix. Through factorization, a low-dimensional embedding is learned for each node. Latent features focus more on global properties and long range effects, because the network's matrix is treated as a whole during factorization. Latent features cannot capture structural similarities between nodes <cit.>, and usually need an extremely large dimension to express some simple heuristics <cit.>. 
Latent features are also transductive. They cannot be transferred to new nodes or new networks. They are also less interpretable than graph structure features.

Network embedding methods <cit.> have gained great popularity recently. They learn low-dimensional representations for nodes too.
Recently, it is shown  that network embedding methods (including DeepWalk <cit.>, LINE <cit.>, and node2vec <cit.>) implicitly factorize some matrix representation of a network <cit.>. For example, DeepWalk approximately factorizes log (vol(G)(1/T∑_r=1^T(D^-1A)^r)D^-1) - log(b), where A is the adjacency matrix of the network G, D is the diagonal degree matrix, T is skip-gram's window size, and b is the number of negative samples. For LINE and node2vec, there also exist such matrices. Since network embedding methods also factorize matrix representations of networks, we may regard them as learning more expressive latent features through factorizing some more informative matrices.

Explicit features are often given by continuous or discrete node attribute vectors. In principle, any side information about the network other than its structure can be seen as explicit features. For example, in citation networks, word distributions are explicit features of document nodes. In social networks, a user's profile information is also explicit feature (however, their friendship information belongs to graph structure features).

These three types of features are largely orthogonal to each other. Many papers have considered using them together for link prediction <cit.> to improve the performance of single-feature-based methods.




















§ MORE DISCUSSION ABOUT NODE LABELING

The necessity of structural node labels for enclosing subgraphs is because, unlike ordinary graphs, enclosing subgraphs intrinsically have a directionality. The center of an enclosing subgraph are two nodes x and y between which the target link is located. Outward from the center, other nodes have larger and larger distance to x and y. Node labeling is to mark such structural differences thus providing additional structural information to facilitate GNN training.

When designing a node labeling for enclosing subgraphs, we always want to ensure that the target nodes x and y have a distinct label so that GNN can distinguish the target link to predict from other edges. Secondly, we want the node labels to reflect nodes' relative positions in their enclosing subgraph. This relative position can be intuitively described by a node i's double-radius with respect to x and y, i.e., (d(i,x),d(i,y)).

We restate our Double-Radius Node Labeling (DRNL) algorithm here. First, assign label 1 to x and y. Then, for any node i with (d(i,x), d(i,y)) = (1,1), assign label f_l(i)=2. Nodes with double-radius (1,2) or (2,1) get label 3. Nodes with double-radius (1,3) or (3,1) get 4. Nodes with (2,2) get 5. Nodes with (1,4) or (4,1) get 6. Nodes with (2,3) or (3,2) get 7. So on and so forth. Our DRNL not only satisfies the above criteria, but also attains the additional benefits that for nodes i and j:

1) if d(i,x) + d(i,y) ≠ d(j,x) + d(j,y), then d(i,x) + d(i,y) < d(j,x) + d(j,y) ⇔  f_l(i) < f_l(j);

2) if d(i,x) + d(i,y) = d(j,x) + d(j,y), then d(i,x)d(i,y) < d(j,x)d(j,y) ⇔  f_l(i) < f_l(j).  

That is, the magnitude of node labels also reflects their distance to the center. Nodes with smaller arithmetic mean distance to the target nodes get smaller labels. If two nodes have the same arithmetic mean distance, the node with a smaller geometric mean distance to the target nodes gets a smaller label. Note that these additional benefits will not be available under one-hot encoding of node labels, since the magnitude information will be lost after one-hot encoding. However, such a labeling is potentially useful when node labels are directly used for training, or used to rank the nodes. Furthermore, our node labeling has a perfect hashing (<ref>) which allows closed-form computation.

We present a lookup table for DRNL and an example labeled subgraph in Figure <ref>. Note that when calculating d(i,x), we temporally remove y from the subgraph, and vice versa. This is because we aim to use the pure distance between i and x without the influence of y. If we do not remove y, d(i,x) will be upper bounded by d(i,y) + d(x,y), obscuring the “true distance” between i and x. 



Our node labeling algorithm is different from the Weisfeiler-Lehman algorithm used in WLNM <cit.>. In WLNM, node labeling is for defining a node order in adjacency matrices – the labels are not really input to machine learning models. To rank nodes with least ties, the node labels should be as fine as possible in WLNM. In comparison, the node labels in SEAL need not be very fine, as their purpose is for indicating nodes' different roles within the enclosing subgraph, not for ranking nodes. 
In addition, node labels in SEAL are encoded into node information matrices and input to machine learning models.







§ DATASET DETAILS


USAir <cit.> is a network of US Air lines with 332 nodes and 2,126 edges. The average node degree is 12.81. NS <cit.> is a collaboration network of researchers in network science with 1,589 nodes and 2,742 edges. The average node degree is 3.45. PB <cit.> is a network of US political blogs with 1,222 nodes and 16,714 edges. The average node degree is 27.36. Yeast <cit.> is a protein-protein interaction network in yeast with 2,375 nodes and 11,693 edges. The average node degree is 9.85. C.ele <cit.> is a neural network of C. elegans with 297 nodes and 2,148 edges. The average node degree is 14.46. Power <cit.> is an electrical grid of western US with 4,941 nodes and 6,594 edges. The average node degree is 2.67. Router <cit.> is a router-level Internet with 5,022 nodes and 6,258 edges. The average node degree is 2.49. E.coli <cit.> is a pairwise reaction network of metabolites in E. coli with 1,805 nodes and 14,660 edges. The average node degree is 12.55. 




§ ADDITIONAL DETAILS ABOUT BASELINES
 
Hyperparameters of heuristic and latent feature methods  
Most hyperparameters are inherited from the original paper of each method. For Katz, we set the damping factor β to 0.001. For PageRank, we set the damping factor α to 0.85. For SimRank, we set γ to 0.8.
For stochastic block model (SBM), we use the implementation of <cit.> using a latent group number 12. For matrix factorization (MF), we use the libFM <cit.> software with the default parameters. 
For node2vec, LINE, and spectral clustering, we first generate 128-dimensional embeddings from the observed networks with default parameters of each software. Then, we use the Hadamard product of two nodes' embeddings as a link's embedding as suggested in <cit.>, and train a logistic regression model with Liblinear <cit.> using automatic hyperparameter selection. For VGAE, we use its default setting.

WLNM   Weisfeiler-Lehman Neural Machine (WLNM) <cit.> is a recent link prediction method that learns general graph structure features. It achieves state-of-the-art performance on various networks, outperforming all handcrafted heuristics. WLNM has three steps: enclosing subgraph extraction, subgraph pattern encoding, and neural network training. In the enclosing subgraph extraction step: for each node pair (x,y), WLNM iteratively extracts x and y's one-hop neighbors, two-hop neighbors, and so on, until the enclosing subgraph has more than K vertices, where K is a user-defined integer. In the subgraph pattern encoding step, WLNM uses the Weisfeiler-Lehman algorithm to define an order for nodes within each enclosing subgraph, so that the neural network can read different subgraphs' nodes in a consistent order and learn meaningful patterns. To unify the sizes of the enclosing subgraphs, after getting the vertex order, the last few vertices are deleted so that all the truncated enclosing subgraphs have the same size K. These truncated enclosing subgraphs are reordered and their fixed-size adjacency matrices are fed into the fully-connected neural network to train a link prediction model. Due to the truncation, WLNM cannot consistently learn from each link's full h-hop neighborhood. The loss of structural information limits WLNM's performance and restrict it from learning complete h-order graph structure features. Following <cit.>, we use K=10 (the best performing K) in our experiments.




WLK   Weisfeiler-Lehman graph kernel (WLK) <cit.> is a state-of-the-art graph kernel. Graph kernels make kernel machines feasible for graph classification by defining some positive semidefinite graph similarity scores. Most graph kernels measure graph similarity by decomposing graphs into small substructures and adding up the pair-wise similarities between these components. Common types of substructures include walks <cit.>, subgraphs <cit.>, paths <cit.>, and subtrees <cit.>. WLK is based on counting common rooted subtrees between two graphs. In our experiments, we train a SVM on the WL kernel matrix. We feed the same enclosing subgraphs as in SEAL to WLK. We search the subtree depth from {0,1,2,3,4,5} on 10% validation links. WLK does not support continuous node information, but supports integer node labels. Thus, we feed the same structural node labels from (<ref>) to WLK too. 

We compare the characteristics of different link prediction methods in Table <ref>.




§ CONFIGURATION DETAILS OF SEAL


In the experiments, we use Deep Graph Convolutional Neural Network (DGCNN) <cit.> as the default GNN engine of SEAL. DGCNN is a recent GNN architecture for graph classification. It has consistently good performance on various benchmark datasets with a single network architecture (avoid hyperparameter tweaking).
DGCNN is equipped with propagation-based graph convolution layers and a novel graph aggregation layer, called SortPooling. We illustrate the overall architecture of DGCNN in Figure <ref>. Given the adjacency matrix A ∈{1,0}^n× n and the node information matrix X ∈ℝ^n× c of an enclosing subgraph, DGCNN uses the following graph convolution layer:

    Z = f( D̃^-1Ã XW),

where Ã = A + I, D̃ is a diagonal degree matrix with D̃_i,i = ∑_j Ã_i,j, W ∈ℝ^c× c' is a matrix of trainable graph convolution parameters, f is an element-wise nonlinear activation function, and Z ∈ℝ^n× c' are the new node states. The mechanism behind (<ref>) is that the initial node states X are first applied a linear transformation by multiplying W, and then propagated to neighboring nodes through the propagation matrix D̃^-1Ã. After graph convolution, the i^th row of Z becomes:

    Z_i = f( 1/|Γ(i)|+1 [X_i W + ∑_j∈Γ(i) X_j W] ),

which summarizes the node information as well as the first-order structure pattern from i's neighbors. DGCNN stacks multiple graph convolution layers (<ref>) and concatenates each layer's node states as the final node states, in order to extract multi-hop node features.




A graph aggregation layer constructs a graph-level feature vector from individual nodes' final states, which is used for graph classification. The most widely used aggregation operation is summing, i.e., nodes' final states after graph convolutions are summed up as the graph's representation. However, the averaging effect of summing might lose much individual nodes' information as well as the topological information of the graph. DGCNN uses a novel SortPooling layer, which sorts the final node states according to the last graph convolution layer's output to achieve an isomorphism invariant node ordering <cit.>. A max-k pooling operation is then used to unify the sizes of the sorted representations of different graphs, which enables training a traditional 1-D CNN on the node sequence.

We use the default setting of DGCNN, i.e., four graph convolution layers as in (<ref>) with 32,32,32,1 channels, a SortPooling layer (with k such that 60% graphs have nodes less than k), two 1-D convolution layers (16 and 32 output channels) and a dense layer (128 neurons), see <cit.>. We train DGCNN on enclosing subgraphs for 50 epochs, and select the model with the smallest loss on the 10% validation data to predict the testing links.

Note that, in any positive training link's enclosing subgraph, we should always remove the edge between the two target nodes before feeding it into a graph classification model. This is because this edge will contain the link existence information, which is not available in any testing link's enclosing subgraph.










§ ADDITIONAL RESULTS


In this section, we show the additional experimental results. We first use 90% observed links as training links and 10% as testing links following the main paper's experiments. The average precision (AP) comparison results with heuristic methods are shown in Table <ref>. The AP comparison results with latent feature methods are shown in Table <ref>. We can see that our proposed SEAL shows great performance improvement over all baselines in both AUC and AP.










To evaluate SEAL's scalability, we show its single-GPU inference time performance in Table <ref>. As we can see, SEAL has good scalability. For networks with over 1E7 potential links, SEAL took less than an hour to make all the predictions. One possible way to further scale SEAL to social networks with millions of users is to first use some simple heuristics such as common neighbors to filter out most unlikely links and then use SEAL to make further recommendations. Another way is to restrict the candidate friend recommendations to be those who are at most 2 or 3 hops away from the target user, which will vastly reduce the number of candidate links to infer for each user and thus further increase the scalability.























Next, we redo the comparisons under 50%–50% train/test split. We randomly remove 50% existing links as positive testing links and use the remaining 50% existing links as positive training links. The same number of negative training and testing links are sampled from the nonexistent links as well. The AUC results are shown in Table <ref> and <ref>. The AP results are shown in Table <ref> and <ref>.

The results are consistent with the 90%–10% split setting.
As we can see, SEAL is still the best among all methods in general. The performance gains over heuristic methods are even larger compared to the 90%-10% split. This indicates that SEAL is able to learn good heuristics even when the network is very incomplete. SEAL also shows more clear advantages over WLNM. On the other hand, we observe that VGAE becomes a strong baseline when network is sparser by achieving the best AUC results on 3 out of 8 datasets. It is thus interesting to study whether replacing the node2vec embeddings in SEAL with the VGAE embeddings can further improve the performance. We leave it to future work.






We further conduct experiments with the setting of the node2vec paper <cit.> on five networks: arXiv (18,722 nodes and 198,110 edges) <cit.>, Facebook (4,039 nodes and 88,234 edges) <cit.>, BlogCatalog (10,312 nodes, 333,983 edges and 39 attributes) <cit.>, Wikipedia (4,777 nodes, 184,812 edges and 40 attributes) <cit.>, and Protein-Protein Interactions (PPI) (3,890 nodes, 76,584 edges and 50 attributes) <cit.>. For each network, 50% of random links are removed and used as testing data, while keeping the remaining network connected. For Facebook and arXiv, all remained links are used as positive training data. For PPI, BlogCatalog and Wikipedia, we sample 10,000 remained links as positive training data. We compare SEAL (h=1, 10 training epochs) with node2vec, LINE, SPC, VGAE, and WLNM (K=10). For node2vec, we use the parameters provided in <cit.> if available. For SEAL and VGAE, the node attributes are used since only these two methods support explicit features. 

Table <ref> shows the results. As we can see, SEAL consistently outperforms all embedding methods. Especially on the last three networks, SEAL (with node2vec embeddings) outperforms pure node2vec by large margins. These results indicate that in many cases, embedding methods alone cannot capture the most useful link prediction information, while effectively combining the power of different types of features results in much better performance. SEAL also consistently outperforms WLNM.









