Chapter 2: Graph Convolutional Networks (GCN)

1. The Mechanics of Message Passing

Message Passing is a processing framework in which each node in a graph receives information (messages) from its neighboring nodes, aggregates them, and uses the result to update its own features. Most graph neural networks (GCN, GraphSAGE, GAT, etc.) are built on top of this message passing framework.

Message passing consists of the following three steps.

• Message (Message Generation): For each edge, compute a "message" from the source node's features (and, if applicable, the edge features)
• Aggregate (Aggregation): Combine the multiple messages received by each node into a single vector using an operation such as sum, mean, or max
• Update: Compute new features using the aggregated messages together with the node's own current features

Expressing this framework mathematically, the update for node \(i\) at layer \(k\) is as follows.

$$\mathbf{m}_{j \to i}^{(k)} = \text{MESSAGE}^{(k)}\left(\mathbf{x}_i^{(k-1)}, \mathbf{x}_j^{(k-1)}\right)$$

$$\mathbf{a}_i^{(k)} = \text{AGGREGATE}^{(k)}\left(\left\{\mathbf{m}_{j \to i}^{(k)} \mid j \in \mathcal{N}(i)\right\}\right)$$

$$\mathbf{x}_i^{(k)} = \text{UPDATE}^{(k)}\left(\mathbf{x}_i^{(k-1)}, \mathbf{a}_i^{(k)}\right)$$

Here \(\mathcal{N}(i)\) is the set of neighboring nodes of node \(i\). By repeating these three steps (stacking layers), each node can indirectly incorporate information from increasingly distant nodes.

Experiencing Message Passing by Hand

Let's first implement the three steps of message passing in plain PyTorch, without using any torch_geometric layers.

import torch

# Graph with 4 nodes (directed edges: 0→1, 1→2, 0→2, 2→3, 1→3)
edge_index = torch.tensor([
    [0, 1, 0, 2, 1],  # source
    [1, 2, 2, 3, 3]   # destination
], dtype=torch.long)

# Initial features for each node (2-dimensional)
x = torch.tensor([
    [1.0, 0.0],
    [0.0, 1.0],
    [1.0, 1.0],
    [0.5, 0.5]
], dtype=torch.float)

num_nodes = x.size(0)
src, dst = edge_index[0], edge_index[1]

# Step 1: Message (here, the source node's features are used as the message as-is)
messages = x[src]  # Shape: [num_edges, num_features]

# Step 2: Aggregate (sum the messages received by each destination node)
aggregated = torch.zeros(num_nodes, x.size(1))
aggregated.index_add_(0, dst, messages)

# Step 3: Update (for simplicity, use the aggregated result directly as the new features)
x_new = aggregated

print("Features before aggregation:\n", x)
print("\nFeatures after aggregation (having received neighbor information):\n", x_new)

This code uses index_add_ to sum the messages for each edge's destination node. In an actual GCN, Step 1 (message generation) incorporates a linear transformation and normalization coefficients, while Step 3 (update) incorporates an activation function.

Implementation Using the MessagePassing Base Class

PyTorch Geometric provides this 3-step framework as a base class called torch_geometric.nn.MessagePassing. By simply defining a message() method and (if needed) an update() method, the propagate() method automatically handles the aggregation process for you.

import torch
from torch_geometric.nn import MessagePassing
from torch_geometric.utils import add_self_loops, degree

class SimpleGCNLayer(MessagePassing):
    def __init__(self, in_channels, out_channels):
        super().__init__(aggr='add')  # Aggregation method: sum
        self.lin = torch.nn.Linear(in_channels, out_channels, bias=False)

    def forward(self, x, edge_index):
        # Add self-loops (include each node's own information in the aggregation)
        edge_index, _ = add_self_loops(edge_index, num_nodes=x.size(0))

        # Linear transformation
        x = self.lin(x)

        # Compute normalization coefficients (inverse square root of degree)
        row, col = edge_index
        deg = degree(col, x.size(0), dtype=x.dtype)
        deg_inv_sqrt = deg.pow(-0.5)
        deg_inv_sqrt[deg_inv_sqrt == float('inf')] = 0
        norm = deg_inv_sqrt[row] * deg_inv_sqrt[col]

        # propagate() automatically calls message → aggregate → update in order
        return self.propagate(edge_index, x=x, norm=norm)

    def message(self, x_j, norm):
        # x_j: features of the source (neighbor) node. PyG collects these automatically
        return norm.view(-1, 1) * x_j

