Graph Neural Networks & Spectral/Spatial Message Passing

Standard deep learning frameworks operate over Euclidean structured data: 1D sequences (audio, natural language) or 2D grids (image pixels). In contrast, many real-world domains, including social networks, molecular chemical bonds, citation graphs, and knowledge graphs, exist in irregular, non-Euclidean topological domains.

Graph Neural Networks (GNNs) extend deep representations to arbitrary graphs by combining local topological connectivity with node and edge attribute features.

1. Graph Data Foundations & Matrix Formalism#

A graph $\mathcal{G} = (\mathcal{V}, \mathcal{E})$ is defined by:

• $\mathcal{V} = \{v_1, v_2, \dots, v_N\}$: A set of $N$ nodes (vertices).
• $\mathcal{E} \subseteq \mathcal{V} \times \mathcal{V}$: A set of $M$ edges connecting node pairs.
• $X \in \mathbb{R}^{N \times d}$: The node feature matrix where row $i$ represents the $d$-dimensional feature vector $x_i$ of node $v_i$.

Matrix Representations#

1. Adjacency Matrix $A \in \{0, 1\}^{N \times N}$:

$A_{i,j} = 1$ if an edge exists between $v_i$ and $v_j$, and $0$ otherwise. For undirected graphs without self-loops, $A$ is symmetric with zeros on the diagonal.

1. Degree Matrix $D \in \mathbb{R}^{N \times N}$:

A diagonal matrix recording the degree (number of connected incident edges) of each node:

2. Spectral Graph Theory & The Graph Laplacian#

In continuous Euclidean domains, the Laplacian operator $\Delta f = \nabla^2 f = \sum_i \frac{\partial^2 f}{\partial x_i^2}$ measures the divergence of the gradient of a function, reflecting how much the value of a function at a point differs from the average of its infinitesimal neighborhood.

In discrete graph theory, the Unnormalized Graph Laplacian $L \in \mathbb{R}^{N \times N}$ is defined as:

Properties of the Graph Laplacian#

For any node signal vector $x \in \mathbb{R}^N$:

Rewriting using $D_{i,i} = \sum_j A_{i,j}$:

• $L$ is symmetric ($L = L^T$) and positive semi-definite ($x^T L x \ge 0$).
• All eigenvalues of $L$ are non-negative real numbers: $0 = \lambda_1 \le \lambda_2 \le \dots \le \lambda_N$.
• The eigenvector corresponding to $\lambda_1 = 0$ is the constant vector $\mathbf{1}$, because $L \mathbf{1} = (D - A) \mathbf{1} = D \mathbf{1} - A \mathbf{1} = \mathbf{0}$.

Normalized Symmetric Laplacian#

To eliminate sensitivity to differences in node degrees:

The elements of $\mathcal{L}_{\text{sym}}$ are:

All eigenvalues of $\mathcal{L}_{\text{sym}}$ lie strictly within the compact interval $[0, 2]$.

3. Derivation of the Kipf & Welling GCN Layer#

Spectral graph convolution is defined as the multiplication of a signal $x$ by a filter $g_\theta$ in the graph Fourier domain:

where $U$ is the matrix of eigenvectors of $\mathcal{L}_{\text{sym}}$ and $\Lambda$ is the diagonal matrix of eigenvalues.

Because eigendecomposition requires $O(N^3)$ computational complexity, Thomas Kipf and Max Welling (ICLR 2017) derived a scalable localized first-order approximation:

Step 1: Truncated Chebyshev Expansion#

Using a first-order Chebyshev polynomial expansion with maximum eigenvalue $\lambda_{\max} \approx 2$:

Constraining $\theta = \theta_0 = -\theta_1$ to constrain parameters and prevent overfitting:

Step 2: The Renormalization Trick#

The matrix operator $\left( I_N + D^{-1/2} A D^{-1/2} \right)$ has eigenvalues in the range $[0, 2]$. Stacking multiple layers causes repeated multiplication by values greater than 1, leading directly to exploding gradients and numerical instability.

To resolve this, Kipf and Welling introduced the Renormalization Trick:

Add self-loops to the adjacency matrix prior to computing the degree matrix:

Define the corresponding degree matrix $\tilde{D}$:

The symmetrically normalized operator becomes:

The entry between connected nodes $i$ and $j$ evaluates to:

This guarantees that eigenvalues of $\tilde{A}_{\text{norm}}$ are bounded within $[-1, 1]$, preventing exponential gradient growth.

Step 3: Multi-Channel Layer Propagation#

Generalizing to a feature matrix $H^{(l)} \in \mathbb{R}^{N \times d_l}$ with learnable weight matrix $W^{(l)} \in \mathbb{R}^{d_l \times d_{l+1}}$ and activation function $\sigma$:

4. The Spatial Message Passing Framework#

Justin Gilmer et al. (2017) demonstrated that most modern GNNs (GCN, GAT, GraphSAGE) can be unified under the Message Passing Neural Network (MPNN) paradigm.

