Understanding Convolutions on Graphs (2024)

This article is one of two Distill publications about graph neural networks. Take a look at A Gentle Introduction to Graph Neural Networks for a companion view on many things graph and neural network related.

Many systems and interactions - social networks, molecules, organizations, citations, physical models, transactions - can be represented quite naturally as graphs. How can we reason about and make predictions within these systems?

One idea is to look at tools that have worked well in other domains: neural networks have shown immense predictive power in a variety of learning tasks. However, neural networks have been traditionally used to operate on fixed-size and/or regular-structured inputs (such as sentences, images and video). This makes them unable to elegantly process graph-structured data.

Graph neural networks (GNNs) are a family of neural networks that can operate naturally on graph-structured data. By extracting and utilizing features from the underlying graph, GNNs can make more informed predictions about entities in these interactions, as compared to models that consider individual entities in isolation.

GNNs are not the only tools available to model graph-structured data: graph kernels and random-walk methods were some of the most popular ones. Today, however, GNNs have largely replaced these techniques because of their inherent flexibility to model the underlying systems better.

In this article, we will illustrate the challenges of computing over graphs, describe the origin and design of graph neural networks, and explore the most popular GNN variants in recent times. Particularly, we will see that many of these variants are composed of similar building blocks.

First, let’s discuss some of the complications that graphs come with.

The Challenges of Computation on Graphs

Lack of Consistent Structure

Graphs are extremely flexible mathematical models; but this means they lack consistent structure across instances. Consider the task of predicting whether a given chemical molecule is toxic :

Looking at a few examples, the following issues quickly become apparent:

• Molecules may have different numbers of atoms.
• The atoms in a molecule may be of different types.
• Each of these atoms may have different number of connections.
• These connections can have different strengths.

Representing graphs in a format that can be computed over is non-trivial, and the final representation chosen often depends significantly on the actual problem.

Node-Order Equivariance

Extending the point above: graphs often have no inherent ordering present amongst the nodes. Compare this to images, where every pixel is uniquely determined by its absolute position within the image!

As a result, we would like our algorithms to be node-order equivariant: they should not depend on the ordering of the nodes of the graph. If we permute the nodes in some way, the resulting representations of the nodes as computed by our algorithms should also be permuted in the same way.

Scalability

Graphs can be really large! Think about social networks like Facebook and Twitter, which have over a billion users. Operating on data this large is not easy.

Luckily, most naturally occuring graphs are ‘sparse’: they tend to have their number of edges linear in their number of vertices. We will see that this allows the use of clever methods to efficiently compute representations of nodes within the graph. Further, the methods that we look at here will have significantly fewer parameters in comparison to the size of the graphs they operate on.

Problem Setting and Notation

There are many useful problems that can be formulated over graphs:

• Node Classification: Classifying individual nodes.
• Graph Classification: Classifying entire graphs.
• Node Clustering: Grouping together similar nodes based on connectivity.
• Link Prediction: Predicting missing links.
• Influence Maximization: Identifying influential nodes.

A common precursor in solving many of these problems is node representation learning: learning to map individual nodes to fixed-size real-valued vectors (called ‘representations’ or ‘embeddings’).

In Learning GNN Parameters, we will see how the learnt embeddings can be used for these tasks.

Different GNN variants are distinguished by the way these representations are computed. Generally, however, GNNs compute node representations in an iterative process. We will use the notation $h_v^{(k)}$hv(k)​ to indicate the representation of node $v$v after the $k^{\text{th}}$kth iteration. Each iteration can be thought of as the equivalent of a ‘layer’ in standard neural networks.

We will define a graph $G$G as a set of nodes, $V$V, with a set of edges $E$E connecting them. Nodes can have individual features as part of the input: we will denote by $x_v$xv​ the individual feature for node $v \in V$v∈V. For example, the ‘node features’ for a pixel in a color image would be the red, green and blue channel (RGB) values at that pixel.

For ease of exposition, we will assume $G$G is undirected, and all nodes are of the same type. These kinds of graphs are called ‘hom*ogeneous’. Many of the same ideas we will see here apply to other kinds of graphs: we will discuss this later in Different Kinds of Graphs.