# Verify behavior
edge_index = torch.tensor([[0, 1, 1, 2],
                           [1, 0, 2, 1]], dtype=torch.long)
x = torch.tensor([[1.0, 2.0], [3.0, 4.0], [5.0, 6.0]], dtype=torch.float)

layer = SimpleGCNLayer(in_channels=2, out_channels=4)
out = layer(x, edge_index)
print("Output shape:", out.shape)
print("Output:\n", out)

This code closely reproduces what the GCNConv layer (introduced later) does internally. The actual GCNConv layer adds further caching and numerical stability enhancements, but the basic principle is the same.

2. Mathematical Background of Graph Convolution

GCN is a method proposed by Kipf & Welling (2017) that formulates a computationally efficient, practical layer by taking a first-order approximation of spectral graph convolution from graph signal processing. Here, we focus on the concept of "symmetric normalization," which is essential for implementation.

Self-Loops and the Normalized Adjacency Matrix

Let \(N\) be the number of nodes and \(A \in \mathbb{R}^{N \times N}\) be the graph's adjacency matrix. Since GCN includes each node's own features in the aggregation, it uses an adjacency matrix with self-loops (Adjacency Matrix with Self-loops), obtained by adding the identity matrix \(I_N\).

$$\tilde{A} = A + I_N$$

Next, define the degree matrix \(\tilde{D}\) as a diagonal matrix.

$$\tilde{D}_{ii} = \sum_{j} \tilde{A}_{ij}$$

Then, to prevent nodes with large degree from having an excessive influence, symmetric normalization (Symmetric Normalization) is applied.

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

Using this \(\hat{A}\), the layer-wise propagation rule for GCN can be written as follows.

$$H^{(l+1)} = \sigma\left(\hat{A} H^{(l)} W^{(l)}\right)$$

Where:

• \(H^{(l)} \in \mathbb{R}^{N \times d_l}\): The feature matrix of all nodes at layer \(l\)
• \(W^{(l)}\): The learnable weight matrix at layer \(l\)
• \(\sigma\): Activation function (e.g., ReLU)

The \(\hat{A} H^{(l)}\) part of this equation corresponds precisely to the "Aggregate" step in the message passing framework from the previous section. With a simple sum that doesn't use symmetric normalization, nodes with large degree would produce excessively large values, destabilizing training; normalization by \(\tilde{D}^{-1/2} \tilde{A} \tilde{D}^{-1/2}\) mitigates this.

import numpy as np

# 4-node undirected graph (no self-loops)
A = np.array([
    [0, 1, 1, 0],
    [1, 0, 1, 1],
    [1, 1, 0, 1],
    [0, 1, 1, 0]
], dtype=float)

# Add self-loops: A_tilde = A + I
I = np.eye(4)
A_tilde = A + I

# Degree matrix D_tilde
D_tilde = np.diag(A_tilde.sum(axis=1))

# Symmetric normalization: D_tilde^(-1/2) A_tilde D_tilde^(-1/2)
D_inv_sqrt = np.diag(1.0 / np.sqrt(np.diag(D_tilde)))
A_norm = D_inv_sqrt @ A_tilde @ D_inv_sqrt

print("Adjacency matrix with self-loops A_tilde:\n", A_tilde)
print("\nDegree matrix D_tilde:\n", D_tilde)
print("\nNormalized adjacency matrix A_norm:\n", np.round(A_norm, 3))

Each element of the output A_norm has the same value as the normalization coefficients used internally by GCNConv. The GCNConv layer we'll look at in the next section performs this sparse normalization computation efficiently (without explicitly constructing the adjacency matrix, working directly with edge indices).

3. Details of the GCNConv Layer

PyTorch Geometric provides the torch_geometric.nn.GCNConv layer, which directly implements the formulas from the previous section. Since it internally handles adding self-loops, symmetric normalization, and the linear transformation, we can build a GCN simply by stacking these layers.

Main Parameters

