A Practical Tutorial on Graph Neural Networks

The key contributions of this tutorial article are as follows:

The structure and taxonomy of this article is outlined in Table 1.

Graphs are formally defined by a set of vertices and the set of edges between these vertices: Put formally, G = G (V, E). Fundamentally, graphs are just a way to encode data, and in that way, every property of a graph represents some real element, or concept in the data. Understanding how graphs can be used to represent complex concepts is key in appreciating their expressiveness and generality as an encoding device (see Figure 1 for examples of this domain agnostic expressiveness).

Vertices represent items, entities, or objects, which can naturally be described by quantifiable attributes and their relationships to other items, entities, or objects. We refer to a set of \(\vert {\mbox{V}}\vert\) vertices as V and the ith single vertex in the set as v \(_{i}\) . Note that there is no requirement for all vertices to be homogenous in their construction.

Edges represent and characterize the relationships that exist between items, entities, or objects. Formally, a single edge can be defined with respect to two (not necessarily unique) vertices. We refer to a set of \(\vert {\mbox{E}} \vert\) edges as E and a single edge between the ith and jth vertices as e \(_{ij}\) .

Neighborhoods are subgraphs within a graph and represent distinct groups of vertices and edges. Most commonly, the neighborhood \(\mathcal {N}_{{\it v_{i}}}\) centered around a vertex v \(_{i}\) comprises of v \(_{i}\) , its adjoining edges (where e \(_{ij}\) \(= 1\) ), and the vertices that are directly connected to it. Neighborhoods can be iteratively grown from a single vertex by considering the vertices attached (via edges) to the current neighborhood. Note that a neighborhood can be defined subject to certain vertex and edge feature criteria (i.e., all vertices within two hops of the central vertex, rather than one hop).

Features are quantifiable attributes that characterize a phenomenon that is under study. In the graph domain, features can be used to further characterize vertices and edges. Extending our social network example, we might have features for each person (vertex) that quantifies the person’s age, popularity, and social media usage. Similarly, we might have a feature for each relationship (edge) that quantifies how well two people know each other, or the type of relationship they have (familial, colleague, etc.). In practice there might be many different features to consider for each vertex and edge, so they are represented by numeric feature vectors referred to as v \(_{i}^{\,F}\) and e \(_{ij}^{\,F}\) , respectively.

Embeddings are compressed feature representations. If we reduce large feature vectors associated with vertices and edges into low dimensional embeddings, then it becomes possible to classify them with low-order models (i.e., if we can make a dataset linearly separable). A key measure of an embedding’s quality is if the points in the original space retain the same similarity in the embedding space. Embeddings can be created (or learned) for vertices, edges, neighborhoods, or graphs. Embeddings are also referred to as representations, encodings, latent vectors, or high-level feature vectors, depending on the context.

Output types change depending on the problem domain. A GNN’s forward pass can be thought of as two key processes: converting input graphs into useful embeddings, performing some downstream task (e.g., classification) on the embeddings, which converts the embeddings into some useful output. We define three commonly observed output types as follows:

Transductive learning methods are exposed to all of the training and testing data before making predictions. For example: Our dataset might consist of a single large graph (e.g., Facebook’s social network graph) and the set of vertices is only partially labelled. The training set consists of the labelled vertices, and the testing set consists of both a small set of labelled vertices (for benchmarking) and the remaining unlabelled vertices. In this case, our learning methods should be exposed to the entire graph during training (including the test vertices), because the additional information (e.g., structural patterns) will be useful to learn from. Transductive learning methods are useful in such cases where it is challenging to separate the training and testing data without introducing biases.

Inductive learning methods reserve separate training and testing datasets. The learning process ingests the training data, and then the learned model is tested using the testing data, which it has not observed before in any capacity.

RGNNs compute embeddings at each vertex in the input graph using a deterministic, shared function called the transition function. It is named the transition function, as it can be interpreted as calculating the next representation of a neighborhood from the neighborhood’s current representation. This transition function can be applied symmetrically at any vertex, even though the size of a vertex’s neighborhood may be variable. This process is illustrated in Figure 2, where the transition function \({f}\:\) calculates an embedding at each vertex for the surrounding neighborhood.

As such, the kth embedding h \(_{i}^{\,k}\) for any given vertex v \(_{i}\) is dependent on the following quantities:

To recurrently apply this learned transition function to compute successive embeddings, \({f}\:\) must have a fixed number of input and output variables. How then can it be dependent on the immediate neighborhood, which might vary in size, depending on where we are in the graph? There are two simple solutions, the first of which is to set a “maximum neighborhood size” and use null vectors when dealing with vertices that have non-existing neighbors [74]. The second approach is to aggregate all neighborhood features in some permutation invariant manner [25], thus ensuring that any neighborhood in the graph is represented by a fixed size feature vector. While both approaches are viable, the first approach does not scale well to “scale-free graphs,” which have degree distributions that follow a power law. Since many real world graphs (e.g., social networks) are scale-free [72], we will use the second solution here. Mathematically, this can be formulated as in Equation (1) [74].

