Graph neural networks: message passing, GCN and expressivity

Before this lesson: lesson 01 (the message and update functions are MLPs) and lesson 09 (self-attention is permutation-equivariant). From the linear-algebra track: lesson 05 (eigenvalues, power iteration) and lesson 19 (the graph Laplacian, the GCN layer derived from spectral filters, the eigenvalues of $\hat A$). Time: about 85 minutes, plus the notebook.

---

Why this matters #

Many objects in science and industry are graphs. A molecule is a graph of atoms and bonds. A crystal is a graph of atoms within a cutoff distance. A weather or finite-element mesh, a road network, a social network and a network of seismic stations are graphs too. GNNs are the core of machine-learned interatomic potentials (SchNet, NequIP, MACE), of learned mesh simulators (MeshGraphNets), of the weather model GraphCast, and of many recommender systems. Interviews ask three standard questions about GNNs: why equivariance, how expressive they are (the WL test), and why deep GNNs fail.

1. Predict first #

2. The picture #

One round of message passing #

Take a path of three nodes, 1 — 2 — 3, with one number per node: $h = (6, 0, 0)$. In one round, every node takes the mean of its own number and its neighbours' numbers:

• node 1 sees $\{6, 0\}$ and gets $3$;
• node 2 sees $\{6, 0, 0\}$ and gets $2$;
• node 3 sees $\{0, 0\}$ and gets $0$.

After round 1, $h = (3, 2, 0)$. Node 3 has not seen node 1 yet: they are 2 edges apart. Round 2 gives $h = (2.5,\ 5/3,\ 1)$. Now node 3 depends on node 1. In general, after $K$ rounds a node depends on its $K$-hop neighbourhood: all nodes at most $K$ edges away. This is the receptive field of the node, as in a CNN (lesson 07). A real GNN layer also multiplies by learned weights and applies a nonlinearity, but information still moves one edge per layer.

Interactive figure: loading…

Compared with a CNN and with attention #

A $3\times3$ convolution is message passing on a pixel grid. There, each neighbour has a fixed position (up, left, …) and its own weight. On a general graph the neighbours have no order, and their number changes from node to node. So the layer must treat them as a multiset: a set in which an element can appear more than once. It must combine them with a symmetric function, that is, a function that gives the same result for every order of its inputs: sum, mean or max. This one constraint shapes the whole lesson.

Self-attention (lesson 09) is message passing on the complete graph: every token sends a message to every other token, with learned weights. A GNN sends messages only along the edges that come with the data. Graph attention (GAT) is attention masked by the adjacency matrix.

Notation #

A graph has $n$ nodes. The adjacency matrix $A \in \{0,1\}^{n\times n}$ has $A_{ij} = 1$ when there is an edge from $j$ to $i$. For an undirected graph, $A$ is symmetric. The neighbours of node $i$ are $\mathcal N(i) = \{j : A_{ij} = 1\}$, and its degree is $d_i = |\mathcal N(i)|$. Node features form a matrix $H \in \mathbb R^{n\times C}$: row $h_i$ belongs to node $i$, and there are $C$ channels (one row per node, as one row per example in a mini-batch).

In code, a sparse graph is stored as edge_index: a $2\times E$ integer array with one column (source $j$, target $i$) per directed edge. An undirected edge appears twice, once in each direction. For the path 1 — 2 — 3 with nodes numbered 0, 1, 2:

edge_index = [[0, 1, 1, 2],    # sources j
              [1, 0, 2, 1]]    # targets i

3. Definitions and main results #

Message passing #

One layer of a message-passing neural network (MPNN) (Gilmer et al. 2017) computes, for every node $i$:

• $\psi$ is the message function, usually an MLP. $e_{ji}$ is an optional edge feature (a bond type, a distance).
• $\bigoplus$ is the aggregation: a function of the multiset of messages that does not depend on their order (sum, mean, max, or an attention-weighted sum).
• $\phi$ is the update function, an MLP or a GRU.

The same $\psi$ and $\phi$ act at every node and every edge. So the number of parameters does not depend on the size of the graph, and one layer costs $O(E)$ message evaluations, where $E$ is the number of edges.

How it is computed. Gather $h_i$ and $h_j$ for every edge, compute all messages in one batch, then add each message into the row of its target node. In PyTorch this last step is out.index_add_(0, target, messages); in NumPy it is np.add.at(out, target, messages). PyTorch Geometric's MessagePassing class has the same three parts (message, aggregate, update).

Two kinds of tasks. For a node-level task (classify each paper), put a small head on each final state $h_i^{(K)}$. For a graph-level task (predict the energy of a molecule), first apply a readout that turns all node states into one vector, such as $y_G = \rho\big(\sum_i h_i^{(K)}\big)$.

The GCN layer #

The graph convolutional network (GCN) layer of Kipf and Welling (2017) is

with a weight matrix $W \in \mathbb R^{C\times F}$ and an elementwise nonlinearity $\sigma$. Row by row:

So GCN is message passing with a linear message $h_jW$ scaled by the degrees, a sum over the neighbours and the node itself, and the update $\phi = \sigma$.