import torch
from torch_geometric.nn import GCNConv

conv = GCNConv(in_channels=16, out_channels=8)

print("Weight matrix shape:", conv.lin.weight.shape)
print("Bias shape:", conv.bias.shape)
print("Number of learnable parameters:", sum(p.numel() for p in conv.parameters()))

# Verify behavior on a simple graph
x = torch.randn(5, 16)  # 5 nodes, 16-dimensional features
edge_index = torch.tensor([[0, 1, 2, 3, 4, 0],
                           [1, 2, 3, 4, 0, 2]], dtype=torch.long)

out = conv(x, edge_index)
print("\nOutput shape:", out.shape)  # [5, 8]

Checking the Effect of the improved Parameter

Setting improved=True uses \(\tilde{A} = A + 2I\), which doubles the weight of self-loops. This is an effective setting when you want to more strongly retain a node's own information.

import torch
from torch_geometric.nn import GCNConv

x = torch.tensor([[1.0], [2.0], [3.0]], dtype=torch.float)
edge_index = torch.tensor([[0, 1, 1, 2],
                           [1, 0, 2, 1]], dtype=torch.long)

torch.manual_seed(42)
conv_normal = GCNConv(1, 1, improved=False, add_self_loops=True)
torch.manual_seed(42)
conv_improved = GCNConv(1, 1, improved=True, add_self_loops=True)

out_normal = conv_normal(x, edge_index)
out_improved = conv_improved(x, edge_index)

print("Output of standard GCNConv (improved=False):\n", out_normal)
print("\nOutput with improved=True (self-loop weight doubled):\n", out_improved)

4. Implementing the Node Classification Task (Cora Dataset)

From here on, we'll use the Cora dataset, introduced in Chapter 1, to implement node classification with GCN. Cora is a citation network with 2,708 papers as nodes and citation relationships as edges; each paper has a 1,433-dimensional Bag-of-Words feature vector and is classified into one of 7 research fields.

Defining a 3-Layer GCN Model

In Chapter 1, we implemented a simple 2-layer GCN, but here we'll implement a 3-layer GCN with two hidden layers, inserting Dropout (explained later) between each layer. Note that the final layer does not pass through an activation function and outputs raw logits directly. This is because CrossEntropyLoss (explained later) performs a softmax-equivalent computation internally, so there is no need for the model itself to produce probabilities.

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

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

class GCN(nn.Module):
    def __init__(self, num_features, hidden_channels, num_classes, dropout=0.5):
        super().__init__()
        self.conv1 = GCNConv(num_features, hidden_channels[0])
        self.conv2 = GCNConv(hidden_channels[0], hidden_channels[1])
        self.conv3 = GCNConv(hidden_channels[1], num_classes)
        self.dropout = dropout

    def forward(self, x, edge_index):
        x = self.conv1(x, edge_index)
        x = F.relu(x)
        x = F.dropout(x, p=self.dropout, training=self.training)

        x = self.conv2(x, edge_index)
        x = F.relu(x)
        x = F.dropout(x, p=self.dropout, training=self.training)

        x = self.conv3(x, edge_index)  # The final layer outputs raw logits
        return x

model = GCN(num_features=dataset.num_features,
            hidden_channels=[32, 16],
            num_classes=dataset.num_classes)

print(model)
print(f"\nNumber of learnable parameters: {sum(p.numel() for p in model.parameters()):,}")

This model progressively compresses the features in the flow 1,433 dimensions → 32 dimensions → 16 dimensions → 7 classes. With each additional layer, the receptive field described in the previous section expands from 1-hop to 2-hop to 3-hop.

5. The Training Loop and Evaluation Metrics

Training with the Adam Optimizer and CrossEntropyLoss

For training, we use the Adam optimizer and torch.nn.CrossEntropyLoss. Since CrossEntropyLoss internally combines softmax with negative log-likelihood loss in its computation, it's fine for the model's output to remain as unnormalized logits.

import torch
import torch.nn as nn
from torch_geometric.datasets import Planetoid

# Prepare the data and model (using the GCN class defined earlier)
dataset = Planetoid(root='/tmp/Cora', name='Cora')
data = dataset[0]