Sometimes we will need to denote a graph property by a matrix $M$M, where each row $M_v$Mv​ represents a property corresponding to a particular vertex $v$v.

Extending Convolutions to Graphs

Convolutional Neural Networks have been seen to be quite powerful in extracting features from images. However, images themselves can be seen as graphs with a very regular grid-like structure, where the individual pixels are nodes, and the RGB channel values at each pixel as the node features.

A natural idea, then, is to consider generalizing convolutions to arbitrary graphs. Recall, however, the challenges listed out in the previous section: in particular, ordinary convolutions are not node-order invariant, because they depend on the absolute positions of pixels. It is initially unclear as how to generalize convolutions over grids to convolutions over general graphs, where the neighbourhood structure differs from node to node. The curious reader may wonder if performing some sort of padding and ordering could be done to ensure the consistency of neighbourhood structure across nodes. This has been attempted with some success , but the techniques we will look at here are more general and powerful.

We begin by introducing the idea of constructing polynomial filters over node neighbourhoods, much like how CNNs compute localized filters over neighbouring pixels. Then, we will see how more recent approaches extend on this idea with more powerful mechanisms. Finally, we will discuss alternative methods that can use ‘global’ graph-level information for computing node representations.

Polynomial Filters on Graphs

The Graph Laplacian

Given a graph $G$G, let us fix an arbitrary ordering of the $n$n nodes of $G$G. We denote the $0-1$0−1 adjacency matrix of $G$G by $A$A, we can construct the diagonal degree matrix $D$D of $G$G as:

where $A_{vu}$Avu​ denotes the entry in the row corresponding to $v$v and the column corresponding to $u$u in the matrix $A$A. We will use this notation throughout this section.

Then, the graph Laplacian $L$L is the square $n \times n$n×n matrix defined as: $L = D - A.$L=D−A.

The graph Laplacian gets its name from being the discrete analog of the Laplacian operator from calculus.

Although it encodes precisely the same information as the adjacency matrix $A$A In the sense that given either of the matrices $A$A or $L$L, you can construct the other. , the graph Laplacian has many interesting properties of its own. The graph Laplacian shows up in many mathematical problems involving graphs: random walks, spectral clustering, and diffusion, to name a few. We will see some of these properties in a later section, but will instead point readers to this tutorial for greater insight into the graph Laplacian.

Polynomials of the Laplacian

Now that we have understood what the graph Laplacian is, we can build polynomials of the form: $p_w(L) = w_0 I_n + w_1 L + w_2 L^2 + \ldots + w_d L^d = \sum_{i = 0}^d w_i L^i.$pw​(L)=w0​In​+w1​L+w2​L2+…+wd​Ld=i=0∑d​wi​Li. Each polynomial of this form can alternately be represented by its vector of coefficients $w = [w_0, \ldots, w_d]$w=[w0​,…,wd​]. Note that for every $w$w, $p_w(L)$pw​(L) is an $n \times n$n×n matrix, just like $L$L.

These polynomials can be thought of as the equivalent of ‘filters’ in CNNs, and the coefficients $w$w as the weights of the ‘filters’.

For ease of exposition, we will focus on the case where nodes have one-dimensional features: each of the $x_v$xv​ for $v \in V$v∈V is just a real number. The same ideas hold when each of the $x_v$xv​ are higher-dimensional vectors, as well.

Using the previously chosen ordering of the nodes, we can stack all of the node features $x_v$xv​ to get a vector $x \in \mathbb{R}^n$x∈Rn.

