Generative Graph Neural Networks for Link Prediction

Owing to their potential in modeling the complex structures of non-Euclidean graphs, GNNs have achieved state-of-the-art performance on almost all graph-based tasks, such as node classification, graph classification, link prediction. Based on different theories and perspectives, a plethora of different GNNs have been proposed over the years. Generally, GNNs can be divided into two categories: spectral-based and spatial-based methods. Of these, spectral-based GNNs are types of GNNs that design graph convolution operators in the spectral domain using Fourier transform. The involved convolution operation is defined as follows:

where $\odot$ denotes an element-wise product. The spectral filter is defined as ${\bf{g}}={{\bf{U}}^{\rm{T}}}{f_{1}}$, and the node signal ${\bf{X}}$ can be processed as follows:

where $\bf{U}$ denotes a matrix of eigenvectors of the normalized Laplacian graph ${\bf{L}}={\bf{I}}-{{\bf{D}}^{-\frac{1}{2}}}{\bf{A}}{{\bf{D}}^{-\frac{1}{2}}}={\bf{U}}\Lambda{{\bf{U}}^{\rm{T}}}$ [11]. Assuming that feature representation of node should be affected only by its $k$-hop neighborhood, [16] proposed a Chebyshev polynomial based $k$-localized convolution and developed a convolutional neural network, ChebNet, which eliminated the need to compute the eigenvectors of the Laplacian. Subsequently, [50] simplified the Chebshev polynomial filter using its first-order approximation and proposed the popular spectral-based method called Graph Convolutional Networks (GCNs). Notably, spatial-based GNNs define graph convolution operator based on graph topology wherein the feature vectors of node’s neighbors are aggregated via a permutation-invariant function. Specifically, [22] proposed a GraphSAGE approach that sampled fixed size neighborhood nodes and used max pooling, mean pooling, and LSTM pooling scheme to aggregate neighbor information. Considering the different weights of node’s neighbors, [46] proposed a Graph Attention Network (GAT) algorithm to calculate attention coefficient and then aggregated the neighborhood information. Other related models include PATCHY-SAN [34], DCNN [4], and further details on GNNs can be found in the review [60].

Following heuristic methods, matrix completion-based methods and network embedding-based methods, neural networks have been gradually applied to link prediction problem and have achieved state-of-the-art results. Specifically, [56] proposed a link prediction method called Weisfeiler-Lehman Neural Machine (WLNM), which labeled nodes using the Weisfeiler-Lehman algorithm and encoded subgraphs to construct a feedforward neural network-based classification model. Next, from the perspective of subgraph classification, [57] proposed a novel GNN-based link prediction framework, SEAL, to learn subgraph structures and node features from local enclosing subgraphs. Along this line, to directly leverage the topology features of local subgraphs, [36] proposed a new random-walk-based pooling scheme, WalkPool, and built features for subgraph classification. Moreover, [18] proposed a neural network-based link prediction method with only one-hop neighborhood information, which demonstrated almost equivalent performance to the WLNM and SEAL. Instead of subgraph classification, [8] converted the original graph into a corresponding line graph and solved the node classification problem for link prediction. To perform link prediction for general directed or undirected complex networks, [48] represented the adjacency matrices of networks as binary images and developed a generative adversarial networks (GANs)-based method. In addition, because existing GNN-based methods do not scale appropriately to large graphs, [30] extracted sparse enclosing subgraphs based on multiple random walks and presented a scalable link prediction solution, called ScaLed. To reduce the time required to determine the distances between two nodes, [27] defined an anchor-based distance and proposed a new distance-enhanced GNN method for link prediction.

Among all existing methods for link prediction, the work closest to the one condidered in this study is the GANs-based method [48]. However, this method predicts potential links via image processing within the GANs framework, whereas the proposed method conducts link prediction via GNNs-based network reconstruction.