We can see that under this formulation Equation (1), \({f}\:\) is well defined. It accepts four feature vectors that all have a defined length, regardless of which vertex in the graph is being considered, regardless of the iteration. This means that the transition function can be applied iteratively, until a stable embedding is reached for all vertices in the input graph. This expression can be interpreted as passing “messages,” or features, throughout the graph; in every iteration, the embedding h \(_{i}^{\,k}\) is dependant on the features and embeddings of its neighbors. This means that with enough recurrent iterations, information will propagate throughout the whole graph: After the first iteration, any vertex’s embedding encodes the features of the neighborhood within a range of a single edge.

In the second iteration, any vertex’s embedding is an encoding of the features of the neighborhood within a range of two edges away, and so on. The iterative passing of “messages” to generate an encoding of the graph is what gives this message passing framework its name.¹

Note that it is typical to explicitly add the identity matrix I \(_{n}\) to the adjacency matrix A, thus ensuring that all vertices become trivially connected to themselves, meaning that a vertex v \(_{i}\) \(\in\) \(\mathcal {N}_{{\it v_{i}}}\) \(\forall i \in {{\mbox{V}}}\) . Moreover, this allows us to directly access the neighborhood by iterating through a single row of the adjacency matrix. This modified adjacency matrix is usually normalized to prevent unwanted scaling of embeddings.

Once we have useful embeddings centered around each vertex in the graph, the goal is to then inference meaningful outputs based on these values (i.e., to perform a downstream task). The output function \({g}\:\) is responsible for taking the converged embeddings of a graph G (V, E) and creating said output. In practice, the output function \({g}\:\) , much like the transition function \({f}\:\) , is implemented by a feed-forward neural network, though other means of returning a single value have been used, including mean operations, dummy super nodes, and attention sums [110].

Intuitively, the combined process of recurrently computing embeddings and subsequently computing downstream inputs can be interpreted as a sequential process of repeated NN computation blocks—or a finite computation graph (see Figure 2). In a supervised setting, a loss signal can be calculated that quantifies the error between the predicted output and a labelled ground truth. Both \({f}\:\) and \({g}\:\) can then be trained via backpropagation of errors, throughout the “unrolled” computation graph. For more detail on this process, see the calculations in [74].

When discussing recurrence thus far, we have referred mainly to computing techniques that are iteratively applied to neighborhoods in a graph to produce embeddings that are dependent on information propagated throughout the graph. However, recurrent techniques may also refer to computing processes over sequential data, e.g., time series data. In the graph domain, sequential data refers to instances that can be interpreted as graphs with features that change over time. These include spatiotemporal graphs [98]. For example, Figure 1(b) illustrates how a graph can represent a skeletal structure in a single image of a hand, however, if we were to create such a graph for every frame of a contiguous video of a moving hand, then we would have a data structure that could be interpreted as a sequence of individual graphs, or a single graph with sequential features, and such data could be used for classifying hand actions in video.

As is the case with traditional sequential data, when processing each state of the sequence, we want to consider not only the current state but also information from the previous states, as outlined in Figure 6(a). A simple solution to this challenge might be to simply concatenate the graph emeddings of previous states to the features of the current state (as in Figure 6(b)), but such approaches do not capture long-term dependencies in the data. In this section, we outline how existing solutions from traditional DL—such as Long Short-Term Memory Networks (LSTMs) and Gated Recurrent Units (GRUs) (outlined in Figure 6)—can be extended to the graph domain.

Gated Recurrent Units (GRUs) provide a less computationally expensive alternative to GLSTMs by removing the need to calculate a cell state in each cell. As such, GRUs have three learnable weight matrices (as illustrated in Figure 6(d)) that serve similar functions to the four learnable weight matrices in GLSTMs. Again, GRUs require some definition of graph concatenation.

Graph LSTMs (GLSTMs) make use of LSTM cells that have been adapted to operate on graph-based data. Whereas the aforementioned recurrent modules (Figure 6(b)) employ a simple concatenation strategy, GLSTMs ensure that long-term dependencies can be encoded in the LSTM’s “cell state” (Figure 6(c)). This alleviates the vanishing gradient problem where long-term dependency signals are exponentially reduced when backpropagated throughout the network [35, 36].

GLSTM cells achieves this through four key processing elements that learn to calculate useful quantities based on the previous state’s embedding and the input from the current state (as illustrated in Figure 6(c)).

To use GLSTMs, we need to define all the operators in Figure 6(e). Since a graph G (V, E) can be thought of as a variably sized set of vertices and edges, we can define graph concatenation as the separate concatenation of vertex features and edge features, where some null padding is used to ensure that the resultant tensor is of a fixed-size. This can be achieved by defining some “max number” of vertices for the input graphs. If the input signal for the GLSTM cell has a fixed size, then all other operators can be interpreted as traditional tensor operations, and the entire process is differentiable when it comes to backpropagation.

In this independent example, we investigate social networks, which represent a rich source of graph data. Due to the popularity of social networking applications, accurate user and community classifications have become exceedingly important for the purpose of analysis, marketing, and influencing. In this example, we look at how the recurrent application of a transition function aids in making predictions on the graph domain, namely, in graph classification.

We will be using the GitHub Stargazer’s dataset [67] (availablehere). GitHub is a code sharing platform with social network elements. Each of the \(12{,}725\) graphs is defined by a group of users (vertices) and their mutual following relationships (undirected edges). Each graph is classified as either a web development group or a machine learning development group. There are no vertex or edge features—all predictions are made entirely from the structure of the graph.