• Self-loop ($A + I$): without it, $h_i'$ would not depend on $h_i$, only on the neighbours.
• Normalization. A plain sum $\tilde AH$ grows with the degree. Nodes with many neighbours (hubs) get large features, and stacked layers blow up. The eigenvalues of $\hat A$ lie in $(-1, 1]$ (linear-algebra lesson 19), so repeated propagation stays bounded. The row-normalized $\tilde D^{-1}\tilde A$ (the mean over the node and its neighbours) has the same eigenvalues, because $\tilde D^{-1}\tilde A = \tilde D^{-1/2}\hat A\tilde D^{1/2}$ is similar to $\hat A$.
• Low-pass filter. $\hat A = I - \tilde L_{\text{sym}}$, where $\tilde L_{\text{sym}}$ is the normalized Laplacian of the graph with self-loops. A graph-Fourier component with Laplacian eigenvalue $\tilde\lambda \in [0, 2)$ is multiplied by $1 - \tilde\lambda$. Smooth components (small $\tilde\lambda$) pass. All others shrink, and the roughest ones ($\tilde\lambda > 1$) also flip sign. Linear-algebra lesson 19 derives the layer as a first-order Chebyshev filter.
• Cost. With a sparse $\hat A$, one layer costs $O(nCF + EF)$, and it has $CF$ weights for a graph of any size.

Permutation equivariance and invariance #

The numbering of the nodes is arbitrary. Node 1 could be called node 5. A good model should not depend on this choice. Renumber the nodes and the node outputs should only be renumbered in the same way (equivariance); a single graph-level output should not change at all (invariance).

Interactive figure: loading…

A relabelling is a permutation $\pi$: old node $i$ becomes new node $\pi(i)$. Its permutation matrix $P$ has $P_{\pi(i),i} = 1$ and zeros elsewhere. Then $(PH)_{\pi(i)} = h_i$ (the rows are moved), and the new adjacency matrix $A' = PAP^\top$ has $A'_{\pi(i)\pi(j)} = A_{ij}$ (the edges are moved with the nodes).

A map $F$ is permutation-equivariant if $F(PH, PAP^\top) = P\,F(H, A)$ for every permutation matrix $P$. A map $R$ to a single vector is permutation-invariant if $R(PH, PAP^\top) = R(H, A)$.

Theorem (equivariance of message passing). Let $F$ be a message-passing layer whose aggregation depends only on the multiset of messages (no edge features, for simplicity). Then:

• (a) $F(PH, PAP^\top) = P\,F(H, A)$ for every permutation matrix $P$;
• (b) a stack of such layers is equivariant;
• (c) the sum readout $R(H) = H^\top\mathbf 1$ satisfies $R(PH) = R(H)$, so "readout after stack" is invariant;
• (d) for GCN: $\hat A' = P\hat AP^\top$, and $\sigma(\hat A'\,PH\,W) = P\,\sigma(\hat AHW)$.

Why this is true, in one sentence: the layer only ever asks "who are my neighbours and what do they send?", and a relabelling changes neither answer, only the names.

Proof.

1. In the new graph, node $\pi(i)$ carries $h_i$, and its neighbours are $\pi(\mathcal N(i))$, because $A'_{\pi(i)\pi(j)} = A_{ij}$.
2. The message from $\pi(j)$ to $\pi(i)$ is $\psi(h_i, h_j) = m_{j\to i}$. So node $\pi(i)$ receives exactly the multiset $\{\{m_{j\to i} : j \in \mathcal N(i)\}\}$, maybe in another order.
3. The aggregation sees only the multiset. So the new aggregate at $\pi(i)$ is $a_i$, and the new output there is $\phi(h_i, a_i) = h_i'$. Row $\pi(i)$ of the new output equals row $i$ of the old output. This is (a).
4. (b): take two layers $F_1$ and $F_2$. By (a) for $F_1$, $F_1(PH, PAP^\top) = PF_1(H, A)$. Then (a) for $F_2$, with input $F_1(H, A)$, gives $F_2(PF_1(H, A), PAP^\top) = PF_2(F_1(H, A), A)$. Induction over the layers gives any depth.
5. (c): $(PH)^\top\mathbf 1 = H^\top P^\top\mathbf 1 = H^\top\mathbf 1$, because every row of $P^\top$ has exactly one 1, so $P^\top\mathbf 1 = \mathbf 1$. Mean and max readouts are invariant for the same reason: they see only the multiset of rows.
6. (d): $\tilde A' = P\tilde AP^\top$ (because $PIP^\top = I$). The degrees are permuted: $\tilde d' = P\tilde AP^\top\mathbf 1 = P\tilde d$, so $\tilde D'^{-1/2} = P\tilde D^{-1/2}P^\top$. Hence $\hat A' = P\hat AP^\top$, and $\sigma(P\hat AP^\top PHW) = P\sigma(\hat AHW)$, because $P^\top P = I$ and an elementwise $\sigma$ commutes with reordering rows. $\blacksquare$