Real-world graphs, also known as complex networks, are abstract representation of complex systems and have been extensively studied in the field of network science. Consequently, numerous studies have revealed that complex networks exhibit rich and diverse connectivity patterns. [32] augmented that the organization of real networks usually embodies both regularities and irregularities, where the former can be modeled and decides the extent to which the formation of a network can be explained. Notably, link predictability reflects the structural regularities in real-world networks and denotes the inherent difficulty of link prediction. [53] proposed a self-representation network model-based method, called NetSRE, for measuring and regulating link predictability of networks. [54] proposed a deep linear coding-based link prediction adversarial attack method by disturbing the underlying structural pattern of networks, which proved that links play global structural roles in network organization. Moreover, [7] suggested that high-order connectivity patterns are essential for understanding the fundamental structures of networks and developed a framework that identified clusters of network motifs. [41] claimed that hierarchical structure plays an important role in complex systems. To prove the existence of hierarchical organization, an unsupervised method for extracting the hierarchical organization of complex networks was introduced and validated.

Although real-world graphs exhibit various structural patterns, most existing neural networks-based link prediction methods simply assume that they are flattened and locally isolated, and these methods judge the existence of links explicitly based only on local enclosing subgraphs. With the exception of local structural features, this study focuses on integrating global and hierarchical structural patterns into neural networks for link prediction.

Notations. Let $\mathbcal{G=(V,E)}$ denote an undirected and unweighted graph, where $\mathbcal{V}=\{{\mathbcal{v}_{1}},\cdots,{\mathbcal{v}_{N}}\}$ denotes the set of nodes and $\mathbcal{E}=\{{\mathbcal{e}_{1}},\cdots,{\mathbcal{e}_{M}}\}$ denotes the set of edges. The adjacency matrix of graph $\mathbcal{G}$ is denoted as ${\bf{A}}\in{\{0,1\}^{N\times N}}$, where ${{\bf{A}}_{ij}}=1$ if nodes $i$ and $j$ are connected and ${{\bf{A}}_{ij}}=0$ otherwise. Each edge $\mathbcal{e}$ can be represented as a node pair $({\mathbcal{u}},{\mathbcal{v}},)$, where $\mathbcal{u},\mathbcal{v}\in\mathbcal{V}$. Let $\mathscr{N}(u)$ denote the neighbors of node $u$, $\mathscr{N}(u)=\mathbcal{\{v|(u,v)\in E\}}$.

Link Prediction. Given an observed graph $\mathbcal{G_{o}=(V,E_{o})}$ that corresponds to the original graph $\mathbcal{G=(V,E)}$, link prediction aims to infer the presence or absence of an edge between a pair of target nodes based on $\mathbcal{G_{o}}$, thereby generating a recovered graph $\mathbcal{G^{*}}$ to approximate the original graph $\mathbcal{G}$. In particular, the prediction problem involves identifying a function that generates a likelihood score for a pair of nodes $(\mathbcal{u},\mathbcal{v})\notin\mathbcal{E}$ to infer the missing link $(\mathbcal{u},\mathbcal{v})$, or to produce a likelihood score for an existing edge $(\mathbcal{u},\mathbcal{v})\in\mathbcal{E}$ to identify spurious links. Thus, the link prediction problem can be formulated as ${s_{\mathbcal{uv}}}=f(\mathbcal{u},\mathbcal{v},{\bf{A}}|\theta)$, where $\theta$ denotes the parameter of link prediction model. In this work, ${\mathbcal{E}_{m}}$ and ${\mathbcal{E}_{s}}$ denote the identified missing and spurious links, respectively.

Note that data augmentation is a set of techniques that increases the amount and diversity of data by creating reasonable virtual data points from existing data, such that better machine learning models can be constructed based on them. According to [59], this study considers graph data augmentation and adopts a random mapping mechanism to produce augmented graph set $\mathbcal{D}$ based on the observed graph $\mathbcal{G_{o}=(V,E_{o})}$. Specifically, the set of all possible edges in the graph $\mathbcal{G_{o}}$ is denoted as $\Omega$, the existing edge set is denoted as $\mathbcal{E}_{o}$, and the non-existing edge set is denoted as $\mathbcal{E}_{non}=\Omega-\mathbcal{E}_{o}$. Thus, the candidate sets for random mapping are defined as follows:

Thereafter, samples are randomly produced from the candidate sets to obtain the edge sets $\mathbcal{E}_{del}$ and $\mathbcal{E}_{add}$. Finally, a new augmented graph is generated by modifying the graph $\mathbcal{G_{o}}$ based on $\mathbcal{E}_{del}$ and $\mathbcal{E}_{add}$:

