Theory of Graph Neural Networks: Representation and Learning

Introduction

Solving these tasks demands a sufficiently rich embedding of the graph or of each node that captures structural properties as well as the attribute information. While graph embeddings have been a widely studied topic, including spectral embeddings and graph kernels, recently, Graph Neural Networks (GNNs) have emerged as an empirically broadly successful model class that, as opposed to, e.g., spectral embeddings, allows to adapt the embedding to the task at hand, generalizes to other graphs of the same input type, and incorporates attributes. Due to space limits, this survey focuses on the popular message passing (spatial) GNNs, formally defined below, and a selection of their rich mathematical connections, with an excursion into higher-order GNNs.

When analyzing learning and risk, three main questions become important:

1. Representational power (Section 2). Which target functions $g$ can be approximated well by a GNN model class $\mathcal{F}$? Answers to this question relate to graph isomorphism testing, approximation theory for neural networks, local algorithms and representing invariance/equivariance under permutations.

2. Generalization (Section 3). Even with sufficient approximation power, we can only estimate a function $\widehat{F}\in\mathcal{F}$ from $\mathcal{D}$: since $\mathcal{P}$ and hence $\mathcal{R}(F)$ is not accessible, the common learning or training procedure is to minimize the empirical risk $\widehat{\mathcal{R}}(F)$:

Generalization asks how well $\widehat{F}$ is performing according to the population risk, i.e., $\mathcal{R}(\widehat{F})$, as a function of $N$ and model properties. Good generalization may demand explicit (e.g., via penalties) or implicit regularization (e.g., via the optimization algorithm, typically variants of stochastic gradient descent). Hence, generalization analyses involve the complexity of the model class $\mathcal{F}$, the target function, the data and the optimization procedure.

3. Generalization under distribution shifts (Section 4). In practice, a learned model $\widehat{F}$ is often deployed on data from a distribution $\mathcal{Q}\neq\mathcal{P}$, e.g., graphs of different size, degree or attribute ranges, for instance, $\operatorname{supp}(Q)\supset\operatorname{supp}(P)$. In which cases can we expect successful extrapolation to $\mathcal{Q}$? This depends on the structure of the graphs and the task, formalizable via graph limits, local structures and algorithmic structures, e.g., dynamic programming.

Beyond these topics, GNNs have close connections to graph signal processing as learnable filters, geometric learning, and probabilistic inference. The first and third connection are covered in .

Lastly, a disclaimer: GNNs are a rapidly evolving research area. Hence, it is almost impossible to include all possible works, and any survey necessarily misses some. The author apologizes in advance for any references that were not covered.

The final node representation $f(v)=h_{v}^{(T)},\forall v\in V$ is the last iterate, possibly concatenated with a linear classifier. Throughout this paper, $\mathcal{N}(v)\subset V$ denotes the neighborhood of $v\in V$, and $\{\!\!\{\cdot\}\!\!\}$ a multiset. Here, $h^{(t)}_{v}$ encodes the $t$-hop neighborhood of node $v$, i.e., the subgraph of all nodes reachable from $v$ within $t$ steps. The number of iterations $T$ is also termed the GNN depth, and one iteration may be viewed as a layer.

The sum may also be replaced by an average, degree-normalized sum or coordinate-wise min or max. In the most general form, the functions $\phi_{1},\phi_{2}$ are implemented as multi-layer perceptrons (MLPs), neural networks that alternate linear transformations and coordinate-wise nonlinear activations such as the ReLU ($\sigma(a)=\max\{a,0\}$) or sigmoid function ($\sigma(a)=(1+\exp(-a))^{-1}$):

The learnable parameters $\theta$ of the MLP are the weight matrices $W^{(j)}$ and bias vectors $b^{(j)}$. The update $f_{\text{Up}}$ in Equation (1.5) is typically a weighted combination with learnable weight matrices:

Finally, if a graph-level prediction is desired, all node representations can be aggregated by a permutation invariant readout function

Here, we assume the readout has the form (1.6) or is a simple sum or average. Typically, all parameters are learned jointly via stochastic gradient descent minimizing the empirical risk.

Throughout this article, $n=|V|$ denotes the number of nodes and $N$ the number of training data points.

Spectral GNNs. Besides message passing GNNs, other architectures have been devised. One important example are spectral graph neural networks , which learn a function of the graph Laplacian $L$, i.e., $F(G)=p(L)X=Vp(\Lambda)V^{\top}X$, where $V$ and $\Lambda$ are the matrices of eigenvectors and eigenvalues (diagonal matrix), respectively, and $p$ is a polynomial. In the sequel, we will focus on message passing GNNs.

Representational power of GNNs

For functions on graphs, representational power has mainly been studied in terms of graph isomorphism: which graphs a GNN can distinguish. Via variations of the Stone-Weierstrass theorem, these results yield universal approximation results. Other works bound the ability of GNNs to compute specific polynomials of the adjacency matrix and to distinguish graphons . Observed limitations of MPNNs have inspired higher-order GNNs (Section 2.3). Moreover, if all node attributes are unique, then analogies to local algorithms yield algorithmic approximation results and lower bounds (Section 2.2).