Rather than use a true RGNN that applies a transition function to hidden states until some convergence criteria is reached, we will instead experiment with limited applications of the transition function. The transition function is a simple message passing aggregator that applies a learned set of weights to create size 16 hidden vector representations. We will see how the prediction task is affected by applying this transition function \(1, 2, 4,\) and 8 times before feeding the hidden representations to an output function for graph classification. We train on 8,096 graphs for 16 epochs and test on 2,048 graphs for each architecture.

As expected, successive transition functions result in more discriminative features being calculated, thus resulting in a more discriminative final representation of the graph (analogous to more convolutional layers in a CNN).

In fact, we can see that the final hidden representations become more linearly separable (see TSNE visualizations in Figure 5), thus, when they are fed to the output function—a linear classifier—the predicted classifications are more often correct. This is a difficult task, since there are no vertex or edge features. State-of-the-art approaches achieved the following mean AUC values averaged over 100 random train/test splits for the same dataset and task: GL2Vec [10]—0.551, Graph2Vec [62]—0.585, SF [16]—0.558, FGSD [89]—0.656. Here, our transition function \({f}\:\) was a “feedforward NN” with just one layer, so more advanced NNs (or other) implementations of \({f}\:\) might result in more performant RGNNs. As more rounds of transition function were applied to our hidden states, the performance—and required computation—increased. Ensuring a consistent number of transition function applications is key in developing simplified GNN architectures and in reducing the amount of computation required in the transition stage. We will explore how this improved concept is realized through CGNNs in Section 4.2.

GRUs are well suited to sequential data when repeating patterns are less frequent, whereas LSTM cells perform well in cases where more frequent pattern information needs to be captured. LSTMs also have a tendency to overfit when compared to GRUs, and as such GRUs outperform LSTM cells when the sample size is low [27].

In this section, we have explained the forward pass that allows an RGNN to produce useful predictions over graph input data. During the forward pass, a transition function \({f}\:\) is recursively applied to an input graph to create high-level features for each neighborhood. The repeated application of \({f}\:\) ensures that at iteration \(k\) , an embedding h \(_{i}^{\,k}\) includes information from vertices \(k\) edges away from \({\it v_{i}}\) . These high-level features can be fed to an output function to solve downstream tasks. During the backward pass, the parameters for the NNs \({f}\:\) and \({g}\:\) are updated with respect to a loss that is backpropagated through the computation graph defined in the forward pass. Recurrent processing units can also refer to approaches for handling graph-based sequential data, which include graph-based extensions to LSTMs and GRUs.

In actuality, the formulation for calculating embeddings provided in Equation (1) represents only one approach to calculating embeddings. This approach will be contextualised in Section 4, where a broader perspective on calculating useful embeddings will be introduced.

While RGNNs offer a simple approach to working with generic graphs, they have a number of shortcomings. Namely, the shared transition function \({f}\:\) means that the same weights are being used to extract features in successive iterations, which may not be ideal for deep learning scenarios where the relationships between low-level features (earlier in the network) are different to the relationships between high-level features (later in the network). Moreover, since RGNNs iterate until convergence, they have variable length encoding networks, which can add implementation complexities. In the next section, we will discuss how these issues can be alleviated by developing formal definitions of convolution in the graph domain.

We define convolution generally as an operation whereby an output is derived from twogiven inputs by integration or summation, which expresses how the one is modified by the other.

Convolution in CNNs involves two matrix inputs: one is the previous layer of activations, and the other is a matrix \(W \times H\) of learned weights, which is “slid” across the activation matrix, aggregating each \(W \times H\) region using a simple linear combination (see Figure 7(a)). In the spatial graph domain, it seems that this type of convolution is not well defined [78]; the convolution of a rigid matrix of learned weights must occur on a rigid structure of activation. How do we reconcile convolutions on unstructured inputs such as graphs?

Note that at no point during our general definition of convolution was the structure of the given inputs alluded to. In fact, convolutional operations can be applied to continuous functions (e.g., audio recordings and other signals), N-dimensional discrete tensors (e.g., semantic vectors in 1D, and images in 2D), and so on. During convolution, one input is typically interpreted as a filter (or kernel) being applied to the other input, and we will adopt this language throughout this section. Specific filters can be utilized to perform specific tasks: In the case of audio recordings, high pass filters can be used to filter out low-frequency signals, and in the case of images, certain filters can be used to increase contrast, sharpen, or blur images. In our previous example of CNNs, filters are learned rather than designed.

One might consider the early RGNNs described in Section 3 as using convolutional operations. In fact, these methods meet the criteria of locality, scalability, and interpretability. First, Equation (1) only operates over the neighborhood of the central vertex v \(_{i}\) , and can be applied on any neighborhood in the graph due to its invariance to permutation and neighborhood size. Second, the NN \({f}\:\) is dependent on a fixed number of weights and has a fixed input and output that is independent of \(\vert {\mbox{V}}\vert\) . Finally, the convolution operation is immediately interpretable as a generalization of image-based convolution: In image-based convolution, neighboring pixel values are aggregated to produce embeddings; in graph-based spatial convolution, neighboring vertex features are aggregated to produce embeddings (see Figure 7). This type of graph convolution is referred to as the spatial graph convolutional operation, since spatial connectivity is used to retrieve the neighborhoods in this process.