Each input graph can be viewed as an instance for link prediction, owing to the generative learning scheme of the models considered in this work. Thus, the dataset containing a series of augmented graphs can be denoted as $\mathbcal{D}=\{{\mathbcal{G}_{i}}|i=1,...,t\}$ and split to yield disjoint training and validation sets. These can be denoted as ${\mathbcal{D}_{train}}$ and ${\mathbcal{D}_{val}}$ respectively, wherein the missing and spurious links of the validation set are guaranteed not to appear in the training set. The observed graph $\mathbcal{G_{o}}$ used to generate the augmented graphs is defined as test set ${\mathbcal{D}_{test}}$.

GCNs are a class of neural networks designed to generalize traditional convolution operator for non-euclidean graph-structured data. In essence, GCNs aim to learn new feature representations of nodes in graphs by exploiting their structural information. Let adjacency matrix ${\bf{A}}\in{\{0,1\}^{N\times N}}$ denote the structural information of the graph $\mathbcal{G}$, and ${\bf{X}}\in{\mathbb{R}^{N\times F}}$ denote the feature matrix of all graph nodes. Mathematically, using the output of the $l$-th layer as the input for the next layer, each neural network layer can be formulated as a nonlinear function:

where ${{\bf{H}}^{(l)}}$ corresponds to the feature matrix of the $l$-th layer, and ${{\bf{H}}^{(0)}}=\bf{X}$ is the input feature matrix of the first layer. Specific GCNs models differ only in the manner in which the nonlinear function $f(\cdot)$ is instantiated. A simple example of $f(\cdot)$ is as follows:

where $\sigma(\cdot)$ denotes a nonlinear activation function, such as a Rectified Linear Unit (ReLU), and ${{\bf{W}}^{(l)}}$ denotes a trainable weight matrix for the $l$-th layer. With this propagation rule, the neighbour’s features are aggregated to represent each node at every layer, and the features become increasingly abstract by stacking layers on top of each other. However, there exist two limitations: the propagation rule simply aggregates the features of neighboring nodes but not the node itself, and the multiplication with ${\bf{A}}$ expected to change the scale of the feature vectors. That is, the nodes with a high degree will have a larger value, and the nodes with a low degree may have smaller values. To address the problems, a new propagation function, $f(\cdot)$, is presented as follows:

where $\hat{\bf{A}}$ is obtained by adding an identity matrix ${\bf{I}}$ to the adjacency matrix $\hat{\bf{A}}={\bf{A}}+{\bf{I}}$, ${{\hat{\bf{D}}}}$ denotes the diagonal node degree matrix of $\hat{\bf{A}}$, and ${{\hat{\bf{D}}}^{-\frac{1}{2}}}\hat{\bf{A}}{{\hat{\bf{D}}}^{-\frac{1}{2}}}$ denotes symmetric normalization.

Traditionally, Principal component analysis (PCA) was proposed to determine a low-dimensional representation of data while retaining as much information as possible. However, the PCA is particularly effective when dealing with Gaussian noise, which is independent and identically distributed with respect to the original data. Hence, the Robust Principal Component Analysis (RPCA) [9] has been proposed to eliminate the effect of erratic noise (outliers). PCA and RPCA methods implicitly assume that the underlying data structure is a single low-rank subspace; however, real-world data may be drawn from a union of multiple subspaces, and therefore, modeling may be inaccurate. To this end, Low-Rank Representation (LRR) [28] has been proposed.

Considering the correlation between the connectivity patterns of nodes in real-world graphs, the adjacency matrix of the graphs should be low-rank. In other words, the rows or columns of the adjacency matrix must not be linearly independent. Thus, assuming that hidden non-zero entries representing missing links can be recovered according to the adjacency matrix, [37] proposed an RPCA-based link prediction method, which is formulated as the following optimization problem:

where ${\rm{rank}}({{\bf{X}}^{*}})$ denotes the rank of matrix ${{\bf{X}}^{*}}$, $||\cdot|{|_{0}}$ is the ${\ell_{0}}-$norm, and $\gamma$ denots the balancing parameter. The method searches for ${{\bf{X}}^{*}}$ with a low rank as low as possible and ${\bf{E}}$ as sparse as possible from $\bf{A}$. Moreover, by representing a network structure with as few representative subgraphs as possible, [53] proposed an LRR-based link prediction method, wherein networks could be modeled via a low-rank and sparse representation, as follows:

where ${\bf{Z}}$ denotes the representation matrix reflecting the organization principle of the network, and $||\cdot|{|_{2,1}}$ is the ${\ell_{2,1}}-$norm.

The notations used in this study are listed in Table 1.

[32] suggested that link formation in real-world graphs is usually driven by both regular and irregular factors, and the former can be explained based on the mixture of multiple mechanisms, such as homophily, triadic closure, preferential attachment. Meanwhile, assuming that high-dimensional data are a mixture of simple data and are drawn from a union of multiple low-dimensional linear subspaces, the LRR has been proposed to represent the data ${\bf{A}}=[{a_{1}},{a_{2}},...,{a_{N}}]$ as a linear combination of the basis in a ”dictionary” ${\bf{D}}=[{d_{1}},{d_{2}},...,{d_{M}}]$:

Thus, the optimal representation matrix ${{\bf{Z}}^{*}}$ uncovers the underlying subspaces in the data. By using each subspace to model a homogeneous subset of the data, multiple subspaces in LRR can capture heterogeneous structures within the data. Therefore, considering the above ideas, the regular structure of real-world graphs can be described appropriately by the LRR model, wherein the generation mechanisms of graph organization essentially corresponds to subspaces and the low rankness constraint captures the global correlation in graphs. Meanwhile, based on the generation mechanisms of graph organization, individual nodes may have similar connection patterns, and substructures that follow the same generation mechanism can be represented by each other, as depicted in Figure 2(b). Therefore, by using the adjacency matrix ${\bf{A}}$ as the dictionary, the real-world graph can be represented by itself, as follows:

In addition to their regular structure, real-world graphs also contain irregular components. Thus, we let matrix ${\bf{E}}$ denote such irregular connections; then, the proposed self-representation model can be modified as ${\bf{A}}={\bf{AZ}}+{\bf{E}}$. According to the LRR, data are considered to be ”sample specific”, and the ${\ell_{21}}-$norm is adopted to constrain the matrix ${\bf{E}}$, i.e., $||{\bf{E}}|{|_{2,1}}$. However, although the proposed method can be used to model real-world graphs, the low-rank model and ${\ell_{21}}-$norm constraints are usually solved using Alternating direction method (ADM), which requires a large number of iterations and has high complexity. Therefore, a reasonable strategy is to relax the constraints with Frobenius norm:

Let $\mathbcal{L}=||{\bf{Z}}||_{\mathbcal}{F}^{2}+\lambda||{\bf{A}}-{\bf{AZ}}||_{\mathbcal}{F}^{2}$ denote the partial derivative of $\mathbcal{L}$ with respect to ${\bf{Z}}$ is ${{\partial\mathbcal{L}}\mathord{\left/{\vphantom{{\partial\mathbcal{L}}{\partial{\bf{Z}}}}}\right.\kern-1.2pt}{\partial{\bf{Z}}}}=2{\bf{Z}}+\lambda(-2{{\bf{A}}^{\rm{T}}}{\bf{A}}+2{{\bf{A}}^{\rm{T}}}{\bf{AZ}})$. By setting ${{\partial L}\mathord{\left/{\vphantom{{\partial L}{\partial{\bf{Z}}}}}\right.\kern-1.2pt}{\partial{\bf{Z}}}}=0$, the optimal representation ${{\bf{Z}}^{*}}$ can be obtained as follows:

where ${\bf{I}}$ denotes the identity matrix. Thus, in the case that the clean data is sufficient enough to represent the graph’s structural patterns and the irregular connections are properly characterized, the structure perturbations can be inferred using ${\bf{A}}{{\bf{Z}}^{*}}$. Hence, the collaborative inference operation is defined as follows:

According to local similarity indices for link prediction, the more the number of paths two nodes possess, the greater the similarity between them. Specifically, two nodes with a high mutual connectivity are more likely to generate a link between them. Thus, n-hop-based $(n\geq 2)$ paths must be explored to characterize the local structural features for link prediction. Using a deep learning framework, the $n$-hop computation can be decomposed into two-hop operations on each neural layer. Hence, a high-order connectivity computation calculates the two hop connectivity of graph nodes in each layer, and the mutual connectivity of two nodes can be estimated by stacking the high-order connectivity computation mechanism. Assuming that the integer powers of the adjacency matrix characterizes the mutual connectivity of graph nodes, that is, ${[{{\bf{A}}^{\rm{n}}}]_{ij}}$ denotes the number of paths with length $n$ connecting nodes $i$ and $j$, the high-order connectivity computation in each neural layer can be defined based on the idea of the second power of adjacency matrix ${\bf{A}}$. From the perspective of graph convolution networks, high-order connectivity computation can be defined as

where the weighted adjacency matrix generated by the proposed collaborative inference operation is viewed as the features of graph nodes. Figure 3 illustrates a high-order connectivity computation. As presented in Equation (15), the global and local structural features can be captured for link prediction at the level of individual nodes and edges. Thus, the nonlinear propagation function can be defined as follows:

Thus, the hierarchical structure of real-world graphs can be characterized by executing the nonlinear propagation function iteratively, in which ${\rm{\mathbcal{HCCA}}}({{\bf{H}}^{(l)}})={{{\bf{\hat{D}}}}^{{\rm{-}}\frac{1}{2}}}{\bf{\hat{A}}}{{{\bf{\hat{D}}}}^{{\rm{-}}\frac{1}{2}}}\mathbcal{CI}({{\bf{H}}^{(l)}})$ represents the high-order connectivity of graph nodes, as depicted in Figure 2(c), and $\mathbcal{CI}({{\bf{H}}^{(l+1)}})$ denotes the collaborative inference.

To estimate the likelihood of potential links, the output of the $(l-1)$-th layer, i.e., ${{\bf{H}}^{(l)}}$, is fed as the input of the $l$-th layer. Based on $\mathbcal{CI}({{\bf{H}}^{(l)}})$ and ${\rm{\mathbcal{HCCA}}}({{\bf{H}}^{(l)}})$, the shallow layers extract the low-order global and local structure features, while the deep layers extract the high-order global and local structure features. Meanwhile, the effective range that the local structure features drawn from increases as the model depth increases. Therefore, the structure features in different range at various order, i.e., ${\rm{\mathbcal{HCCA}}}({{\bf{H}}^{(l)}})$ and $\mathbcal{CI}({{\bf{H}}^{(l)}})$, $0\leq l\leq L$, all contribute to the inference of potential links, although the exact extent of their contribution depends on the graph data.

To overcome the issues mentioned above, in addition to being used as the inputs of the next layer, the outputs of neural network layers are mapped to skip a block of several layers based on residual connections, as illustrated in Figure 2(a). Next, all outputs are concatenated and used as the input of a two-layer Multi-layer Perceptron (MLP), which is defined as:

where ${\bf{O}}$ is a vector containing the probabilities of links between all possible node pairs, and missing and spurious links can be inferred based on it.

To train the proposed model, augmented graphs generated based on the observed graph are used as training data, and the adjacency matrix of the observed graph is flatten as its labels ${\bf{Y}}$, where ${{{\bf{Y}}_{i*N+j}}}$ denotes the existence of the link between nodes $i$ and $j$. Correspondingly, ${\bf{O}}$ represents the prediction results obtained by the proposed model for all possible links. Here, the Binary Cross-Entropy (BCE) is used as the loss function:

The learned model is then deployed on the observed graph to reconstruct the original graph. The training process of GraphLP is outlined in Algorithm 1.

(1) Generalized local similarity indices. The high-order connectivity computation ${\rm{\mathbcal{HCCA}}}({{\bf{H}}^{(l)}})$ in every neural network layer is essentially the second power of the adjacency matrix, and it obtains the connectivity of node pairs within two-hop neighborhood. As the model depth increases, the connectivity of node pairs in a wider range is considered. Thus, GraphLP can degenerate to ${S_{ij}}={{\bf{A}}^{2}}+\alpha{{\bf{A}}^{3}}+\beta{{\bf{A}}^{4}}+\cdot\cdot\cdot$ when collaborative inference and deep learning mechanism are abolished.