A standard characterization of the discriminative power of GNNs is via the hierarchy of the Weisfeiler-Leman (WL) algorithm for graph isomorphism testing, also known as color refinement or vertex classification , which was inspired by the work of Weisfeiler and Leman . The WL algorithm does not entirely solve the graph ismomorphism problem, but its power has been widely studied.

A labeled graph is a graph endowed with a node coloring $l:V(G)\to\Sigma$ for some sufficiently large alphabet $\Sigma$. Given a labeled graph $(G,l)$, the 1-dimensional WL algorithm (1-WL) iteratively computes a node coloring $c_{l}^{(t)}:V(G)\to\Sigma$. Starting with $c_{l}^{(0)}$ in iteration $t=0$, in iteration $t>0$ it sets for all $v\in V$

where Hash is an injective map from the input pair to $\Sigma$, i.e., it assigns a unique color to each neighborhood pattern. To compare two graphs $G,G^{\prime}$, the algorithm compares the multisets $\{\!\!\{c_{l}^{(t)}(v)\mid v\in V(G)\}\!\!\}$ and $\{\!\!\{c_{l}^{(t)}(u)\mid u\in V(G^{\prime})\}\!\!\}$ in each iteration. If the sets differ, then it determines that $G\neq G^{\prime}$. Otherwise, it terminates when the number of colors in iteration $t$ and $t-1$ are the same, which occurs after at most $\max\{|V(G)|,|V(G^{\prime})|\}$ iterations.

The computational analogy between the 1-WL algorithm and MPNNs is obvious. Since the WL algorithm uniquely colors each neighborhood, the coloring $c_{l}^{(t)}(v)$ always refines the coloring $h^{(t)}_{v}$ from a GNN.

If for two graphs $G,G^{\prime}$ a message passing GNN outputs $f_{G}(G)\neq f_{G}(G^{\prime})$, then the 1-WL algorithm will determine that $G\neq G^{\prime}$. For any $t$, there exists an MPGNN such that $c_{l}^{(t)}\equiv h^{(t)}$. A sufficient condition is that the aggregate, update and readout operations are injective multiset functions.

GNNs that use the degree for normalization in the aggregation can be equivalent to the 1-WL agorithm too, but with one more iteration in the WL algorithm .

Theorem 1 demands the neighbor aggregation $f_{\text{Agg}}$ to be an injective multiset function on sets $S$ ($|S|\leq M$). Theorem 2 shows how to universally approximate multiset functions.

The proof idea is to show that there exists an injective function of the form $\sum_{s\in S}\phi(s)$. The above result is an extension of a universal approximation result for set functions , and suggests a neural network model for sets where $\phi_{1},\phi_{2}$ are approximated by MLPs. The Graph Isomorphism Network (GIN) implements this sum decomposition in the aggregation function to ensure the ability of injective operations.

1.2 Implications for graph distinction

Theorem 1 allows to directly transfer any known result for the 1-WL test to MPNNs. For instance, the 1-WL test succeeds to distinguish graphs sampled uniformly from all graphs on $n$ nodes with high probability, and failure probability going to zero as $n\to\infty$ . 1-WL can also distinguish any non-isomorphic pair of trees . It fails for regular graphs, as all node colors will be the same. The graphs that the 1-WL algorithm can distinguish from any non-isomorphic graph can be recognized in quasi-linear time . See also for more detailed results on the expressive power of variants of the WL algorithm.

1.3 Computation trees and structural graph properties

To further illustrate the implications of GNNs’ discriminative power, we look at some specific examples. The maximum information contained in any embedding $h_{v}^{(t)}$ can be characterized by a computation tree $\mathcal{T}(h^{(t)}_{v})$, i.e., an “unrolling” of the message passing procedure. The 1-WL test essentially colors computation trees. The tree $\mathcal{T}(h^{(t)}_{v})$ is constructed recursively: let $\mathcal{T}(h^{(0)}_{v})=x_{v}$ for all $v\in V$. For $t>0$, construct a root with label $x_{v}$ and, for any $u\in\mathcal{N}(v)$ construct a child subtree $\mathcal{T}(h^{(t-1)}_{u})$. Figure 1 illustrates an example.

If for two nodes $u\neq v$, we have $\mathcal{T}(h^{(t)}_{v})=\mathcal{T}(h^{(t)}_{u})$, then $h^{(t)}_{v}=h^{(t)}_{u}$.

Comparing computation trees directly implies that MPNNs cannot distinguish regular graphs. It also shows further limitations with practical impact, as indicated in Figure 1, in particular for learning combinatorial algorithms and for predicting properties of molecules, where functional groups are of key importance. We say a class of models $\mathcal{F}$ decides a graph property if there exists an $F\in\mathcal{F}$ such that for any two $G,G^{\prime}$ that differ in the property, we obtain $F(G)\neq F(G^{\prime})$.

MPNNs cannot decide girth, circumference, diameter, radius, existence of a conjoint cycle, total number of cycles, and existence of a $k$-clique . MPNNs cannot count induced (attributed) subgraphs for any connected pattern of 3 or more nodes, except star-shaped patterns .

Motivated by these limitations, generalizations of GNNs were proposed that provably increase their representational power. Two main directions are to (1) introduce node IDs (Section 2.2), and (2) use higher-order functions that act on tuples of nodes (Section 2.3).

2 Node IDs, local algorithms, combinatorial optimization and lower bounds