device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')
model = GCN(num_features=dataset.num_features,
            hidden_channels=[32, 16],
            num_classes=dataset.num_classes,
            dropout=0.5).to(device)
data = data.to(device)

optimizer = torch.optim.Adam(model.parameters(), lr=0.01, weight_decay=5e-4)
criterion = nn.CrossEntropyLoss()

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

@torch.no_grad()
def evaluate(mask):
    model.eval()
    out = model(data.x, data.edge_index)
    pred = out.argmax(dim=1)
    correct = (pred[mask] == data.y[mask]).sum().item()
    return correct / mask.sum().item()

for epoch in range(1, 201):
    loss = train()
    if epoch % 20 == 0:
        train_acc = evaluate(data.train_mask)
        val_acc = evaluate(data.val_mask)
        print(f"Epoch {epoch:03d} | Loss: {loss:.4f} | Train Acc: {train_acc:.4f} | Val Acc: {val_acc:.4f}")

test_acc = evaluate(data.test_mask)
print(f"\nTest Accuracy: {test_acc:.4f}")

Sample Output:

Epoch 020 | Loss: 1.5893 | Train Acc: 0.7357 | Val Acc: 0.6660
Epoch 040 | Loss: 0.9142 | Train Acc: 0.9214 | Val Acc: 0.7580
Epoch 060 | Loss: 0.5271 | Train Acc: 0.9714 | Val Acc: 0.7740
Epoch 080 | Loss: 0.3520 | Train Acc: 0.9857 | Val Acc: 0.7780
Epoch 100 | Loss: 0.2647 | Train Acc: 0.9929 | Val Acc: 0.7820
Epoch 120 | Loss: 0.2103 | Train Acc: 1.0000 | Val Acc: 0.7800
Epoch 140 | Loss: 0.1782 | Train Acc: 1.0000 | Val Acc: 0.7760
Epoch 160 | Loss: 0.1549 | Train Acc: 1.0000 | Val Acc: 0.7740
Epoch 180 | Loss: 0.1401 | Train Acc: 1.0000 | Val Acc: 0.7720
Epoch 200 | Loss: 0.1298 | Train Acc: 1.0000 | Val Acc: 0.7700

Test Accuracy: 0.7930

While training accuracy (Train Acc) reaches 100%, validation accuracy (Val Acc) peaks around epoch 100 and then begins to gradually decline. This is a typical sign of overfitting, which we'll address in the next section.

What Accuracy Alone Doesn't Show

Accuracy is a basic metric for model evaluation, but it can hide differences in performance across classes. In particular, when the number of samples is imbalanced across classes, the overall accuracy can appear high even if the model gets almost none of the minority-class predictions correct. Therefore, we'll also check precision, recall, and F1 score, which let us evaluate performance on a per-class basis.

• Precision: The proportion of predictions for a class that were actually correct
• Recall: The proportion of actual members of a class that were correctly predicted
• F1 Score: The harmonic mean of precision and recall

from sklearn.metrics import classification_report, confusion_matrix

@torch.no_grad()
def detailed_evaluation(mask, mask_name):
    model.eval()
    out = model(data.x, data.edge_index)
    pred = out.argmax(dim=1)

    y_true = data.y[mask].cpu().numpy()
    y_pred = pred[mask].cpu().numpy()

    print(f"=== Evaluation metrics for {mask_name} ===")
    print(classification_report(y_true, y_pred, digits=3, zero_division=0))
    print("Confusion Matrix:\n", confusion_matrix(y_true, y_pred))

detailed_evaluation(data.test_mask, "test data")

6. Countermeasures for Overfitting (Dropout, Regularization)

As we saw in the learning curve in the previous section, the phenomenon where validation accuracy plateaus (or declines) even as training accuracy reaches 100% is called overfitting. In Cora's standard split, only 140 labeled nodes (about 5% of the total) are available for training, making GCN relatively prone to overfitting in this setting. Here, we'll examine two representative countermeasures.

Dropout and Weight Decay

Dropout is a technique that randomly sets some units' outputs to zero during training, preventing the model from becoming overly dependent on specific combinations of nodes or features. We've already incorporated this into this chapter's GCN model as F.dropout.