(2) Connection to WalkPool. WalkPool [36] first generates node representations based on GNN and encodes them into edge weights of the extracted enclosing subgraphs; following this, it uses the edge weights to compute the transition probabilities of random walk. Next, the method calculates a list of features based on the transition probabilities to classify the subgraphs. However, for an enclosing subgraph $\mathbcal{G}=(\mathbcal{V},\mathbcal{E})$, its variants ${\mathbcal{G}^{+}}=(\mathbcal{V},\mathbcal{E}\cup\{i,j\})$ and ${\mathbcal{G}^{-}}=(\mathbcal{V},\mathbcal{E}\backslash\{i,j\})$ are used as positive and negative samples, respectively. In essence, this method discriminates only those subgraphs that differ by a single edge and is not suitable for practical link prediction scenarios. In contrast, GraphLP can predict any potential links based on graph structure features.

(3) Connection to LFLP. The LFLP [53] constructs an adjacency matrix based on a self-representation model and then combines it with the observed network to identify missing and spurious links. The collaborative inference operation $\mathbcal{CI}({{\bf{H}}^{(l)}})$ of our work is similar to that in the LFLP with respect to modeling the global structure of graphs; however, the difference is that only low-order global structural features are considered in LFLP, whereas multi-order global and local structural features are characterized based on the deep-learning framework in GraphLP.

For 90$\%$ of the observed links, the results about the AUC and AP with standard deviations are presented in Table 4 and 5, which indicate that GraphLP significantly outperforms other state-of-the-art algorithms in terms of both AUC and AP, with exception of the NS and Router datasets. The results demonstrate that the learning of local and global graph structure entirely characterizes the underlying structural patterns; thus, the missing links and spurious links can be better identified. Table 4 indicates that GraphLP significantly improves the AUC on the PB, C.ele, and Router datasets, with approximately 4$\%$, 7$\%$, and 3$\%$ performance improvement, respectively, compared to the WP algorithm. In addition, the proposed method still performs better than other state-of-the-art methods on the USAir, Yeast and NS datasets. Moreover, the results for the AP presented in Table 5 also indicate that GraphLP outperforms state-of-the-art methods on most of datasets, and GraphLP achieves a maximum performance enhancement of approximately 9$\%$ compared to the best performing graph neural network method WP.

For 50$\%$ of the observed links, the results also demonstrate that the proposed model achieves remarkable performance compared to the methods, as described in Table 6 and 7. The results illustrate that, as the amount of structure perturbation increases, GraphLP can still appropriately learn the real graph structure, thus recovering the original graph effectively. Therefore, the values of the AUC and AP decreased to a lower extent. Furthermore, by comparing Table 4 with Table 6 and Table 5 with Table 7, we can infer that the AUC and AP values drop faster for the other state-of-the-art methods than those for GraphLP, which demonstrates that GraphLP can better capture the underlying structural patterns to demonstrate better performance.

To further verify the proposed method, the precision of GraphLP and traditional link prediction methods are calculated based on the following datasets: (1) Macaque [15], cortical networks of the macaque monkey; (2) Mangwet [5], the food web in Mangrove Estuary during the wet season; (3) Jazz [19], a collaboration network of jazz musicians; (4) Metabolic [17], a metabolic network of C.elegans; (5) USAir, (6) C.ele, (7) E.coli and (8) Yeast. Here, six representative traditional link prediction methods are selected for comparison.

(1) The CN [29] metric is among the most widely used methods for link prediction problem. It assumes that two nodes will be more easily connected if they share more common neighbors;

(2) The RA [61] metric is inspired by the physical processes involved in resource allocation, which suppresses the contribution of high-degree common neighbors;

(3) The LP [61] index measures the structural similarity of node pairs within three-hops;

(4) The Non-negative Matrix Factorization (NMF) [26] model is used for structure prediction by learning the latent features of real-world graphs;

(5) The Robust Principal Component Analysis (RPCA) [37] represents a real-world graph based on the sparsity and low rank property of its adjacency matrix and infers potential links based on matrix completion.

(6) The LFLP [53] uses a self-representation model to reconstruct the original graph based on a few representative subgraphs.