For each node $v \in \mathcal{V}$, layer $l$ executes three sequential phases:

1. Message Generation: Each neighbor $u \in \mathcal{N}(v)$ emits a message vector:

1. Permutation-Invariant Aggregation: The node aggregates incoming messages using an invariant operator $\bigoplus \in \{\sum, \text{mean}, \max\}$:

1. Node State Update: The node updates its embedding combining its own prior state with the aggregated neighborhood message:

5. Pathologies in Graph Neural Networks#

Deep GNNs encounter unique architectural failure modes:

6. From-Scratch Implementation: Graph Convolutional Network (GCN)#

We implement a complete Graph Convolutional Network layer from scratch using pure Python and Rust standard libraries. The implementation creates the adjacency matrix, adds self-loops, performs symmetric degree normalization $\tilde{D}^{-1/2} \tilde{A} \tilde{D}^{-1/2}$, and computes forward activations.

Python Implementation (Pure Standard Library)#

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
"""
Graph Convolutional Network (GCN) layer from scratch
using only the Python standard library (math, random).
"""
import math
from typing import List

def matmul(A: List[List[float]], B: List[List[float]]) -> List[List[float]]:
    rows_A, cols_A = len(A), len(A[0])
    rows_B, cols_B = len(B), len(B[0])
    assert cols_A == rows_B, f"Dimension mismatch: {cols_A} != {rows_B}"

    out = [[0.0 for _ in range(cols_B)] for _ in range(rows_A)]
    for i in range(rows_A):
        for k in range(cols_A):
            a_ik = A[i][k]
            if abs(a_ik) < 1e-12:
                continue
            for j in range(cols_B):
                out[i][j] += a_ik * B[k][j]
    return out

class GCNLayer:
    def __init__(self, in_features: int, out_features: int):
        self.in_features = in_features
        self.out_features = out_features

        # Xavier initialization
        limit = math.sqrt(6.0 / (in_features + out_features))
        self.W = [
            [((i * 13 + j * 7) % 100 / 50.0 - 1.0) * limit for j in range(out_features)]
            for i in range(in_features)
        ]

    def compute_normalized_adjacency(self, A: List[List[float]]) -> List[List[float]]:
        """Computes A_norm = D_tilde^(-1/2) * (A + I) * D_tilde^(-1/2)"""
        N = len(A)
        # 1. Add self-loops: A_tilde = A + I
        A_tilde = [[A[i][j] + (1.0 if i == j else 0.0) for j in range(N)] for i in range(N)]

        # 2. Compute degree vector: d_i = sum_j A_tilde_ij
        degrees = [sum(A_tilde[i]) for i in range(N)]

        # 3. Compute D_tilde^(-1/2)
        inv_sqrt_d = [1.0 / math.sqrt(d) if d > 0 else 0.0 for d in degrees]

        # 4. Multiply D^(-1/2) * A_tilde * D^(-1/2)
        A_norm = [[0.0 for _ in range(N)] for _ in range(N)]
        for i in range(N):
            for j in range(N):
                A_norm[i][j] = inv_sqrt_d[i] * A_tilde[i][j] * inv_sqrt_d[j]

        return A_norm

    def forward(self, H: List[List[float]], A: List[List[float]]) -> List[List[float]]:
        """H: [N, in_features], A: [N, N] -> returns H_next: [N, out_features]"""
        A_norm = self.compute_normalized_adjacency(A)

        # 1. Message passing aggregate: H_agg = A_norm * H
        H_agg = matmul(A_norm, H)

        # 2. Feature projection: Z = H_agg * W
        Z = matmul(H_agg, self.W)

        # 3. Non-linear activation: ReLU
        H_next = [[max(0.0, val) for val in row] for row in Z]
        return H_next

if __name__ == "__main__":
    # Graph of 4 nodes: (0-1), (1-2), (2-3), (3-0) ring
    A = [
        [0.0, 1.0, 0.0, 1.0],
        [1.0, 0.0, 1.0, 0.0],
        [0.0, 1.0, 0.0, 1.0],
        [1.0, 0.0, 1.0, 0.0],
    ]

    # 4 nodes with 3 input features each
    H0 = [
        [1.0, 0.0, 0.5],
        [0.2, 1.0, 0.1],
        [0.0, 0.5, 1.0],
        [0.8, 0.2, 0.0],
    ]

    gcn = GCNLayer(in_features=3, out_features=2)
    H1 = gcn.forward(H0, A)

    print("Input Node Features H0:")
    for i, row in enumerate(H0):
        print(f"Node {i}: {[round(x, 2) for x in row]}")

    print("\nOutput Node Embeddings H1 (after GCN conv):")
    for i, row in enumerate(H1):
        print(f"Node {i}: {[round(x, 4) for x in row]}")

Rust Implementation (Pure Standard Library)#

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
// Graph Convolutional Network (GCN) layer from scratch in pure Rust standard library.
// Computes symmetric renormalization D^(-1/2) (A + I) D^(-1/2) and forward propagation.