Once we have constructed the feature vector $x$x, we can define its convolution with a polynomial filter $p_w$pw​ as: $x’ = p_w(L) \ x$x′=pw​(L)x To understand how the coefficients $w$w affect the convolution, let us begin by considering the ‘simplest’ polynomial: when $w_0 = 1$w0​=1 and all of the other coefficients are $0$0. In this case, $x’$x′ is just $x$x: $x’ = p_w(L) \ x = \sum_{i = 0}^d w_i L^ix = w_0 I_n x = x.$x′=pw​(L)x=i=0∑d​wi​Lix=w0​In​x=x. Now, if we increase the degree, and consider the case where instead $w_1 = 1$w1​=1 and and all of the other coefficients are $0$0. Then, $x’$x′ is just $Lx$Lx, and so: $\begin{aligned} x’_v = (Lx)_v &= L_v x \\ &= \sum_{u \in G} L_{vu} x_u \\ &= \sum_{u \in G} (D_{vu} - A_{vu}) x_u \\ &= D_v \ x_v - \sum_{u \in \mathcal{N}(v)} x_u \end{aligned}$xv′​=(Lx)v​​=Lv​x=u∈G∑​Lvu​xu​=u∈G∑​(Dvu​−Avu​)xu​=Dv​xv​−u∈N(v)∑​xu​​ We see that the features at each node $v$v are combined with the features of its immediate neighbours $u \in \mathcal{N}(v)$u∈N(v). For readers familiar with Laplacian filtering of images, this is the exact same idea. When $x$x is an image, $x’ = Lx$x′=Lx is exactly the result of applying a ‘Laplacian filter’ to $x$x.

At this point, a natural question to ask is: How does the degree $d$d of the polynomial influence the behaviour of the convolution? Indeed, it is not too hard to show that: This is Lemma 5.2 from . $\text{dist}_G(v, u) > i \quad \Longrightarrow \quad L_{vu}^i = 0.$distG​(v,u)>i⟹Lvui​=0. This implies, when we convolve $x$x with $p_w(L)$pw​(L) of degree $d$d to get $x’$x′: $\begin{aligned} x’_v = (p_w(L)x)_v &= (p_w(L))_v x \\ &= \sum_{i = 0}^d w_i L_v^i x \\ &= \sum_{i = 0}^d w_i \sum_{u \in G} L_{vu}^i x_u \\ &= \sum_{i = 0}^d w_i \sum_{u \in G \atop \text{dist}_G(v, u) \leq i} L_{vu}^i x_u. \end{aligned}$xv′​=(pw​(L)x)v​​=(pw​(L))v​x=i=0∑d​wi​Lvi​x=i=0∑d​wi​u∈G∑​Lvui​xu​=i=0∑d​wi​distG​(v,u)≤iu∈G​∑​Lvui​xu​.​

Effectively, the convolution at node $v$v occurs only with nodes $u$u which are not more than $d$d hops away. Thus, these polynomial filters are localized. The degree of the localization is governed completely by $d$d.

To help you understand these ‘polynomial-based’ convolutions better, we have created the visualization below. Vary the polynomial coefficients and the input grid $x$x to see how the result $x’$x′ of the convolution changes. The grid under the arrow shows the equivalent convolutional kernel applied at the highlighted pixel in $x$x to get the resulting pixel in $x’$x′. The kernel corresponds to the row of $p_w(L)$pw​(L) for the highlighted pixel. Note that even after adjusting for position, this kernel is different for different pixels, depending on their position within the grid.

ChebNet

What is the motivation behind these choices?

• $L$L is actually positive semi-definite: all of the eigenvalues of $L$L are not lesser than $0$0. If $\lambda_{\max}(L) > 1$λmax​(L)>1, the entries in the powers of $L$L rapidly increase in size. $\tilde{L}$L~ is effectively a scaled-down version of $L$L, with eigenvalues guaranteed to be in the range $[-1, 1]$[−1,1]. This prevents the entries of powers of $\tilde{L}$L~ from blowing up. Indeed, in the visualization above: we restrict the higher-order coefficients when the unnormalized Laplacian $L$L is selected, but allow larger values when the normalized Laplacian $\tilde{L}$L~ is selected, in order to show the result $x’$x′ on the same color scale.
• The Chebyshev polynomials have certain interesting properties that make interpolation more numerically stable. We won’t talk about this in more depth here, but will advise interested readers to take a look at as a definitive resource.

Polynomial Filters are Node-Order Equivariant