Although the RGNN technique meets the definition of spatial convolution, there are numerous improvements in the literature. For example, the choice of aggregation function is not trivial—different aggregation functions can have notable effects on performance and computational cost.

A notable framework that investigated aggregator selection is the GraphSAGE framework [29], which demonstrated that learned aggregators can outperform simpler aggregation functions (such as taking the mean of embeddings) and thus can create more discriminative, powerful vertex embeddings. Regardless of the aggregation function, GraphSAGE works by computing embeddings based on the central vertex and an aggregation of its neighborhood (see Equation (2)). By including the central vertex, it ensures that vertices with near identical neighborhoods have different embeddings. GraphSAGE has since been outperformed on accepted benchmarks [20] by other frameworks [6], but the framework is still competitive and can be used to explore the concept of learned aggregators (see Section 4.3).

Alternatively, Message Passing Neural Networks (MPNNs) compute directional messages between vertices with a message function that is dependent on the source vertex, the destination vertex, and the edge connecting them [25]. Rather than aggregate the neighbor’s features and concatenating them with the central vertex’s features as in GraphSAGE, MPNNs sum the incoming messages, and pass the result to a readout function alongside the central vertex’s features (see Equation (3)). Both the message function and readout function can be implemented with simple NNs in practice. This generalizes the concepts outlined in Equation (1) and allows for more meaningful patterns to be identified by the learned functions.

One of the most popular spatial convolution methods is Graph Convolutional Networks (GCNs), which produce embeddings by summing features extracted from each neighboring vertex and then applying non-linearity [97]. These methods are highly scalable, local, and furthermore, they can be “stacked” to produce layers in a CGNN. Each of these features is normalized based on the relative neighborhood scales of the current and neighbor vertex, thus ensuring that embeddings do not “explode” in scale during the forward pass.

Graph Attention Networks (GATs) extend GCNs: Instead of using the size of the neighborhoods to weight the importance of v \(_{i}\) to v \(_{j}\) , they implicitly calculate this weighting based on the normalized product of an attention mechanism [87]. In this case, the attention mechanism is dependent on the embeddings of two vertices and the edge between them. Vertices are constrained to only be able to attend to neighboring vertices, thus localizing the filters. GATs are stabilized during training using multi-head attention and regularization and are considered less general than MPNNs [88]. Although GATs limit the attention mechanism to the direct neighborhood, the scalability to large graphs is not guaranteed, as attention mechanisms have compute complexities that grow quadratically with the number of vertices being considered.

Interestingly, all of these approaches consider information from the direct neighborhood and the previous embeddings, aggregate this information in some symmetric fashion, apply learned weights to calculate more complex features, and “activate” these results in some way to produce an embedding that captures non-linear relationships.

In this section, we discuss another class of convolution approaches that evolved from the perspective of Graph Signal Processing (GSP) [78, 82]. These methods are attractive, as they are well grounded in a formal definition of convolution and can be directly interpreted as signal processing techniques in the domain of graph structured data.

The path to defining spectral graph convolution is described by the following series of statements.

The GraphSAGE (SAmple and aggreGatE) algorithm [29] emerged in 2017 as a method for not only learning useful vertex embeddings, but also for predicting vertex embeddings on unseen vertices. This allows powerful high-level feature vectors to be produced for vertices that were not seen at train time; enabling us to effectively work with dynamic graphs, or very large graphs (>100,000 vertices).

In this example, we use the Cora dataset (see Figure 8) as provided by the deep learning library DGL [92]. The Cora dataset is oft considered “the MNIST of graph-based learning” and consists of \(2{,}708\) scientific publications (vertices), each classified into one of seven subfields in AI (or classes). Each vertex has a 1,433 element binary feature vector, which indicates if each of the 1,433 designated words appeared in the publication.

GraphSAGE operates on a simple assumption: vertices with similar neighborhoods should have similar embeddings. In this way, when calculating a vertex’s embedding, GraphSAGE considers the vertex’s neighbors’ embeddings. The function that produces the embedding from the neighbors’ embeddings is learned, rather than the embedding being learned directly. Consequently, this method is not transductive, it is inductive, in that it generates general rules that can be applied to unseen vertices, rather than reasoning from specific training cases to specific test cases.

Importantly, the GraphSAGE loss function is unsupervised and uses two distinct terms to ensure that neighboring vertices have similar embeddings and distant or disconnected vertices have embedding vectors that are numerically far apart. This ensures that the calculated vertex embeddings are highly discriminative.

In this worked example, we experiment by changing the aggregator functions used in each GNN and observe how this affects our overall test accuracy. In all experiments, we use 2 hidden GraphSAGE convolution layers, 16 hidden channels (i.e., embedding vectors have 16 elements), and we train for 120 epochs before testing our vertex classification accuracy. We consider the mean, pool, and LSTM (long short-term memory) aggregator functions.

