Graph machine learning: message passing, expressiveness, and scale

Why this subject matters now

A decade ago "machine learning on graphs" meant hand-designed features fed to a linear model, or kernel methods whose cost grew quadratically in the number of graphs. Three things changed. First, the 2014–2017 wave (DeepWalk, node2vec, ChebNet, GCN, GraphSAGE, GAT) showed that differentiable, parameter-shared architectures could learn representations of nodes and graphs end to end, and that a single layer type, aggregate messages from neighbors then update, covers a wide range of problems. Second, the theory caught up. In 2019 two groups independently proved that this entire family is exactly as powerful as a graph isomorphism heuristic from 1968, which turned architecture design from folklore into a question with theorems, and spawned a hierarchy of provably more expressive models. Third, the applications became industrial, with recommendation at Pinterest over three billion nodes, molecular property prediction and machine learning force fields accurate enough to replace quantum-chemistry calls, physics and weather simulation at DeepMind, fraud detection at every large payments company, and combinatorial optimization inside commercial solvers. A practitioner today is expected to know not just how to call a GNN library but why sum aggregation beats mean for structure discrimination, why ten layers is usually worse than three, what a graph transformer buys and costs, and how any of this trains when the graph does not fit in memory. Those are exactly the questions this page works through.

Graphs as data

The objects and the tasks

A graph \( G = (V, E) \) has \( n = |V| \) nodes and \( m = |E| \) edges, an adjacency matrix \( A \in \{0,1\}^{n \times n} \) (weighted variants replace the 1s), a degree matrix \( D = \diag(d_1, \dots, d_n) \) with \( d_i = \sum_j A_{ij} \), and usually a node feature matrix \( X \in \R^{n \times d} \) whose row \( x_i \) describes node \( i \). Edges may carry features too. Prediction problems come at three granularities. Node-level tasks attach a label to each node, such as classifying papers in a citation network, flagging fraudulent accounts, or predicting protein function. These are often transductive, the test nodes are present in the training graph with labels hidden, which is unusual relative to standard supervised learning and matters for evaluation. Edge-level tasks predict whether or what kind of edge exists between a pair, as in link prediction in a social network, recommendation (user-item edges), and knowledge-graph completion. Graph-level tasks attach a label to a whole graph, such as a molecule's solubility, a program's correctness, or a circuit's timing. Each granularity ends in a different readout, respectively a per-node head, a pairwise decoder on two node embeddings, and a pooling over all nodes, but the representation machinery underneath is shared.

Permutation invariance and equivariance, precisely

The defining property of graph data is that node identity is arbitrary. Storing the same graph with nodes numbered differently changes \( A \) and \( X \) but not the object they describe. Formally, let \( P \in \{0,1\}^{n \times n} \) be a permutation matrix (exactly one 1 per row and column, \( P\T P = I \)). Renumbering nodes by the permutation \( \pi \) maps \( A \mapsto P A P\T \) and \( X \mapsto P X \). A graph-level function \( f \) is permutation invariant if

meaning the output does not depend on the ordering at all. A molecule's boiling point cannot change because the atoms were listed in a different order. A node-level function \( F : (A, X) \mapsto \R^{n \times d'} \), producing one output row per node, is permutation equivariant if

meaning relabeling the input relabels the output the same way, and nothing else changes. Node \( i \)'s embedding is a function of the graph as seen from node \( i \), not of the integer \( i \). Equivariance is the right property for the internal layers of any graph network (each layer maps node states to node states), and invariance is what a graph-level readout must add at the end, which is why readouts are sums, means, or maxes over nodes, since symmetric functions are exactly the invariant ones. The composition rule is what makes this a design principle rather than an afterthought. A stack of equivariant layers followed by one invariant pooling is invariant end to end, so it suffices to enforce the symmetry layer by layer. This framing, symmetry first, architecture second, is the organizing idea of the geometric deep learning program of Bronstein, Bruna, Cohen, and Veličković, which recovers CNNs (translation symmetry), GNNs (permutation symmetry), and spherical networks (rotation symmetry) from the same recipe.

Why standard architectures fail on graphs

An MLP applied to the flattened adjacency matrix \( \mathrm{vec}(A) \in \R^{n^2} \) fails all three requirements at once. It is not invariant, because the input dimension assigned to entry \( (i,j) \) has its own private weights, so two orderings of the same graph produce different outputs unless the network happens to learn all \( n! \) symmetries from data, which for \( n = 20 \) is \( 2.4 \times 10^{18} \) equivalent presentations of each training example. It is not size-generalizing, because the weight matrix is bound to a fixed \( n \). And it has no locality prior, since nothing ties the weight for edge \( (1,2) \) to the weight for edge \( (7,9) \), though the statistics of both are identical, so parameters are not shared where the symmetry says they should be. CNNs fail more subtly. A convolution is exactly the equivariant linear map for the translation group on a grid, and a grid is a graph, but a very special one. Every interior node has an identically shaped, ordered neighborhood (up, down, left, right), which is what lets a 3×3 kernel assign a distinct weight to each neighbor position. General graphs have neighborhoods of varying size with no canonical ordering, so there is no consistent way to say which neighbor gets which kernel weight. RNNs impose a linear order that the data does not have. Any serialization of a graph into a sequence breaks permutation symmetry the same way the MLP does. The conclusion, which the rest of the page builds on, is that a graph layer may use a neighbor's identity only through the multiset of neighbor states. Aggregation must be a symmetric function, and all remaining expressive power must come from what is computed before and after the aggregation.

Classical graph features and kernels

Degree, centrality, and local structure

Before learned representations, node features were computed. The degree \( d_i \) is the zeroth-order description. Eigenvector centrality defines importance recursively, a node is important if its neighbors are important, \( c_i = \tfrac{1}{\lambda} \sum_j A_{ij} c_j \), i.e. \( A c = \lambda c \), and taking the leading eigenvector (which Perron-Frobenius guarantees is entrywise nonnegative for a connected graph) gives a well-defined score. PageRank is the damped, stochastic version, the stationary distribution of a walk that with probability \( \beta \) follows a uniformly random out-edge and with probability \( 1 - \beta \) teleports uniformly, \( \pi\T = \beta \pi\T D^{-1} A + (1-\beta) \tfrac{1}{n} \mathbf{1}\T \). Betweenness counts the fraction of shortest paths through a node, closeness is the inverse mean distance to all others, and the local clustering coefficient measures how close a node's neighborhood is to a clique, \( C_i = 2\,T_i / \big(d_i (d_i - 1)\big) \) where \( T_i \) is the number of edges among node \( i \)'s neighbors, equivalently the number of triangles through \( i \). These remain strong baselines. A later section shows that message-passing networks provably cannot compute some of them (triangle counts among them), which is one reason concatenating such features onto \( X \) still helps in practice.

Graphlets

Graphlets generalize the triangle count. Fix all connected graphs on up to \( k \) nodes (there are 2 on three nodes, 6 on four, 21 on five, counting only connected ones), and describe a node by the vector counting how many induced copies of each graphlet touch it, broken down by which position (orbit) the node occupies. The graph-level version counts graphlet occurrences in the whole graph, and the graphlet kernel between two graphs is the inner product of their normalized count vectors. The obstruction is cost, since exact counting of \( k \)-node graphlets is \( O(n^{k}) \) in general, so \( k \) stops at 5 and sampling estimators take over. The kernel view matters historically because it framed the question the right way. A graph similarity is legitimate exactly when it is computed from permutation-invariant statistics, and the field's job is to find statistics that are cheap, discriminative, and (later) learnable.

Color refinement: the Weisfeiler-Lehman test, step by step

The most important classical algorithm for this page is 1-dimensional Weisfeiler-Lehman color refinement (Weisfeiler and Lehman, 1968). Every node starts with the same color (or its discrete label if it has one). Each round, every node builds a signature consisting of its own color plus the multiset of its neighbors' colors, and an injective relabeling (a hash) maps each distinct signature to a fresh color.

The double braces denote a multiset. Multiplicity counts, order does not. The partition of nodes by color can only refine over rounds, so after at most \( n \) rounds it stabilizes. To test two graphs for isomorphism, run refinement on both and compare the color histograms each round. If they ever differ, the graphs are certainly not isomorphic. If they stabilize equal, the test is inconclusive (the graphs may or may not be isomorphic). The WL kernel of Shervashidze et al. (2011) turns this into a graph similarity by running \( T \) rounds and taking the inner product of the concatenated color histograms, which is computable in \( O(Tm) \) and was the state of the art for graph classification for years. The worked problem below runs the algorithm completely on a pair of trees. The expressiveness section returns to the pairs where it fails.

T1:   1 -- 2 -- 3 -- 4 -- 5        T2:   1 -- 2 -- 3 -- 4 -- 5
           |                                  |
           6                                  6
      (leaf 6 attached to node 2)        (leaf 6 attached to node 3)

Shallow node embeddings

The encoder-decoder view

The first learned alternative to hand-built features assigns each node \( v \) a free vector \( z_v \in \R^d \) (the "encoder" is a lookup table \( Z \in \R^{n \times d} \)) and trains \( Z \) so that a simple decoder on pairs, usually the dot product \( z_u\T z_v \), reconstructs some notion of neighborhood similarity. Choosing the similarity to be adjacency itself and the loss to be squared error gives classical matrix factorization, \( \min_Z \| A - Z Z\T \|_F^2 \), whose optimum is the top-\( d \) eigendecomposition of \( A \). Choosing multi-hop or walk-based similarities gives the DeepWalk family. The unifying statement, made precise below, is that essentially all of these methods factorize some fixed matrix \( M(A) \), and differ only in which \( M \) and which loss.

DeepWalk and node2vec

DeepWalk (Perozzi, Al-Rfou, Skiena, 2014) imported the word2vec recipe. Sentences become truncated random walks, words become nodes, and the skip-gram objective asks \( z_u \) to predict the nodes appearing within a window \( T \) of \( u \) on the walks. node2vec (Grover and Leskovec, 2016) kept the objective and changed the walk. Its walk is second-order Markov. To choose the step after arriving at \( v \) from \( t \), each neighbor \( x \) of \( v \) is weighted by

where \( d(t,x) \) is the graph distance between the previous node and the candidate. The return parameter \( p \) controls immediate backtracking. The in-out parameter \( q \) tilts between staying near \( t \) and departing. Small \( q \) (say \( 0.25 \)) makes outward steps cheap, so walks behave DFS-like, ranging far and encoding homophily, nodes in the same community land near each other. Large \( q \) keeps walks in the immediate neighborhood, BFS-like, and the embeddings instead encode structural roles. Two hub nodes in different communities see similar local multisets and get similar vectors. This is measurable on even a toy graph. On a five-node graph, 10-step walks from a fixed start with \( q = 0.25 \) reached an average maximum distance of 2.71 hops from the start, versus 1.36 hops with \( q = 4 \), the DFS/BFS split in one number (measured with the walk implementation in the code section).

The negative-sampling objective, with its gradient

The exact skip-gram loss needs a softmax over all \( n \) nodes per prediction, which is intractable at scale. Negative sampling replaces it with binary discrimination. For each observed (center, context) pair \( (u, v) \) from the walks, and \( k \) "negative" nodes \( n_i \) drawn from a noise distribution \( P_n \) (empirically \( P_n(v) \propto d_v^{3/4} \)), maximize

The gradient is worth doing once by hand. Using \( \tfrac{d}{dx} \log \sigma(x) = \sigma(-x) \) and \( \tfrac{d}{dx} \log \sigma(-x) = -\sigma(x) \),