Equivariance builds the symmetry into the model, as translation equivariance does for a CNN. Node outputs move with the nodes; only a readout is invariant. (Prince's book stores nodes as columns; his Eq. 13.6 is the same statement.)

Expressivity: the 1-WL test #

How powerful can a message-passing GNN be? The answer comes from an old test for graph isomorphism. Two graphs are isomorphic if one is a relabelling of the other.

1-WL colour refinement (Weisfeiler and Leman). Give every node a colour: its input label, or one common colour if there are no features. Then repeat:

with one injective HASH for both graphs ("injective" means different inputs always get different colours). The new colour of a node encodes its old colour and the multiset of its neighbours' colours. Then compare the colour histograms: the number of nodes of each colour in each graph. If at some round the histograms differ, the graphs are not isomorphic. If they stay equal, 1-WL cannot distinguish the graphs; they may or may not be isomorphic.

Interactive figure: loading…

Theorem (Xu et al. 2019, Lemma 2; Morris et al. 2019). If 1-WL gives two graphs the same colour histograms in every round, then every MPNN whose input features are a function of the initial colours gives them the same graph-level output.

Why this is true: an MPNN layer is one round of 1-WL in which HASH is replaced by a learned function that may merge colours but never splits them. So it can never see more than 1-WL.

Proof.

1. Claim: after $t$ layers, $h_i^{(t)} = g_t(c_i^{(t)})$ for some function $g_t$, for every node of both graphs. In words: equal colours imply equal embeddings.
2. For $t = 0$ this holds by the assumption on the input features.
3. Step from $t$ to $t + 1$. Let $c_i^{(t+1)} = c_k^{(t+1)}$ for nodes $i$, $k$ in either graph. HASH is injective, so $c_i^{(t)} = c_k^{(t)}$ and the multisets of neighbour colours are equal. By the claim for $t$, $h_i^{(t)} = h_k^{(t)}$ and the multisets of neighbour embeddings are equal. So the multisets of messages, the aggregates and the outputs of $\phi$ are equal: $h_i^{(t+1)} = h_k^{(t+1)}$.
4. Equal colour histograms then give equal multisets of final embeddings, hence equal sum, mean or max readouts. $\blacksquare$

Example. In two triangles and in a 6-ring, every node has degree 2. Round 1 gives every node of both graphs the colour HASH(0, {{0, 0}}), and no colour ever splits. So every MPNN gives both graphs the same output, although only one is connected and only one has triangles. In chemistry, decalin (two 6-rings that share a bond) and bicyclopentyl (two 5-rings joined by a bond) are also indistinguishable by 1-WL (checked in section 5). Plain MPNNs cannot count cycles in general.

Ways out: add structural or positional features (cycle counts, Laplacian eigenvectors), use higher-order GNNs ($k$-WL, cost $O(n^k)$), add random node identifiers (exact equivariance is lost), or use geometry (3-D positions; see the deeper box at the end of this section).

Sum, mean or max? (GIN) #

Simplest case: one-hot features. Let each neighbour carry the one-hot vector of its type, $a = (1, 0)$ or $b = (0, 1)$. Then:

• the sum is the count vector, which is the multiset: it loses nothing (it is injective);
• the mean is the count vector divided by the degree: it keeps the proportions and loses the size;
• the max is the 0/1 indicator of the types present: it keeps the set and loses the counts.

Interactive figure: loading…

Xu et al. (2019) prove the general version. For countable features, there is a function $f$ such that $X \mapsto \sum_{x\in X}f(x)$ is injective on bounded multisets (Lemma 5). Mean captures only the distribution of the multiset, and max only its underlying set (Corollaries 8 and 9). A GNN with injective aggregation, update and readout is as powerful as 1-WL (Theorem 3). Their graph isomorphism network (GIN) meets these conditions:

with $\varepsilon$ learned or fixed. The factor $1 + \varepsilon$ keeps the node's own state apart from the neighbour sum (Corollary 6). The MLP matters: one linear map followed by a ReLU, as in GCN, is not enough (§5.1 of the paper). GCN and GraphSAGE (mean or max, with one-layer perceptrons) are strictly less powerful than 1-WL. Mean is still a good choice when proportions matter and the features are rich. Sum is the choice when counts matter: the number of hydrogen atoms, a total energy.

Over-smoothing #

Simplest case: drop the weights and the nonlinearity. Then $K$ layers compute $\hat A^Kx$ for a feature column $x \in \mathbb R^n$. Each step is a low-pass filter. What is left after many steps?

Interactive figure: loading…

Theorem (over-smoothing of linear propagation). Let the graph be connected. Write $\hat A = \sum_{i=1}^n\lambda_iu_iu_i^\top$ with orthonormal eigenvectors $u_i$ and $\lambda_1 = 1$. Let $\rho = \max_{i\ge2}|\lambda_i|$. Then $\rho < 1$, $u_1 = \tilde D^{1/2}\mathbf 1/\|\tilde D^{1/2}\mathbf 1\|$, and

In the limit, node $i$ holds $\sqrt{\tilde d_i}\cdot\sum_j\sqrt{\tilde d_j}\,x_j\big/\sum_j\tilde d_j$: one common number times $\sqrt{\tilde d_i}$.

Why this is true: repeated multiplication by $\hat A$ is power iteration (linear-algebra lesson 05). Every eigen-direction is multiplied by $\lambda_i^K$; only the direction with $\lambda_1 = 1$ survives.

Proof (four steps)

1. $\hat A$ is symmetric, so it has real eigenvalues and orthonormal eigenvectors. All eigenvalues lie in $(-1, 1]$ (linear-algebra lesson 19).
2. $\hat A\tilde D^{1/2}\mathbf 1 = \tilde D^{-1/2}\tilde A\mathbf 1 = \tilde D^{-1/2}\tilde d = \tilde D^{1/2}\mathbf 1$. So $\tilde D^{1/2}\mathbf 1$ is an eigenvector with eigenvalue 1.
3. The eigenvalue 1 is simple. With the Laplacian $L = D - A = \tilde D - \tilde A$, we have $I - \hat A = \tilde D^{-1/2}L\tilde D^{-1/2}$. So $\hat Ax = x$ $\iff$ $L(\tilde D^{-1/2}x) = 0$ $\iff$ $\tilde D^{-1/2}x$ is constant (for a connected graph, $\mathbf 1$ spans the null space of $L$) $\iff$ $x \propto \tilde D^{1/2}\mathbf 1$. With step 1, every other eigenvalue has $|\lambda_i| < 1$, so $\rho < 1$.
4. $\hat A^Kx = \sum_i\lambda_i^Ku_i(u_i^\top x)$. The $i = 1$ term is $u_1u_1^\top x$ for every $K$. Every other coefficient shrinks by a factor $|\lambda_i|^K \le \rho^K$, and orthogonality gives $\|\hat A^Kx - u_1u_1^\top x\|^2 = \sum_{i\ge2}\lambda_i^{2K}(u_i^\top x)^2 \le \rho^{2K}\|x - u_1u_1^\top x\|^2$. $\blacksquare$

What it means. Every feature column tends to a multiple of $(\sqrt{\tilde d_i})_i$, so the nodes differ only by their degree. With mean aggregation, every node tends to the same value. Dense, well-mixed graphs (small $\rho$) over-smooth in a few layers.

• Weights do not save you. A linear GCN computes $\hat A^KXW_1\cdots W_K$: $\hat A$ acts on the nodes and the weight matrices $W_l$ act on the channels, so the weights cannot undo the filter over the nodes. With ReLU, Oono and Suzuki (2020) show that the distance to the "degree-only" subspace is $O((s\rho)^K)$, with $s$ the largest singular value of the weights. When $s\rho < 1$, the output keeps only degree (and connected-component) information.
• Smoothing is also the signal. When neighbours usually share the label (homophily), a few rounds remove noise and keep the class. This is why a 2-layer GCN works on citation graphs. Section 5 shows both sides on one graph.
• Fixes: residual connections $H' = H + \sigma(\hat AHW)$, initial residuals (GCNII), normalization (PairNorm), DropEdge, or only 2–3 layers.

A common measure of smoothness is the Dirichlet energy $\mathcal E(H) = \operatorname{tr}\big(H^\top(I - \hat A)H\big) = \tfrac12\sum_{i,j}\tilde A_{ij}\big\|h_i/\sqrt{\tilde d_i} - h_j/\sqrt{\tilde d_j}\big\|^2$. It is 0 exactly on the over-smoothed limit and decays to 0 with depth.

Over-squashing #

After $K$ layers, node $i$ sees its $K$-hop neighbourhood. On trees and small-world graphs this neighbourhood grows exponentially with $K$. All of it must fit into one fixed-size vector $h_i \in \mathbb R^F$, and it must pass through the few edges that join distant regions (the bottlenecks). This is over-squashing (Alon and Yahav 2021).

Topping et al. (2022, Lemma 1) bound it for an MPNN that aggregates with the weights $\hat A_{ij}$. If the derivatives of $\phi$ and $\psi$ are bounded by $\alpha$ and $\beta$, a node $s$ at distance $r + 1$ from $i$ satisfies

Here $x_s$ is the input feature of node $s$. A walk of length $r + 1$ between nodes at distance $r + 1$ must be a shortest path. So $(\hat A^{r+1})_{is}$ is a sum, over the shortest paths, of products of $r + 1$ factors $1/\sqrt{\tilde d_u\tilde d_v}$. The bound is small when few paths join $i$ and $s$. In a tree there is only one path, and the bound decays exponentially in $r$.

	Over-smoothing	Over-squashing
What goes wrong	node states become too similar	distant information cannot get through
Cause	repeated averaging (a low-pass filter)	exponential receptive field, fixed-size vectors, bottleneck edges
Measured by	distance to the smooth subspace, $\rho^K$; Dirichlet energy	the Jacobian $\partial h_i/\partial x_s$
Fixes	residuals, GCNII, PairNorm, fewer layers	rewiring (add edges), a virtual node, a fully connected last layer, global attention, multi-scale meshes

Graph attention #

Graph attention (GAT) (Veličković et al. 2018) replaces the fixed weights $\hat A_{ij}$ by learned attention weights. With $z_i = h_iW$, a score $e_{ij} = \operatorname{LeakyReLU}(a^\top[z_i \,\|\, z_j])$ ($\|$ is concatenation, $a$ a learned vector), weights $\alpha_{ij} = \operatorname{softmax}_{j \in \mathcal N(i)\cup\{i\}}(e_{ij})$, and output $h_i' = \sigma\big(\sum_j \alpha_{ij}z_j\big)$. This is self-attention whose softmax runs only over the neighbours: attention masked by $A + I$. GAT scores pairs with a small additive network, not with a scaled dot product, but the idea is the same. A transformer is the special case of the complete graph, with positional encodings (lesson 10) in place of the graph structure.

E(3)-equivariant GNNs: what geometry adds

Atoms also have positions $x_i \in \mathbb R^3$. Rotating, reflecting or translating the whole system (the group E(3)) leaves the energy unchanged and rotates the forces. A plain MPNN on the bond graph ignores geometry. Raw coordinates as input features break the symmetry, which the model must then learn from data, and only approximately. There are two ways to build the symmetry in:

• Invariant models use only invariant inputs in the messages: distances $\|x_i - x_j\|$ (SchNet), plus angles (DimeNet).
• Equivariant models also carry vector or tensor features that rotate with the input. EGNN (Satorras et al. 2021) uses $m_{ij} = \psi(h_i, h_j, \|x_i - x_j\|^2)$, $x_i' = x_i + C\sum_{j\ne i}(x_i - x_j)\,\phi_x(m_{ij})$ and $h_i' = \phi(h_i, \sum_{j\ne i}m_{ij})$. The vectors $x_i - x_j$ rotate with the system and ignore translations, and $\phi_x(m_{ij})$ is invariant, so $x'$ is equivariant. NequIP and MACE use features labelled by the irreducible representations of SO(3) (spherical harmonics of degree $\ell$), combined by Clebsch–Gordan tensor products.

What this buys: exact symmetry (an invariant energy $E$, and forces $F = -\partial E/\partial X$ from autodiff that are equivariant automatically; proof in linear-algebra lesson 19); data efficiency (NequIP reports better accuracy than earlier models with up to three orders of magnitude less training data); and information that the bond graph lacks (distances, angles, orientation). The ML-for-science track covers these models in depth.

Training a small GCN #

Kipf and Welling's model for semi-supervised node classification is

with the cross-entropy over the labelled nodes only, and one forward pass over the full graph per step (16 hidden units, dropout 0.5, weight decay $5\cdot10^{-4}$, Adam with learning rate 0.01). On Cora (2,708 papers, 7 classes, 20 labels per class) it reaches 81.5% test accuracy, against 55.1% for an MLP on the features alone (their Table 3). The setting is transductive: the unlabelled nodes are in the graph during training, so their features flow through $\hat A$ into the predictions of the labelled nodes. An inductive model must handle new nodes or new graphs at test time.

Answers to section 1. (1) No. Every node of both graphs has degree 2, so 1-WL gives every node the same colour in every round, and the theorem gives equal outputs. For GCN you can also see it directly: $\tilde d_i = 3$ everywhere and $\hat A\mathbf 1 = \mathbf 1$, so a constant feature stays constant. You need extra features (triangle counts), a more expressive model, or geometry. (2) The accuracy drops. Each layer is a low-pass filter, and 32 of them push all node states towards $\tilde D^{1/2}\mathbf 1$ times a common vector (over-smoothing). A deep plain stack is also hard to optimize. Kipf and Welling found 2–3 layers best, and beyond about 7 layers training was difficult without residual connections (their Appendix B).

4. Worked example #

The star graph: a hub with three leaves. Node 1 is the centre, joined to leaves 2, 3 and 4. (A "virtual node" joined to every node is such a hub.)

Step 1: degrees with self-loops. The centre has 3 neighbours plus itself: $\tilde d_1 = 4$. Each leaf has 1 neighbour plus itself: $\tilde d = (4, 2, 2, 2)$.

Step 2: the matrix $\hat A$. Each entry is $\hat A_{ij} = \tilde A_{ij}/\sqrt{\tilde d_i\tilde d_j}$. Centre–centre: $1/\sqrt{4\cdot4} = 1/4$. Centre–leaf: $1/\sqrt{4\cdot2} = 1/\sqrt8$. Leaf–itself: $1/\sqrt{2\cdot2} = 1/2$. Leaf–leaf: 0 (no edge).

Step 3: one propagation of the feature $h = (0, 1, 2, 3)^\top$, with three aggregations:

• GCN: centre $= 0/4 + (1 + 2 + 3)/\sqrt8 = 6/\sqrt8 = 2.121$; leaf 2 $= 0/\sqrt8 + 1/2 = 0.5$; so $\hat Ah = (2.121,\ 0.5,\ 1.0,\ 1.5)$.
• Plain sum with self-loops: $\tilde Ah = (6, 1, 2, 3)$. The centre value grows with its degree.
• Mean with self-loops: $\tilde D^{-1}\tilde Ah = (6/4, 1/2, 2/2, 3/2) = (1.5, 0.5, 1.0, 1.5)$.

Step 4: the spectrum.

• $\lambda = 1$, with eigenvector $\tilde D^{1/2}\mathbf 1 = (2, \sqrt2, \sqrt2, \sqrt2)$. Check the centre row: $2/4 + 3\sqrt2/\sqrt8 = 1/2 + 3/2 = 2$ ✓. A leaf row: $2/\sqrt8 + \sqrt2/2 = \sqrt2$ ✓.
• $\lambda = 1/2$, twice: differences between leaves with the centre at 0, such as $(0, 1, -1, 0)$ and $(0, 1, 1, -2)$. The centre receives $(1 - 1)/\sqrt8 = 0$, and each leaf keeps half its value.
• The last eigenvalue from the trace: $1/4 + 3/2 = 1 + 1/2 + 1/2 + \lambda_4$, so $\lambda_4 = -1/4$, with eigenvector $(3, -\sqrt2, -\sqrt2, -\sqrt2)$.
• So $\rho = \max(1/2, 1/4) = 1/2$.

Step 5: over-smoothing. $u_1 = (2, \sqrt2, \sqrt2, \sqrt2)/\sqrt{10}$, and $u_1^\top h = (0 + \sqrt2 + 2\sqrt2 + 3\sqrt2)/\sqrt{10} = 6\sqrt2/\sqrt{10}$. So

Every leaf ends at 1.2, whatever its own value was. The centre ends at $1.2\sqrt2$ (the ratio $\sqrt{\tilde d_1/\tilde d_2} = \sqrt2$). The rest of $h$ is $(0, -1, 0, 1)$ in the $\lambda = 1/2$ space (squared norm 2), plus a part with squared norm 4.8 in the $\lambda = -1/4$ space. So the error after $K$ steps is $\sqrt{2\cdot0.25^K + 4.8\cdot0.0625^K}$: 0.894 at $K = 1$, 0.044 at $K = 5$, and 0.0014 at $K = 10$. With mean aggregation, all four nodes tend to $\sum_j\tilde d_jh_j/\sum_j\tilde d_j = (0 + 2 + 4 + 6)/10 = 1.2$.

Step 6: multisets. One-hot types $a = (1, 0)$ and $b = (0, 1)$:

neighbours	sum	mean	max
$\{a, b\}$	(1, 1)	(0.5, 0.5)	(1, 1)
$\{a, a, b, b\}$	(2, 2)	(0.5, 0.5)	(1, 1)
$\{a, a, b\}$	(2, 1)	(0.667, 0.333)	(1, 1)

Sum separates all three. Mean confuses the first two (same proportions). Max confuses all three (same set).

Step 7: parameters. Kipf and Welling's Cora model has $1433\cdot16 + 16\cdot7 = 23{,}040$ weights (their Eq. 9 has no biases; 23,063 with biases). The same weights run on a graph of any size.

5. In code #

NumPy checks the worked example, the 1-WL examples, equivariance of masked attention, and shows smoothing at work on a 3-community graph. The notebook for this lesson has these cells and the exercises.

import numpy as np
from sklearn.linear_model import LogisticRegression
from sklearn.preprocessing import StandardScaler, normalize
np.set_printoptions(precision=3, suppress=True)

# 1. Worked example: the star graph (centre 0, leaves 1, 2, 3)
A = np.zeros((4, 4))
A[0, 1:] = A[1:, 0] = 1
A_tilde = A + np.eye(4)                                  # self-loops
d = A_tilde.sum(1)                                       # degrees with self-loops: [4. 2. 2. 2.]
A_hat = A_tilde / np.sqrt(np.outer(d, d))                # D^-1/2 (A + I) D^-1/2
h = np.array([0., 1., 2., 3.])
print(A_hat @ h)                                         # [2.121 0.5   1.    1.5  ]
print(np.linalg.eigvalsh(A_hat))                         # [-0.25  0.5   0.5   1.  ]
u1 = np.sqrt(d) / np.linalg.norm(np.sqrt(d))             # top eigenvector, D^1/2 1 normalized
limit = u1 * (u1 @ h)
print(limit)                                             # [1.697 1.2   1.2   1.2  ]
print([round(float(np.linalg.norm(np.linalg.matrix_power(A_hat, K) @ h - limit)), 4)
       for K in (1, 5, 10)])                             # [0.8944, 0.0442, 0.0014]

# 2. 1-WL colour refinement: equal colour histograms => no MPNN separates the graphs
def wl_histogram(edges, n, rounds=4):
    nbrs = [[] for _ in range(n)]
    for a, b in edges:
        nbrs[a].append(b)
        nbrs[b].append(a)
    colour = [0] * n                                     # no node features: one colour
    for _ in range(rounds):
        colour = [hash((colour[i], tuple(sorted(colour[j] for j in nbrs[i])))) for i in range(n)]
    return sorted(colour)

ring = lambda k, o=0: [(o + i, o + (i + 1) % k) for i in range(k)]
decalin = ring(6) + [(1, 6), (6, 7), (7, 8), (8, 9), (9, 0)]    # two 6-rings sharing the bond 0-1
bicyclopentyl = ring(5) + ring(5, 5) + [(0, 5)]                 # two 5-rings joined by the bond 0-5
print(wl_histogram(ring(3) + ring(3, 3), 6) == wl_histogram(ring(6), 6),
      wl_histogram(decalin, 10) == wl_histogram(bicyclopentyl, 10))     # True True

# 3. Sum, mean and max over one-hot neighbour features a = (1, 0), b = (0, 1)
a, b = np.eye(2)
for name, X in {"{a,b}": [a, b], "{a,a,b,b}": [a, a, b, b], "{a,a,b}": [a, a, b]}.items():
    X = np.array(X)
    print(name, X.sum(0), X.mean(0).round(3), X.max(0))
# {a,b} [1. 1.] [0.5 0.5] [1. 1.]
# {a,a,b,b} [2. 2.] [0.5 0.5] [1. 1.]
# {a,a,b} [2. 1.] [0.667 0.333] [1. 1.]

# 4. Graph attention = self-attention masked by A + I: equivariant if nodes and mask move together
rng = np.random.default_rng(0)
n, C = 8, 4
adj = rng.random((n, n)) < 0.3
adj = adj | adj.T | np.eye(n, dtype=bool)                # undirected, with self-loops
Wq, Wk, Wv = rng.standard_normal((3, C, C))
def masked_attention(X, adj):
    S = (X @ Wq) @ (X @ Wk).T / np.sqrt(C)
    S = np.where(adj, S, -np.inf)                        # a node only attends to its neighbours
    P = np.exp(S - S.max(1, keepdims=True))
    return (P / P.sum(1, keepdims=True)) @ (X @ Wv)
W_mlp = rng.standard_normal((n * C, n * C))              # a dense layer on the flattened node table
mlp = lambda X: (X.reshape(-1) @ W_mlp).reshape(n, C)
X = rng.standard_normal((n, C))
perm = rng.permutation(n)
err_att = np.abs(masked_attention(X[perm], adj[perm][:, perm]) - masked_attention(X, adj)[perm]).max()
err_mlp = np.linalg.norm(mlp(X[perm]) - mlp(X)[perm]) / np.linalg.norm(mlp(X))
print(f"{err_att:.1e}", round(float(err_mlp), 2))       # 4.4e-16 1.33

# 5. Smoothing on a 3-community graph: it helps, then it over-smooths
rng = np.random.default_rng(0)
n, k = 300, 3
y = np.repeat(np.arange(k), n // k)
U = np.triu(rng.random((n, n)) < np.where(y[:, None] == y[None, :], 0.08, 0.005), 1)
A = (U | U.T).astype(float)
X = 0.35 * rng.standard_normal((k, 16))[y] + rng.standard_normal((n, 16))   # weak features
train = np.zeros(n, bool)
for c in range(k):
    train[rng.choice(np.flatnonzero(y == c), 10, replace=False)] = True
P = normalize(A + np.eye(n), norm="l1")                  # mean over the node and its neighbours
for K in (0, 1, 2, 8, 32, 128, 512):
    Z = np.linalg.matrix_power(P, K) @ X
    spread = np.linalg.norm(Z - Z.mean(0)) / np.linalg.norm(Z)
    acc = [LogisticRegression(max_iter=5000).fit(F[train], y[train]).score(F[~train], y[~train])
           for F in (Z, StandardScaler().fit_transform(Z))]
    print(K, f"{spread:.1e}", np.round(acc, 3))
# 0 9.8e-01 [0.504 0.493]
# 1 8.8e-01 [0.93 0.9 ]
# 2 7.6e-01 [0.993 0.996]
# 8 4.1e-01 [0.993 0.993]
# 32 2.5e-02 [0.993 0.993]
# 128 6.8e-07 [0.333 0.989]
# 512 4.1e-15 [0.333 0.333]

Part 4: masked attention is equivariant up to round-off; the dense layer is not. Part 5 prints, for each $K$, the relative spread between nodes and the test accuracy of a logistic regression on the raw and on the standardized features. 89% of the edges join nodes of the same class, so two rounds of averaging lift the accuracy from 0.50 to 0.99. Then the spread shrinks like $0.9^K$ (0.9 is the second eigenvalue of $P$, which comes from the communities). At $K = 128$ the differences between nodes are about $10^{-6}$ of the common part, and the classifier on raw features predicts one class. After standardization the class signal survives until round-off erases it at $K = 512$.

6. Where it shows up #

• Molecules and materials. Interatomic potentials (SchNet, NequIP, MACE) predict the energy as a sum of atomic terms (a sum readout) and the forces by autodiff. They run molecular dynamics far more cheaply than density functional theory (DFT).
• Weather and climate. GraphCast has an encoder (from the 0.25° grid to a mesh), a processor of 16 GNN layers with unshared weights on a multi-mesh, and a decoder (from the mesh back to the grid). The multi-mesh is an icosahedron refined 6 times (40,962 nodes), with the edges of all refinement levels. Coarse edges carry information far in a few hops: a fix for over-squashing. ECMWF's AIFS uses GNN encoders and decoders around a transformer processor.
• PDE surrogates. MeshGraphNets (Pfaff et al. 2021) use encode–process–decode with 15 message-passing steps on simulation meshes.
• Seismology. Van den Ende and Ampuero (2020) locate earthquakes and estimate magnitudes with a GNN over the station network, for any number of stations in any order.
• Recommendation and citation graphs. Node classification and link prediction with GCN, GraphSAGE (sampled neighbours, inductive) and GAT.

7. Common mistakes #

• Wrong edge direction. Messages flow from source to target. Aggregating by the wrong row of edge_index reverses the graph. An undirected edge must appear in both directions.
• No self-loops or no normalization. Without self-loops, a node ignores its own state. Without normalization, features scale with the degree and explode with depth.
• A mean readout for an extensive target. A total energy grows with the number of atoms; use a sum readout.
• Stacking plain GCN layers to see further. Depth over-smooths. Use residuals, GCNII, or long-range edges.
• "The MLP inside makes it universal." The 1-WL limit comes from what the aggregation sees. No MLP removes it.
• Breaking equivariance by accident: ordering neighbours by their index, using node IDs as features, or using Laplacian eigenvectors without handling their sign ambiguity.
• Label leakage. In transductive training, the labels of test nodes must not be used, not even for early stopping.

8. Check yourself #

1. Compute $\hat A$ for a triangle (3 nodes, all joined). What does one GCN propagation do?
   Answer

   $\tilde A$ is the all-ones matrix and $\tilde d_i = 3$, so $\hat A = \frac13\mathbf 1\mathbf 1^\top$. Then $\hat Ah = \bar h\,\mathbf 1$: every node gets the mean of all nodes. The eigenvalues are 1, 0, 0, so $\rho = 0$: one propagation over-smooths completely.
2. Give two neighbour multisets that both mean and max confuse but sum separates. Then give two that max confuses but mean separates.
   Answer

   $\{a\}$ and $\{a, a\}$: the mean is $a$ and the max is $a$ for both, but the sums are $a$ and $2a$. $\{a, b\}$ and $\{a, a, b\}$ with $a \ne b$: the same set, so the same max; the means $(a + b)/2$ and $(2a + b)/3$ differ by $(b - a)/6 \ne 0$.
3. A layer concatenates the features of a node's neighbours in increasing node-index order (padded with zeros up to degree 3) and applies an MLP. Is it permutation-equivariant?
   Answer

   No. After a relabelling, the same neighbours appear in a different order, so the MLP sees a permuted input and in general gives a different output: the layer depends on the numbering. A symmetric aggregation (sum, mean, max, attention) fixes this.
4. The non-top eigenvalues of $\hat A$ have $\rho = 0.95$. How many propagation steps until the non-smooth part of $x$ is guaranteed to be below 1% of its initial size? And for the star graph?
   Answer

   $0.95^K \le 0.01 \iff K \ge \ln 0.01/\ln 0.95 = 89.8$: 90 steps. The star has $\rho = 0.5$: 7 steps ($0.5^7 \approx 0.0078$).
5. You add a virtual node joined to all $n$ nodes. Which depth problem does this address, and what is its cost?
   Answer

   Over-squashing and long-range transport: any two nodes are now at most 2 hops apart. The cost: the virtual node is itself one fixed-size vector for all $n$ nodes (a new bottleneck), and the extra edges speed up mixing, which can make over-smoothing faster.
6. You predict the total energy of molecules of different sizes. Sum or mean readout? And how do you get forces that rotate with the molecule?
   Answer

   Sum: energy is extensive. Two far-apart copies of a molecule have twice the energy; a sum readout doubles, a mean does not. For forces, predict an E(3)-invariant energy and take $F = -\partial E/\partial X$ by autodiff; the chain rule makes $F$ rotate with the molecule.

Recall cardsAnswer in your head first, then reveal. Add the cards to your review deck to see them again just before you would forget.

1. Write the three steps of one message-passing layer.

   Show answer

2. Write the GCN layer.

   Show answer

3. Why does GCN add self-loops and use the symmetric normalization?

   Show answer

4. State permutation equivariance of a GNN layer $F(H, A)$.

   Show answer

5. Why is the sum readout $H^\top\mathbf 1$ permutation-invariant?

   Show answer

6. What is the 1-WL limit of message-passing GNNs?

   Show answer

7. Why is sum aggregation more expressive than mean or max?

   Show answer

8. After $K$ GCN propagations on a connected graph, what does a feature column converge to, and how fast?

   Show answer

9. What is over-squashing?

   Show answer

10. Compute $\hat A$ for the triangle graph.

    Show answer

Add these cards to my reviewOpen the review →

9. Interview angle #

• "Why $\tilde D^{-1/2}(A + I)\tilde D^{-1/2}$?" — Self-loops keep the node's own state. The normalization keeps the eigenvalues in $(-1, 1]$, so stacking is stable. It is a first-order spectral filter.
• "How powerful are GNNs?" — At most as powerful as 1-WL. GIN reaches it with sum aggregation and an MLP. Two triangles vs a 6-ring is the standard counterexample. Fixes: structural or positional features, higher-order GNNs, geometry.
• "Why does depth not help GNNs?" — Over-smoothing ($\rho^K$) and over-squashing (the Jacobian bound, bottlenecks), with the fixes for each.
• "GNN vs transformer?" — A transformer is a GNN on the complete graph with attention weights. Positional encodings stand in for the graph structure. GAT is masked attention.
• "Design a model for molecules, or for weather on the sphere." — An E(3)-invariant energy with autodiff forces and equivariant messages; encode–process–decode on a multi-scale mesh.

10. Practice #

The notebook for this lesson runs in Colab. It uses NumPy for the graph operations and PyTorch for training. Five exercises, each with a check cell and a folded solution:

1. Build the normalized adjacency $\hat A$ from edge_index, checked on the star graph.
2. Message passing with scatter-add (np.add.at), checked against the dense product $\hat AHW$.
3. A permutation-equivariance test for your layer, and a test that must fail for a layer that is not equivariant.
4. 1-WL colour refinement on two graphs: separate a path from a triangle plus an edge, and fail on two triangles vs a hexagon.
5. Train a 2-layer GCN on a synthetic 3-community graph, then show over-smoothing as the depth grows.

11. Watch and read #

Theoretical Foundations of Graph Neural NetworksPetar Veličković · 1 h 12 min

Permutation invariance and equivariance from first principles, then message passing, GCN and GAT; the language of sections 2 and 3.

Graph Convolutional Networks - Oxford Geometric Deep LearningFederico Barbero · 14 min

The GCN layer, self-loops and the symmetric normalization, as in section 3.

How powerful are Graph Neural Networks? - Oxford Geometric Deep LearningFederico Barbero · 9 min

The 1-WL test, the limit of message passing and GIN.

Over-squashing and Over-smoothingWelcomeAIOverlords · 10 min

The two depth problems side by side, with the bottleneck picture.

ICLR 2021 Keynote - "Geometric Deep Learning: The Erlangen Programme of ML" - M BronsteinMichael Bronstein · 38 min

Symmetry as the organizing idea: grids, graphs and groups, and where equivariant GNNs fit.

• Simon Prince, Understanding Deep Learning (free PDF), chapter 13: §13.3–13.4 (tasks, equivariance, parameter sharing, a GCN layer), §13.6–13.8 (transductive training; mean, Kipf normalization, max and attention layers), and the chapter notes on the power of GNNs and over-smoothing.
• Sanchez-Lengeling et al., A Gentle Introduction to Graph Neural Networks (Distill, 2021): interactive figures of message passing, pooling and readouts.
• Xu et al., How Powerful are Graph Neural Networks? (2019): §2 (the WL test), §4 (Lemma 5, Corollary 6, GIN), §5 (one-layer perceptrons, mean and max).
• Kipf and Welling, Semi-Supervised Classification with Graph Convolutional Networks (2017): §2 (the propagation rule), §3.1 (the two-layer model), Table 3 (propagation models), Appendix B (depth and residual connections).
• Topping et al., Understanding over-squashing and bottlenecks on graphs via curvature (2022): Lemma 1 (the Jacobian bound).
• PyTorch Geometric, Creating Message Passing Networks: the MessagePassing class and the GCN layer step by step.