Another representative countermeasure is weight decay, also known as L2 regularization. By adding the sum of squared weights to the loss function, this penalizes overly large weights and keeps the model smooth. In PyTorch, this can easily be specified via the weight_decay argument of the Adam optimizer.

Let's run an experiment to see how the gap between training accuracy and test accuracy (the generalization gap) changes with and without these two settings.

import torch
import torch.nn as nn
from torch_geometric.datasets import Planetoid

def run_experiment(dropout_rate, weight_decay, epochs=200, seed=0):
    torch.manual_seed(seed)
    dataset = Planetoid(root='/tmp/Cora', name='Cora')
    data = dataset[0]

    # Uses the GCN class defined in "4. Implementing the Node Classification Task"
    model = GCN(num_features=dataset.num_features,
                hidden_channels=[32, 16],
                num_classes=dataset.num_classes,
                dropout=dropout_rate)
    optimizer = torch.optim.Adam(model.parameters(), lr=0.01, weight_decay=weight_decay)
    criterion = nn.CrossEntropyLoss()

    for epoch in range(epochs):
        model.train()
        optimizer.zero_grad()
        out = model(data.x, data.edge_index)
        loss = criterion(out[data.train_mask], data.y[data.train_mask])
        loss.backward()
        optimizer.step()

    model.eval()
    with torch.no_grad():
        out = model(data.x, data.edge_index)
        pred = out.argmax(dim=1)
        train_acc = (pred[data.train_mask] == data.y[data.train_mask]).float().mean().item()
        test_acc = (pred[data.test_mask] == data.y[data.test_mask]).float().mean().item()

    return train_acc, test_acc

# Case 1: No regularization (prone to overfitting)
train_acc_no_reg, test_acc_no_reg = run_experiment(dropout_rate=0.0, weight_decay=0.0)

# Case 2: With Dropout + weight_decay
train_acc_reg, test_acc_reg = run_experiment(dropout_rate=0.5, weight_decay=5e-4)

print("Setting                  | Train Acc | Test Acc | Generalization Gap")
print("-" * 62)
print(f"No regularization        | {train_acc_no_reg:.4f}    | {test_acc_no_reg:.4f}   | {train_acc_no_reg - test_acc_no_reg:.4f}")
print(f"Dropout + weight_decay    | {train_acc_reg:.4f}    | {test_acc_reg:.4f}   | {train_acc_reg - test_acc_reg:.4f}")

Sample Output:

Setting                  | Train Acc | Test Acc | Generalization Gap
--------------------------------------------------------------
No regularization        | 1.0000    | 0.7460   | 0.2540
Dropout + weight_decay    | 1.0000    | 0.7930   | 0.2070

We can see that in the no-regularization setting, training accuracy reaches 100% while test accuracy stagnates, resulting in a larger generalization gap. Combining Dropout and L2 regularization narrows this gap to some extent.

Early Stopping

Another practical countermeasure is early stopping. This monitors the validation loss every epoch, and if it fails to improve for a certain number of consecutive epochs (patience), training is stopped and the model from the point with the lowest validation loss so far is adopted.

import copy
import torch
import torch.nn as nn
from torch_geometric.datasets import Planetoid

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

model = GCN(num_features=dataset.num_features,
            hidden_channels=[32, 16],
            num_classes=dataset.num_classes,
            dropout=0.5)
optimizer = torch.optim.Adam(model.parameters(), lr=0.01, weight_decay=5e-4)
criterion = nn.CrossEntropyLoss()

best_val_loss = float('inf')
patience = 20
patience_counter = 0
best_model_state = None

for epoch in range(1, 301):
    model.train()
    optimizer.zero_grad()
    out = model(data.x, data.edge_index)
    loss = criterion(out[data.train_mask], data.y[data.train_mask])
    loss.backward()
    optimizer.step()

    model.eval()
    with torch.no_grad():
        out = model(data.x, data.edge_index)
        val_loss = criterion(out[data.val_mask], data.y[data.val_mask]).item()

    if val_loss < best_val_loss:
        best_val_loss = val_loss
        best_model_state = copy.deepcopy(model.state_dict())
        patience_counter = 0
    else:
        patience_counter += 1

    if patience_counter >= patience:
        print(f"Epoch {epoch}: Stopping training because validation loss did not improve for {patience} epochs")
        break