The major weaknesses of MPNNs arise from their inability to identify nodes as the origin of specific messages. Hence, MPNNs can be strengthened by making nodes more distinguishable. The gained representational power follows from connections with local algorithms, where the input graph defines both the computational problem and the network topology of a distributed system: there, each node $v\in V$ is a local machine and generates a local output, and all nodes execute the same algorithm, without faults.

Approximation Algorithms. Sato et al. achieve a partial node distinction by transferring the idea of a port numbering from local algorithms. Edges incident to each node are numbered as outgoing ports. In each round, each node simultaneously sends a message to each port, but the messages can differ across ports:

Permutation invariance, though, is not immediate. This corresponds to the vector-vector consistent (VVC) model for local algorithms . The VVC analogy allows to transfer results on representing approximation algorithms. CPNGNN is a specific VVC GNN model.

There exists a CPNGNN that can compute a $(\deg(G)+1)$-approximation for the minimum dominating set problem, a CPNGNN that can compute a 2-approximation for the minimum vertex cover problem, but no CPNGNN can do better. No CPNGNN can compute a constant-factor approximation for the maximum matching problem.

Adding a weak vertex 2-coloring leads to further results. Despite the increased power compared to MPNN, CPNGNNs retain most limitations of Proposition 2 .

A more powerful alternative is to endow nodes with fully unique identifiers . E.g., augmenting the GIN model (the most expressive MPNN) with random node identifiers yields a model that can decide subgraphs that MPNN and CPNGNN cannot . This model can further achieve better approximation results for minimum dominating set ($H(\deg(G)+1)+\epsilon$), where $H$ is the harmonic number) and maximum matching ($1+\epsilon$).

Turing completeness. Analogies to local algorithms further imply that MPNNs with unique node IDs are Turing complete, i.e., they can compute any function that a Turing machine can compute, including graph isomorphism. In particular, the proof shows an equivalence to the Turing universal LOCAL model from distributed computing .

If $f_{\text{Up}}$ and $f_{\text{Agg}}$ are Turing complete functions and the message passing GNN gets unique node IDs, then the classes GNN and LOCAL are equivalent. For any MPNN $F$ there exists a local algorithm $\mathcal{A}$ of the same depth, such that $F(G)=\mathcal{A}(G)$, and vice versa.

Lower bounds. For bounded size, GNNs lose computational power. Via analogies to the CONGEST model , which bounds message sizes, one can transfer results on decision, optimization and estimation problems on graphs. These lead to lower bounds on the product of depth and width of the GNN if the nodes do not have access to a random generator. Here, as before, width of a GNN refers to the dimensionality of the embeddings $h^{(t)}_{v}$.

If a problem cannot be solved in less than $T$ rounds in CONGEST using messages of at most $b$ bits, then it cannot be solved by an MPNN of width $w\leq(b-\log_{2}n)/p=O(b/\log n)$ and depth $T$, where $p=\Theta(n)$.

Theorem 5 directly implies lower bounds for solving combinatorial problems, e.g., $Tw=\Omega(n/\log n)$ for cycle detection and computing diameter, and $T\sqrt{w}=\Omega(\sqrt{n}/\log n)$ for minimum spanning tree, minimum cut and shortest path .

Random node IDs. While unique node IDs are powerful in theory, in many practical examples the input graphs do not have unique IDs. An alternative is to assign random node IDs . This can still yield GNNs that are essentially permutation invariant: while their outputs are random, the outputs for different graphs are still sufficiently separated . This leads to a probabilistic universal approximation result:

The proof builds on a result by that states that any logical sentence in FOC2 can be expressed by the addressed GNN. The logic considered here is a fragment of first-order (FO) predicate logic that allows to incorporate counting quantifiers of the form $\exists^{\geq k}x\psi(x)$, i.e., there are at least $k$ elements $x$ satisfying $\psi$, but is restricted to two variables. FOC2 is tightly linked with the 1-WL test: for any nodes $u,v\in V$ in any graph, 1-WL colors $u$ and $v$ the same if and only if they are classified the same by all FOC2 classifiers .

Another successful idea is to augment node attribute vectors with attributes that contain structural or topological information .

3 Higher-order GNNs

Instead of adding unique node IDs, one may increase the expressive power of GNNs by encoding subsets of $V$ that are larger than the single nodes used in MPNNs. Three such directions are: (1) neural network versions of higher-dimensional WL algorithms, (2) (non)linear equivariant operations, and (3) recursion. Other strategies that could not be covered here use, e.g., simplicial and cell complexes .

Extending analogies between MPNNs and the 1-WL algorithm , the first class of higher-order GNNs imitates versions of the $k$-dimensional WL algorithm. The $k$-WL algorithms are defined on $k$-tuples of nodes, and different versions differ in their aggregation and definition of neighborhood. In iteration , the $k$-WL algorithm labels each $k$-tuple $s\in V^{k}$ by a unique ID for its isomorphism type. Then it aggregates over neighborhoods $\mathcal{N}^{\text{WL}}_{i}(s)=\{(s_{1},s_{2},\ldots s_{i-1},v,s_{i+1},\ldots s_{k})\mid\forall v\in V\}$ for $1\leq i\leq k$:

For two graphs $G,G^{\prime}$ the $k$-WL algorithm then decides “not isomorphic” if $\{\!\!\{c^{(t)}(s)\mid s\in V(G)^{k}\}\!\!\}\neq\{\!\!\{c^{(t)}(s^{\prime})\mid s^{\prime}\in V(G^{\prime})^{k}\}\!\!\}$ for some $t$, and returns “maybe isomorphic” otherwise. Like the 1-WL test, the $k$-WL test decides “not isomorphic” only if $G\ncong G^{\prime}$. The Folklore $k$-WL algorithm ($k$-FWL) differs in its update rule, which “swaps” the order of the aggregation steps :

The 1-WL and 2-WL tests are equivalent, and for $k\geq 2$, the $(k+1)$-WL test can distinguish strictly more graphs than the $k$-WL test . The $k$-FWL algorithm is as powerful as the $(k+1)$-WL algorithm for $k\geq 2$ .

Set-WL GNN. Since computations on $k$-tuples are expensive, Morris et al. consider a GNN that corresponds to a set version of a $k$-WL algorithm. For any set $S\subseteq V$ with $|S|=k$, let $\mathcal{N}^{\text{set}}(S)=\{T\subset V,|T|=k\mid|S\cap T|=k-1\}$. The set-based WL ($k$-SWL) algorithm then updates as

its GNN analogue uses the aggregation and update (cf. Eqns. (1.6),(1.8))

where $\sigma$ is a coordinatewise nonlinearity (e.g., sigmoid or ReLU). This family of GNNs is equivalent in power to the $k$-SWL test (Theorem 8). For computational efficiency, a local version restricts the neighborhood of $S$ to sets $T$ such that the nodes $\{u,v\}=S\Delta T$ in the symmetric difference are connected in the graph. This local version is weaker .

Folklore WL GNN. In analogy to the $k$-FWL algorithm, Maron et al. define $k$-FGNNs with aggregations

The family of $k$-FGNNs is a class of nonlinear equivariant networks, and is equivalent in power to the $k$-FWL test and the $(k+1)$-WL test (Theorem 8).

3.2 Linear equivariant layers.

The associated GNN model uses one parameter (coefficient) for each basis tensor. Importantly, the number of parameters is independent of the number of nodes. The proof for identifying the basis tensors sets up a fixed point equation with Kronecker products of any permutation matrix that any equivariant tensor must satisfy. The solutions to these equations are defined by equivalence classes of multi-indices in $[n]^{k}$. Each equivalence class is represented by a partition $\gamma$ of $[k]$, e.g., $\gamma=\{\{1\},\{2,3\}\}$ includes all multi-indices $(i_{1},i_{2},i_{3})$ where $i_{1}\neq i_{2},i_{3}$ and $i_{2}=i_{3}$. The basis tensors $B^{\gamma}\in\{0,1\}^{n^{k}}$ are then such that $B^{\gamma}_{s}=1$ if and only if $s\in\gamma$.

Linear equivariant GNNs of order $k$ ($k$-LEGNNs) parameterized with the full basis are as discriminative as the $k$-WL algorithm (Theorem 8). To achieve this discriminative power, each entry $H^{(0)}_{s,:}$ in the input tensor encodes an initial coloring of the isomorphism type of the subgraph indexed by the $k$-tuple $s$.

3.3 Summary of Representational Power via WL

The following theorem summarizes equivalence results between the GNNs discussed so far and variants of the WL test. Following , we here use equivalence relations, as they suffice for universal approximation in Section 2.4. For a set $\mathcal{F}$ of functions defined on $\mathcal{G}$, define an equivalence relation $\rho$ via the joint discriminative power of all functions $F\in\mathcal{F}$, i.e., for any two graphs $G,G^{\prime}\in\mathcal{G}$:

The above GNN families have the following equivalences:

Analogous results hold for equivariant models (for node representations), with the exception of equality (2.15), which becomes an inclusion: \rho(\text{k-LGNN}_{E})\subseteq\rho(\text{k-WL}_{E}) .

3.4 Relational Pooling.

where $X_{\pi}$ is $X$ with permuted rows, and $H^{(0)}$ is the tensor combining adjacency matrix and node attributes. Here, $g$ is any permutation-sensitive function, and may be modeled via various nonlinear function approximators, e.g. neural networks such as fully connected networks (MLPs), recurrent neural networks or a combination of a convolutional network applied to $A$ and an MLP applied to $X$. In particular, this model allows to implement graph isomorphism testing via node IDs (cf. Section 2.2) if $g$ is a universal approximator . For instance, node IDs may be permuted over nodes and concatenated with the node attributes:

The drawback of the Relational Pooling (2.17) is its computational intractability. Various approximations have been considered, e.g., defining canonical orders, stochastic approximations, and applying $g$ to all possible $k$-subsets of $V$. In the latter case, increasing $k$ strictly increases the expressive power. Local Relational Pooling is a variant that applies relational pooling to the $k$-hop subgraphs centered at each node, and then aggregates the results. This operation provably allows to identify and count subgraphs of size up to $k$ .

3.5 Recursion

A general strategy for encoding a graph is to encode a collection of subgraphs and then aggregate these encodings. The question of what graphs $G,G^{\prime}$ this process allows to distinguish ($F(G)\neq F(G^{\prime})$), depends on the collection of subgraphs used, the subgraph encoding function and the aggregation function. As a special case, this process includes the reconstruction hypothesis , i.e., the question whether any graph $G$ can be reconstructed from the collection of its subgraphs $G\setminus\{v\}$, for all $v$ in $G$. One challenge with the reconstruction hypothesis is that no alignment of the $G\setminus\{v\}$ is available. Indeed, node correspondences across subgraphs provide important information .