The mean aggregator function sums the neighborhood’s vertex embeddings and then divides the result by the number of vertices considered. The pool aggregator function is actually a single fully connected layer with a non-linear activation function that then has its output element-wise max pooled. The layer weights are learned, thus allowing the most important features to be selected. The LSTM aggregator function is an LSTM cell. Since LSTMs consider input sequence order, this means that different orders of neighbor embedding produce different vertex embeddings. To minimize this effect, the order of the input embeddings is randomized. This introduces the idea of aggregator symmetry; an aggregator function should produce a constant result, invariant to the order of the input embeddings.

The mean, pool, and LSTM aggregators score test accuracies of \(\mathbf {66.0\%}\) , \(\mathbf {74.4\%}\) , and \(\mathbf {68.3\%}\) , respectively. As expected, the learned pool and LSTM aggregators are more effective than the simple mean operation, though they incur significant training overheads and may not be suitable for smaller training graphs or graph datasets. Indeed, in the original GraphSAGE paper [29], it was found that the LSTM and pool methods generally outperformed the mean and GCN aggregation methods across a range of datasets.

At the time of publication, GraphSAGE outperformed the state-of-the-art on a variety of graph-based tasks on common benchmark datasets. Since that time, a number of inductive learning variants of GraphSAGE have been developed, and their performance on benchmark datasets is regularly updated.²

The Laplacian is a second order differential operator that is calculated as the divergence of the gradient of a function in Euclidean space. The Laplacian occurs naturally in equations that model physical interactions, including but not limited to electromagnetism, heat diffusion, celestial mechanics, and pixel interactions in computer vision applications. Similarly, it arises naturally in the graph domain, where we are interested in the “diffusion of information” throughout a graph structure.

More formally, if we define flux as the quantity passing outward through a surface, then the Laplacian represents the density of the flux of the gradient flow of a given function. A step-by-step visualization of the Laplacian’s calculation is provided in Figure 9. Note that the definition of the Laplacian is dependant on three things: functions, the gradient of a function, and the divergence of the gradient. Since we are seeking to define the Laplacian in the graph domain, we need to define how these constructs operate in the graph domain.

Functions in the graph domain (referred to as graph signals in GSP) are a mapping from every vertex in a graph to a scalar value: \(f({\it G}\,({\it V, {\mbox{E}} \hspace{-0.85358pt}})): {{\mbox{V}}} \mapsto \hspace{1.42262pt}\mathbb {R}\) . Multiple graph functions can be defined over a given graph, and we can interpret a single graph function as a single feature vector defined over the vertices in a graph. See Figure 10 for an example of a graph with two graph functions.

The gradient of a function in the graph domain describes the the direction and the rate of fastest increase of graph signals. In a graph structure, when we refer to “direction,” we are referring to the edges of the graph; the avenues by which a graph function can change. For example, in Figure 10, the graph functions are 8-dimensional vectors (defined over the vertices), but the gradients of the functions for this graph are 12-dimensional vectors (defined over the edges) and are calculated as in Equations (7). Refer to Table 3 for a formal definition of the incident matrix M.

In Equations (7), the gradient vectors describe the difference in graph function value across the vertices/along the edges. Specifically, note that the largest magnitude value is 20,866 and corresponds to e \(_{12}\) , the edge between Hobart and Melbourne in Figure 10. In other words, the greatest magnitude edge is between the city with the least cases and the city with the most cases. Similarly, the lowest magnitude edge is e \(_{2}\) ; the edge between Perth and Adelaide, which has the least difference in cases.

The divergence of a gradient function in the graph domain describes the outward flux of the gradient function at every vertex. To continue with our example, we could interpret the divergence of the gradient function \(f_{cases}\) as the outgoing “flow” of infectious disease cases from each capital city. Whereas the gradient function was defined over the graph’s edges, the divergence of a gradient function is defined over the graph’s vertices and is calculated as in Equation (8).

The maximum value in the divergence vector for the infectious disease graph signal is 134,671, corresponding to Melbourne (the 7th vertex). Again, this can be interpreted as the magnitude of the “source” of infectious disease cases from Melbourne. Contrastively, the minimum value is \(-\) 281,123, corresponding to Canberra, the largest “sink” of infections disease.

Note as well that the dimensionality of the original graph function is 8—corresponding to the vector space, its gradient’s dimensionality is 12—corresponding to its edge space, and the Laplacian’s dimensionality is again 8—corresponding to the vertex space. This mimics the calculation of the Laplacian in Figure 9, where the original scalar field (representing the magnitude at each point) is converted to a vector field (representing direction) and then back to a scalar field (representing how each point acts as a source).

The graph Laplacian appears naturally in these calculations as a \(\vert {\mbox{V}}\vert \times \vert {\mbox{V}}\vert\) matrix operator in the form \({\bf L} ={\bf M} {\bf M} ^T\) (see Equation (8)). This corresponds to the formulation provided in Table 3, as shown in Equation (9), and this formulation is referred to as the combinatorial definition \({\bf L} = {\bf D}- {\bf A_{W}}\) (the normalized definition is defined as L \(\!^{sn}\) [17]). The graph Laplacian is pervasive in the fields of GSP [14].