fn matmul(a: &[Vec<f64>], b: &[Vec<f64>]) -> Vec<Vec<f64>> {
    let rows_a = a.len();
    let cols_a = a[0].len();
    let rows_b = b.len();
    let cols_b = b[0].len();
    assert_eq!(cols_a, rows_b, "Matrix inner dimensions must match");

    let mut out = vec![vec![0.0; cols_b]; rows_a];
    for i in 0..rows_a {
        for k in 0..cols_a {
            let a_ik = a[i][k];
            if a_ik.abs() < 1e-12 {
                continue;
            }
            for j in 0..cols_b {
                out[i][j] += a_ik * b[k][j];
            }
        }
    }
    out
}

struct GCNLayer {
    in_features: usize,
    out_features: usize,
    weights: Vec<Vec<f64>>, // [in_features x out_features]
}

impl GCNLayer {
    fn new(in_features: usize, out_features: usize) -> Self {
        let limit = (6.0 / (in_features + out_features) as f64).sqrt();
        let mut weights = vec![vec![0.0; out_features]; in_features];
        for i in 0..in_features {
            for j in 0..out_features {
                weights[i][j] = (((i * 13 + j * 7) % 100) as f64 / 50.0 - 1.0) * limit;
            }
        }
        Self { in_features, out_features, weights }
    }

    fn compute_normalized_adjacency(&self, a: &[Vec<f64>]) -> Vec<Vec<f64>> {
        let n = a.len();
        // A_tilde = A + I
        let mut a_tilde = vec![vec![0.0; n]; n];
        let mut degrees = vec![0.0; n];

        for i in 0..n {
            for j in 0..n {
                a_tilde[i][j] = a[i][j] + if i == j { 1.0 } else { 0.0 };
                degrees[i] += a_tilde[i][j];
            }
        }

        let inv_sqrt_d: Vec<f64> = degrees.iter()
            .map(|&d| if d > 0.0 { 1.0 / d.sqrt() } else { 0.0 })
            .collect();

        let mut a_norm = vec![vec![0.0; n]; n];
        for i in 0..n {
            for j in 0..n {
                a_norm[i][j] = inv_sqrt_d[i] * a_tilde[i][j] * inv_sqrt_d[j];
            }
        }
        a_norm
    }

    fn forward(&self, h: &[Vec<f64>], a: &[Vec<f64>]) -> Vec<Vec<f64>> {
        let a_norm = self.compute_normalized_adjacency(a);
        let h_agg = matmul(&a_norm, h);
        let z = matmul(&h_agg, &self.weights);

        // Apply ReLU activation
        z.into_iter()
            .map(|row| row.into_iter().map(|val| if val > 0.0 { val } else { 0.0 }).collect())
            .collect()
    }
}

fn main() {
    let a = vec![
        vec![0.0, 1.0, 0.0, 1.0],
        vec![1.0, 0.0, 1.0, 0.0],
        vec![0.0, 1.0, 0.0, 1.0],
        vec![1.0, 0.0, 1.0, 0.0],
    ];

    let h0 = vec![
        vec![1.0, 0.0, 0.5],
        vec![0.2, 1.0, 0.1],
        vec![0.0, 0.5, 1.0],
        vec![0.8, 0.2, 0.0],
    ];

    let gcn = GCNLayer::new(3, 2);
    let h1 = gcn.forward(&h0, &a);

    println!("GCN Layer 1 Output (Rust std):");
    for (i, row) in h1.iter().enumerate() {
        println!("Node {}: [{:.4}, {:.4}]", i, row[0], row[1]);
    }
}

7. Knowledge Check#

question: "Why does Kipf and Welling's GCN apply the renormalization trick (A + I) instead of directly using (I + D^(-1/2) A D^(-1/2))?"

options:

• "(Recommended) Because (I + D^(-1/2) A D^(-1/2)) has eigenvalues in [0, 2], causing numerical instability and exploding gradients in multi-layer deep networks"
• "Because adding identity matrices converts graphs from undirected to directed"
• "Because it completely eliminates the need for computing node degrees"
• "Because it guarantees that all graph cycles are removed"

explanation: "The operator (I + D^(-1/2) A D^(-1/2)) has eigenvalues in [0, 2]. Multiplying repeatedly by values up to 2 causes numerical explosion. Adding self-loops before degree normalization bounds the eigenvalues strictly within [-1, 1]."

:::

question: "What causes the oversmoothing phenomenon in deep Graph Neural Networks?"

options:

• "(Recommended) Repeatedly applying normalized adjacency matrix multiplication acts as a low-pass Laplacian smoothing filter, causing node features to converge to uniform values"
• "The learning rate decaying to zero too early"
• "The degree matrix possessing zero values on the diagonal"
• "Using cross-entropy loss instead of mean squared error"

explanation: "As layer depth increases, GCN message passing repeatedly averages neighbor representations. Mathematically, this acts as repeated low-pass filtering, causing all node embeddings across connected components to collapse into stationary vectors proportional only to node degree."

:::