Indeed, the expressive power of a model based on subgraph encodings depends on the set of subgraphs, the type of subgraph encodings and the aggregation. Tahmasebi et al. use recursion as a powerful tool: instead of iterative message passing or layering, a recursive application of the above subgraph embedding step, even with a simple set aggregation like (1.6), can enable a GNN that can count any bounded-size subgraphs, as opposed to MPNNs (Prop. 2).

Let $\mathcal{N}_{r}(v)$ be the $r$-hop neighborhood of $v$ in $G$. Recursive neighborhood pooling (RNP) encodes intersections of such neighborhoods of different radii. Given an input graph $G$ with node attributes $(h^{\text{in}}_{u})_{u\in V(G)}$ and a sequence $(r_{1},\ldots,r_{t})$ of radii, RNP-GNN recursively encodes the node-deleted $r_{1}$-neighborhoods $G_{v}=\mathcal{N}_{r_{1}}(v)\setminus\{v\}$ of all nodes $v\in V$ after marking the deletion in augmented representations $h^{\text{aug}}_{u}$, $u\in V$. It then combines the results, and returns node representations of all nodes. Concretely, for each $v\in V$, it computes $G_{v}$ and

If the sequence of radii is empty (base case), then the algorithm returns the input attributes $h^{\text{in}}_{u}$. In contrast to iterative message passing, the encoded subgraphs here correspond to intersections of local neighborhoods. Together with the node deletions and markings that retain node correspondences, this maintains more structural information. Formally, if the sequence of radii dominates a covering sequence for a subgraph $H$ of interest, then, with appropriate parameters, RNP can count the induced and non-induced subgraphs of $G$ isomorphic to $H$ . The computational cost is $O(n^{k})$ for recursion depth $k$, and better for very sparse graphs, in line with computational lower bounds.

4 Universal approximation

Distinguishing given graphs is closely tied to approximating continuous functions on graphs. In early work, Scarselli et al. take a fixed point view and show a universal approximation result for infinite-depth MPNNs whose layers are contraction operators, for functions on equivalence classes defined by computation trees. Dehmamy et al. analyze the ability of GNNs to compute polynomials of the adjacency matrix.

Later works derive universal approximation results for graph and permutation-equivariant functions from graph discrimination results via extensions of the Stone-Weierstrass theorem . Maron et al. argue that $H$-invariant networks (for a permutation group $H$) can universally approximate $H$-invariant polynomials, which in turn can universally approximate any invariant function . Keriven and Peyré do not fix the size of the graph and show that shallow equivariant networks can, with a single set of parameters, well approximate a function on graphs of varying size. Both constructions involve very large tensors.

Generalization

Beyond approximation power, a second important question in machine learning is generalization. Generalization asks how well the estimated function $\widehat{F}$ is performing according to the population risk, i.e., $\mathcal{R}(\widehat{F})$, as a function of the number of data points $N$ and model properties. Good generalization may demand explicit (e.g., via a penalty term) or implicit regularization (e.g., via the optimization algorithm). Hence, generalization analyses involve aspects of the complexity of the model class $\mathcal{F}$, the target function we aim to learn, the data and the optimization procedure. This is particularly challenging for neural networks, due to the nested functional form and the non-convexity of the empirical risk.

A classic learning theoretic perspective bounds the generalization gap $\mathcal{R}(\widehat{F})-\widehat{\mathcal{R}}(\widehat{F})$ via the complexity of the model class $\mathcal{F}$ (Section 3.1). These approaches do not take into account possible implicit regularization via the optimization procedure. One possibility to do so is via the Neural Tangent Kernel approximation (Section 3.2). Finally, for more complex, structured target functions, e.g., algorithms or physics simulations, one may want to also consider the structure of the target task. One such option is Algorithmic Alignment (Section 3.3). Another strategy for obtaining generalization bounds is via algorithmic stability, the condition that, if one data point is replaced, the outcome of the learning algorithm does not change much. This strategy led to some early bounds for spectral GNNs .

Vapnik-Chervonenkis dimension. The first GNN generalization bound was based on bounding the Vapnik-Chervonenkis (VC) dimension of the GNN function class $\mathcal{F}$. The VC dimension of $\mathcal{F}$ expresses the maximum cardinality of a set of data points such that for any binary labeling of the data, some GNN in $\mathcal{F}$ can perfectly fit, i.e., shatter, the set. The VC dimension directly leads to a bound on the generalization gap. Here, we only state the results for sigmoid activation functions.

The VC dimension of GNNs with $p$ parameters, $H$ hidden neurons (in the MLP) and input graphs of size $n$ is $O(p^{2}H^{2}n^{2})$.

Strictly speaking, Theorem 10 is for node classification with one hidden layer in the aggregation function MLPs. The VC dimension directly yields a bound on the generalization gap: for a class $\mathcal{F}$ with VC dimension $D$, with probability $1-\delta$, it holds that

