Graph Convolutional Networks

Reading time: 55 min | Interview relevance: Very High - GCN is the canonical GNN; every GNN interview starts here | Target roles: ML Engineer, AI Engineer, Research Engineer, Graph ML Engineer

---

The Real Interview Moment

"Derive the GCN layer from first principles," the interviewer says. "I want to understand where the $\tilde{D}^{-1/2} \tilde{A} \tilde{D}^{-1/2}$ comes from."

You are at a research-oriented ML company. The question is standard for graph ML roles. The interviewer wants to see whether you understand the derivation or just memorized the formula. These are very different things.

The formula is three steps: spectral graph convolution (requires eigendecomposition, $O(n^3)$) → Chebyshev polynomial approximation (no eigendecomposition, $O(K|E|)$) → first-order simplification with renormalization trick (one matrix multiply, stable across deep networks). At the end you get $H^{(l+1)} = \sigma(\tilde{D}^{-1/2}\tilde{A}\tilde{D}^{-1/2}H^{(l)}W^{(l)})$ - the GCN layer. Understanding this derivation tells you not just what GCN computes, but why those specific matrix operations and not others. It also explains why GCN has two layers in practice rather than twenty.

---

Why Graph Neural Networks?

Before GNNs, handling graph-structured data meant either flattening it (losing the relational structure) or using hand-crafted graph features (requiring domain expertise). Consider three problems:

Citation network classification: predict the research topic of a paper given its features (word frequencies) and its citation relationships to other papers. Standard MLP ignores the citations. A paper about deep learning that cites mostly deep learning papers is almost certainly about deep learning - the graph structure encodes crucial information.

Molecular property prediction: predict whether a drug molecule will bind to a protein target. The molecule is naturally a graph (atoms = nodes, bonds = edges). Flattening the molecule to a feature vector destroys the spatial and topological relationships that determine binding affinity.

Social network recommendation: predict whether two users will connect. The existing network of connections is the primary signal. Embedding each user independently ignores the rich context of their mutual connections.

In all these cases, the graph structure is not incidental - it is central to the prediction task. GNNs are the natural tool because they process graph-structured data directly, propagating information along edges.

---

Why This Exists: The Convolution Problem on Graphs

CNNs work because convolution on regular grids (images, audio) is efficient and captures spatial locality. The Convolution Theorem says: convolution in the spatial domain equals pointwise multiplication in the frequency domain. On a 2D grid, this means you can efficiently learn which spatial patterns (edges, textures, objects) predict the output.

Can we do the same for graphs? The challenge: graphs lack the regular structure of grids. There is no canonical "left neighbor" or "right neighbor" for a graph node. Nodes have variable-degree neighborhoods. The graph has no inherent spatial coordinate system.

The answer from signal processing: use the graph's own structure to define its "frequency domain." The eigenvectors of the normalized graph Laplacian form an orthogonal basis for functions on the graph - they are the graph's natural frequency basis, analogous to Fourier modes on a circle.

---

Step 1: Spectral Graph Convolution

The normalized graph Laplacian is:

$\mathcal{L} = I - D^{-1/2} A D^{-1/2} = U \Lambda U^\top$

where:

• $A \in \mathbb{R}^{n \times n}$ - adjacency matrix (1 if edge, 0 otherwise)
• $D$ - degree matrix: $D_{ii} = \sum_j A_{ij}$
• $U = [\mathbf{u}_1, \ldots, \mathbf{u}_n]$ - eigenvectors (graph Fourier basis)
• $\Lambda = \text{diag}(\lambda_1, \ldots, \lambda_n)$ - eigenvalues $\in [0, 2]$

Spectral graph convolution of signal $\mathbf{x}$ with filter $g_\theta$:

$g_\theta \star \mathbf{x} = U\, g_\theta(\Lambda)\, U^\top \mathbf{x}$

where $g_\theta(\Lambda) = \text{diag}(g_\theta(\lambda_1), \ldots, g_\theta(\lambda_n))$ is a learnable filter applied in the spectral domain. The filter $g_\theta$ has one parameter per eigenvalue - $n$ parameters total.

The intractable problem: this requires:

1. Computing eigenvectors $U$: $O(n^3)$ - for millions of nodes, completely infeasible
2. Storing $U$: $O(n^2)$ - also infeasible for large graphs
3. The filter has $n$ parameters - one per eigenvalue - does not generalize across graphs of different sizes

---

Step 2: ChebNet - Chebyshev Polynomial Approximation

Defferrard et al. (2016) published "Convolutional Neural Networks on Graphs with Fast Localized Spectral Filtering" (NeurIPS 2016). Their key insight: approximate the spectral filter $g_\theta(\Lambda)$ with a polynomial of the eigenvalues:

$g_\theta(\Lambda) \approx \sum_{k=0}^{K} \theta_k T_k(\tilde{\Lambda})$

where $\tilde{\Lambda} = \frac{2}{\lambda_{\max}} \Lambda - I$ rescales eigenvalues to $[-1, 1]$, and $T_k$ is the Chebyshev polynomial of degree $k$:

$T_0(x) = 1, \quad T_1(x) = x, \quad T_k(x) = 2x\, T_{k-1}(x) - T_{k-2}(x)$

Why Chebyshev polynomials? They form an orthogonal basis on $[-1, 1]$ and can be evaluated recursively. Crucially, computing $T_k(\tilde{\mathcal{L}})\mathbf{x}$ does not require eigendecomposition - it is computed recursively via sparse matrix-vector products:

$\mathbf{z}^{(0)} = \mathbf{x}, \quad \mathbf{z}^{(1)} = \tilde{\mathcal{L}}\mathbf{x}, \quad \mathbf{z}^{(k)} = 2\tilde{\mathcal{L}} \mathbf{z}^{(k-1)} - \mathbf{z}^{(k-2)}$

The full spectral convolution then becomes $O(K|E|)$ - linear in the number of edges times the filter order $K$. This is tractable for large sparse graphs.

$K$-locality property: a Chebyshev polynomial of order $K$ is $K$-localized. That is, the filter $g_\theta$ applied with order $K$ only depends on nodes within $K$ hops of the center node. This is the graph analog of a $K$-wide spatial filter in image convolutions - the "receptive field" is the $K$-hop neighborhood.

---

Step 3: Kipf & Welling (2017) - The GCN Simplification

Kipf & Welling (2017) introduced "Semi-Supervised Classification with Graph Convolutional Networks" (ICLR 2017), one of the most cited papers in graph ML. Their simplification: take $K = 1$ (first-order Chebyshev) and assume $\lambda_{\max} \approx 2$:

$g_\theta \star \mathbf{x} \approx \theta_0 \mathbf{x} + \theta_1 \left(\mathcal{L} - I\right)\mathbf{x} = \theta_0 \mathbf{x} - \theta_1 D^{-1/2} A D^{-1/2} \mathbf{x}$

With a single shared parameter $\theta = \theta_0 = -\theta_1$ (reducing from two parameters to one):

$g_\theta \star \mathbf{x} \approx \theta \left(I + D^{-1/2} A D^{-1/2}\right) \mathbf{x}$

This is a mean of the node's own features and its symmetrically-normalized neighborhood features.

---

The Renormalization Trick

$I + D^{-1/2} A D^{-1/2}$ has eigenvalues in $[0, 2]$ (both 0 and 2 can be achieved). When you stack $L$ layers of this operator, the operator raised to the $L$-th power has eigenvalues up to $2^L$. For $L = 10$, that is 1024× amplification - numerical explosion.

The fix: instead of $A$, use $\tilde{A} = A + I$ (add self-loops to every node). Recompute the degree matrix for $\tilde{A}$:

$\tilde{D}_{ii} = \sum_j \tilde{A}_{ij} = \deg(i) + 1$

Replace the unstable operator with the renormalized version:

$\hat{A} = \tilde{D}^{-1/2} \tilde{A} \tilde{D}^{-1/2}$

This operator has eigenvalues in $[0, 1]$. Raising it to any power keeps eigenvalues in $[0, 1]$ - numerically stable. The self-loop term $+I$ also has a semantic interpretation: every node's own representation is mixed into its neighborhood aggregate, preventing the node from "losing itself" in the aggregation.

Why symmetric normalization $\tilde{D}^{-1/2}\tilde{A}\tilde{D}^{-1/2}$ rather than row normalization $\tilde{D}^{-1}\tilde{A}$?

Row normalization gives: node $i$'s new representation = simple average of its neighbors' representations. This is a random walk normalization - the row sum is 1. Symmetric normalization gives: each edge's contribution is normalized by the geometric mean of endpoint degrees. Node $u$'s contribution to node $v$ scales as $1/\sqrt{\tilde{d}_u \tilde{d}_v}$. High-degree hub nodes spread their signal more thinly (large $\tilde{d}_u$) while nodes receiving from hubs also discount more (large $\tilde{d}_v$). In practice, symmetric normalization outperforms row normalization on most graph benchmarks.

---

The GCN Layer Formula

For a graph with $n$ nodes and node feature matrix $H^{(l)} \in \mathbb{R}^{n \times d_l}$ at layer $l$:

$\boxed{H^{(l+1)} = \sigma\!\left(\tilde{D}^{-1/2} \tilde{A}\, \tilde{D}^{-1/2} H^{(l)} W^{(l)}\right)}$

where:

• $\tilde{A} = A + I_n$ - adjacency matrix with added self-loops
• $\tilde{D}_{ii} = \sum_j \tilde{A}_{ij} = \deg(i) + 1$ - degree matrix of $\tilde{A}$
• $W^{(l)} \in \mathbb{R}^{d_l \times d_{l+1}}$ - learnable weight matrix
• $\sigma$ - nonlinear activation (ReLU for hidden, softmax for output)

Per-Node Interpretation

For node $v$, the update is:

$\mathbf{h}_v^{(l+1)} = \sigma\!\left(W^{(l)} \sum_{u \in \mathcal{N}(v) \cup \{v\}} \frac{1}{\sqrt{\tilde{d}_v \tilde{d}_u}} \mathbf{h}_u^{(l)}\right)$

where $\tilde{d}_v = \deg(v) + 1$. This is a degree-normalized aggregation of neighboring features plus the node's own features.

Node $v$ with 5 neighbors and degree 5 ($\tilde{d}_v = 6$) aggregates:

• Its own feature: weight $= 1/\sqrt{6 \times 6} = 1/6$
• Each neighbor $u$ with degree $\tilde{d}_u$: weight $= 1/\sqrt{6 \times \tilde{d}_u}$

A high-degree hub neighbor ($\tilde{d}_u = 100$): weight $= 1/\sqrt{6 \times 100} = 1/24.5$ - small contribution. A leaf neighbor ($\tilde{d}_u = 2$): weight $= 1/\sqrt{6 \times 2} = 1/3.5$ - larger contribution.

The normalization makes the hub node's influence proportional to how many other nodes it is also influencing - a natural fairness property.

---

Semi-Supervised Node Classification

Kipf & Welling's paper was specifically about semi-supervised learning on graphs. The Cora dataset split: 2708 papers, 5429 citations, 7 classes, but only 140 labeled nodes (about 5% of the total). GCN achieves ~81% test accuracy on 1000 test nodes despite seeing only 140 labels.

How? Three mechanisms work together:

1. Label propagation through structure: labeled nodes' class information propagates along citation edges. After 2 GCN layers, each node has aggregated information from all 2-hop neighbors. Papers that cite the same 3 papers as a labeled paper tend to belong to the same class.

2. Semi-supervised setting: all node features and all edges are visible during training, even for unlabeled nodes. The model learns a good feature space from the full graph, not just from labeled examples.

3. Homophily inductive bias: GCN's low-pass filtering enforces the assumption that connected nodes have similar representations - reasonable for citation networks where papers citing each other tend to have similar topics.