The results of missing link prediction with respect to Precision are shown in Table 8. For each network, the bold number in the corresponding column indicates the highest accuracy. The results presented in Table 8 demonstrate that the proposed model GraphLP model performs the best among the methods. Furthermore, the link prediction accuracy of the proposed model is far higher than that of the other methods, which can be at least three times higher than that of the best-performing method. For spurious links prediction, the results measured by Precision are listed in Table 9. For all networks, GraphLP performs the best among the methods and is remarkably better than the second best algorithm. The results presented in Table 8 and 9 demonstrate that GraphLP has stronger ability to learn structural features, and can recover the structure of the original network more accurately. Based on Table 2, it can be observed that the Precision of our proposed model performs best, despite the large differences between the ACC and AD across all the datasets; thus indicates that the proposed model performs well for heterogeneous graph structures.

To verify the effectiveness of the proposed model of missing and spurious links inference, the topology of the recovered graphs in the model training process on Club dataset is visually compared, as depicted in Figure 4 and 5. The top half depicts the topology of the graphs, wherein the red links denote the missing links, the blue links denote the spurious links, and the gray links denote the original links. The bottom half depicts the likelihood scores of missing and spurious links. Based on the results, it can be concluded that with an increase in epoch times, the likelihood scores of missing links increases gradually, and the likelihood scores of spurious links decline gradually. The widths of the lines indicate the following process: when the number of epochs reaches 100, the likelihood scores of missing links approach 1.0, and the likelihood scores of spurious links approach 0.0. This proves that the proposed GraphLP model can distinguish between missing links and spurious links and infer them effectively. Moreover, to further prove the effectiveness of the proposed model, the topology of the recovered graphs in the model training process is visualized when 20% of the links of Club are perturbed, as depicted in Figure 5. Compared to Figure 4, we determined that the likelihood of missing and spurious links is weakened with an increase in the structure perturbation ratio. However, with an increase in the training epochs, the model is still able to distinguish and infer the missing and spurious links according to the structural patterns. For instance, when the training epoch reaches 140, the likelihood scores of missing links are greater than 0.5, and those of spurious links are less than 0.3, indicating that the proposed model can still predict missing and spurious links with high accuracy.

Next, the performance of GraphLP is explored at various model depths. As depicted in Figure 6, the performance of the model in terms of the AUC and AP trained with one-layer neural network is poor; however, its performance improves significantly with an increase in the number of layers. In particular, when the model depth is equal to two, a significant performance improvement is noted. The primary reason for this is that the model with two-layers multi-order global and local structural features is integrated adaptively based on the MLP component, which considerably improves the performance of the model. Subsequently, as the layer number increases, a slight improvement in model performance is still noted. When the number of layers is four, the accuracy of GraphLP declines significantly on NS and fluctuates on other datasets. A possible reason for this is that the model with four layers becomes more complex, thereby requiring more training iterations or an appropriate learning rate [14]. In general, the performance of the proposed model is optimal when the depth is three, and a deep architecture is necessary.

To examine the sensitivity of the proposed model to the trade-off parameter, the AUC and AP values of link prediction methods with different $\lambda$ are presented in Figure 7. Based on the results, it can be concluded that the performance of the proposed model is not sensitive to $\lambda$ for most datasets. In Figure 7(a), for the USAir and NS datsets, the AUC value varies significantly under different $\lambda$, but the performance is still better than other of the other algorithms. In Figure 7(b), the AP value remains stable for different $\lambda$ values, indicating that the proposed model is insensitive to different $\lambda$. Overall, our proposed algorithm exhibited satisfactory performance on most datasets with various $\lambda$.

Generally, GraphLP converges to optimal values after approximitely 200 epochs on most datasets. In particular, Figure 8 plots the learning curves of GraphLP on the USAir and C.ele datasets, including the training loss, validation AUC, validation AP, and validation loss. The results indicate that the AUC and AP values increase rapidly with the decrease in training loss and validation loss, and these values converge to the optimal value when the validation loss approaches a minimum value. Additionally, we discover that validation loss is lower than training loss, and the difference between them remains relatively stable. A possible reason for this is that the dropout manipulation is only applied to the training process.