Since \({\bf L} = {\bf D}- {\bf A_{W}}\) , the graph Laplacian must be a real ( \({\bf L} _{ij} \in \hspace{1.42262pt}\mathbb {R}, \hspace{2.84526pt} \forall \hspace{2.84526pt} 0 \le i,j \lt \vert {\mbox{V}}\vert\) ) and symmetric ( \({\bf L} = {\bf L} ^T\) ) matrix. As such, it will have an eigensystem composed of a set of \(\vert {\mbox{V}}\vert\) orthonormal eigenvectors, each associated with a single real eigenvalue [78]. We denote the ith eigenvector with \(u_i\) , and the associated eigenvalue with \(\lambda _i\) , each satisfying \({\bf L} u_i = \lambda _i u_i\) , where the eigenvectors \(u_i\) are the \(\vert {\mbox{V}}\vert\) -dimensional columns in the matrix (Fourier basis) U. The Laplacian can be factored as three matrices such that \({\bf L} = {\bf U} \mathbf {\Lambda } {\bf U} ^T\) through a process known as eigenvector decomposition. A variety of algorithms exist for solving this kind of eigendecomposition problem (e.g., the QR algorithm and Singular Value Decomposition).

These eigenvectors form a basis in \(\hspace{1.42262pt}\mathbb {R}^{\vert {\mbox{V}}\vert }\) , and as such, we can express any discrete graph function as a linear combination these eigenvectors. We define the graph Fourier transform of any graph function/signal as \(\hat{f} = {\bf U} ^T f \in \hspace{1.42262pt}\mathbb {R}^{\vert {\mbox{V}}\vert }\) , and its inverse as \(f = {\bf U} \hat{f} \in \hspace{1.42262pt}\mathbb {R}^{\vert {\mbox{V}}\vert }\) . To complete our goal of performing convolution in the spectral domain, we now complete the following steps:

This simple formulation has a number of downsides. Foremost, the approach is not localized—it has global support—meaning that the filter is applied to all vertices (i.e., the entirety of the graph function). This means that useful filters are not shared, and that the locality of graph structures is not being exploited. Second, it is not scalable: The number of learnable parameters grows with the size of the graph (not the scale of the filter) [17], the \(O(N^2)\) cost of matrix multiplication scales poorly to large graphs, and the \(O(N^3)\) time complexity of QR-based eigendecomposition [81] is prohibitive on large graphs. Moreover, directly computing this transform requires the diagonalization of the Laplacian and is infeasible for large graphs (where \(\vert {\mbox{V}}\vert\) exceeds a few thousand vertices) [31]. Finally, since the structure of the graph dictates the values of the Laplacian, graphs with dynamic topologies cannot use this method of convolution.

To alleviate the locality issue, Reference [21] noted that the smoothing of filters in the frequency space would result in localization in the vertex space. Instead of learning the filter directly, they formed the filter as a combination of smooth polynomial functions and instead learned the coefficients to these polynomials. Since the Laplacian is a local operator affecting only direct neighbors of any given vertex, then a polynomial of degree \(r\) affects vertices \(r\) -hops away. By approximating the spectral filter in this way (instead of directly learning it), spatial localization is thus guaranteed [77]. Furthermore, this improved scalability; learning \(K\) coefficients of the predefined smooth polynomial functions meant that the number of learnable parameters was no longer dependent on the size of the input graph. Additionally, the learned model could be applied to other graphs, too, as opposed to spectral filter coefficients that are basis-dependant. Since then, multiple potential polynomials have been used for specialized effects (e.g., Chebyshev polynomials, Cayley polynomials [51]).

Equation (11) outlines this approach. The learnable parameters are \(\theta _k\) —vectors of Chebyshev polynomial coefficients—and \(T_k(\widetilde{\mathbf {\Lambda }})\) is the Chevyshev polynomial of order \(k\) (dependent on the normalized diagonal matrix of scaled eigenvalues \(\widetilde{\mathbf {\Lambda }}\) ). Chevyshev polynomials can be computed recursively with a stable recurrence relation and form an orthogonal basis [83]. We recommend Reference [65] for a full treatment on Chebyshev polynomials.

Interestingly, these approximate approaches demonstrate an equivalence between spatial and spectral techniques. Both are spatially localized and allow for a single filter to extract repeating patterns of interest throughout a graph, both have a number of learnable parameters that are independent of the input graph size, and each have meaningful and intuitive interpretations from a spatial (Figure 7) and spectral (Figure 10) perspective. In fact, GCNs can be viewed as a first-order approximation of Chebyshev polynomials [56]. For an in-depth treatment on the topic of GSP, we recommend References [82] and [78].

As noted, spatial techniques and spectral techniques each have their advantages and disadvantages. The most popular spatial techniques are localized, generally scalable, and easily interpretable as methods for extracting features of interest from neighborhoods within graphs. As a downside, some of the more popular approaches (i.e., GATs) require expensive computations that scale in compute complexity quadratically with the size of their inputs, making them unsuitable for large graphs. However, spectral techniques can also be localized, scalable, and physically interpreted, but in some cases require rigorous computations to calculate the graph Laplacian. In general, eigendecomposition-based techniques cannot be used for graph inputs that have dynamic topologies and are computationally prohibitive for large graphs.