The polynomial filters we considered here are actually independent of the ordering of the nodes. This is particularly easy to see when the degree of the polynomial $p_w$pw​ is $1$1: where each node’s feature is aggregated with the sum of its neighbour’s features. Clearly, this sum does not depend on the order of the neighbours. A similar proof follows for higher degree polynomials: the entries in the powers of $L$L are equivariant to the ordering of the nodes.

Embedding Computation

We now describe how we can build a graph neural network by stacking ChebNet (or any polynomial filter) layers one after the other with non-linearities, much like a standard CNN. In particular, if we have $K$K different polynomial filter layers, the $k^{\text{th}}$kth of which has its own learnable weights $w^{(k)}$w(k), we would perform the following computation:

Note that these networks reuse the same filter weights across different nodes, exactly mimicking weight-sharing in Convolutional Neural Networks (CNNs) which reuse weights for convolutional filters across a grid.

Modern Graph Neural Networks

ChebNet was a breakthrough in learning localized filters over graphs, and it motivated many to think of graph convolutions from a different perspective.

As we noted before, this is a $1$1-hop localized convolution. But more importantly, we can think of this convolution as arising of two steps:

• Aggregating over immediate neighbour features $x_u$xu​.
• Combining with the node’s own feature $x_v$xv​.

Key Idea: What if we consider different kinds of ‘aggregation’ and ‘combination’ steps, beyond what are possible using polynomial filters?

By ensuring that the aggregation is node-order equivariant, the overall convolution becomes node-order equivariant.

These convolutions can be thought of as ‘message-passing’ between adjacent nodes: after each step, every node receives some ‘information’ from its neighbours.

By iteratively repeating the $1$1-hop localized convolutions $K$K times (i.e., repeatedly ‘passing messages’), the receptive field of the convolution effectively includes all nodes upto $K$K hops away.

Embedding Computation

Message-passing forms the backbone of many GNN architectures today. We describe the most popular ones in depth below:

• Graph Convolutional Networks (GCN)
• Graph Attention Networks (GAT)
• Graph Sample and Aggregate (GraphSAGE)
• Graph Isomorphism Network (GIN)

Thoughts

An interesting point is to assess different aggregation functions: are some better and others worse? demonstrates that aggregation functions indeed can be compared on how well they can uniquely preserve node neighbourhood features; we recommend the interested reader take a look at the detailed theoretical analysis there.

Here, we’ve talk about GNNs where the computation only occurs at the nodes. More recent GNN models such as Message-Passing Neural Networks and Graph Networks perform computation over the edges as well; they compute edge embeddings together with node embeddings. This is an even more general framework - but the same ‘message passing’ ideas from this section apply.

Interactive Graph Neural Networks

Below is an interactive visualization of these GNN models on small graphs. For clarity, the node features are just real numbers here, shown inside the squares next to each node, but the same equations hold when the node features are vectors.

In practice, each iteration above is generally thought of as a single ‘neural network layer’. This ideology is followed by many popular Graph Neural Network libraries, For example: PyTorch Geometric and StellarGraph. allowing one to compose different types of graph convolutions in the same model.

From Local to Global Convolutions

The methods we’ve seen so far perform ‘local’ convolutions: every node’s feature is updated using a function of its local neighbours’ features.

While performing enough steps of message-passing will eventually ensure that information from all nodes in the graph is passed, one may wonder if there are more direct ways to perform ‘global’ convolutions.

The answer is yes; we will now describe an approach that was actually first put forward in the context of neural networks by , much before any of the GNN models we looked at above.

Spectral Convolutions

As before, we will focus on the case where nodes have one-dimensional features. After choosing an arbitrary node-order, we can stack all of the node features to get a ‘feature vector’ $x \in \mathbb{R}^n$x∈Rn.

Key Idea: Given a feature vector $x$x, the Laplacian $L$L allows us to quantify how smooth $x$x is, with respect to $G$G.

How?