# Restore the model from the point with the lowest validation loss
model.load_state_dict(best_model_state)

model.eval()
with torch.no_grad():
    out = model(data.x, data.edge_index)
    pred = out.argmax(dim=1)
    test_acc = (pred[data.test_mask] == data.y[data.test_mask]).float().mean().item()

print(f"Final test accuracy with Early Stopping: {test_acc:.4f}")

Exercises

import torch
import torch.nn as nn
import torch.nn.functional as F
from torch_geometric.data import Data
from torch_geometric.nn import GCNConv

# 5-node synthetic graph (undirected, ring structure)
edge_index = torch.tensor([
    [0, 1, 1, 2, 2, 3, 3, 4, 4, 0],
    [1, 0, 2, 1, 3, 2, 4, 3, 0, 4]
], dtype=torch.long)

torch.manual_seed(0)
x = torch.randn(5, 8)  # 5 nodes, 8-dimensional features
y = torch.tensor([0, 1, 0, 1, 0], dtype=torch.long)  # 2-class classification

data = Data(x=x, edge_index=edge_index, y=y)

class SmallGCN(nn.Module):
    def __init__(self, in_channels, hidden_channels, num_classes):
        super().__init__()
        self.conv1 = GCNConv(in_channels, hidden_channels)
        self.conv2 = GCNConv(hidden_channels, num_classes)

    def forward(self, x, edge_index):
        x = F.relu(self.conv1(x, edge_index))
        x = F.dropout(x, p=0.3, training=self.training)
        x = self.conv2(x, edge_index)
        return x

model = SmallGCN(in_channels=8, hidden_channels=4, num_classes=2)
optimizer = torch.optim.Adam(model.parameters(), lr=0.05, weight_decay=1e-4)
criterion = nn.CrossEntropyLoss()

for epoch in range(100):
    model.train()
    optimizer.zero_grad()
    out = model(data.x, data.edge_index)
    loss = criterion(out, data.y)
    loss.backward()
    optimizer.step()

model.eval()
with torch.no_grad():
    pred = model(data.x, data.edge_index).argmax(dim=1)
    acc = (pred == data.y).float().mean().item()

print(f"Final Loss: {loss.item():.4f}")
print(f"Accuracy: {acc:.4f}")

dropout_rates = [0.0, 0.3, 0.5, 0.7]
results = {}

for rate in dropout_rates:
    train_acc, test_acc = run_experiment(dropout_rate=rate, weight_decay=5e-4)
    results[rate] = (train_acc, test_acc)
    print(f"dropout={rate}: Train Acc={train_acc:.4f}, Test Acc={test_acc:.4f}")

Summary

In this chapter, we learned about the mechanics of message passing, the mathematical background of GCN, and the actual implementation of a node classification task. Let's review the learning objectives set out at the beginning.

• ✅ Explain the mechanics of message passing: We learned the 3 steps — Message, Aggregate, Update — and how to implement them using the MessagePassing base class
• ✅ Understand the mathematical background of graph convolution: We confirmed the adjacency matrix with self-loops, symmetric normalization, and the layer-wise propagation rule \(H^{(l+1)} = \sigma(\hat{A} H^{(l)} W^{(l)})\) that uses them
• ✅ Make effective use of PyTorch Geometric's GCNConv layer: We understood the meaning and appropriate use of its main parameters (improved, cached, normalize, etc.)
• ✅ Implement a node classification task on the Cora dataset: We defined a 3-layer GCN model and actually trained it on the Cora dataset
• ✅ Evaluate a model with appropriate evaluation metrics: In addition to accuracy, we learned a multifaceted evaluation approach using precision, recall, F1 score, and the confusion matrix
• ✅ Implement countermeasures for overfitting: We implemented three representative countermeasures — Dropout, weight_decay (L2 regularization), and early stopping — and compared their effects

Reference Resources

• PyTorch Geometric Official Documentation
• PyTorch Geometric Official Tutorial: Creating Message Passing Networks
• PyTorch Geometric GitHub
• GCN Paper: Semi-Supervised Classification with Graph Convolutional Networks (Kipf & Welling, 2017)