Transductive vs inductive: GCN is trained transductively. Test nodes are present during training (their features are used to compute other nodes' embeddings), but their labels are withheld. This is different from the traditional "train then test on unseen data" paradigm.

Training: use full graph G (all nodes, all edges, all features)
          compute forward pass for ALL nodes
          compute loss only on train_mask nodes

Inference: test nodes already have embeddings from training
           just read their logits (no new computation needed)

---

Over-Smoothing: Why Deeper ≠ Better

Each GCN layer applies a low-pass filter $\hat{A} = \tilde{D}^{-1/2}\tilde{A}\tilde{D}^{-1/2}$. This filter smooths node features by mixing them with neighbors'. After $k$ layers, the operation is $\hat{A}^k H^{(0)} W$. As $k \to \infty$:

$\hat{A}^k \to \mathbf{v}_1 \mathbf{v}_1^\top$

where $\mathbf{v}_1$ is the principal eigenvector of $\hat{A}$ (corresponding to eigenvalue 1). All node representations are projected onto this single vector - every node in the same connected component converges to the same representation. Distinction between nodes vanishes.

Measuring Over-Smoothing: MAD

Chen et al. (2020) quantified over-smoothing with the Mean Average Distance (MAD) metric:

$\text{MAD} = \frac{1}{|V|} \sum_{v \in V} \frac{\sum_{u \in V} \mathbf{1}[u \neq v] \cdot d(\mathbf{h}_v, \mathbf{h}_u)}{\sum_{u \in V} \mathbf{1}[u \neq v]}$

where $d(\mathbf{h}_v, \mathbf{h}_u) = 1 - \frac{\mathbf{h}_v \cdot \mathbf{h}_u}{\|\mathbf{h}_v\| \|\mathbf{h}_u\|}$ is the cosine distance. High MAD = diverse representations. Low MAD (near 0) = all nodes similar = over-smoothed.

Empirically:

• Depth 2: MAD ≈ 3.2, test accuracy ≈ 81%
• Depth 4: MAD ≈ 2.1, test accuracy ≈ 80%
• Depth 8: MAD ≈ 0.9, test accuracy ≈ 76%
• Depth 16: MAD ≈ 0.1, test accuracy ≈ 64%
• Depth 32: MAD ≈ 0.02, test accuracy ≈ 40%

The curve is clear: after 2–3 layers, deeper networks rapidly lose performance as node representations collapse.

---

Full Code: Manual GCN, PyG GCNConv, Over-Smoothing

import torch
import torch.nn as nn
import torch.nn.functional as F
from torch_geometric.datasets import Planetoid
from torch_geometric.utils import add_self_loops, degree
from torch_geometric.nn import GCNConv
import numpy as np

# ─── LOAD CORA ─────────────────────────────────────────────────────────────────

dataset = Planetoid(root='/tmp/cora', name='Cora')
data = dataset[0]

print(f"Nodes: {data.num_nodes}")         # 2708
print(f"Edges: {data.num_edges}")         # 10556 (undirected, stored both ways)
print(f"Features: {data.num_features}")   # 1433 (bag-of-words)
print(f"Classes: {dataset.num_classes}")  # 7
print(f"Train nodes: {data.train_mask.sum().item()}")   # 140
print(f"Val nodes:   {data.val_mask.sum().item()}")     # 500
print(f"Test nodes:  {data.test_mask.sum().item()}")    # 1000

# ─── MANUAL GCN LAYER ─────────────────────────────────────────────────────────

class GCNLayerManual(nn.Module):
    """
    Manual GCN layer implementing:
    H' = D̃^(-1/2) Ã D̃^(-1/2) H W

    Steps:
    1. Add self-loops: Ã = A + I
    2. Compute degree of Ã: D̃_{ii} = deg(i) + 1
    3. Compute normalization: D̃^(-1/2) per node
    4. For each edge (u,v): weight = 1/sqrt(d̃_u * d̃_v)
    5. Aggregate: new_h_v = W * sum_{u in N(v)∪{v}} weight(u,v) * h_u
    """
    def __init__(self, in_dim: int, out_dim: int):
        super().__init__()
        self.W = nn.Linear(in_dim, out_dim, bias=False)
        nn.init.xavier_uniform_(self.W.weight)

    def forward(self, x, edge_index, num_nodes):
        # Step 1: Add self-loops
        edge_index, _ = add_self_loops(edge_index, num_nodes=num_nodes)

        row, col = edge_index
        # Step 2: Degree of à (includes self-loop)
        deg = degree(row, num_nodes=num_nodes)  # D̃ diagonal

        # Step 3: D̃^(-1/2)
        deg_inv_sqrt = deg.pow(-0.5)
        deg_inv_sqrt[deg_inv_sqrt == float('inf')] = 0.0

        # Step 4: Edge weights = 1/sqrt(d̃_u * d̃_v)
        norm = deg_inv_sqrt[row] * deg_inv_sqrt[col]  # (num_edges,)

        # Step 5: Transform features then aggregate
        x = self.W(x)  # (n, out_dim)

        # Weighted message passing: for each edge (src→dst), add norm * x[src] to x[dst]
        out = torch.zeros_like(x)
        out.scatter_add_(
            0,
            col.unsqueeze(1).expand(-1, x.shape[1]),
            norm.unsqueeze(1) * x[row]
        )
        return out

# ─── TWO-LAYER GCN WITH PyG ───────────────────────────────────────────────────

class GCN(nn.Module):
    def __init__(self, num_features: int, num_classes: int, hidden: int = 64):
        super().__init__()
        self.conv1 = GCNConv(num_features, hidden)
        self.conv2 = GCNConv(hidden, num_classes)
        self.dropout = nn.Dropout(p=0.5)

    def forward(self, x, edge_index):
        x = self.conv1(x, edge_index)
        x = F.relu(x)
        x = self.dropout(x)
        x = self.conv2(x, edge_index)
        return x   # (num_nodes, num_classes)

    def embed(self, x, edge_index):
        """Return penultimate layer embeddings (before final projection)."""
        with torch.no_grad():
            x = F.relu(self.conv1(x, edge_index))
        return x

model = GCN(dataset.num_features, dataset.num_classes)
optimizer = torch.optim.Adam(model.parameters(), lr=0.01, weight_decay=5e-4)

# ─── TRAINING LOOP ────────────────────────────────────────────────────────────

def train(model, data, optimizer):
    model.train()
    optimizer.zero_grad()
    out = model(data.x, data.edge_index)
    loss = F.cross_entropy(out[data.train_mask], data.y[data.train_mask])
    loss.backward()
    optimizer.step()
    return loss.item()

@torch.no_grad()
def evaluate(model, data):
    model.eval()
    out = model(data.x, data.edge_index)
    pred = out.argmax(dim=1)
    results = {}
    for split in ['train', 'val', 'test']:
        mask = getattr(data, f'{split}_mask')
        correct = (pred[mask] == data.y[mask]).sum()
        results[split] = correct.item() / mask.sum().item()
    return results

print("\nTraining 2-layer GCN on Cora...")
best_val = 0
best_test = 0
for epoch in range(200):
    loss = train(model, data, optimizer)
    if epoch % 20 == 0:
        acc = evaluate(model, data)
        if acc['val'] > best_val:
            best_val = acc['val']
            best_test = acc['test']
        print(f"Epoch {epoch:03d} | Loss: {loss:.4f} | "
              f"Train: {acc['train']:.3f} | Val: {acc['val']:.3f}")

print(f"\nBest Val: {best_val:.3f} | Test: {best_test:.3f}")
# Expected: Test accuracy ~0.814

# ─── OVER-SMOOTHING MEASUREMENT ───────────────────────────────────────────────

def mean_average_distance(h: torch.Tensor) -> float:
    """
    MAD: average cosine distance between all pairs of node embeddings.
    Near 0 = all nodes have similar representations (over-smoothed).
    Near 1 = all nodes are maximally different.
    """
    h_norm = F.normalize(h, p=2, dim=1)  # (n, d), L2 normalized
    # Pairwise cosine similarity
    sim = h_norm @ h_norm.T              # (n, n)
    dist = 1 - sim                        # cosine distance
    # Average off-diagonal elements
    n = h.shape[0]
    mask = ~torch.eye(n, dtype=torch.bool)
    return dist[mask].mean().item()

# Build GCNs of increasing depth and measure MAD
print("\nOver-smoothing: MAD vs depth")
print(f"{'Depth':>6} {'Test Acc':>10} {'MAD':>10}")
print("-" * 30)

for depth in [1, 2, 4, 8, 16]:
    convs = nn.ModuleList()
    convs.append(GCNConv(dataset.num_features, 64))
    for _ in range(depth - 2):
        convs.append(GCNConv(64, 64))
    if depth > 1:
        convs.append(GCNConv(64, dataset.num_classes))
    else:
        convs.append(GCNConv(dataset.num_features, dataset.num_classes))

    deep_model = nn.Sequential()  # placeholder for deep GCN
    # (Training each deep GCN would take too long for demo - show expected results)
    # Expected results from literature:
    expected = {1: (0.71, 3.5), 2: (0.814, 3.2), 4: (0.80, 2.1),
                8: (0.76, 0.9), 16: (0.64, 0.1)}
    test_acc, mad = expected[depth]
    print(f"{depth:>6} {test_acc:>10.3f} {mad:>10.3f}")

# ─── MITIGATION 1: RESIDUAL CONNECTIONS ─────────────────────────────────────

class ResGCN(nn.Module):
    """GCN with residual connections to mitigate over-smoothing."""
    def __init__(self, num_features, num_classes, hidden=64, num_layers=4):
        super().__init__()
        self.input_proj = nn.Linear(num_features, hidden)
        self.convs = nn.ModuleList([
            GCNConv(hidden, hidden) for _ in range(num_layers)
        ])
        self.bns = nn.ModuleList([
            nn.BatchNorm1d(hidden) for _ in range(num_layers)
        ])
        self.output = nn.Linear(hidden, num_classes)

    def forward(self, x, edge_index):
        x = F.relu(self.input_proj(x))
        for conv, bn in zip(self.convs, self.bns):
            # Residual: new_x = σ(BN(conv(x))) + x
            x = F.relu(bn(conv(x, edge_index))) + x
        return self.output(x)

# ─── MITIGATION 2: DROPEDGE ──────────────────────────────────────────────────

class DropEdgeGCN(nn.Module):
    """GCN with DropEdge (Rong et al. 2020) - randomly remove edges during training."""
    def __init__(self, num_features, num_classes, hidden=64, drop_edge_rate=0.3):
        super().__init__()
        self.conv1 = GCNConv(num_features, hidden)
        self.conv2 = GCNConv(hidden, num_classes)
        self.drop_edge_rate = drop_edge_rate

    def forward(self, x, edge_index):
        from torch_geometric.utils import dropout_edge
        # During training: randomly drop 30% of edges
        # This acts as data augmentation and prevents over-smoothing
        if self.training:
            edge_index, _ = dropout_edge(
                edge_index, p=self.drop_edge_rate, training=True
            )
        x = F.relu(self.conv1(x, edge_index))
        x = F.dropout(x, p=0.5, training=self.training)
        return self.conv2(x, edge_index)

# ─── MITIGATION 3: JUMPING KNOWLEDGE NETWORKS ────────────────────────────────

class JKNet(nn.Module):
    """
    Jumping Knowledge Network (Xu et al. 2018).
    Aggregates representations from ALL layers, not just the last.
    Nodes with different optimal receptive fields use different depths.
    """
    def __init__(self, in_channels, hidden, out_channels, num_layers=4, mode='cat'):
        super().__init__()
        self.convs = nn.ModuleList()
        self.convs.append(GCNConv(in_channels, hidden))
        for _ in range(num_layers - 1):
            self.convs.append(GCNConv(hidden, hidden))

        self.mode = mode
        if mode == 'cat':
            # Concatenate all layer outputs
            self.out = nn.Linear(hidden * num_layers, out_channels)
        else:
            # Max pooling across layer outputs
            self.out = nn.Linear(hidden, out_channels)

    def forward(self, x, edge_index):
        representations = []
        for conv in self.convs:
            x = F.relu(conv(x, edge_index))
            representations.append(x)

        if self.mode == 'cat':
            # Each node concatenates its representation from all layers
            x = torch.cat(representations, dim=-1)
        else:
            # Max across layers (elementwise)
            x = torch.stack(representations, dim=0).max(dim=0)[0]

        return self.out(x)

# ─── VISUALIZING LEARNED EMBEDDINGS WITH T-SNE ────────────────────────────────

def visualize_embeddings(model, data, title="GCN Embeddings"):
    """
    Extract penultimate-layer embeddings and visualize with t-SNE.
    Correctly learned representations: nodes of same class cluster together.
    Over-smoothed representations: all nodes collapse to one cluster.
    """
    model.eval()
    with torch.no_grad():
        # Get embeddings from penultimate layer
        embeddings = model.embed(data.x, data.edge_index)
        embeddings_np = embeddings.numpy()
        labels = data.y.numpy()

    try:
        from sklearn.manifold import TSNE
        tsne = TSNE(n_components=2, random_state=42, perplexity=30)
        emb_2d = tsne.fit_transform(embeddings_np)

        print(f"\n{title} - t-SNE visualization")
        print("(In a Jupyter notebook, use matplotlib to show scatter plot)")
        print(f"Embedding shape: {embeddings_np.shape}")
        print(f"Number of classes: {len(set(labels.tolist()))}")

        # Measure cluster quality: intra-class vs inter-class distance
        class_embeddings = {}
        for cls in range(dataset.num_classes):
            mask = labels == cls
            class_embeddings[cls] = emb_2d[mask]

        intra_class_dists = []
        for cls, embs in class_embeddings.items():
            if len(embs) > 1:
                diffs = embs[:, None] - embs[None, :]
                dists = np.sqrt((diffs ** 2).sum(axis=-1))
                intra_class_dists.append(dists[np.triu_indices(len(embs), k=1)].mean())

        print(f"Mean intra-class distance: {np.mean(intra_class_dists):.2f}")
        print("(Smaller = tighter clusters = better class separation)")

    except ImportError:
        print("Install sklearn for t-SNE visualization")

visualize_embeddings(model, data, "2-Layer GCN on Cora")

---

Node Classification vs Link Prediction vs Graph Classification

GCN can be applied to three different prediction tasks with different output heads:

Node classification: predict a class for each node. Output head: linear layer on each node's final embedding. Loss: cross-entropy on labeled nodes. GCN propagates information from labeled to unlabeled nodes through graph structure.

Link prediction: predict whether an edge exists between two nodes. Output head: dot product or MLP on pairs of node embeddings. Loss: binary cross-entropy with negative sampling. GCN produces embeddings where adjacent nodes should be similar.

Graph classification: predict a label for the entire graph. Output head: global pooling (mean, max, or sum of node embeddings) followed by a linear layer. Loss: cross-entropy across graphs. Requires mini-batch training on multiple graphs.

---

Scalability: Full-Graph vs Mini-Batch

Standard GCN training loads the entire graph into GPU memory:

• Node features: $n \times d$ float tensor
• Adjacency: sparse $n \times n$ matrix

For Cora (2708 nodes, 1433 features): $2708 \times 1433 \times 4$ bytes ≈ 15MB - trivial. For Reddit (232,965 nodes, 602 features): $232,965 \times 602 \times 4$ bytes ≈ 560MB - manageable. For OGB-Proteins (132,534 nodes, 8 features, 39M edges): adjacency alone ≈ 600MB - tight. For Pinterest (3B nodes): entirely infeasible.

The solution for large graphs: GraphSAGE (next lesson) introduces mini-batch training with neighbor sampling, eliminating the full-graph memory requirement. For GCN specifically, Cluster-GCN (Chiang et al. 2019) partitions the graph into clusters and trains one cluster per mini-batch.

---

GCN Hyperparameter Guide

Hyperparameter	Recommended	Effect
Number of layers	2 (rarely 3)	More layers → over-smoothing
Hidden dimension	64–256	Larger = more capacity
Dropout rate	0.3–0.6	Prevents overfit on small labeled sets
Weight decay	5e-4	L2 regularization on $W$
Learning rate	0.01 (Adam)	Adam is robust here
Activation	ReLU	SELU and ELU also work
Self-loop	Always add	The renormalization trick requires it

---

YouTube Resources

Resource	Creator	Focus
Graph Convolutional Networks Explained	Aleksa Gordić	GCN derivation and PyG implementation
Semi-Supervised Classification with GCNs	Yannic Kilcher	Kipf & Welling 2017 paper review
PyTorch Geometric Tutorial	Antonio Longa	Hands-on PyG tutorial on Cora
Over-Smoothing in Graph Neural Networks	Standford CS224W	Over-smoothing analysis and fixes
Spectral Graph Theory for GNNs	Standford CS224W	Laplacian eigenvalues and Chebyshev approx

---

Common Mistakes

:::danger Common Mistake 1: Not adding self-loops before applying GCNConv Without self-loops, a node's representation in layer $l+1$ is computed only from its neighbors, not from itself. In a 2-layer GCN, a node's final embedding does not include any direct contribution from its own initial features - the signal gets mixed away. Always add self-loops: $\tilde{A} = A + I$. PyG's GCNConv adds self-loops by default, but if you are implementing manually, this is the most common oversight. :::

:::danger Common Mistake 2: Using more than 2-3 GCN layers without mitigations A 10-layer GCN on Cora will almost certainly perform worse than a 2-layer GCN due to over-smoothing. If you need more depth for a large receptive field, use residual connections (ResGCN), DropEdge, or Jumping Knowledge Networks. The node classification error curve peaks around depth 2 and degrades rapidly after depth 4. :::

:::warning Common Mistake 3: Confusing transductive and inductive evaluation In the standard Cora split, test nodes are present in the graph during training (their features are used for message passing even though their labels are withheld). This is transductive evaluation. If you split the graph and hide test nodes entirely, performance will be lower because the model cannot use their structural context. When reporting GCN results, specify whether the evaluation is transductive (standard) or inductive (more realistic for production). :::

:::warning Common Mistake 4: Using GCN for graphs with heterophily GCN's low-pass filtering assumes that connected nodes have similar representations (homophily). On heterophilic graphs (where connected nodes have different labels - e.g., some social networks, some bioinformatics graphs), this assumption is wrong and GCN performs poorly. For heterophilic settings, GAT or H2GCN are better choices. :::

---

Interview Q&A

Q1: Derive the GCN layer from spectral graph convolution. What are the three approximations?

Spectral graph convolution is $g_\theta \star x = U g_\theta(\Lambda) U^\top x$, requiring eigendecomposition of the Laplacian: $O(n^3)$ and $O(n^2)$ memory. Three approximations reduce this to the GCN layer. Approximation 1 (Defferrard 2016, ChebNet): replace $g_\theta(\Lambda)$ with a $K$-th order Chebyshev polynomial $\sum_{k=0}^K \theta_k T_k(\tilde{\Lambda})$. This avoids eigendecomposition - Chebyshev polynomials are evaluated recursively via sparse matrix products. Cost drops from $O(n^3)$ to $O(K|E|)$. Approximation 2 (Kipf 2017): take $K=1$ and assume $\lambda_{\max} \approx 2$. Simplifies to $\theta(I + D^{-1/2}AD^{-1/2})x$ with one parameter. Approximation 3 - Renormalization: $I + D^{-1/2}AD^{-1/2}$ has eigenvalues in $[0,2]$; stacking layers causes explosion. Add self-loops: $\tilde{A} = A+I$, $\tilde{D}_{ii} = \deg(i)+1$, use $\tilde{D}^{-1/2}\tilde{A}\tilde{D}^{-1/2}$ which has eigenvalues in $[0,1]$. Result: $H^{(l+1)} = \sigma(\tilde{D}^{-1/2}\tilde{A}\tilde{D}^{-1/2}H^{(l)}W^{(l)})$.

Q2: Explain symmetric normalization $\tilde{D}^{-1/2}\tilde{A}\tilde{D}^{-1/2}$. Why not just row normalization $\tilde{D}^{-1}\tilde{A}$?

Without normalization ($AH$): high-degree nodes accumulate much larger feature vectors than low-degree nodes, causing instability - a node with 1000 neighbors gets 1000× the signal of a leaf node. Row normalization ($\tilde{D}^{-1}\tilde{A}H$): each node computes a simple average of its neighbors. This is stable but asymmetric - a hub contributing to 1000 nodes is treated the same as a leaf. Symmetric normalization ($\tilde{D}^{-1/2}\tilde{A}\tilde{D}^{-1/2}H$): each edge $(u,v)$ has weight $1/\sqrt{\tilde{d}_u \tilde{d}_v}$. Hub nodes have large $\tilde{d}_u$, so their contribution to each neighbor is proportionally smaller - they spread their signal more thinly. Nodes receiving from hubs also discount more (large $\tilde{d}_v$). This is more balanced: hub nodes do not dominate just by having many connections. In practice, symmetric normalization consistently outperforms row normalization on node classification benchmarks.

Q3: What is over-smoothing? How do you detect it and what are the mitigations?

Over-smoothing is the convergence of node representations to the same value as GCN depth increases. Each layer applies a low-pass filter $\hat{A}^k$; as $k \to \infty$, $\hat{A}^k$ projects all nodes onto the same principal eigenvector direction - all representations become identical within connected components. Detection: use the MAD metric - average cosine distance between node pairs. MAD near 0 means over-smoothed. Mitigations: (1) Residual connections - add $x$ directly to $\text{conv}(x)$ at each layer; the residual preserves the node's original representation. (2) DropEdge - randomly remove a fraction of edges during training; reduces the mixing effect and acts as data augmentation. (3) Jumping Knowledge Networks - concatenate or take max of representations from all intermediate layers; nodes can use the depth with the right receptive field for their position. (4) PairNorm - normalize representations to have fixed variance, preventing collapse. In practice: use 2 layers as default, add residuals if depth is needed.

Q4: GCN achieves ~81% test accuracy on Cora with only 140 training labels. How?

Three mechanisms: First, message passing propagates label information through citation edges. A paper citing 5 highly-labeled "deep learning" papers strongly suggests it is also a deep learning paper. After 2 GCN layers, each node aggregates information from all nodes within 2 hops - this captures broader neighborhood context. Second, semi-supervised learning: all node features and edges are visible during training (even unlabeled test nodes participate in message passing). The model learns a global feature space from the entire graph structure, not just from 140 examples. Third, homophily inductive bias: GCN's low-pass filter encourages similar representations for connected nodes. In Cora, citation connection strongly correlates with topic similarity (papers in the same area cite each other), so this bias is well-calibrated. The 140 labeled nodes effectively "anchor" 7 regions of the graph to their respective classes, and the GCN propagates these anchors outward through the citation structure.

Q5: What is the difference between node classification, link prediction, and graph classification in GCN?

Node classification: predict a label for each node. GCN outputs a $(n, C)$ matrix of logits, one row per node. Loss is cross-entropy on labeled rows only. Standard for Cora, CiteSeer. Link prediction: predict whether an edge exists between a pair of nodes. GCN produces node embeddings; link probability is computed from pairs - typically dot product $\sigma(z_u^T z_v)$ or MLP on $[z_u; z_v]$. Loss is binary cross-entropy with negative sampling (random pairs without edges as negatives). Standard for recommendation systems, knowledge graph completion. Graph classification: predict a label for an entire graph (not individual nodes). Requires a readout function to aggregate node embeddings into a graph-level representation - sum, mean, or max pooling across all nodes. Then a linear head predicts the graph label. Used in molecular property prediction (each molecule is a graph), program analysis. Requires batching multiple graphs per training step, which PyG handles with a batch pointer tensor.

---

Historical Context: From Spectral to Spatial

The history of GCNs is a story of simplification. Bruna et al. (2014) introduced the first spectral graph convolution, but it required storing and computing full eigenvectors - completely intractable for large graphs. Henaff et al. (2015) proposed learning smooth spectral multipliers but retained the eigendecomposition requirement.

Defferrard et al. (2016, ChebNet) made the breakthrough: approximate spectral filters with Chebyshev polynomials, eliminating eigendecomposition. ChebNet could handle graphs with millions of nodes. The parameter count dropped from $O(n)$ (one per eigenvalue) to $O(K)$ (one per polynomial degree).

Kipf & Welling (2017) simplified ChebNet further to $K=1$, trading expressiveness for computational efficiency and eliminating the need to tune polynomial order. The resulting model is simple enough to reason about analytically and fast enough to train on CPUs. Its paper accumulated over 10,000 citations within five years.

The spatial interpretation of GCN emerged post-hoc: the $\tilde{D}^{-1/2}\tilde{A}\tilde{D}^{-1/2}$ operation is equivalent to mean aggregation with degree normalization - a message passing operation that happens to emerge from spectral theory. This connection between spectral and spatial graph convolution unified two previously separate research threads.

The Graph Isomorphism Network (GIN, Xu et al., 2019) later showed that sum aggregation (not mean or max) is necessary for a GNN to be as expressive as the Weisfeiler-Lehman graph isomorphism test. GCN's mean aggregation is strictly less expressive than GIN - important context for tasks requiring structural discrimination.

---

Key Takeaways

GCN is the canonical entry point for graph neural networks. Its derivation - from spectral graph convolution through Chebyshev approximation to first-order simplification with renormalization - tells you why each design choice exists. The formula $H^{(l+1)} = \sigma(\tilde{D}^{-1/2}\tilde{A}\tilde{D}^{-1/2}H^{(l)}W^{(l)})$ is not arbitrary; it is the most efficient stable approximation of a localized spectral filter.

In practice: use 2 layers, add dropout (0.5), add weight decay (5e-4), and expect ~81% on Cora. For deeper models, add residual connections or DropEdge. For large graphs, switch to GraphSAGE with mini-batch training. For heterophilic graphs or tasks requiring attention, switch to GAT or GATv2. GCN is the baseline - understand it deeply, then extend it appropriately for your specific constraints.

---

Choosing Between GCN, GAT, and GraphSAGE

The three major GNN architectures serve different problem profiles. This decision framework helps select the right tool.

Property	GCN	GAT	GraphSAGE
Neighbor weighting	Fixed (degree-normalized)	Learned (attention)	Fixed (mean/max/LSTM)
Scalability	Full-graph only	Full-graph (dense)	Mini-batch, large graphs
Inductive capability	No (transductive)	Yes	Yes
Heterophily	Poor	Good	Moderate
Interpretability	Low	High (attention weights)	Low
Default aggregation	Normalized mean	Attention-weighted sum	Concatenate + transform
Benchmark (Cora)	~81%	~83%	~81%
Benchmark (PPI)	~59%	~97%	~97%

The practical workflow: always start with GCN as baseline. If performance is insufficient and the graph is heterophilic, move to GATv2. If the graph is large or you need to embed new nodes at inference time, move to GraphSAGE. If you need the best of both - attention-weighted and scalable - use GraphSAGE with attention aggregators (Graph Attention Sampling Aggregator, GASA).

---

GCN for Molecular Property Prediction

Node classification on citation networks is the canonical demo, but GCN's most commercially valuable application is molecular property prediction. Molecules are naturally represented as graphs: atoms are nodes (features: atomic number, charge, hybridization), bonds are edges (features: bond type, aromaticity).

from torch_geometric.datasets import MoleculeNet
from torch_geometric.data import DataLoader
from torch_geometric.nn import GCNConv, global_mean_pool

dataset = MoleculeNet(root='data/ESOL', name='ESOL')  # solubility prediction
print(f"Molecules: {len(dataset)}, Features: {dataset.num_features}, "
      f"Targets: {dataset.num_tasks}")

class MolGCN(torch.nn.Module):
    def __init__(self, in_channels, hidden, out_channels=1):
        super().__init__()
        self.conv1 = GCNConv(in_channels, hidden)
        self.conv2 = GCNConv(hidden, hidden)
        self.conv3 = GCNConv(hidden, hidden)
        self.head = torch.nn.Linear(hidden, out_channels)

    def forward(self, x, edge_index, batch):
        # 3 GCN layers with residual
        h = F.relu(self.conv1(x, edge_index))
        h = F.relu(self.conv2(h, edge_index)) + h  # residual
        h = F.relu(self.conv3(h, edge_index)) + h  # residual
        # Graph-level readout: mean over all atoms in each molecule
        h = global_mean_pool(h, batch)
        return self.head(h)  # regression: predict log solubility

model = MolGCN(in_channels=dataset.num_features, hidden=128)
loader = DataLoader(dataset, batch_size=32, shuffle=True)

# Training loop for molecular regression
optimizer = torch.optim.Adam(model.parameters(), lr=1e-3)
for epoch in range(100):
    model.train()
    for batch in loader:
        out = model(batch.x.float(), batch.edge_index, batch.batch)
        loss = F.mse_loss(out.squeeze(), batch.y.squeeze())
        optimizer.zero_grad()
        loss.backward()
        optimizer.step()

The key difference from node classification: the global_mean_pool aggregates over all atoms in a molecule to produce a fixed-size graph-level embedding. This transforms GCN from a node-level to a graph-level predictor. RMSE on ESOL solubility typically reaches ~0.6 log mol/L with this architecture - competitive with state-of-the-art deep learning methods from 2020.

---

Scalability: Full-Batch vs Mini-Batch for GCN

Standard GCN requires the full adjacency matrix and all node features in GPU memory. For graphs beyond ~100K nodes this becomes infeasible. The transition options:

Approach	Memory	Accuracy	When to Use
Full-batch GCN	$O(n \cdot d)$	Best	Graphs with <100K nodes
Cluster-GCN	$O(\text{cluster})$	≈Full-batch	Medium graphs, partition-friendly
GraphSAGE NeighborLoader	$O(S^K \cdot B)$	-1 to 2%	Graphs with millions of nodes
Graph-SAINT	$O(\text{subgraph})$	≈Full-batch	Graphs where random walks are cheap

Cluster-GCN (Chiang et al., 2019) partitions the graph into subgraphs using METIS, then trains on one subgraph per batch. Within-cluster edges are preserved (no gradient noise from cross-cluster edge removal), but between-cluster information is lost for that batch. Works well when the graph has natural cluster structure (social networks, citation networks).

For production systems: if your graph fits in GPU memory, use standard PyG full-batch GCN. If it does not, the choice between Cluster-GCN and NeighborLoader depends on whether your graph has natural partitioning structure (favor Cluster-GCN) or highly irregular degree distributions (favor NeighborLoader).

---

GCN Expressiveness and the WL Test

A critical question for any GNN: when do two structurally different graphs produce the same embedding? This is the graph isomorphism problem, and the answer determines how much structural information a GNN can capture.

The Weisfeiler-Lehman (WL) graph isomorphism test (1968) is a classical algorithm for determining whether two graphs are isomorphic. It iteratively assigns color labels to nodes based on their neighborhood colors, much like GNN message passing. WL can distinguish most pairs of non-isomorphic graphs but fails on some symmetric structures (e.g., two regular graphs with the same degree).

Xu et al. (2019) proved that GNNs using summation aggregation are at most as powerful as the WL test. Mean and max aggregation - used by GCN and GraphSAGE respectively - are strictly less expressive. Two graphs that WL distinguishes will always be distinguished by GIN (sum aggregation); some pairs that WL distinguishes will be mapped to the same embedding by GCN (mean aggregation).

Example where GCN fails to distinguish:
Graph A: node v has 3 neighbors each with feature [1]
Graph B: node v has 1 neighbor with feature [3]

GCN mean aggregation: both produce AGG = [1] - identical
GIN sum aggregation:  A produces AGG = [3], B produces AGG = [3] - also identical here

In practice, the WL expressiveness limitation matters primarily for graph classification tasks where the global structure must be discriminated (e.g., distinguishing cubic graphs, detecting cycles). For node classification with feature-rich nodes (Cora's 1433-dimensional bag-of-words), the feature content dominates and mean aggregation is usually sufficient. The distinction matters most in molecular property prediction, where graph topology is the information (no node features beyond atom type).

---

Advanced Interview Q&A

Q6: How does GCN handle directed graphs, and what changes in the normalization?

Standard GCN assumes undirected graphs. The renormalization trick produces a symmetric normalized adjacency $\tilde{D}^{-1/2}\tilde{A}\tilde{D}^{-1/2}$, which requires $\tilde{A}$ to be symmetric (i.e., if edge $(i,j)$ exists, edge $(j,i)$ also exists). For directed graphs, two options: (1) Treat as undirected: add both directions for every edge, losing directional information. Simplest approach, works surprisingly well for citation networks (citations are directional but mutual citation is common). (2) Use asymmetric normalization: $\tilde{D}^{-1}\tilde{A}$ (row normalization) or $\tilde{A}\tilde{D}^{-1}$ (column normalization). Row normalization: each node averages its out-neighbors, so information flows in the direction of edges. Column normalization: each node's information is weighted by its in-degree, so popular nodes contribute less per edge. For knowledge graphs and heterogeneous directed graphs, specialized models like RGCN (with separate weights per edge direction) are preferred. PyG's GCNConv has an add_self_loops and normalize parameter; for directed graphs, pass normalize=False and compute your own asymmetric normalization.

Q7: What is the relationship between GCN and spectral clustering?

Both GCN and spectral clustering use the graph Laplacian $L = D - A$. Spectral clustering computes the $k$ smallest eigenvectors of $L$ and uses them as node features for k-means clustering. GCN approximates a function of the graph Laplacian's eigenvalues (the spectral filter), parameterized by learnable weights $\theta$. The connection is tight: the GCN propagation rule $\tilde{D}^{-1/2}\tilde{A}\tilde{D}^{-1/2}$ is equivalent to one step of a normalized random walk, which is closely related to the normalized Laplacian $L_{sym} = I - D^{-1/2}AD^{-1/2}$. In the limit of infinite layers and zero weight matrix, GCN converges to the leading eigenvectors of $\tilde{D}^{-1/2}\tilde{A}\tilde{D}^{-1/2}$ - exactly what spectral clustering computes. GCN thus generalizes spectral clustering by: (1) making the spectral filter learnable, (2) incorporating node features, and (3) scaling linearly in edges rather than cubically in nodes.

---

GCN as a Low-Pass Graph Filter

The spectral interpretation of GCN reveals something non-obvious: GCN is fundamentally a low-pass filter on graph signals.

In signal processing, a low-pass filter suppresses high-frequency components (rapid oscillations) and preserves low-frequency components (smooth, slowly varying signals). On a graph, "frequency" is defined in terms of the graph Laplacian eigenvalues: low-frequency eigenvectors are smooth (connected nodes have similar values), high-frequency eigenvectors oscillate rapidly across edges (connected nodes have opposite values).

GCN's propagation rule $\tilde{D}^{-1/2}\tilde{A}\tilde{D}^{-1/2}$ applies a spectral filter with eigenvalues $\lambda_k \in [0, 2]$. The GCN filter maps eigenvalue $\lambda_k$ to weight $1 - \lambda_k/2$ - eigenvalues near 0 (smooth) are amplified, eigenvalues near 2 (oscillatory) are damped. This is exactly a low-pass filter.

The practical implication: GCN is well-suited for homophilic graphs (where connected nodes tend to have the same label), because low-pass filtering encourages similar representations for connected nodes. For heterophilic graphs (connected nodes tend to have different labels), the low-pass bias works against you - it smooths away exactly the information you need.

import torch
import numpy as np
from torch_geometric.utils import to_dense_adj, get_laplacian

def analyze_gcn_as_filter(data):
    """
    Analyze the GCN propagation matrix as a spectral filter.
    Shows which frequencies are amplified vs attenuated.
    """
    # Compute normalized adjacency (the GCN propagation matrix)
    n = data.num_nodes
    edge_index = data.edge_index

    # Add self-loops
    from torch_geometric.utils import add_self_loops
    edge_index_sl, _ = add_self_loops(edge_index, num_nodes=n)

    # Compute degree matrix
    deg = torch.zeros(n)
    deg.scatter_add_(0, edge_index_sl[0], torch.ones(edge_index_sl.shape[1]))

    # Normalized adjacency: D^{-1/2} A D^{-1/2}
    norm = deg.pow(-0.5)
    norm[norm == float('inf')] = 0
    A_dense = to_dense_adj(edge_index_sl, max_num_nodes=n).squeeze(0)
    norm_mat = torch.diag(norm)
    prop_matrix = norm_mat @ A_dense @ norm_mat

    # Eigenvalue decomposition
    eigenvalues, _ = torch.linalg.eigh(prop_matrix)
    # For normalized adjacency with self-loops, eigenvalues in ~[0, 1]
    # GCN filter response: eigenvalue directly (passes low freq, attenuates high)

    print("GCN spectral filter response:")
    print(f"  Min eigenvalue: {eigenvalues.min():.4f} (most attenuated frequency)")
    print(f"  Max eigenvalue: {eigenvalues.max():.4f} (least attenuated frequency)")
    print(f"  Low-freq components (>0.9): {(eigenvalues > 0.9).sum().item()} eigenvalues")
    print(f"  High-freq components (<0.1): {(eigenvalues < 0.1).sum().item()} eigenvalues")

    return eigenvalues.numpy()

This low-pass interpretation also explains over-smoothing: stacking $L$ GCN layers applies the filter $L$ times. Each application shrinks the high-frequency components further toward zero, until after enough layers, every node's embedding is proportional to the dominant low-frequency eigenvector - the same vector for all nodes. Over-smoothing is not a bug in GCN's implementation; it is a fundamental consequence of repeated low-pass filtering.

This insight motivates high-pass and band-pass graph filters (GPRGNN, ChebNetII) for heterophilic graphs, and multi-scale methods (JK-Networks, APPNP) that mix different filtering depths to avoid the smoothing collapse.

---

GCN Debugging Checklist

When GCN does not converge or produces poor accuracy, work through this checklist systematically.

Symptom	Likely Cause	Fix
Loss doesn't decrease	Learning rate too high or too low	Try lr=0.01 (Adam), 0.1 (SGD)
Train acc high, val acc low	Overfitting	Add dropout (0.5), weight decay (5e-4)
Val acc plateaus at 50–60%	Over-smoothing from too many layers	Reduce to 2 layers; add residuals
NaN loss	Numerical instability in normalization	Check for isolated nodes (degree=0); add self-loops
Poor performance on heterophilic data	Low-pass filter bias	Switch to GAT or H2GCN
OOM on large graph	Full-batch training	Switch to NeighborLoader (GraphSAGE)
Slow convergence	No feature normalization	Add NormalizeFeatures() transform

# Quick diagnostic: check graph properties before choosing GCN
from torch_geometric.utils import homophily

h = homophily(data.edge_index, data.y, method='node')
print(f"Node homophily: {h:.3f}")
# h > 0.7: GCN will work well (high homophily)
# 0.4 < h < 0.7: try GCN, but GAT may outperform
# h < 0.4: low homophily - use GAT, H2GCN, or GPRGNN

# Check for isolated nodes (can cause NaN gradients)
from torch_geometric.utils import degree
deg = degree(data.edge_index[0], num_nodes=data.num_nodes)
isolated = (deg == 0).sum().item()
print(f"Isolated nodes: {isolated}")
if isolated > 0:
    print("Warning: add self-loops or handle isolated nodes before GCN training")

:::tip 🎮 Interactive Playground

Visualize this concept: Try the GNN Message Passing demo on the EngineersOfAI Playground - no code required.

:::