Interestingly, in these bounds, GNNs are a generalization of recurrent neural networks . The VC dimension bounds for GNNs are the same as for recurrent neural networks ; the bounds for fully connected MLPs are missing the factor $n^{2}$ .

Rademacher Complexity. Bounds that are in many cases tighter can be obtained via Rademacher complexity. The empirical Rademacher complexity $\widehat{\mathfrak{R}}_{S}(\mathcal{F})$ of a function class $\mathcal{F}$ measures how well it can fit “noise” in the form of uniform random variables $\sigma=(\sigma_{1},\ldots,\sigma_{N})$ in $\{-1,+1\}$:

Let $\mathcal{C}$ be the product of the Lipschitz constants of $\phi,\rho,g$ and $B$; $T$ the number of GNN iterations; $w$ the dimension of the embeddings $h^{(t)}_{v}$, and $d$ the maximum branching factor in the computation tree. Then the generalization gap of the GNN can be bounded as: $\widetilde{O}(\frac{wd}{\sqrt{N}\gamma})$ for $\mathcal{C}<1/d$, $\widetilde{O}(\frac{wdT}{\sqrt{N}\gamma})$ for $\mathcal{C}=1/d$ and $\widetilde{O}(\frac{wd\sqrt{wT}}{\sqrt{N}\gamma})$ for $\mathcal{C}>1/d$.

2 Generalization bounds via the Neural Tangent Kernel

In contrast to the results in Section 3.1, the complexity measure $\mathbf{y}^{\top}\widetilde{H}^{-1}\mathbf{y}$ of the target function is data-dependent. If the target function to be learned follows a simple GNN structure with a polynomial, then this bound can be polynomial:

Let $\bar{h}_{v}=c_{v}\sum_{u\in\mathcal{N}(v)\cup\{v\}}h_{u}$. If the labels $y_{i}$, $1\leq i\leq N$, satisfy

3 Generalization via Algorithmic Alignment

The Graph NTK analysis shows a polynomial sample complexity if the function to be learned is close to the computational structure of the GNN, in a simple way. While this applies to mainly simpler learning tasks, the idea of an “alignment” of computational structure carries further. Recently, there has been growing interest in learning scientific tasks, e.g., given a set of particles or planets along with their location, mass and velocity, predict the next state of the system , and in “algorithmic reasoning”, e.g., learning to solve combinatorial optimization problems in particular over graphs . In such cases, the target function corresponds to an algorithm, e.g., a dynamic program.

While many neural network architectures have the power to represent such tasks, empirically, they do not learn them equally well from data. In particular, GNNs perform well here, i.e., their architecture encodes suitable inductive biases . As a concrete example, consider the Shortest Path problem. The computational structure of MPNNs matches that of the Bellman-Ford (BF) algorithm very well: both “algorithms” iterate, and in each iteration $t$, update the state as a function of the neighboring nodes and edge weights $w(u,v)$:

Hence, the GNN can simulate the BF algorithm if it uses sufficiently many iterations, and if the aggregation function approximates the BF state update (relaxation step). Intuitively, this update is a much simpler function to learn than the full algorithm as a black box, i.e., the GNN encodes much of the algorithmic structure, sparsity and invariances in the architecture. More generally, MPNNs match the structure of many dynamic programs in an analogous way , as long as the updates are permutation invariant or sufficient node identification is provided as input, in light of the results in Section 2. Dudzik and Veličković refine and generalize the relations between GNNs and dynamic programming by using category theory.

The NTK results formalize simplicity by a small function norm in the RKHS associated with the Graph NTK; this can become complicated with more complex tasks and multiple layers. To quantify structural match, Xu et al. define algorithmic alignment by viewing a neural network as a structured arrangement of learnable modules – in a GNN, the (MLPs in the) aggregation functions – and define complexity via sample complexity of those modules in a PAC-learning framework. Sample complexity in PAC learning is defined as follows: We are given a data sample $\{(x_{i},y_{i})\}_{i=1}^{N}$ drawn i.i.d. from a distribution $\mathcal{P}$ that satisfies $y_{i}=g(x_{i})$ for an underlying target function $g$. Let $f=\mathcal{A}(\{x_{i},y_{i}\}_{i=1}^{N})$ be the function output by a learning algorithm $\mathcal{A}$. For a fixed error $\epsilon$ and failure probability $1-\delta$, the function $g$ is $(N,\epsilon,\delta)$-PAC learnable with $\mathcal{A}$ if

The sample complexity $\mathcal{C}_{\mathcal{A}}(g,\epsilon,\delta)$ is the smallest $N$ so that $g$ is $(N,\epsilon,\delta)$-learnable with $\mathcal{A}$.

Let $g$ be a target function and $\mathcal{N}$ a neural network with $M$ modules $\mathcal{N}_{i}$. The module functions $f_{1},...,f_{M}$ generate $g$ for $\mathcal{N}$ if, by replacing $\mathcal{N}_{i}$ with $f_{i}$, the network $\mathcal{N}$ simulates $g$. Then $\mathcal{N}$ $(N,\epsilon,\delta)$-algorithmically aligns with $g$ if (1) $f_{1},...,f_{M}$ generate $g$ and (2) there are learning algorithms $\mathcal{A}_{i}$ for learning $f_{i}$ with $\mathcal{N}_{i}$, with sample complexity $M\cdot\max_{i}C_{\mathcal{A}_{i}}(f_{i},\epsilon,\delta)\leq N$.