AEs work in a two-step fashion, first encoding the input data into a latent space and then decoding this compressed representation to reconstruct the original input data, as depicted in Figure 11 (though in some cases higher dimensionality latent space representations have been used). The AE is then trained to minimize the reconstruction loss, which is calculated using only the input data and can therefore be trained in an unsupervised manner. In its simplest form, such a loss is defined as \({\mbox{Loss}}_{AE}=\Vert X-\hat{X}\Vert ^{2}\) , where we have an input instance \(X\) and the reconstructed input \(\hat{X}\) .

The difference between AEs to GAEs is illustrated in Figure 11 and requires the definition of encoders and decoders that take in and put out graph structures, respectively. One of the most common methods for doing this is to replace the encoder with a CGNN and replace the decoder with a method that can reconstruct the graph structure of the input [46].

With a well-defined loss function, we can perform end-to-end learning across this network to optimize the encoding and decoding to strike a balance between both sensitivity to inputs and generalizability—we do not want the network to overfit and “memorize” all training inputs. Rather, the goal is for the encoder network to represent repeating patterns within the input data in a more compressed and efficient format.

Once trained, GAEs (like AEs), can be split into their component networks to perform specific tasks. A popular use case for the encoder is to generating robust embeddings for supervised downstream tasks (e.g., classification, visualization, regression, clustering), and a use for the decoder is to generate new graph instances that include properties from the original dataset. This allows the generation of large synthetic datasets.

Rather than representing inputs with single points in latent space, variational autoencoders (VAEs) learn to encode inputs as probability distributions in latent space. Figure 13 shows a VGAE that predicts a multivariate Guassian-like distribution \(q(Z|A,X)\) for a given input.

In this example, we will implement GAEs and VGAEs to perform unsupervised learning. After training, the learned embeddings can be used for both vertex classification and edge prediction tasks, even though the model was not trained to perform these tasks initially, thus demonstrating that the embeddings are meaningful representations of vertices. We will focus on edge prediction in citation networks, though these models can be easily applied to many task contexts and problem domains (e.g., vertex and graph classification).

To investigate GAEs, we use the Citeseer, Cora, and PubMed datasets, which are accessible via PyTorchGeometric [24]. In each of these graphs, the vertex features are word vectors indicating the presence or absence of predefined keywords (see Section 4.3 for another example of the Cora dataset being used).

We first implement a GAE. This model uses a single GCN to encode input features into a latent space. An inner product decoder is then applied to reconstruct the input from the latent space embedding (as described in Section 5.2). During training, we optimize the network by reducing the reconstruction loss. We then apply the model to an edge prediction task to test whether the embeddings can be used in performing downstream machine learning tasks and not just in reconstructing the inputs.

We then implement the VGAE as first described in Reference [46]. Unlike GAEs, there are now two GCNs in the encoder model (one each for the mean and variance of a probability distribution). The loss is also changed to Kullback–Leibler (KL) divergence to optimize for an accurate probability distribution. From here, we follow the same steps as for the GAE method: An inner product decoder is applied to the embeddings to perform input reconstruction. Again, we will test the efficacy of these learned embeddings on downstream machine learning tasks.

To test GAEs and VGAEs on each graph, we average the results from 10 experiments. In each experiment, the model is trained for 200 iterations to learn a 16-dimensional embedding.

In alignment with Reference [46], we see the VGAE outperform GAE on the Cora and Citeseer graphs, while the GAE outperforms the VGAE on the PubMed graph. Performance of both algorithms was significantly higher on the PubMed graph, likely owing to PubMed’s larger number of vertices ( \(\vert {\mbox{V}}\vert = 19,\!717\) ), and therefore more training examples, than Citeseer ( \(\vert {\mbox{V}}\vert = 3,\!327\) ) or Cora ( \(\vert {\mbox{V}}\vert = 2,\!708\) ). Moreover, while Cora and Citeseer vertex features are simple binary word vectors (of sizes 1,433 and 3,703, respectively), PubMed uses the more descriptive TF-IDF word vector, which accounts for the frequency of terms. This feature may be more discriminative and thus more useful when learning vertex embeddings.

This creates a “smoother” latent space that covers the full spectrum of inputs, rather than leaving “gaps,” where an unseen latent space vector would be decoded into a meaningless output. This has the effect of increasing generalization to unseen inputs and regularizing the model to avoid overfitting. Ultimately, this approach transforms the GAE into a more suitable generative model.

Unlike in GAEs—where the loss is simple the mean squared error between the input and the reconstructed input—a VGAE’s loss imposes an additional penalty that ensures that the latent distributions are normalized. More specifically, this term regularizes the latent space distributions by ensuring that they do not diverge significantly from some prior distribution with desirable properties. In our example, we use the normal distribution (denoted as \(N(0, 1)\) ). This divergence is quantified in our case using Kulback-Leibler divergence (KL), though other similarity metrics (e.g., Wassertein space distance or ranking loss) can be used successfully. Without this loss penalty, the VGAE encoder might generate distributions with small variances or high magnitude means; both of which would make it harder to sample from the distribution effectively.

Graph Adversarial Techniques (GAdvTs) use adversarial learning methods, whereby an AI model acts as an adversary to another during training to mutually improve the performance of both models in tandem. Due to the adversarial nature of GAdvT’s, developments in this area have been described as an “arms race between attackers and defenders” [12]. As with traditional adversarial techniques, common goals for GAdvTs include:

The field of GAdvTs is broad, with multiple different kinds of attacks, training regimes, and use cases. In this tutorial, we will look at how GAdvTs can be used to extend VGAEs to create robust encoding networks, and well regularized generative mechanisms. Figure 14 describes a typical architecture for adversarially training a VGAE.

To ensure that the sampling operation is differentiable, VGAEs leverage a “reparameterization trick,” where a random Guassian sample is generated from \(N(0,1)\) outside the forward pass of the network. The sample is then transformed by the parameterization of the generated distribution \(q(z)\) , rather than having the sample be generated directly from \(q(z)\) [19]. Since this approach is entirely differentiable, it allows for end-to-end training via backpropagation of an unsupervised loss signal.

In this section, we have explained the mechanics behind traditional and graph autoencoders. Analogous to how AEs use NN to perform encoding, GAEs use CGNNs to perform encoding and create embeddings [46]. Similarly, an unsupervised reconstruction error term is defined; in the case of GAEs, this is between the original adjacency matrix and the predicted adjacency matrix (produced by the decoder). GAEs and VGAEs represent a simple method for performing unsupervised training, which allows us to learn powerful embeddings in the graph domain without any labelled data, but requires regularization techniques to smooth their latent space representations and reparameterization tricks to ensure differentiability.

Before the seminal work in Reference [46] on VGAEs, a number of deep GAEs had been developed for unsupervised training on graph-structured data, including Deep Neural Graph Representations (DNGR) [8] and Structure Deep Network Embeddings (SDNE) [91]. These methods operate on only the adjacency matrix, so information about both the entire graph and the local neighborhoods is lost. More recent work mitigates this by using an encoder that aggregates information from a vertex’s local neighborhood to learn latent vector representations. For example, Reference [70] proposes a linear encoder that uses a single weights matrix to aggregate information from each vertex’s one-step local neighborhood, showing competitive performance on numerous benchmarks. Despite this, typical GAEs use more complex encoders—primarily CGNNs—to capture nonlinear relationships in the input data and larger local neighborhoods [5, 8, 46, 64, 84, 85, 91, 102].

Recent advancements in deep learning have allowed deeper NNs to be developed and applied throughout the field of AI. As the mechanics that drive predictions (and thus decisions) become more complex, the path by which those decisions are reached becomes more obfuscated. Explainable AI (XAI) promises to address this issue.

Explainability in the graph domain promises much of the same benefits as it does across AI, including more interpretable outputs, clearer relationships between inputs and outputs, more interpretable models, and in general, more trust between AI and human operators across problem domains (e.g., digital pathology [40], knowledge graphs [95]).

While the suite of available XAI algorithms has been consistently growing over the recent years (e.g., LIME, SHAP), graph specific XAI algorithms are relatively few and far between [103]. A key reason for this might be the requirement for graph explanations to incorporate not just the relationships among the input features, but also the relationships surrounding the input’s structural/topological information. In particular, the exploration of instance-level explainers—including high-fidelity perturbative and gradient-based methods—may provide good approximations of input importance in graph prediction tasks. Further techniques, especially those that assign quantitative importances to a graph’s structure and its features will give a more holistic view of explainability in the graph domain.

In traditional deep learning, a common technique for dealing with extremely large datasets and AI models is to distribute and parallelize computations where possible. In the graph domain, the non-rigid structure of graphs presents additional challenges. For example, how can a graph be uniformly partitioned across multiple devices? How can message-passing frameworks be efficiently implemented in a distributed system? These questions are especially pertinent for extremely large graphs. Recent developments suggest that sampling-based approaches may provide appropriate solutions in the near future [109], though such solutions are non-trivial, especially when graphs are stored on distributed systems [76].

Moreover, the scalability of GNN modules themselves may be improved by further directed research. For example, popular GNN variants such as MPNNs can in practice only be applied to small graphs due to the large computational overheads associated with the message-passing framework. Methods such as GATs show promising results regarding scalability, but attentional mechanisms still incur a quadratic time complexity, which may be prohibitive for graphs with large neighborhoods (on average). An exciting further avenue of research regarding GATs is their equivalence to Transformer networks [41, 47, 86, 101]. Further directed research in this area may contribute not only to the development of exciting new graph-based techniques, but also the understanding of Transformer networks as a whole. Breakthroughs in this area may address challenges specific to Transformers, such as the design of efficient positional encodings, effective warm-up strategies, and the quantification of inductive biases.

Self-supervised Learning (SSL) techniques have recently been suggested as the “next step” in AI training paradigms, as they close the gap—and even outperform—fully supervised approaches in visual tasks [2, 9]. Contemporary approaches to SSL include using models that learn from one another [9, 26, 104]. In related work, recent research suggests that contrastive objectives can be designed in the graph domain by selecting views of a single graph instance, thus permitting the capture of universal structural properties across graph without the need for large labelled datasets [42, 66, 111, 112]. These initial investigations demonstrate that the graph domain is well suited to the application of advanced learning paradigms techniques. Further research in this area may produce more general pretrained GNNs, allow the leveraging of large unlabelled graph datasets, and yield further insight into nature unsupervised/weakly supervised learning in the development of intelligence.