After normalizing $x$x such that $\sum_{i = 1}^n x_i^2 = 1$∑i=1n​xi2​=1, if we look at the following quantity involving $L$L: $R_L$RL​ is formally called the Rayleigh quotient. $R_L(x) = \frac{x^T L x}{x^T x} = \frac{\sum_{(i, j) \in E} (x_i - x_j)^2}{\sum_i x_i^2} = \sum_{(i, j) \in E} (x_i - x_j)^2.$RL​(x)=xTxxTLx​=∑i​xi2​∑(i,j)∈E​(xi​−xj​)2​=(i,j)∈E∑​(xi​−xj​)2. we immediately see that feature vectors $x$x that assign similar values to adjacent nodes in $G$G (hence, are smooth) would have smaller values of $R_L(x)$RL​(x).

$L$L is a real, symmetric matrix, which means it has all real eigenvalues $\lambda_1 \leq \ldots \leq \lambda_{n}$λ1​≤…≤λn​. An eigenvalue $\lambda$λ of a matrix $A$A is a value satisfying the equation $Au = \lambda u$Au=λu for a certain vector $u$u, called an eigenvector. For a friendly introduction to eigenvectors, please see this tutorial. Further, the corresponding eigenvectors $u_1, \ldots, u_{n}$u1​,…,un​ can be taken to be orthonormal: $u_{k_1}^T u_{k_2} = \begin{cases} 1 \quad \text{ if } {k_1} = {k_2}. \\ 0 \quad \text{ if } {k_1} \neq {k_2}. \end{cases}$uk1​T​uk2​​={1ifk1​=k2​.0ifk1​≠k2​.​ It turns out that these eigenvectors of $L$L are successively less smooth, as $R_L$RL​ indicates: This is the min-max theorem for eigenvalues. $\underset{x, \ x \perp \{u_1, \ldots, u_{i - 1}\}}{\text{argmin}} R_L(x) = u_i. \qquad \qquad \qquad \min_{x, \ x \perp \{u_1, \ldots, u_{i - 1}\}} R_L(x) = \lambda_i.$x,x⊥{u1​,…,ui−1​}argmin​RL​(x)=ui​.x,x⊥{u1​,…,ui−1​}min​RL​(x)=λi​. The set of eigenvalues of $L$L are called its ‘spectrum’, hence the name! We denote the ‘spectral’ decomposition of $L$L as: $L = U \Lambda U^T.$L=UΛUT. where $\Lambda$Λ is the diagonal matrix of sorted eigenvalues, and $U$U denotes the matrix of the eigenvectors (sorted corresponding to increasing eigenvalues): $\Lambda = \begin{bmatrix} \lambda_{1} & & \\ & \ddots & \\ & & \lambda_{n} \end{bmatrix} \qquad \qquad \qquad \qquad U = \begin{bmatrix} \\ u_1 \ \cdots \ u_n \\ \end{bmatrix}.$Λ=⎣⎡​λ1​​⋱​λn​​⎦⎤​U=⎣⎡​u1​⋯un​​⎦⎤​. The orthonormality condition between eigenvectors gives us that $U^T U = I$UTU=I, the identity matrix. As these $n$n eigenvectors form a basis for $\mathbb{R}^n$Rn, any feature vector $x$x can be represented as a linear combination of these eigenvectors: $x = \sum_{i = 1}^n \hat{x_i} u_i = U \hat{x}.$x=i=1∑n​xi​^​ui​=Ux^. where $\hat{x}$x^ is the vector of coefficients $[x_0, \ldots x_n]$[x0​,…xn​]. We call $\hat{x}$x^ as the spectral representation of the feature vector $x$x. The orthonormality condition allows us to state: $x = U \hat{x} \quad \Longleftrightarrow \quad U^T x = \hat{x}.$x=Ux^⟺UTx=x^. This pair of equations allows us to interconvert between the ‘natural’ representation $x$x and the ‘spectral’ representation $\hat{x}$x^ for any vector $x \in \mathbb{R}^n$x∈Rn.

Spectral Representations of Natural Images

As discussed before, we can consider any image as a grid graph, where each pixel is a node, connected by edges to adjacent pixels. Thus, a pixel can have either $3, 5,$3,5, or $8$8 neighbours, depending on its location within the image grid. Each pixel gets a value as part of the image. If the image is grayscale, each value will be a single real number indicating how dark the pixel is. If the image is colored, each value will be a $3$3-dimensional vector, indicating the values for the red, green and blue (RGB) channels. We use the alpha channel as well in the visualization below, so this is actually RGBA.

This construction allows us to compute the graph Laplacian and the eigenvector matrix $U$U. Given an image, we can then investigate what its spectral representation looks like.

To shed some light on what the spectral representation actually encodes, we perform the following experiment over each channel of the image independently:

• We first collect all pixel values across a channel into a feature vector $x$x.
• Then, we obtain its spectral representation $\hat{x}$x^. $\hat{x} = U^T x$x^=UTx
• We truncate this to the first $m$m components to get $\hat{x}_m$x^m​. By truncation, we mean zeroing out all of the remaining $n - m$n−m components of $\hat{x}$x^. This truncation is equivalent to using only the first $m$m eigenvectors to compute the spectral representation. $\hat{x}_m = \text{Truncate}_m(\hat{x})$x^m​=Truncatem​(x^)
• Then, we convert this truncated representation $\hat{x}_m$x^m​ back to the natural basis to get $x_m$xm​. $x_m = U \hat{x}_m$xm​=Ux^m​

Finally, we stack the resulting channels back together to get back an image. We can now see how the resulting image changes with choices of $m$m. Note that when $m = n$m=n, the resulting image is identical to the original image, as we can reconstruct each channel exactly.

As $m$m decreases, we see that the output image $x_m$xm​ gets blurrier. If we decrease $m$m to $1$1, the output image $x_m$xm​ is entirely the same color throughout. We see that we do not need to keep all $n$n components; we can retain a lot of the information in the image with significantly fewer components. We can relate this to the Fourier decomposition of images: the more eigenvectors we use, the higher frequencies we can represent on the grid.

To complement the visualization above, we additionally visualize the first few eigenvectors on a smaller $8 \times 8$8×8 grid below. We change the coefficients of the first $10$10 out of $64$64 eigenvectors in the spectral representation and see how the resulting image changes:

These visualizations should convince you that the first eigenvectors are indeed smooth, and the smoothness correspondingly decreases as we consider later eigenvectors.

For any image $x$x, we can think of the initial entries of the spectral representation $\hat{x}$x^ as capturing ‘global’ image-wide trends, which are the low-frequency components, while the later entries as capturing ‘local’ details, which are the high-frequency components.

Embedding Computation

We now have the background to understand spectral convolutions and how they can be used to compute embeddings/feature representations of nodes.

As before, the model we describe below has $K$K layers: each layer $k$k has learnable parameters $\hat{w}^{(k)}$w^(k), called the ‘filter weights’. These weights will be convolved with the spectral representations of the node features. As a result, the number of weights needed in each layer is equal to $m$m, the number of eigenvectors used to compute the spectral representations. We had shown in the previous section that we can take $m \ll n$m≪n and still not lose out on significant amounts of information.

Thus, convolution in the spectral domain enables the use of significantly fewer parameters than just direct convolution in the natural domain. Further, by virtue of the smoothness of the Laplacian eigenvectors across the graph, using spectral representations automatically enforces an inductive bias for neighbouring nodes to get similar representations.

Assuming one-dimensional node features for now, the output of each layer is a vector of node representations $h^{(k)}$h(k), where each node’s representation corresponds to a row of the vector.

We fix an ordering of the nodes in $G$G. This gives us the adjacency matrix $A$A and the graph Laplacian $L$L, allowing us to compute $U_m$Um​. Finally, we can describe the computation that the layers perform, one after the other:

The method above generalizes easily to the case where each $h^{(k)} \in \mathbb{R}^{d_k}$h(k)∈Rdk​, as well: see for details.

With the insights from the previous section, we see that convolution in the spectral-domain of graphs can be thought of as the generalization of convolution in the frequency-domain of images.

Spectral Convolutions are Node-Order Equivariant

We can show spectral convolutions are node-order equivariant using a similar approach as for Laplacian polynomial filters.