Algorithmic alignment resembles Kolmogorov complexity . Thus, it can be hard to obtain the optimal alignment between a neural network and an algorithm. But, any algorithmic alignment yields a bound, and any with acceptable sample complexity may suffice. The complexity of the MLP modules in GNNs may be measured with a variety of techniques. One option is the NTK framework. The module-based bounds then resemble the polynomial bound in Theorem 13, since both are extensions of . However, here, the bounds are applied at a module level, and not for the entire GNN as a unit. Theorem 14 translates these bounds, in a simplified setting, into sample complexity bounds for the full network.

Fix $\epsilon$ and $\delta$. Suppose $\left\{G_{i},y_{i}\right\}_{i=1}^{N}\sim\mathcal{P}$, where $\lvert V(G_{i})\rvert\leq n$, and $y_{i}=g(G_{i})$ for some $g$. Suppose $\mathcal{N}_{1},...,\mathcal{N}_{M}$ are network $\mathcal{N}$’s MLP modules in sequential order of processing. Suppose $\mathcal{N}$ and $g$ $(N,\epsilon,\delta)$-algorithmically align via functions $f_{1},...,f_{M}$ for a constant $M$. Under the following assumptions, $g$ is $(N,O(\epsilon),O(\delta))$-learnable by $\mathcal{N}$. a) Sequential learning. We train $\mathcal{N}_{i}$’s sequentially: $\mathcal{N}_{1}$ has input samples $\{\hat{x}_{i}^{(1)},f_{1}(\hat{x}_{i}^{(1)})\}_{i=1}^{N}$, with $\hat{x}_{i}^{(1)}$ obtained from $G_{i}$. For $j>1$, the input $\hat{x}_{i}^{(j)}$ for $\mathcal{N}_{j}$ are the outputs of the previous modules, but labels are generated by the correct functions $f_{j-1},...,f_{1}$ on $\hat{x}_{i}^{(1)}$. b) Algorithm stability. Let $\mathcal{A}$ be the learning algorithm for the $\mathcal{N}_{i}$’s. Suppose $f=\mathcal{A}(\left\{x_{i},y_{i}\right\}_{i=1}^{N})$, and $\hat{f}=\mathcal{A}(\left\{\hat{x}_{i},y_{i}\right\}_{i=1}^{N})$. For any $x$, $\|f(x)-\hat{f}(x)\|\leq L_{0}\cdot\max_{i}\|x_{i}-\hat{x}_{i}\|$, for some $L_{0}<\infty$. c) Lipschitzness. The learned functions $\hat{f}_{j}$ satisfy $\|\hat{f}_{j}(x)-\hat{f}_{j}(\hat{x})\|\leq L_{1}\|x-\hat{x}\|$, for some $L_{1}<\infty$.

The big $O$ notation here hides factors including the Lipschitz constants, number of modules and graph size. When measuring module complexity via the NTK, Theorem 14 indeed yields a gap in upper bounds between fully connected networks and GNNs in simple cases , supporting empirical results. While some works use sequential training in experiments , empirically, better alignment improves learning and generalization even with more common “end-to-end” training, i.e., optimizing all parameters simultaneously .

At a general level, these alignment results indicate how incorporating expert knowledge, e.g. in terms of algorithmic techniques or physics, into the design of the learning method can improve sample efficiency.

Extrapolation

Structural similarity of graphs. One possibility to guarantee successful extrapolation to larger graphs is to assume sufficient structural similarity between the graphs in $\mathcal{P}$ and $\mathcal{Q}$, in particular, structural properties that matter for the GNN family under consideration. For spectral GNNs, this assumption has been formalized as the graphs arising from the same underlying topological space, manifold or graphon. Under such conditions, spectral GNNs – with conditions on the employed filters – can generalize to larger graphs . The underlying structure also ensures similar local structure of the graphs.

For spatial message passing GNNs, whose representations rely on computation trees as local structures (Section 2.1), an agreement in the distributions of the computation trees in the graphs sampled from $\mathcal{P}$ and $\mathcal{Q}$ is necessary . This is violated, for instance, if the degree distribution is a function of the graph size, as is the case for random graphs under the Erdős-Rényi or Preferential Attachment models. The computation tree of depth $t$ rooted at a node $v$ corresponds to the color $c^{(t)}(v)$ assigned by the 1-WL algorithm.

Let $\mathcal{P}$ and $\mathcal{Q}$ be finitely supported distributions of graphs. Let $\mathcal{P}^{t}$ be the distribution of colors $c^{(t)}(v)$ over $\mathcal{P}$ and similarly $\mathcal{Q}^{t}$ for $\mathcal{Q}$. Assume that any graph in $\mathcal{Q}$ contains a node with a color in $\mathcal{Q}^{t}\setminus\mathcal{P}^{t}$. Then, for any graph regression task solvable by a GNN with depth $t$ there exists a GNN with depth at most $t+3$ that perfectly solves the task on $\mathcal{P}$ and predicts an answer with arbitrarily large error on all graphs from $\mathcal{Q}$.