The positive term pulls \( z_u \) toward \( z_v \) with strength equal to how unconfident the model still is, and each negative pushes \( z_u \) away with strength equal to how wrongly confident it is. Both forces vanish as the classification saturates correctly. Concretely, if \( z_u\T z_v = 1 \) and one negative has \( z_u\T z_n = 0.5 \), the loss is \( -\log \sigma(1) - \log \sigma(-0.5) = 0.313 + 0.974 = 1.287 \), and the pull coefficient on \( z_v \) is \( \sigma(-1) = 0.269 \) while the push coefficient on \( z_n \) is \( \sigma(0.5) = 0.622 \). The wrong negative currently gets the larger correction, which is the behavior you want.

What matrix is being factorized

Levy and Goldberg showed in 2014 that skip-gram with negative sampling has a closed-form pointwise optimum. Fix a pair \( (u, v) \), let \( x = z_u\T z_v \), let \( \#(u,v) \) be the number of times the pair occurs in the walk corpus \( D \), and \( \#(u), \#(v) \) the marginal counts. Collecting all terms of the objective that involve this pair (the pair appears \( \#(u,v) \) times as a positive, and \( v \) is drawn as a negative for \( u \) about \( k \, \#(u) \tfrac{\#(v)}{|D|} \) times) gives

Setting \( \tfrac{d\L}{dx} = \#(u,v)\, \sigma(-x) - k \tfrac{\#(u)\#(v)}{|D|} \sigma(x) = 0 \) and using \( \sigma(x)/\sigma(-x) = e^{x} \) gives

the pointwise mutual information of the pair, shifted by \( \log k \). So the embedding is implicitly factorizing the PMI matrix of walk co-occurrences. Qiu et al. (2018, NetMF, from Tsinghua) carried the computation one step further by evaluating the co-occurrence statistics analytically for random walks. On a connected undirected graph the probability that a length-\( T \) window contains the pair at offset \( r \) involves the \( r \)-step transition matrix \( P^r \) with \( P = D^{-1} A \), and the limit matrix DeepWalk factorizes is

Every ingredient is a polynomial in the degree-normalized adjacency, which ties the walk methods directly to the spectral story of the next section. DeepWalk is a low-pass filter of the normalized adjacency, learned by sampling instead of by eigendecomposition. The limits of the whole shallow family follow from the encoder being a lookup table. Parameters grow as \( O(nd) \) with no sharing, a node unseen at training time has no embedding (transductive only), and node features \( X \) are ignored entirely. Fixing those three limits is what graph neural networks are for.

Spectral graph theory and the road to GCN

The Laplacian and its variants

For an undirected graph, define the combinatorial Laplacian, the symmetric normalized Laplacian, and the random-walk Laplacian.

\( L \) and \( L_{\mathrm{sym}} \) are symmetric positive semidefinite, so they have an orthonormal eigenbasis \( U = [u_1, \dots, u_n] \) with real eigenvalues \( 0 = \lambda_1 \le \lambda_2 \le \dots \le \lambda_n \). For \( L_{\mathrm{sym}} \) there is additionally \( \lambda_n \le 2 \), with equality exactly when the graph has a bipartite component (Chung's book is the standard reference for these facts). The constant vector is in the null space of \( L \) (each row sums to zero), and the multiplicity of eigenvalue 0 equals the number of connected components. \( \lambda_2 \), the Fiedler value, measures how hard the graph is to disconnect and controls mixing times and, as shown later, over-smoothing rates.

Dirichlet energy and the Rayleigh quotient, derived

The reason eigenvalues mean "frequencies" is a two-line identity. For any signal \( x \in \R^n \) (one number per node),

Rewrite the first term using \( \sum_i d_i x_i^2 = \sum_{i,j} A_{ij} x_i^2 = \tfrac{1}{2} \sum_{i,j} A_{ij} (x_i^2 + x_j^2) \) (symmetry of \( A \)), so

This is the Dirichlet energy, the total squared disagreement across edges, a discrete integral of the squared gradient. The Rayleigh quotient \( R(x) = \tfrac{x\T L x}{x\T x} \) therefore measures smoothness per unit norm, and by the Courant-Fischer theorem the eigenvectors are exactly the stationary points of \( R \), where \( u_1 \) (constant) is the smoothest possible signal with \( R = 0 \), \( u_2 \) is the smoothest signal orthogonal to constants, and so on up to \( u_n \), the most oscillatory signal the graph supports. A low eigenvalue means smooth, or "low frequency", and a high eigenvalue means rapidly alternating across edges, or "high frequency". Problem 5 verifies this numerically on a 4-cycle, where the energy computed edge by edge matches \( \sum_i \lambda_i \hat{x}_i^2 \) exactly.

The graph Fourier transform and spectral convolution

On the discretized circle, the Laplacian eigenvectors are the discrete Fourier modes, so it is natural to define the graph Fourier transform by analogy, with analysis \( \hat{x} = U\T x \) and synthesis \( x = U \hat{x} \). The classical convolution theorem says convolution in the original domain is pointwise multiplication in frequency, and since general graphs have no translation operator to define convolution directly, the theorem is promoted to a definition. A spectral filter with frequency response \( \hat{g} : \R \to \R \) acts as

Bruna et al. (2014) trained \( \hat{g}(\lambda_i) \) as free parameters, one per eigenvalue. This has three problems. The eigendecomposition costs \( O(n^3) \) once and each filtering costs \( O(n^2) \). The filter is tied to one graph's basis, so nothing transfers across graphs. And a generic spectral filter has no locality, since changing \( x \) at one node can change the output everywhere.

ChebNet: polynomial filters are local

Defferrard, Bresson, and Vandergheynst (2016, EPFL) fixed all three problems with one move, restricting \( \hat{g} \) to a degree-\( K \) polynomial. If \( \hat{g}(\lambda) = \sum_{k=0}^{K} \theta_k \lambda^k \), then

because \( U \Lambda^k U\T = (U \Lambda U\T)^k = L^k \). No eigendecomposition is needed, the parameters \( \theta_0, \dots, \theta_K \) are graph-independent, and since \( (L^k)_{ij} \ne 0 \) only when \( j \) is within \( k \) hops of \( i \), the filter is exactly \( K \)-localized. For numerical conditioning ChebNet uses the Chebyshev basis \( T_0(x) = 1 \), \( T_1(x) = x \), \( T_k(x) = 2x\,T_{k-1}(x) - T_{k-2}(x) \) on the rescaled operator \( \tilde{L} = \tfrac{2}{\lambda_{\max}} L_{\mathrm{sym}} - I \) (mapping the spectrum into \( [-1, 1] \), where the Chebyshev polynomials are bounded by 1),

computed by the recurrence with \( K \) sparse matrix-vector products, \( O(Km) \) total.

From ChebNet to GCN: the first-order approximation and the renormalization trick

Kipf and Welling (2017, Amsterdam) took the crudest useful version of ChebNet and turned it into the default graph layer. Set \( K = 1 \) and approximate \( \lambda_{\max} \approx 2 \) (its upper bound for \( L_{\mathrm{sym}} \)), so \( \tilde{L} = L_{\mathrm{sym}} - I = -D^{-1/2} A D^{-1/2} \). Then

Two free parameters per filter invite overfitting and doubling of work, so constrain \( \theta = \theta_0 = -\theta_1 \), which gives

Now a subtlety with real consequences. The operator \( I + D^{-1/2} A D^{-1/2} \) has eigenvalues in \( [0, 2] \) (it is \( 2I - L_{\mathrm{sym}} \)), so stacking \( \ell \) layers applies an operator with spectral radius up to \( 2^{\ell} \), and deep stacks explode or vanish. The renormalization trick replaces the sum of identity and normalized adjacency with a single normalized operator on the self-looped graph.

\( \hat{S} \) is symmetric with eigenvalues in \( (-1, 1] \) (the self-loops also pull the most negative eigenvalue strictly above \( -1 \)), so repeated application is non-expansive. Applying the filter to every feature channel with a weight matrix \( W \in \R^{d \times d'} \) mixing channels, and adding a nonlinearity, gives the GCN layer.

Entrywise, node \( i \) computes \( \sum_{j \in \mathcal{N}(i) \cup \{i\}} \tfrac{1}{\sqrt{(d_i + 1)(d_j + 1)}}\, W\T h_j \), a degree-normalized average of neighbor features, linearly transformed. Every derivation choice is visible in the final behavior. The polynomial restriction is why GCN is local, the \( K = 1 \) truncation is why one layer sees one hop, the shared \( \theta \) is why it has so few parameters, and the renormalization is why it can be stacked at all. The same operator \( \hat{S} \) reappears below as the exact cause of over-smoothing.

The message-passing framework

The general aggregate-update form

Gilmer et al. (2017, Google) observed that GCN, and essentially every graph network proposed before or since, is an instance of one template. Each layer, every node receives a message from each neighbor, reduces the messages with a permutation-invariant aggregator, and updates its state.

where \( \phi \) is the message function (it may read the edge feature \( e_{uv} \)), \( \bigoplus \) is a symmetric reduction (sum, mean, max, or an attention-weighted sum), and \( \psi \) is the update. Because \( \bigoplus \) is symmetric in the neighbors and \( \phi, \psi \) are shared across all nodes, each layer is permutation equivariant by construction, which is the entire symmetry argument from the first section discharged in one line. After \( L \) layers, \( h_v^{(L)} \) is a function of the \( L \)-hop neighborhood of \( v \). A graph-level readout \( h_G = \mathrm{READOUT}(\{\!\{ h_v^{(L)} \}\!\}) \) with a symmetric READOUT makes the composition invariant. The named architectures are choices of \( \phi, \bigoplus, \psi \).

GCN as message passing

GCN sets \( \phi(h_v, h_u) = \tfrac{1}{\sqrt{(d_v+1)(d_u+1)}} W\T h_u \), \( \bigoplus = \) sum over \( \mathcal{N}(v) \cup \{v\} \), and \( \psi = \sigma \). The implementation-relevant reading is that with the graph stored as an edge list \( (\mathrm{src}, \mathrm{dst}) \), the layer is a dense matmul \( H W \), a gather \( (HW)[\mathrm{src}] \), a per-edge scale by the normalization coefficient, and a scatter-add into the destination rows. That four-op decomposition is exactly what the code section implements, and it is the computational signature of every message-passing layer, one dense GEMM plus irregular gather/scatter traffic proportional to \( m \).

GraphSAGE: sampling and concatenation

Hamilton, Ying, and Leskovec (2017) designed GraphSAGE for the inductive, large-graph setting. Two changes from GCN. First, the update concatenates the node's own transformed state with the aggregated neighborhood instead of averaging itself in with the neighbors.

optionally followed by \( \ell_2 \) normalization. Keeping the self-channel separate matters when a node should not be smoothed into its neighborhood, an early nod to heterophily. Second, \( \mathcal{N}(v) \) is not the full neighborhood but a fixed-size uniform sample of it, redrawn each minibatch, which bounds the cost of a layer regardless of degree. The scaling section quantifies what this buys. The paper studies three aggregators. Mean averages neighbor vectors. Max-pooling applies a learned pointwise MLP then an elementwise max, \( \max_{u} \sigma(W_{\mathrm{pool}} h_u + b) \), which detects the presence of feature patterns in the neighborhood regardless of how many neighbors express them. LSTM runs a sequence model over a random permutation of the neighbors. Strictly it is not permutation invariant, and the paper accepts the broken symmetry for the extra capacity, randomizing the order so the model cannot rely on it. The now-standard judgment is that the LSTM aggregator was a dead end (symmetric aggregators with more expressive \( \phi \) do better), but mean and max remain defaults, and the fixed-fanout sampling idea became the backbone of all industrial-scale GNN training.

GAT: learned edge weights via attention

GCN's edge weights are fixed by degrees. Veličković et al. (2018) let the network decide how much each neighbor matters. For each edge \( (j \to i) \) (self-loops included), compute a scalar compatibility from the transformed endpoint features, then normalize over each node's in-neighborhood with a softmax.

Here \( a \in \R^{2d'} \) is a learned vector and the softmax runs over a different support per node, its incident edges, which is why implementations need a segment-softmax (per-destination max-subtraction, exponentiation, per-destination sum, divide), not a dense one. Multi-head attention runs \( K \) independent copies and concatenates (or averages, at the final layer). Note where this sits in the design space. Attention reweights the aggregation, but the support is still the one-hop neighborhood, so GAT is message passing with a learned \( \bigoplus \), not a transformer. Brody et al. (2022) identified one structural weakness. Since \( a\T [W h_i \| W h_j] = a_1\T W h_i + a_2\T W h_j \), the ranking of \( e_{ij} \) over \( j \) is determined by \( a_2\T W h_j \) alone, identical for every query node \( i \) up to an additive constant that the softmax cancels. The original GAT can only express this "static" attention. GATv2 moves the nonlinearity inside, \( e_{ij} = a\T \mathrm{LeakyReLU}( W [h_i \| h_j] ) \), and recovers query-dependent rankings at the same cost.

GIN: why sum aggregation and an MLP

Xu, Hu, Leskovec, and Jegelka (2019, MIT and Stanford) asked which aggregate-update choices lose information, and the answer clarified the whole design space. The layer's job, viewed abstractly, is to map the pair (own state, multiset of neighbor states) to a new state. If two different pairs ever map to the same output, the network has irreversibly merged two distinguishable neighborhoods. So maximal discriminative power requires the layer to be injective on (state, multiset) pairs. Check the standard aggregators against multisets over a feature vocabulary, say \( a = (1, 0) \) and \( b = (0, 1) \).

• Mean collapses multiplicity scaling, since \( \{\!\{a, b\}\!\} \) and \( \{\!\{a, a, b, b\}\!\} \) both average to \( (0.5, 0.5) \). Mean recovers the distribution of neighbor features, never the counts.
• Max collapses both multiplicity and proportion, since any multiset containing at least one \( a \) and one \( b \) maxes to \( (1, 1) \). Max recovers only the underlying set.
• Sum gives \( (1, 1) \) versus \( (2, 2) \), which are distinct. Sum preserves the full multiset when features are one-hot, because the sum of one-hots is exactly the histogram.

The positive result is a lemma. If the feature space \( \mathcal{X} \) is countable and multiset sizes are bounded, there exists \( f : \mathcal{X} \to \R^n \) such that \( h(S) = \sum_{x \in S} f(x) \) is unique for every multiset \( S \). The construction is digit encoding. Enumerate \( \mathcal{X} \) by \( z : \mathcal{X} \to \N \) and set \( f(x) = N^{-z(x)} \) where \( N \) bounds the multiset size. Then \( h(S) \) writes the multiplicities of \( S \) in base \( N \), and distinct multisets give distinct digits. Moreover \( (c, S) \mapsto (1 + \epsilon) f(c) + \sum_{x \in S} f(x) \) is injective on pairs for irrational \( \epsilon \), because the irrational weighting prevents the center's contribution from being confused with any integer combination of neighbor contributions. Any function on multisets then factors as \( g(h(S)) \), and since \( f \) and \( g \) are unknown, both are approximated by MLPs (this is where universal approximation enters, and why a linear layer after the sum is not enough, as a 1-layer perceptron cannot represent the required \( g \circ \) digit-decode). The layer that implements the argument is the Graph Isomorphism Network.

with \( \epsilon \) learned or fixed. GIN is provably the most expressive architecture in the message-passing class. It matches the 1-WL ceiling that the next section establishes for the whole family. In practice, GIN-style sum aggregation dominates on graph-level tasks where structure matters (molecules), while mean/attention aggregation often wins on noisy node-level tasks where the distribution of neighbor features is the signal and degree is a nuisance variable. Expressiveness is necessary for the former and nearly irrelevant to the latter.

Expressiveness: message passing meets Weisfeiler-Lehman

The equivalence, at the level of the argument

The 2019 results of Xu et al. and, independently, Morris et al. (whose paper also introduced higher-order k-GNNs) say that anonymous message-passing networks are at most as powerful as 1-WL at distinguishing graphs, and can match it with injective components. Both directions are inductions worth internalizing.

Upper bound. The claim is that if 1-WL assigns two nodes (possibly in different graphs) the same color at iteration \( t \), then any message-passing network assigns them equal features at layer \( t \). Induction on \( t \). For the base case, with identical (or no) initial labels, both statements hold at \( t = 0 \). For the step, suppose it holds at \( t \). Equal WL colors at \( t + 1 \) means equal colors at \( t \) and equal neighbor color multisets at \( t \) (the hash is injective, so equal outputs force equal signatures). By the induction hypothesis, the two nodes then have equal features at layer \( t \) and their neighborhoods have equal feature multisets at layer \( t \). But a message-passing layer computes the new feature as a function only of (own feature, neighbor feature multiset), the same function for every node, so the two nodes get equal features at layer \( t + 1 \). The contrapositive is the useful reading. A GNN can only separate what WL separates, because the GNN's per-layer information is a possibly lossy function of the WL signature, while WL's hash is lossless. Summing over nodes for a graph readout preserves the statement at graph level. Equal WL color histograms force equal pooled embeddings.

Matching the bound. If every layer is injective on (feature, multiset) pairs, and the readout is injective on the final multiset of node features, then by the same induction run in reverse, distinct WL signatures produce distinct features at every step, so the network distinguishes exactly the graphs WL distinguishes after the same number of rounds. GIN's sum-MLP construction supplies the injectivity, hence "GIN is as powerful as 1-WL". The theorem's fine print is part of the lesson. It concerns anonymous networks (features carry no identity), \( L \) layers match \( L \) WL rounds, and "distinguishing" is about equality of representations, not about what a downstream classifier can learn from finite data.

What 1-WL cannot see

The ceiling is low in concrete, checkable ways. The refinement update depends only on color multisets, so on any \( r \)-regular graph with uniform initial colors, every node has the identical signature \( (c, \{\!\{ c^{(r)} \}\!\}) \) every round, refinement never gets off the ground, and any two \( r \)-regular graphs on \( n \) nodes are indistinguishable. Two canonical pairs follow.

C6 (one 6-cycle)          2 x C3 (two triangles)
    1 -- 2                    1        4
   /      \                  / \      / \
  6        3                2---3    5---6
   \      /
    5 -- 4
Both 2-regular, 6 nodes, 6 edges. 1-WL: identical color histograms
at every round. Connectivity differs; WL cannot tell.

K3,3 (complete bipartite)     Prism (C3 x K2)
  1   2   3                     1 ------- 4
  |\ /|\ /|                    /|         |\
  | X | X |    vs             2-+---------+-5
  |/ \|/ \|                    \|         |/
  4   5   6                     3 ------- 6
(every top connected          (two triangles 1-2-3, 4-5-6
 to every bottom)              joined by a perfect matching)
Both 3-regular, 6 nodes, 9 edges. K3,3 is triangle-free and
bipartite; the prism has two triangles. 1-WL: identical.

Running the color-refinement implementation from the code section confirms both failures. The color histograms agree at every round for C6 versus the two triangles, and for K3,3 versus the prism. The corollary that matters for applications is that since a message-passing network computes a function of the WL colors, no such network can compute any quantity that differs between a WL-equivalent pair. Triangle counts differ between K3,3 (zero) and the prism (two), so no anonymous MPNN can count triangles, or cycles generally, or distinguish many molecular substructures (decalin versus bicyclopentyl is the standard chemistry example of a WL-equivalent pair). Girth, diameter, and connectivity are similarly out of reach. C6 is connected, the two triangles are not, and WL cannot tell.

The k-WL hierarchy

The classical fix is to color tuples instead of nodes. \( k \)-WL maintains a color for every \( k \)-tuple of nodes, initialized by the tuple's isomorphism type (which nodes are equal, which pairs are adjacent), and refines it using the colors of the \( n \) neighbors obtained by substituting each position of the tuple. The hierarchy is strict for every \( k \ge 2 \), since \( (k{+}1) \)-WL distinguishes graphs that \( k \)-WL cannot (Cai, Fürer, and Immerman constructed the separating families in 1992), and \( k \)-WL decides isomorphism only for graphs of bounded treewidth. Morris et al.'s k-GNNs implement neural analogues and inherit the power, at the cost that defines the whole area, \( O(n^k) \) states and up to \( O(n^{k+1}) \) work per iteration. 3-WL (which can count triangles and distinguishes K3,3 from the prism) already costs \( O(n^3) \) memory, impractical beyond a few thousand nodes. Provably-powerful architectures at 3-WL strength with \( O(n^2) \) memory exist (the PPGN line of Maron et al., built from matrix products of equivariant tensors), but the quadratic object is still the whole \( n \times n \) relation, so none of this scales to large sparse graphs.

Practical routes beyond 1-WL, and their price

Four families trade the clean hierarchy for deployability. Positional and structural encodings augment \( X \) with features that break WL-equivalence before message passing starts, such as Laplacian eigenvectors, random-walk return probabilities \( (\hat{S}^r)_{vv} \) for \( r = 1..R \) (these count closed walks, so they see triangles), distances to anchor sets, or explicit subgraph counts. The runtime cost is near zero, with some preprocessing. The price is that "positions" are not canonical (the eigenvector sign problem, treated with graph transformers below). Random features (Sato et al. 2021, Abboud et al. 2021, Oxford) assign each node a random identifier. Two WL-equivalent nodes now differ almost surely, and with enough width such networks are universal in expectation. The price is that the function is no longer deterministic or exactly equivariant, and generalization relies on averaging over draws. Subgraph GNNs (ESAN of Bevilacqua et al., node-marking variants) run a base GNN on a bag of perturbed copies of the graph, one per deleted or marked node, and aggregate. The result is strictly more powerful than 1-WL, bounded by 3-WL. The price is a factor-\( n \) blowup in compute per graph. Higher-order networks (k-GNN, PPGN) buy provable power at the polynomial memory costs above. The honest summary is that on molecule-sized graphs (tens of nodes), subgraph and higher-order methods are affordable and measurably better. On million-node graphs, structural encodings are the only entry on the menu, and in many industrial tasks the 1-WL ceiling is not the binding constraint anyway, feature quality and sampling noise are.

Training pathologies

Over-smoothing, derived from the spectrum

Strip a deep GCN to its propagation skeleton by ignoring weights and nonlinearities, leaving \( H^{(L)} = \hat{S}^L X \) with \( \hat{S} = \tilde{D}^{-1/2} \tilde{A} \tilde{D}^{-1/2} \). Since \( \hat S \) is symmetric, expand in its orthonormal eigenbasis, \( \hat{S} = \sum_i \mu_i v_i v_i\T \) with \( 1 = \mu_1 > \mu_2 \ge \dots \ge \mu_n > -1 \) for a connected non-bipartite graph (self-loops kill bipartiteness, so the strict inequalities hold for any connected graph after renormalization). First verify the claimed top eigenpair directly.

using \( \tilde{A} \mathbf{1} = \tilde{D} \mathbf{1} \) (row sums are degrees). So \( v_1 \propto \tilde{D}^{1/2} \mathbf{1} \), the vector with entries \( \sqrt{d_i + 1} \), has eigenvalue exactly 1. Powering the expansion gives

Every column of the feature matrix converges geometrically to a multiple of \( v_1 \). Node \( i \)'s limiting feature is proportional to \( \sqrt{d_i + 1} \) times a graph-wide constant per channel. All information except degree is destroyed, and the rate is governed by the spectral gap. Well-connected graphs (large gap, large \( \lambda_2 \) of the Laplacian) smooth fastest. This is the content of Li, Han, and Wu (2018), who first framed GCN as Laplacian smoothing. Oono and Suzuki (2020, Tokyo) extended the argument to the full nonlinear network. With ReLU (which is 1-Lipschitz and, being nonnegative, plays well with the Perron direction) and weight matrices of maximum singular value \( s \), the distance from \( H^{(L)} \) to the degree-scaled subspace contracts by at least \( s \bar\mu \) per layer, so whenever \( s < 1 / \bar\mu \), expressive power decays exponentially in depth. Nonlinearities do not save you.

The effect is easy to measure. On the Cora citation graph (2,708 nodes, 10,556 directed edge entries), repeatedly applying \( \hat{S} \) to the raw features and tracking the mean pairwise cosine similarity between node feature vectors (1,000-node sample) gives, from the run in the code section, 0.057 at depth 0, 0.148 after 1 step, 0.245 after 2, 0.344 after 4, 0.491 after 8, 0.621 after 16, 0.739 after 32, 0.813 after 64. With random Xavier-initialized weight matrices and ReLU between propagations, collapse is much faster, reaching 0.399 after 1 layer, 0.623 after 2, 0.821 after 4, 0.928 after 8, 0.977 after 16, and 0.9985 after 64, at which point all nodes are essentially the same vector up to scale. The pure-propagation curve converges at rate \( \bar\mu^L \) toward the degree profile. The nonlinear one collapses faster because untrained weights also lose rank.

Over-squashing and bottlenecks

Over-smoothing is about signals becoming too similar. Over-squashing, named by Alon and Yahav (2021, Technion), is about too much signal forced through too small a pipe. For a task with radius-\( L \) interactions, node \( v \)'s receptive field after \( L \) layers contains every node within \( L \) hops, a set that grows exponentially in \( L \) on anything tree-like (on a tree with branching factor \( b \) it is \( \Theta(b^L) \)), yet everything must be compressed into one fixed-width vector \( h_v \in \R^d \). The quantitative handle is the input-output Jacobian. For a GCN-style network, an application of the chain rule across layers bounds

with \( c \) collecting Lipschitz constants of the layer maps. Sensitivity to a distant input decays with the \( (v, u) \) entry of the powered propagation matrix. When the only paths from \( u \) to \( v \) squeeze through a narrow cut, that entry is exponentially small, and gradients cannot carry the dependency, regardless of the loss. Topping et al. (2022, Oxford) localized the blame in geometry. They define a discrete balanced Forman curvature on edges (positive inside dense cliques where neighborhoods overlap, negative on bridge-like edges whose endpoints share few neighbors, the discrete analogue of negative Ricci curvature), prove that strongly negatively curved edges cause the exponentially small Jacobian entries, and propose rewiring. Their SDRF algorithm adds support edges around the most negatively curved edge and removes the most positively curved ones, surgically widening bottlenecks before training. The same goal drives simpler tricks. Alon and Yahav's fix was to make the last layer fully adjacent (every node talks to every node once), and later work connects graphs to expanders or adds virtual global nodes. All of these trade a little structural fidelity for a usable gradient path, and all foreshadow the graph transformer, which deletes the bottleneck by attending globally at every layer.

Depth, residuals, jumping knowledge, and normalization

The practical consequence of both pathologies is that GNN depth does not behave like CNN depth. On homophilous node classification, 2 to 4 layers is almost always optimal, and going to 16 or 64 requires machinery. Residual connections \( h^{(\ell+1)} = h^{(\ell)} + f(h^{(\ell)}) \) make identity the default and slow the spectral collapse (the operator becomes \( I + \hat{S} \), whose eigenvalues \( 1 + \mu_i \) no longer shrink the top of the spectrum relative to the rest as aggressively). The initial residual of GCNII (Chen et al., 2020) mixes a fraction of \( H^{(0)} \) into every layer, which provably prevents convergence to the degree profile since the limit must retain an \( H^{(0)} \) component. Jumping knowledge (Xu et al., 2018) feeds the readout a combination (concatenation, max, or attention) of all intermediate layers \( h_v^{(1)}, \dots, h_v^{(L)} \), so each node effectively selects its own receptive-field radius rather than inheriting the deepest one. DropEdge randomly deletes edges each epoch, thinning the propagation operator and delaying smoothing. Normalization needs graph-specific care. BatchNorm across a node minibatch mixes statistics across graphs of different sizes. PairNorm (Zhao and Akoglu, 2020) instead recenters and rescales so the mean pairwise feature distance is held constant across layers, directly counteracting the collapse measured above, and GraphNorm (Cai et al., 2021) normalizes within each graph with a learnable shift, which stabilizes graph-level training. None of these make depth free. They raise the useful maximum from about 3 to about 10 to 60 layers depending on task, and the honest default remains to use the shallowest network whose receptive field covers the interaction radius of the problem.

Graph transformers

Full attention over nodes

A graph transformer treats the node set as an unordered sequence and runs standard dense self-attention. Every node attends to every node, and the graph enters only through added encodings and optional attention biases. Two immediate consequences. Cost is \( O(n^2) \) per layer independent of sparsity, so the approach is natural for molecule-sized graphs (\( n \le 10^3 \)) and needs linear-attention surgery beyond that. And the bottleneck geometry disappears. Every pair of nodes has a length-1 interaction path, so over-squashing in the structural sense cannot occur, which is exactly the property message passing lacks. The price is that with no encodings, the architecture is blind to the graph entirely, attention over a bare node set is permutation equivariant but sees no edges, so all structural information must be injected explicitly. This makes the positional-encoding problem the heart of the area, not an implementation detail.

The positional-encoding problem

Sequences have a canonical coordinate, position \( 1..T \). Graphs do not, and every candidate coordinate has a defect to manage. Laplacian eigenvector encodings (Dwivedi and Bresson, 2021) take the \( k \) eigenvectors of \( L_{\mathrm{sym}} \) with smallest nonzero eigenvalues and append \( (u_2(v), \dots, u_{k+1}(v)) \) to each node's features, a soft coordinate system in which nearby nodes get similar values, the graph analogue of sinusoidal encodings (which are exactly the circle's Laplacian eigenvectors). The defect is sign ambiguity. If \( u \) is a unit eigenvector so is \( -u \), and the eigensolver's choice is arbitrary, so the same graph can arrive with any of \( 2^k \) sign patterns. With repeated eigenvalues the ambiguity is a full orthogonal transform within each eigenspace. Training-time random sign flipping teaches approximate invariance. SignNet (Lim et al., 2022, MIT) enforces it exactly by processing each eigenvector through \( \phi(u) + \phi(-u) \), an even function by construction. Random-walk structural encodings (RWSE) sidestep signs entirely. Node \( v \) gets the return probabilities \( \big( (\hat S)_{vv}, (\hat S^2)_{vv}, \dots, (\hat S^R)_{vv} \big) \), which are canonical, permutation equivariant, and encode local cycle structure (a nonzero 3-step return requires a triangle), though as purely local statistics they provide structure rather than position. Relative encodings bias attention logits per pair instead of per node. Graphormer (Ying et al., 2021, Microsoft Research) adds a learned scalar bias indexed by shortest-path distance \( b_{\mathrm{SPD}(i,j)} \), plus an edge-feature term averaged along a shortest path, plus degree embeddings on nodes, and this recipe won the OGB-LSC quantum chemistry challenge. The SPD bias is a graph analogue of relative position and is unaffected by sign issues, but costs an all-pairs BFS at preprocessing.

Hybrids, and when attention actually wins

GraphGPS (Rampášek et al., 2022) is the assembly kit that made the comparison fair. Each layer runs a local message-passing block and a global attention block in parallel on the same representation and sums them, with LapPE/SignNet or RWSE injected at the input, and linear-attention options (Performer) for large \( n \). The empirical picture from GraphGPS and the Long Range Graph Benchmark has two halves. On molecular property tasks with long-range effects, on tasks whose labels depend on global geometry, and wherever the graph has severe bottlenecks, the attention path pays for itself. On classic homophilous node classification and on tasks decided by local substructure, well-tuned MPNNs match transformers at a fraction of the cost, and the hybrid's MPNN block does most of the work. A clarifying way to state it is that a graph transformer is message passing on the complete graph with learned edge weights, so the choice is not a new paradigm but a choice of support, and dense support helps exactly when the task's dependency structure is denser than the input graph.

Scaling to large graphs

The neighbor-explosion problem, with the arithmetic

Minibatch training on graphs is not like minibatch training on images, because the loss at one node depends on its \( L \)-hop neighborhood. With mean degree \( \bar d \), the expected receptive field of a single seed node after \( L \) layers is roughly \( \sum_{\ell=0}^{L} \bar d^{\,\ell} \) (ignoring overlap, which helps only in dense regions). Problem 4 works the numbers for a real product graph. The conclusion is that a 3-layer full-neighborhood computation for a batch of 1,024 seeds can require gathering features for tens of millions of node visits, more than 50 GB of feature traffic per batch, which no accelerator holds. Every scalable method is a different answer to "which part of the receptive field do we refuse to compute".

Node-wise sampling and its variance

GraphSAGE's estimator replaces the sum over \( \mathcal{N}(v) \) with a sum over a uniform sample \( S \subset \mathcal{N}(v) \), \( |S| = f \), scaled to be unbiased for the mean aggregator. The estimate \( \tfrac{1}{f} \sum_{u \in S} h_u \) has expectation \( \tfrac{1}{d_v} \sum_{u \in \mathcal{N}(v)} h_u \) and variance proportional to \( \tfrac{1}{f}\big(1 - \tfrac{f-1}{d_v - 1}\big) \Var_u[h_u] \). Noise shrinks as \( 1/f \) but never reaches zero while \( f < d_v \), and crucially the noise is inside a nonlinearity, so the layer output is biased even when the pre-activation estimate is not (\( \E[\sigma(\hat m)] \ne \sigma(\E[\hat m]) \)). Depth compounds both the bias and the variance since each layer samples again. Historical-activation methods (VR-GCN of Chen et al., 2018) reduce the variance by using stale cached activations as control variates, sampling only the change since the cache. GNNAutoScale (Fey et al., 2021) pushed this to its logical end, keeping full historical embeddings for out-of-batch neighbors in CPU memory so the in-GPU computation touches only the batch, with staleness as the sole approximation.

Layer-wise and subgraph sampling

Node-wise sampling re-explodes per layer because each sampled node samples its own children. Layer-wise methods (FastGCN, LADIES) instead sample one shared set of nodes per layer, with importance weights, holding the per-layer cost fixed. The difficulty is keeping enough connectivity between consecutive layer samples. Subgraph sampling abandons the layered structure. Draw one subgraph, run a full GNN on it, and let the minibatch be the subgraph. Cluster-GCN (Chiang et al., 2019, Google) partitions the graph with METIS into dense clusters, trains on one cluster (or a random union of a few, which restores some of the cut edges) at a time. The bias is that all edges across partition boundaries are dropped from gradient computation, which METIS minimizes by construction. GraphSAINT (Zeng et al., 2020) samples subgraphs by node, edge, or short random-walk samplers and then removes the sampling bias analytically. Each node's loss contribution is weighted by \( 1/\P(v \in \mathcal{G}_s) \) and each edge's message by \( 1/\P((u,v) \in \mathcal{G}_s) \), with the inclusion probabilities estimated by pre-sampling, so the minibatch gradient is unbiased for the full-graph gradient while every batch is a small dense graph with no neighbor explosion at all. In practice GraphSAINT-style random-walk sampling and Cluster-GCN batching are the workhorses for datasets like the 100-million node OGB papers graph.

Distributed training and the systems constraints

Past a single machine, the graph itself is partitioned (DistDGL, PyG's distributed stack, Alibaba's AliGraph). Each worker owns a partition plus a halo of remote neighbors, and every minibatch triggers RPCs to fetch remote features. The dominant cost in production GNN training is usually not the GEMMs but feature gathering. Sampled neighborhoods touch essentially random rows of a feature matrix that lives in CPU RAM or across the network, so throughput is set by random-access bandwidth and the effectiveness of caching high-degree nodes (which appear in a large fraction of batches, so an LRU or degree-based cache captures most traffic). On the accelerator side, scatter-add aggregation is memory-bound. It does \( O(1) \) arithmetic per loaded value, so its ceiling is DRAM bandwidth, about 3 TB/s on an H100 80GB (measured 2,992 GB/s float32 copy on this machine's H100), while the dense transform \( HW \) can run at 745 TFLOPS in bf16 (measured, 4,096-dim square matmul). The design pressure this creates is visible in every serious system. Fuse the normalization into the gather, reorder nodes for locality (RCM or METIS orderings raise cache hit rates), use CSR with segment reductions instead of atomic scatter when degrees are skewed, and overlap sampling (CPU), feature fetch (PCIe/NVLink), and compute (GPU) in a producer-consumer pipeline, which is precisely the architecture PinSage described in 2018 and every framework since has reimplemented.

Heterogeneous graphs and knowledge graphs

Relational GCNs and metapaths

Real graphs are typed. A citation network has authors, papers, and venues, and a medical knowledge graph has drugs, targets, and diseases with dozens of edge types. Running one shared aggregator over all of it conflates semantically different relations. The relational GCN (Schlichtkrull et al., 2018, Amsterdam) gives each relation \( r \) its own weight matrix and normalizer.

which with hundreds of relations would need hundreds of full matrices, so R-GCN regularizes by basis decomposition, \( W_r = \sum_{b=1}^{B} a_{rb} V_b \), with all relations sharing \( B \) basis transforms and differing only in the mixing coefficients \( a_{rb} \), cutting parameters from \( |\mathcal{R}| d d' \) to \( B d d' + |\mathcal{R}| B \) and, as importantly, letting rare relations borrow statistical strength from common ones. The complementary idea for typed graphs is the metapath. A type sequence such as author \( \to \) paper \( \to \) author defines a derived "co-authorship" adjacency, and models like metapath2vec (walks constrained to follow a metapath) and HAN (attention over several metapath-derived graphs) push the type semantics into the connectivity rather than the weights. Metapaths encode domain knowledge and are interpretable. R-GCN learns the mixing but pays in parameters. Production systems commonly use both, a small set of curated metapaths to densify the graph, then a relational GNN on top.

Knowledge-graph embeddings and their relation-pattern algebra

Knowledge-graph completion scores triples \( (h, r, t) \), head entity, relation, tail entity, each embedded, with a scoring function \( f_r(h, t) \) trained so true triples score above corrupted ones. The classic models, with entity embeddings \( \mathbf{h}, \mathbf{t} \) and relation embedding \( \mathbf{r} \), are

where \( \circ \) is elementwise product. RotatE's constraint \( |r_i| = 1 \) makes each relation coordinate a rotation \( e^{i\theta_i} \) of the complex plane. What separates these models is which relation patterns they can represent exactly, namely symmetry (\( r(x,y) \Rightarrow r(y,x) \), e.g. "married to"), antisymmetry ("parent of"), inversion (\( r_2 = r_1^{-1} \), "advisor of" vs "advisee of"), and composition (\( r_3 = r_2 \circ r_1 \), "mother's husband is father"). Problem 6 proves the key entries, and the table below summarizes them.

None of the four handles one-to-many relations gracefully (TransE forces all tails of a fixed head-relation pair to the same point), which motivated the TransH/TransR projection line and, later, replacing the decoder entirely with a GNN encoder (R-GCN with a DistMult decoder was an early instance). At industrial scale, embedding tables dominate. A 100-million entity graph at \( d = 400 \) float32 is 160 GB of parameters before the model starts, which is why systems like DGL-KE shard entity embeddings across machines with asynchronous sparse updates.

Evaluating link prediction correctly

Link-prediction evaluation is a ranking exercise with two standard traps. The protocol, for each test triple \( (h, r, t) \), is to corrupt one side at a time, score \( (h, r, t') \) for every candidate entity \( t' \), rank the true tail among them, and likewise for heads. Report MRR \( = \tfrac{1}{|T|} \sum 1/\mathrm{rank} \) (reciprocal rank, dominated by the top of the list, unlike mean rank which one pathological example can ruin) and Hits@k, the fraction of cases ranked in the top \( k \). The first trap is the filtered setting. Many corruptions are themselves true triples present in train, valid, or test (a person has several children). Ranking the test answer below another correct answer is not an error, so all known true triples except the one being tested are removed from the candidate list before ranking. Unfiltered ("raw") numbers systematically understate good models. The second trap is ties. If the model gives many candidates identical scores (degenerate but common with ReLU outputs), placing the true triple first among ties inflates Hits@1. Sun et al. (2020) showed several published models lost most of their reported gains under tie-averaged ranking, and the accepted convention is now the average rank over ties (or random tie-breaking). A third, dataset-level trap is that FB15k and WN18 contain near-duplicate inverse relations, so a rule "if \( (t, r^{-1}, h) \) in train, predict \( t \)" scores near state of the art. FB15k-237 and WN18RR exist because the originals leak, and results quoted on the originals should be discounted.

Temporal and dynamic graphs

Real graphs change. Transactions arrive with timestamps, social edges form and dissolve, road conditions evolve. Two modeling regimes. Discrete-time methods take graph snapshots \( G_1, \dots, G_T \) and combine a GNN per snapshot with a sequence model across them, either on the node states (GCN then GRU per node) or on the parameters themselves (EvolveGCN runs an RNN whose hidden state is the GCN weight matrix, useful when nodes churn so per-node recurrence is ill-defined). Continuous-time methods consume an event stream \( (u_i, v_i, t_i, e_i) \) directly. The Temporal Graph Network of Rossi et al. (2020) is the reference design. Each node carries a memory vector updated by a recurrent cell whenever an event touches it. An embedding module computes usable representations by attending over recent temporal neighbors. Time enters through a learned encoding of the elapsed interval \( \Delta t \) (in TGAT, via Bochner's theorem, random-feature cosines of \( \omega \Delta t \), the temporal analogue of positional encodings). Training uses temporal batching with the constraint that a node's memory may only reflect events strictly before the one being predicted, otherwise the model reads its own answer. Evaluation carries its own leakage taxonomy. Negatives must be sampled respecting time (an edge that appears next week is not a valid negative today), transductive and inductive splits (new nodes at test time) must be reported separately, and the standard random-negative protocol saturates so quickly that recent benchmarks (TGB) moved to harder historical negatives, under which method rankings visibly reshuffle.

Generative models for graphs

Generating graphs inherits every difficulty of this page plus one more, that the likelihood itself is symmetry-afflicted. A distribution over graphs assigns mass to isomorphism classes, but a model over adjacency matrices assigns mass to labeled matrices, and one \( n \)-node graph corresponds to up to \( n! \) matrices (fewer if the graph has automorphisms). An order-dependent model defines \( p(G) = \sum_{\pi} p(A_\pi) \), a sum with \( n! \) terms that cannot be evaluated, so reported "likelihoods" from autoregressive graph models are likelihoods of one ordering, a lower bound of unquantified looseness, and comparing models by such numbers is meaningless unless the ordering distribution is identical. The three families manage the symmetry differently. Autoregressive. GraphRNN (You et al., 2018) generates nodes one at a time, emitting each new node's connections to existing nodes with an edge-level RNN inside a node-level RNN. The ordering problem is tamed (not solved) by training on BFS orderings only, which collapses many permutations into one canonical family and bounds the lookback window, since in a BFS order a new node can only connect to a recent frontier. One-shot latent-variable. GraphVAE-style models decode a whole adjacency matrix (all \( n^2 \) entries) from a latent vector. The decoder output must be compared to the target graph under unknown correspondence, which requires approximate graph matching in the loss, cubic-cost and brittle, which is why one-shot VAEs stalled at small molecules. Diffusion. This family is the current state of the art. GDSS diffuses continuous noise on node features and adjacency jointly. DiGress (Vignac et al., 2023, EPFL) runs a discrete diffusion directly on categorical node and edge types, with a graph-transformer denoiser. Diffusion fits the symmetry unusually well. If the noise process is exchangeable and the denoiser is permutation equivariant, the generated distribution is permutation invariant by construction, no ordering ever enters, and the model trains by denoising score matching rather than likelihood, sidestepping the \( n! \) sum. The residual costs are the \( O(n^2) \) edge-matrix denoiser and hundreds of sampling steps. SPECTRE (Martinkus et al., 2022, ETH Zurich) generates dominant Laplacian eigenvalues and eigenvectors first and conditions the rest of the graph on them, an interesting middle road that controls global structure directly.

Applications, with their specifics

Molecules and equivariant geometric GNNs

Molecular property prediction was the first domain where graph-shaped learning beat feature engineering end to end (Gilmer et al.'s MPNN paper was written for QM9), and it is where the symmetry story acquires a second group. A molecule is not only a graph. It has a conformation \( \{ \mathbf{x}_i \in \R^3 \} \), and physical quantities transform under the Euclidean group \( E(3) \) of rotations, reflections, and translations in specific ways. The energy is invariant, \( E(R\mathbf{x} + b) = E(\mathbf{x}) \), while forces are equivariant, transforming as vectors, \( F(R\mathbf{x} + b) = R\, F(\mathbf{x}) \). An architecture that does not respect this either wastes capacity learning the symmetry from augmentation or, worse, makes unphysical predictions that depend on the arbitrary lab frame. The design lineage tracks which geometric information enters and how. SchNet (Schütt et al., 2017) uses only interatomic distances, expanded in radial basis functions and turned into continuous convolution filters. This is rigorously invariant, but distances alone cannot see angles, and two different local geometries with equal distance multisets collide, the geometric cousin of the WL failure. DimeNet (Gasteiger et al., 2020) adds angles via directional message passing, messages live on edges and interact through the angle between adjacent edges, strictly increasing geometric expressiveness while staying invariant. EGNN (Satorras, Hoogeboom, Welling, 2021, Amsterdam) achieves \( E(n) \)-equivariance with almost no machinery. Alongside invariant features \( h_i \), it updates coordinates by

a sum of relative vectors scaled by invariant coefficients. Rotating the input rotates every \( \mathbf{x}_i - \mathbf{x}_j \), hence the update, hence the output, and equivariance follows in one line with no spherical harmonics. NequIP (Batzner et al., 2022, Harvard) goes the full representation-theoretic route. Features are geometric tensors labeled by rotation order \( \ell \) (scalars, vectors, higher tensors), messages combine neighbor features with spherical harmonics of the edge direction through Clebsch-Gordan tensor products (the e3nn library implements the algebra), and the payoff is data efficiency that changed the field's economics. NequIP reported matching or beating prior force fields with roughly a thousand DFT training configurations where invariant models needed one to two orders of magnitude more, which matters because the training data is quantum-chemistry compute. The through-line for interviews is that equivariance is a constraint that shrinks the hypothesis class to the physically consistent subset, and the gains show up precisely where data is expensive.

Drug discovery, recommendation, fraud

In drug discovery the best-known result is Stokes et al. (2020, MIT), where a directed-MPNN trained on about 2,300 measured compounds, screening the Drug Repurposing Hub in silico, surfaced halicin, an antibiotic with a novel mechanism validated in mice. The honest reading is that the GNN was a good ranker over a well-chosen library with tight experimental follow-up, not an oracle. Property predictors (absorption, toxicity, solubility) built on MPNN or geometric backbones are now routine components of screening funnels, with the caveat from the benchmarking literature that on several public property datasets, fingerprint-based models remain competitive. In recommendation, PinSage (Ying et al., 2018) remains the canonical industrial deployment, 3 billion nodes and 18 billion edges at Pinterest, where the innovations are systems and sampling rather than layers. Neighborhoods are defined not by adjacency but by short random walks, taking the \( T \) most-visited nodes as the neighborhood with visit counts as importance pooling weights, which handles power-law degrees and gives every node a fixed-size, relevance-weighted neighborhood. Training uses hard negatives scheduled by a curriculum (progressively harder impostors ranked by PageRank proximity), and inference over billions of nodes runs as a MapReduce pipeline with a producer-consumer GPU feed. The architecture is a two-layer GraphSAGE variant. The paper's lasting lesson is that at industrial scale the neighborhood definition, negative mining, and pipeline are the model. Fraud detection inverts a core assumption of most of this page. Fraudsters connect to victims and to legitimate hubs precisely to look normal, so the graphs are heterophilous and adversarial, camouflage in the literature's term. Low-pass aggregation then actively hurts, smoothing fraud signal into the benign majority, and deployed systems use relation-aware neighbor filtering (CARE-GNN and successors), separate aggregation of similar and dissimilar neighbors, and features from the temporal transaction stream, an instance of the general rule that the homophily assumption should be tested, not presumed.

Simulation and combinatorial optimization

The simulation line treats a mesh or particle system as a graph and learns the simulator. MeshGraphNets (Pfaff et al., 2021, DeepMind) is the reference, using encode-process-decode with 10 to 15 message-passing steps on the simulation mesh, augmented with world-space edges connecting nodes that are close in space though far on the mesh (how collisions and contact are seen), plus learned adaptive remeshing. It runs one to two orders of magnitude faster than the classical solvers it imitates on cloth, structural mechanics, and fluid benchmarks and generalizes to larger meshes than it trained on, because the learned local update is resolution-agnostic, the same parameter-sharing argument that makes GNNs size-generalize at all. The same encode-process-decode template with distance-based particle graphs (GNS, Sanchez-Gonzalez et al., 2020) handles sand, water, and goop, and the traffic instance of the pattern (diffusion-convolution plus a sequence model, DCRNN) has run inside production ETA systems. DeepMind's GraphCast pushed a related mesh-GNN onto a global icosahedral weather mesh and outperformed the leading operational forecast on most medium-range variables in 2023. In combinatorial optimization the sober summary is that GNNs work best as learned heuristics inside exact frameworks. Gasse et al. (2019) encode a MILP's variable-constraint bipartite graph and imitate strong branching, replacing an expensive expert computation inside branch and bound, and neural local-search and neural-guided sampling variants follow the same pattern. End-to-end neural solvers that emit tours or assignments directly still trail tuned classical solvers (LKH, Concorde, Gurobi) at realistic sizes and generalize poorly across instance scale, and the neural algorithmic reasoning program (Veličković and Blundell) is the research response, training networks to imitate algorithm traces, aligning the architecture with the computation rather than the answer.

An honest assessment

The reproducibility literature earned its place in this subject and should be part of how results are read. Shchur et al. (2018, TU Munich) showed that on the standard citation benchmarks, the ranking of GNN architectures reshuffles under different random train/test splits and equal hyperparameter budgets. Many claimed architectural advances were split luck. Errica et al. (ICLR 2020) ran a fair protocol over graph-classification datasets and found that on several chemical benchmarks a structure-agnostic baseline (an MLP over aggregated node features, no edges used) matched or beat published GNNs, and that reported numbers were often not comparable due to inconsistent selection protocols. Dwivedi et al.'s benchmarking suite made the same point constructively, with fixed budgets and splits on harder tasks. On node classification, Huang et al. (2021, Cornell) showed that "Correct and Smooth", a linear model or MLP followed by two rounds of classical label propagation, matched or beat large GNNs on several OGB leaderboards with orders of magnitude fewer parameters. Much of what message passing earns on homophilous graphs is recoverable by propagating labels rather than features. The field's structural response was OGB (Hu et al., 2020), with large graphs, fixed and meaningful splits (scaffold splits for molecules, temporal splits for citation and product graphs, so test distributions actually shift), standardized evaluators, and public leaderboards with required code. The current equilibrium is worth stating plainly. GNNs are clearly the right tool when the target is a function of structure (molecules, simulation, combinatorial problems) or when relational signal is the only signal (cold-start recommendation, fraud rings). On feature-rich homophilous node classification their edge over "good features plus propagation" is real but modest. And any single-benchmark claim of architectural superiority without budgeted, multi-split evaluation should be treated as unproven.

Worked problems

Problems 1 through 4 appear in their sections above, covering WL refinement on a tree pair, aggregator injectivity, the spectrum of the renormalized propagation matrix, and the neighbor-explosion arithmetic. Two more round out the set.

Implementation

Everything below was run on this machine (PyTorch 2.7, JAX 0.6, PyTorch Geometric 2.8, an H100 80GB available for the larger runs). The design decision that makes all of it work is the edge-list representation. A graph is two integer arrays src and dst of length \( m \), a message-passing layer is gather, per-edge compute, scatter-reduce, and every named architecture is a few lines of difference in the middle step.

Message passing from scratch: the GCN layer

The first block implements the sparse GCN layer exactly as derived. Append self-loops, compute the symmetric normalization \( 1/\sqrt{\tilde d_{\mathrm{src}} \tilde d_{\mathrm{dst}}} \) per edge, transform, gather from source nodes, scale, scatter-add into destination nodes. PyTorch uses index_add_ / scatter_add_. JAX uses segment_sum, which XLA compiles to the same scatter pattern. Both handle directed edge lists. For an undirected graph, store each edge in both directions.

import torch

def gcn_norm(edge_index, num_nodes):
    """Self-loops + symmetric normalization for the renormalized S-hat.

    edge_index: (2, E) int64, messages flow src=edge_index[0] -> dst=edge_index[1]
    returns edge_index with loops (2, E+N) and edge weights (E+N,)
    """
    loop = torch.arange(num_nodes, device=edge_index.device)
    src = torch.cat([edge_index[0], loop])            # (E+N,)
    dst = torch.cat([edge_index[1], loop])            # (E+N,)
    ones = torch.ones(src.numel(), device=src.device)
    deg = torch.zeros(num_nodes, device=src.device)   # d_i + 1 (loop included)
    deg.scatter_add_(0, dst, ones)                    # in-degree; equals out-degree
    dinv = deg.pow(-0.5)                              # for an undirected edge list
    w = dinv[dst] * dinv[src]                         # (E+N,) 1/sqrt(d~_i d~_j)
    return torch.stack([src, dst]), w

def gcn_layer(x, edge_index, weight, W, b=None):
    """One GCN layer: sigma is left to the caller.

    x: (N, d_in)   W: (d_in, d_out)   returns (N, d_out)
    """
    src, dst = edge_index
    h = x @ W                                         # (N, d_out) dense GEMM
    msg = h[src] * weight.unsqueeze(-1)               # (E+N, d_out) gather+scale
    out = torch.zeros_like(h)
    out.index_add_(0, dst, msg)                       # scatter-add reduce
    return out if b is None else out + b

import jax
import jax.numpy as jnp

def gcn_layer(x, edge_index, W, num_nodes, b=None):
    """One GCN layer with segment_sum as the scatter-add.

    x: (N, d_in)  edge_index: (2, E)  W: (d_in, d_out)  returns (N, d_out)
    """
    loop = jnp.arange(num_nodes)
    src = jnp.concatenate([edge_index[0], loop])      # (E+N,)
    dst = jnp.concatenate([edge_index[1], loop])      # (E+N,)
    ones = jnp.ones_like(src, dtype=x.dtype)
    deg_out = jax.ops.segment_sum(ones, src, num_nodes)   # d~ at source side
    deg_in  = jax.ops.segment_sum(ones, dst, num_nodes)   # d~ at target side
    w = deg_in[dst] ** -0.5 * deg_out[src] ** -0.5    # (E+N,)
    h = x @ W                                         # (N, d_out)
    out = jax.ops.segment_sum(h[src] * w[:, None], dst, num_nodes)
    return out if b is None else out + b

# num_nodes must be static under jit; close over it or mark it static:
gcn_layer_jit = jax.jit(gcn_layer, static_argnums=3)

Cross-checked against each other on a random directed graph, the two implementations agree to float32 roundoff (max abs difference 1.5e-4 across a 20-node, 4-channel output whose entries are order 1, a discrepancy from reduction-order nondeterminism, not logic).

GAT and GIN layers

GAT adds one genuinely new primitive, a segment softmax, a numerically stable softmax computed independently over each destination node's incoming edges (max-subtract per segment, exponentiate, normalize by the per-segment sum). GIN is the shortest layer in the family, a plain scatter-add with no normalization, followed by an MLP. The absence of degree normalization is the point, as derived in the injectivity section.

import torch
import torch.nn.functional as F

def gat_layer(x, edge_index, W, a_src, a_dst, neg_slope=0.2):
    """Single-head GAT. Caller supplies edge_index WITH self-loops.

    x: (N, d_in)  W: (d_in, d_out)  a_src, a_dst: (d_out,)
    e_ij = LeakyReLU(a_src . Wh_j + a_dst . Wh_i), softmax over j into i.
    """
    n = x.size(0)
    src, dst = edge_index
    h = x @ W                                          # (N, d_out)
    e = (h[src] * a_src).sum(-1) + (h[dst] * a_dst).sum(-1)   # (E,)
    e = F.leaky_relu(e, neg_slope)
    # segment softmax over incoming edges of each dst node
    e_max = torch.full((n,), float("-inf"), device=x.device)
    e_max = e_max.scatter_reduce(0, dst, e, reduce="amax")    # (N,)
    p = (e - e_max[dst]).exp()                         # (E,) stable exp
    denom = torch.zeros(n, device=x.device).scatter_add_(0, dst, p)
    alpha = p / denom[dst]                             # (E,) sums to 1 per dst
    out = torch.zeros_like(h)
    out.index_add_(0, dst, h[src] * alpha.unsqueeze(-1))
    return out

def gin_layer(x, edge_index, mlp, eps=0.0):
    """GIN: h_v = MLP((1+eps) x_v + sum_{u in N(v)} x_u). No normalization."""
    src, dst = edge_index
    agg = torch.zeros_like(x).index_add_(0, dst, x[src])   # (N, d) sum-aggregate
    return mlp((1.0 + eps) * x + agg)

import jax
import jax.numpy as jnp

def gat_layer(x, edge_index, W, a_src, a_dst, num_nodes, neg_slope=0.2):
    """Single-head GAT; edge_index already contains self-loops.

    x: (N, d_in)  W: (d_in, d_out)  a_src, a_dst: (d_out,)
    """
    src, dst = edge_index
    h = x @ W                                          # (N, d_out)
    e = (h[src] * a_src).sum(-1) + (h[dst] * a_dst).sum(-1)   # (E,)
    e = jnp.where(e > 0, e, neg_slope * e)             # LeakyReLU
    e_max = jax.ops.segment_max(e, dst, num_nodes)     # (N,)
    p = jnp.exp(e - e_max[dst])                        # (E,)
    denom = jax.ops.segment_sum(p, dst, num_nodes)     # (N,)
    alpha = p / denom[dst]                             # (E,)
    return jax.ops.segment_sum(h[src] * alpha[:, None], dst, num_nodes)

def gin_layer(x, edge_index, mlp_apply, mlp_params, num_nodes, eps=0.0):
    """GIN with a caller-supplied MLP (e.g. flax/haiku apply fn)."""
    src, dst = edge_index
    agg = jax.ops.segment_sum(x[src], dst, num_nodes)  # (N, d) sum-aggregate
    return mlp_apply(mlp_params, (1.0 + eps) * x + agg)

Both scratch layers were checked against the PyTorch Geometric references with copied parameters on a random 20-node graph (self-loops and duplicate edges removed first, since GATConv re-canonicalizes them). GAT gave maximum absolute difference 0.0 versus GATConv, and GIN gave maximum absolute difference 0.0 versus GINConv.

WL color refinement and node2vec walks

Two classical algorithms, both pure Python. The refinement uses sorted-tuple signatures and a dictionary as the injective hash. It is the code that produced the histograms quoted in Problem 1 and the C6 / K3,3 failure claims. The walker implements the exact \( (p, q) \) bias from the node2vec derivation.

from collections import Counter
import random

def wl_refine(adj, iters=None):
    """1-WL color refinement. adj: {node: [neighbors]}.
    Returns the stable color histogram (a Counter)."""
    nodes = sorted(adj)
    colors = {v: 0 for v in nodes}
    for _ in range(iters or len(nodes)):
        sigs = {v: (colors[v], tuple(sorted(colors[u] for u in adj[v])))
                for v in nodes}
        table = {}                       # injective hash: signature -> new color
        new = {}
        for v in nodes:
            if sigs[v] not in table:
                table[sigs[v]] = len(table)
            new[v] = table[sigs[v]]
        if all(new[v] == colors[v] for v in nodes):   # partition stabilized
            break
        colors = new
    return Counter(colors.values())

# Failure pairs from the expressiveness section, verified:
c6     = {0:[1,5],1:[0,2],2:[1,3],3:[2,4],4:[3,5],5:[4,0]}
two_c3 = {0:[1,2],1:[0,2],2:[0,1],3:[4,5],4:[3,5],5:[3,4]}
k33    = {0:[3,4,5],1:[3,4,5],2:[3,4,5],3:[0,1,2],4:[0,1,2],5:[0,1,2]}
prism  = {0:[1,2,3],1:[0,2,4],2:[0,1,5],3:[0,4,5],4:[1,3,5],5:[2,3,4]}
assert wl_refine(c6) == wl_refine(two_c3)     # cannot distinguish
assert wl_refine(k33) == wl_refine(prism)     # cannot distinguish

def node2vec_walk(adj, start, length, p, q, rng):
    """One biased second-order walk. alpha = 1/p (back), 1 (dist 1), 1/q (away)."""
    walk = [start]
    if not adj[start]:
        return walk
    walk.append(rng.choice(adj[start]))
    while len(walk) < length:
        prev, cur = walk[-2], walk[-1]
        nbrs = adj[cur]
        prev_nbrs = set(adj[prev])
        wts = [(1.0 / p) if x == prev else
               1.0 if x in prev_nbrs else
               (1.0 / q) for x in nbrs]
        walk.append(rng.choices(nbrs, weights=wts, k=1)[0])
    return walk

# Measured on a 5-node graph, 2000 walks of length 10 from node 0:
#   q = 0.25 (DFS-like): mean max distance reached = 2.71 hops
#   q = 4.0  (BFS-like): mean max distance reached = 1.36 hops

Verification against PyTorch Geometric on Cora

Two levels of verification on the Cora citation graph (2,708 nodes, 10,556 directed edge entries, 1,433-dimensional bag-of-words features, 7 classes, the standard public split). First, numerical. Copy the weights out of a randomly initialized GCNConv into the scratch layer and compare outputs on the full graph. The result was maximum absolute difference 0.0, mean absolute difference 0.0, bitwise identical, because both reduce with the same deterministic scatter kernel over the same edge ordering. Second, end-to-end. Train the standard 2-layer, 16-hidden-unit GCN (Adam, learning rate 0.01, weight decay 5e-4, dropout 0.5, 200 epochs) three seeds each. The scratch implementation reached test accuracies 0.788, 0.800, 0.805 (mean 0.798), and PyTorch Geometric reached 0.808, 0.806, 0.804 (mean 0.806). Both sit at the published GCN number for this split (81.5% with early stopping in Kipf and Welling). The small gap between the two is seed noise in initialization order, not a logic difference, as the layer-level check proves.

import torch
import torch.nn.functional as F
from torch_geometric.datasets import Planetoid
from torch_geometric.nn import GCNConv

data = Planetoid(root="./cora", name="Cora")[0]
x, edge_index, y = data.x, data.edge_index, data.y     # (2708,1433) (2,10556)
n = x.size(0)

# --- layer-level agreement, identical parameters ---
conv = GCNConv(1433, 16, bias=True)
W = conv.lin.weight.data.t().clone()                   # PyG stores (out,in)
b = conv.bias.data.clone()
ei, w = gcn_norm(edge_index, n)
ours, theirs = gcn_layer(x, ei, w, W, b), conv(x, edge_index)
print((ours - theirs).abs().max())    # -> 0.0  (bitwise identical)

# --- end-to-end training, scratch layer only ---
class ScratchGCN(torch.nn.Module):
    def __init__(self, d_in, d_h, d_c):
        super().__init__()
        self.W1 = torch.nn.Parameter(torch.empty(d_in, d_h))
        self.W2 = torch.nn.Parameter(torch.empty(d_h, d_c))
        torch.nn.init.xavier_uniform_(self.W1)
        torch.nn.init.xavier_uniform_(self.W2)
        self.b1 = torch.nn.Parameter(torch.zeros(d_h))
        self.b2 = torch.nn.Parameter(torch.zeros(d_c))
    def forward(self, x, ei, w):
        h = F.relu(gcn_layer(x, ei, w, self.W1, self.b1))
        h = F.dropout(h, 0.5, self.training)
        return gcn_layer(h, ei, w, self.W2, self.b2)

model = ScratchGCN(1433, 16, 7)
opt = torch.optim.Adam(model.parameters(), lr=0.01, weight_decay=5e-4)
for epoch in range(200):
    model.train(); opt.zero_grad()
    out = model(x, ei, w)
    F.cross_entropy(out[data.train_mask], y[data.train_mask]).backward()
    opt.step()
model.eval()
acc = (model(x, ei, w).argmax(1)[data.test_mask]
       == y[data.test_mask]).float().mean()
# 3 seeds: 0.788 / 0.800 / 0.805  (PyG same recipe: 0.808 / 0.806 / 0.804)

Measuring over-smoothing

The experiment behind the numbers in the over-smoothing section applies the propagation operator repeatedly to Cora's features and tracks the mean pairwise cosine similarity over a fixed 1,000-node sample, once for pure propagation \( \hat S^L X \) and once interleaving random Xavier weight matrices and ReLU (an untrained deep GCN). If the derivation is right, similarity must rise toward 1 geometrically, and it does.

import torch
import torch.nn.functional as F

def mean_pairwise_cos(h, sample=1000, seed=0):
    g = torch.Generator().manual_seed(seed)
    idx = torch.randperm(h.size(0), generator=g)[:sample]
    hs = F.normalize(h[idx], dim=1)                    # (S, d) unit rows
    sim = hs @ hs.t()                                  # (S, S)
    s = sim.size(0)
    return ((sim.sum() - s) / (s * (s - 1))).item()    # off-diagonal mean

ei, w = gcn_norm(edge_index, n)
src, dst = ei

h = x.clone()
for depth in range(1, 65):
    msg = h[src] * w.unsqueeze(-1)                     # propagate: h <- S-hat h
    h = torch.zeros_like(h).index_add_(0, dst, msg)
    if depth in (1, 2, 4, 8, 16, 32, 64):
        print(depth, round(mean_pairwise_cos(h), 4))

# Measured (Cora, this machine):
# pure propagation S^L X:
#   depth  0: 0.057   1: 0.148   2: 0.245   4: 0.344
#   depth  8: 0.491  16: 0.621  32: 0.739  64: 0.813
# untrained GCN (Xavier weights + ReLU each step):
#   depth  1: 0.399   2: 0.623   4: 0.821   8: 0.928
#   depth 16: 0.977  32: 0.994  64: 0.9985

Reading the numbers against the theory, pure propagation climbs steadily toward the degree profile at the rate set by \( \bar\mu \), while the untrained network with nonlinearities reaches 0.93 mean cosine similarity by depth 8, at which point a linear probe on the features is already close to useless. This is why the residual, initial-residual, and PairNorm interventions from the pathologies section exist, and why a 2-layer GCN is the default and not a compromise.

How it is done in practice

The gap between the derivations and a production system is mostly systems engineering, and it concentrates in four places. The sparsity format. Research code uses COO edge lists with atomic scatter-adds. Production kernels convert to CSR and run segmented reductions, which replace contended atomics with per-row loops and vectorized loads. The workload is memory-bound either way. A scatter-add performs one multiply-add per 8 bytes loaded, so its ceiling is DRAM bandwidth, about 3 TB/s on an H100 80GB (2,992 GB/s measured float32 copy on this machine), while the dense feature transform of the same layer can run at 745 TFLOPS in bf16 (measured at 4,096-square). A GNN layer is therefore a bandwidth-bound half glued to a compute-bound half, and the standard optimizations, fusing normalization into the gather, reordering nodes with METIS or RCM so neighbors share cache lines, batching small graphs into one block diagonal graph, all attack the bandwidth half. The sampling pipeline. In industrial training (PinSage through today's DGL and PyG stacks), GPU utilization is set by whether CPU-side neighbor sampling and feature gathering keep up. The standard shape is a producer-consumer design with pinned-memory prefetch, feature caches keyed by degree, and, increasingly, GPU-side sampling on unified memory. Inference. Full-graph layer-by-layer inference (materialize all nodes' layer-1 outputs, then layer 2, and so on) avoids sampling bias and costs \( O(Lm) \) total rather than per-node neighborhood explosion. Embeddings are then served from a vector index, and the GNN reruns only on cache misses or drifted neighborhoods. Monitoring. The failure modes worth alerting on are graph-specific. They include degree-distribution drift (a new spam campaign changes the fan-out arithmetic), feature staleness on historical embeddings, and homophily drift, since a model trained on a homophilous graph degrades quietly as adversaries rewire around it.

The current research frontier

Four active fronts, with the groups pushing them. First, expressiveness with a budget, meaning subgraph GNNs and their unifying theory (Bevilacqua, Frasca, Maron and collaborators across Technion and Imperial), homomorphism- and substructure-counting views of expressive power, and the question of whether any of it transfers to large sparse graphs. In parallel, the graph neural tangent kernel line begun at CMU connects infinite-width GNNs to kernel regression and gives clean theory testbeds. Second, transformers and foundation-model ambitions, spanning GraphGPS-style hybrids, linear-attention scaling, Graphormer descendants in molecular property services, and a wave of attempts at pretrained "graph foundation models", none of which has yet produced the cross-domain transfer that pretraining produced for text and images. The open question is whether graphs are one modality or many. Third, geometric learning for science, currently the field's clearest wins. These include NequIP and MACE-style equivariant force fields (Harvard, Cambridge) adopted by materials groups, DeepMind's GNoME materials-discovery and GraphCast weather results, and protein-interface work descending from AlphaFold's use of attention over residue graphs. The contested engineering question is equivariant tensor-product networks versus cheaper invariant-plus-frames constructions, with strong results on both sides (EGNN's minimalism from Amsterdam versus the e3nn lineage). Fourth, evaluation and data, with OGB-LSC at hundred-million-node scale, the Long Range Graph Benchmark probing over-squashing claims, temporal benchmarks with hard negatives, and a growing reproducibility literature (Oxford, TU Munich, and others) whose recurring finding, that tuned baselines close most published gaps, is the healthiest pressure the field has.

Open source to read

Ordered roughly by how much a careful read teaches.

• pyg-team/pytorch_geometric, the reference PyTorch GNN library. Open torch_geometric/nn/conv/message_passing.py first. The MessagePassing base class is the aggregate-update template of this page turned into code, and every layer in the library is a subclass overriding message and update.
• dmlc/dgl, the other major framework, stronger on heterogeneous graphs and distributed training. Open python/dgl/nn/pytorch/conv/graphconv.py to see the same GCN expressed in DGL's message/reduce idiom.
• google-deepmind/jraph, minimal JAX graph nets. Open jraph/_src/models.py. The GraphNetwork function is the Battaglia et al. relational template (edge, node, global blocks) in about a hundred lines of segment-sum code, the cleanest implementation of the abstraction anywhere.
• rusty1s/pytorch_scatter, the scatter/segment kernels underneath PyG. Open torch_scatter/scatter.py for the dispatch logic, then the CUDA sources to see how segmented reductions avoid atomic contention. This is where the memory-bound analysis above becomes concrete.
• snap-stanford/ogb, the Open Graph Benchmark datasets, splits, and evaluators. Open ogb/nodeproppred/dataset_pyg.py to see how the standardized splits are loaded and enforced, and the Evaluator classes to see exactly which metric each leaderboard computes.
• awslabs/dgl-ke, knowledge-graph embeddings at scale (TransE, DistMult, ComplEx, RotatE with sharded async training). Open python/dglke/models/general_models.py for the scoring functions side by side, then the partitioning code for how a 160 GB embedding table trains on commodity machines.
• mir-group/nequip, the E(3)-equivariant force field. Open nequip/nn/_convnetlayer.py to see an equivariant interaction block assembled from e3nn primitives, with the irreps bookkeeping explicit.
• e3nn/e3nn, the equivariance algebra itself. Open e3nn/o3/_tensor_product/_tensor_product.py. The Clebsch-Gordan tensor product that makes "multiply two geometric tensors and stay equivariant" a library call is defined here, and reading it demystifies every equivariant architecture built on top.

Common misconceptions

"More GNN layers means a bigger receptive field means better." The over-smoothing derivation says otherwise. Repeated application of \( \hat S \) converges geometrically to the degree profile, and the measured collapse (0.93 mean cosine similarity by depth 8 on Cora, untrained) arrives long before large receptive fields do. Depth in GNNs buys hops, not hierarchy, and hops beyond the task's interaction radius are pure damage without residuals, jumping knowledge, or normalization.

"GAT is a graph transformer." GAT normalizes attention over the one-hop neighborhood. Its support is the input graph, so it inherits every message-passing pathology including over-squashing, and its original scoring function cannot even express query-dependent attention rankings (the GATv2 static attention result). A graph transformer attends over all pairs and must re-inject structure through encodings. They are different designs with different failure modes.

"Sum, mean, max aggregation are interchangeable details." They have provably different discriminative power. Mean sees the neighbor-feature distribution, max sees the support, and only sum sees the full multiset, which is why GIN needs sum to reach the 1-WL ceiling. Problem 2's four multisets separate the three aggregators with two-dimensional arithmetic.

"GNNs learn graph structure, so they can count substructures like triangles." Anonymous message passing is bounded by 1-WL, which cannot distinguish K3,3 (zero triangles) from the triangular prism (two). No amount of training data fixes an architectural ceiling. Counting requires positional features, subgraph methods, or higher-order networks.

"Node embeddings like node2vec are a learned, superior replacement for spectral methods." The Levy-Goldberg/NetMF analysis shows skip-gram with negative sampling implicitly factorizes a log-PMI matrix built from powers of \( D^{-1} A \), a spectral object approximated by sampling. The methods differ in estimator and loss, not in kind, and inherit the same transductive limits.

"Beating a GCN baseline on Cora-style benchmarks demonstrates a better architecture." The reproducibility literature (Shchur et al., Errica et al., Dwivedi et al.) found ranking reshuffles across splits, structure-agnostic baselines matching GNNs on several graph-classification datasets, and Correct-and-Smooth matching large models on OGB leaderboards. Claims need fixed splits, equal tuning budgets, and strong simple baselines, which is exactly what OGB institutionalized.

"Knowledge-graph embedding quality is whatever MRR the paper reports." Raw versus filtered ranking, tie-breaking policy, and inverse-relation leakage (FB15k, WN18) each move reported numbers by large margins. Several published gains vanished under tie-averaged filtered evaluation. The protocol is part of the result.

"Equivariance is mathematical pedantry, and augmentation learns the same thing." Augmentation teaches approximate invariance at extra sample cost and provides no guarantee off-distribution. An equivariant architecture satisfies the symmetry exactly for every input, and the practical consequence is measured, not aesthetic. NequIP-class force fields reach target accuracy with orders of magnitude fewer expensive quantum-chemistry labels than unconstrained models.

Self-check

References

1. Hamilton, W. L. Graph Representation Learning. Morgan & Claypool, 2020. cs.mcgill.ca/~wlh/grl_book
2. Bronstein, M., Bruna, J., Cohen, T., Veličković, P. Geometric Deep Learning: Grids, Groups, Graphs, Geodesics, and Gauges. 2021. arXiv:2104.13478
3. Chung, F. Spectral Graph Theory. AMS CBMS Regional Conference Series 92, 1997.
4. Perozzi, B., Al-Rfou, R., Skiena, S. "DeepWalk: Online Learning of Social Representations." KDD 2014. arXiv:1403.6652
5. Grover, A., Leskovec, J. "node2vec: Scalable Feature Learning for Networks." KDD 2016. arXiv:1607.00653
6. Qiu, J., Dong, Y., Ma, H., Li, J., Wang, K., Tang, J. "Network Embedding as Matrix Factorization: Unifying DeepWalk, LINE, PTE, and node2vec." WSDM 2018. arXiv:1710.02971
7. Bruna, J., Zaremba, W., Szlam, A., LeCun, Y. "Spectral Networks and Deep Locally Connected Networks on Graphs." ICLR 2014. arXiv:1312.6203
8. Defferrard, M., Bresson, X., Vandergheynst, P. "Convolutional Neural Networks on Graphs with Fast Localized Spectral Filtering." NeurIPS 2016. arXiv:1606.09375
9. Kipf, T. N., Welling, M. "Semi-Supervised Classification with Graph Convolutional Networks." ICLR 2017. arXiv:1609.02907
10. Hamilton, W. L., Ying, R., Leskovec, J. "Inductive Representation Learning on Large Graphs." NeurIPS 2017. arXiv:1706.02216
11. Veličković, P., Cucurull, G., Casanova, A., Romero, A., Liò, P., Bengio, Y. "Graph Attention Networks." ICLR 2018. arXiv:1710.10903
12. Xu, K., Hu, W., Leskovec, J., Jegelka, S. "How Powerful are Graph Neural Networks?" ICLR 2019. arXiv:1810.00826
13. Morris, C., Ritzert, M., Fey, M., Hamilton, W. L., Lenssen, J. E., Rattan, G., Grohe, M. "Weisfeiler and Leman Go Neural: Higher-Order Graph Neural Networks." AAAI 2019. arXiv:1810.02244
14. Gilmer, J., Schoenholz, S., Riley, P., Vinyals, O., Dahl, G. "Neural Message Passing for Quantum Chemistry." ICML 2017. arXiv:1704.01212
15. Li, Q., Han, Z., Wu, X.-M. "Deeper Insights into Graph Convolutional Networks for Semi-Supervised Learning." AAAI 2018. arXiv:1801.07606
16. Oono, K., Suzuki, T. "Graph Neural Networks Exponentially Lose Expressive Power for Node Classification." ICLR 2020. arXiv:1905.10947
17. Alon, U., Yahav, E. "On the Bottleneck of Graph Neural Networks and its Practical Implications." ICLR 2021. arXiv:2006.05205
18. Topping, J., Di Giovanni, F., Chamberlain, B. P., Dong, X., Bronstein, M. "Understanding Over-Squashing and Bottlenecks on Graphs via Curvature." ICLR 2022. arXiv:2111.14522
19. Dwivedi, V. P., Bresson, X. "A Generalization of Transformer Networks to Graphs." 2021. arXiv:2012.09699
20. Ying, C., Cai, T., Luo, S., Zheng, S., Ke, G., He, D., Shen, Y., Liu, T.-Y. "Do Transformers Really Perform Bad for Graph Representation?" (Graphormer). NeurIPS 2021. arXiv:2106.05234
21. Rampášek, L., Galkin, M., Dwivedi, V. P., Luu, A. T., Wolf, G., Beaini, D. "Recipe for a General, Powerful, Scalable Graph Transformer." NeurIPS 2022. arXiv:2205.12454
22. Chiang, W.-L., Liu, X., Si, S., Li, Y., Bengio, S., Hsieh, C.-J. "Cluster-GCN: An Efficient Algorithm for Training Deep and Large Graph Convolutional Networks." KDD 2019. arXiv:1905.07953
23. Zeng, H., Zhou, H., Srivastava, A., Kannan, R., Prasanna, V. "GraphSAINT: Graph Sampling Based Inductive Learning Method." ICLR 2020. arXiv:1907.04931
24. Ying, R., He, R., Chen, K., Eksombatchai, P., Hamilton, W. L., Leskovec, J. "Graph Convolutional Neural Networks for Web-Scale Recommender Systems." (PinSage). KDD 2018. arXiv:1806.01973
25. Bordes, A., Usunier, N., Garcia-Durán, A., Weston, J., Yakhnenko, O. "Translating Embeddings for Modeling Multi-relational Data." (TransE). NeurIPS 2013.
26. Sun, Z., Deng, Z.-H., Nie, J.-Y., Tang, J. "RotatE: Knowledge Graph Embedding by Relational Rotation in Complex Space." ICLR 2019. arXiv:1902.10197
27. Schlichtkrull, M., Kipf, T. N., Bloem, P., van den Berg, R., Titov, I., Welling, M. "Modeling Relational Data with Graph Convolutional Networks." (R-GCN). ESWC 2018. arXiv:1703.06103
28. Satorras, V. G., Hoogeboom, E., Welling, M. "E(n) Equivariant Graph Neural Networks." ICML 2021. arXiv:2102.09844
29. Batzner, S., Musaelian, A., Sun, L., Geiger, M., Mailoa, J. P., Kornbluth, M., Molinari, N., Smidt, T., Kozinsky, B. "E(3)-Equivariant Graph Neural Networks for Data-Efficient and Accurate Interatomic Potentials." (NequIP). Nature Communications 2022. arXiv:2101.03164
30. Pfaff, T., Fortunato, M., Sanchez-Gonzalez, A., Battaglia, P. "Learning Mesh-Based Simulation with Graph Networks." (MeshGraphNets). ICLR 2021. arXiv:2010.03409
31. You, J., Ying, R., Ren, X., Hamilton, W. L., Leskovec, J. "GraphRNN: Generating Realistic Graphs with Deep Auto-regressive Models." ICML 2018. arXiv:1802.08773
32. Hu, W., Fey, M., Zitnik, M., Dong, Y., Ren, H., Liu, B., Catasta, M., Leskovec, J. "Open Graph Benchmark: Datasets for Machine Learning on Graphs." NeurIPS 2020. arXiv:2005.00687
33. Dwivedi, V. P., Joshi, C. K., Luu, A. T., Laurent, T., Bengio, Y., Bresson, X. "Benchmarking Graph Neural Networks." 2020 (JMLR 2023). arXiv:2003.00982
34. Errica, F., Podda, M., Bacciu, D., Micheli, A. "A Fair Comparison of Graph Neural Networks for Graph Classification." ICLR 2020. arXiv:1912.09893