The proof exploits the fact that GNN predictions on nodes only depend on the associated computation tree and that a sufficiently flexible GNN can assign arbitrary target labels to any computation tree . I.e., the available information allows for multiple local minima of the empirical risk. A similar result can be shown for node prediction tasks. (A “sufficiently large” GNN here means depth at least $t+2$ layers and width $\max\{(\max\deg(G)+1)^{t}\cdot|C|,2\sqrt{|P|}\}$, where the max degree refers to any graph in the support, $|C|$ is the finite number of possible input node attributes and $P$ the set of colors encountered in graphs in the support.)

Conditions on the GNN. If sufficient structural similarity of the input graphs cannot be guaranteed, then further restrictions on the GNN can enable extrapolation to different graph sizes, structures and ranges of input node attributes. If there are no training observations in a certain range of attributes or set of local structures, then the predictions of the learned model depend on the inductive biases induced by the model architecture, loss function and training algorithm. Which prediction function, out of multiple fitting functions, a model will choose, depends on these biases.

Xu et al. analyze such biases to obtain conditions on the GNN for extrapolation. Taking the perspective of algorithmic alignment (Section 3.3), they first analyze how individual module functions, i.e., the MLPs in the aggregation function of a GNN, extrapolate, and then transfer this to the entire GNN. The aggregation functions enter the extrapolation regime, e.g., if the node attributes, node degrees or computation trees are different under $\mathcal{Q}$ compared to $\mathcal{P}$, as they determine the inputs to the aggregations. The following theorem states that, sufficiently far away from $\operatorname{supp}(\mathcal{P})$, MLPs implement directionally linear functions.

The linear function and the constant terms in the convergence rate depend on the training data and the direction $v$. The proof of Theorem 16 relies on the fact that a neural network in the NTK regime learns a minimum-norm interpolation function . Although Theorem 16 uses a simplified setting of a wide 2-layer network, similar results hold empirically for more general MLPs .

To appreciate the implications of this result in the context of GNNs, consider the example of Shortest Path in Equation (3.4). For the aggregation function to mimic the Bellman-Ford algorithm, the MLP must approximate a nonlinear function. But, in the extrapolation regime, it implements a linear function and therefore is expected to not approximate Bellman Ford well any more. Indeed, empirical works that successfully extrapolate GNNs for Shortest Path use a different aggregation function of the form

Here, the nonlinear parts do not need to be learned, allowing to extrapolate with a linear learned MLP. More generally, the directionally linear extrapolation suggests that (1) the architecture or (2) the input encoding should be set up such that the target function can be approximated when MLPs learn linear functions (linear algorithmic alignment). An example for (2) may be found in forecasting physical systems, e.g., predicting the evolution of $n$ objects in a gravitational system, and the node (object) attributes are mass, location and velocity at time $t$. The position of an object at time $t+1$ is a nonlinear function of the attributes of the other objects. When encoding the nonlinear function as transformed edge attributes, the function to be learned becomes linear. Many empirical works that successfully extrapolate implement the idea of linear algorithmic alignment .

Finally, the geometry of the training data also plays an important role. show empirical results and initial theoretical results for learning max-degree, that, even with linear algorithmic alignment, sufficient diversity in the training data is needed to identify the correct linear functions. These data conditions are weaker than those implied by Theorem 15, due to the linear algorithmic alignment assumption.

For the case when the target test distribution $\mathcal{Q}$ is known, Yehudai et al. propose approaches for combining elements of $\mathcal{P}$ and $\mathcal{Q}$ to enhance the range of the data seen by the GNN.

Conclusion

This survey summarized three main topics in theoretically understanding GNNs: representation and approximation, generalization, and extrapolation. As GNNs are an active research area, many results could not be included. E.g., we focused on MPNNs and main ideas for higher-order GNNs, but neglected spectral GNNs, which closely relate to ideas in graph signal processing. Other emergent topics include adversarial robustness, optimization behavior of the empirical risk and its improvements, and computational scalability and approximations. Overall, GNNs have a rich set of mathematical connections, a selection of which was covered here.

Many questions remain. Regarding approximation capabilities, the limitations of MPNNs have motivated powerful higher-order GNNs. However, these are still computationally expensive. What efficiency is theoretically possible? Moreover, most applications may not require full graph isomorphism power, or $k$-WL power for large $k$. What other measures of representational power make sense? Do they allow better and sharper complexity results? Initial works consider subgraph counting as a benchmark task .

The generalization results so far need to use simplifications in the analysis, similar to most theoretical analyses of deep learning. To what extent can they be relaxed? Do more specific tasks or graph classes allow sharper results? Which modifications of GNNs would allow them to generalize better, and how do higher-order GNNs generalize? Similar questions pertain to extrapolation and reliability under distribution shifts, a topic that has been studied even less than GNN generalization.

In general, revealing further mathematical connections may enable the design of richer models and enable a more thorough understanding of GNNs’ learning abilities and limitations, and eventual improvements.

The author would like to thank Keyulu Xu, Derek Lim, Behrooz Tahmasebi, Vikas Garg, Tommi Jaakkola, Andreas Loukas, Jingling Li, Mozhi Zhang, Simon Du, Ken-ichi Kawarabayashi, Weihua Hu, Jure Leskovec, Joan Bruna and Yusu Wang for discussions on the theory of GNNs, collaborations and pointers.

This work was partially supported by NSF CAREER award 1553284, NSF SCALE MoDL award 2134108, and NSF CCF-2112665 (TILOS AI Research Institute).

References