import os
import torch
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams.update({'font.size': 16})
N=20
from torch_geometric.nn import GCNConv, ChebConv
from torch_geometric.datasets import Planetoid
import torch.nn.functional as FSpectral Graph Convolutions
In a non-graph setting, we would apply a CNN (Convolutional Neural Network) to a 2- or 3-D matrix. Those would be Spatial Methods (Grids).
Applying the CNN grid methods is not possible on the Graphs by default; therefore, we turn to the Spectral Methods, which deal with frequencies and waves.
The convolution theorem The convolution operation (\(\mathcal{v} * \mathcal{w}\)) from the time domain can be transformed to a multiplication operation in the frequency domain. We can then invert it back to the original time domain, resulting in the convolution output. Based on this, we use: 1. Laplacian to transform the graph to the frequency domain. 2. The Laplacian of the graph has multiple components; we separate the wave components, which are the eigenvectors and use them for transforming the graph features to the frequency domain. 3. This proves to be a slow operation, so for computational efficiency, we use Chebyshev polynomials that allow us to perform an approximate equivalent of convolution directly in the spatial domain without ever calculating \(U\), with just the Laplacian.
Fourier Transform
It is a mathematical function \(\mathcal{F}\), that decomposes a function into its constituent frequencies.
Transformation \(\mathcal{F}: \mathbb{R}^N \rightarrow \mathbb{R}^N\). It has the following property that makes it desirable to us:
\(\mathcal{F}^{-1}(\mathcal{F(\nu)}) = \mathcal{v}\), i.e. \(\mathcal{F}\) is an Invertible function.
\(\mathcal{F}(\mathcal{v} * \mathcal{w}) = \mathcal{F}(\mathcal{v}) \cdot \mathcal{F}(\mathcal{w})\), i.e. It transforms convolution function into matrix multiplication function.
This implies \[\begin{align*} v * w &= \mathcal F^{-1}(\mathcal F (v * w))\\ &= \mathcal F^{-1}(\mathcal F(v) \cdot \mathcal F(w)) \end{align*}\]
Computing Fourier Transform:
\[\mathcal{F}(\mathcal{v}) = U \cdot\mathcal{v},\;\; \mathcal{F}^{-1} = \frac{1}{N} U^{H} \cdot \mathcal{v} \tag{i}\]
where the \(N\times N\) matrix \(U\) is defined as \[ \begin{align*} \\ U = \begin{bmatrix} u_0(0) & u_1(0) & \dots & u_{N-1}(0)\\ u_0(1) & u_1(1) & \dots & u_{N-1}(1)\\ \vdots & \vdots& & \vdots\\ u_0(N-1) & u_1(N-1) & \dots & u_{N-1}(N-1)\\ \end{bmatrix} \end{align*} \tag{ii} \]
And \(u_n\) is calculated as:
\[u_n(x)= \cos\left(2 \pi \frac{n}{N} x\right) - i \sin\left(2 \pi \frac{n}{N} x\right), \tag{iii}\]
# implementation of u_{n} and the function to generate U i.e. eq ii
def matrix_U(N):
# generates tensors for x in 0 to N-1 for eq (iii),
u = lambda n, N: np.cos(2 * np.pi/N * np.arange(N)) - 1j * np.sin(2 * np.pi/N * n * np.arange(N))
U = np.empty((N, 0))
for n in range(N):
U = np.c_[U, u(n, N)] #stack each u as a column in U
return U
# implementing
def fourier_transform(v):
N = v.shape[0]
U = matrix_U(N)
return U @ v
def inverse_fourier_transform(v):
N = v.shape[0]
U = matrix_U(N)
return (U.conj().transpose() @ v) / NGraph Laplacian
It is computed as: \[L = D - A\]
where, \(D\): The Degree Matrix, \(A\): Adjacency Matrix
Laplacian is a measure of “smoothness”. When applied to a graph signal, it informs how much a node differs from the average of its neighbours. - Low Frequency: The signal changes slowly across the graph. (neighbors have very similar values) - High Frequency: The signal changes rapidly. (neighbors have very different values)
Convolution on a Graph
We using the graph’s “waves” (the matrix \(U\)), we can define convolution as: \[ \mathcal{v} * \mathcal{w} = U ((U^{T}\mathcal{w}) \odot (U^{T}\mathcal{v}))\]
The Major Problem: Calculating \(U\) (eigenvectors) requires decomposing a huge matrix \((N×N)\). For a graph with millions of nodes, this is computationally impossible.
Chebyshev Filters
\[\mathcal{v} * \mathcal{w} \approx \sum_{k=0}^{K} \theta_k \cdot T_k(\tilde{L})x\]
\[T_{k+1}(x) + T_{k-1}(x) = 2x \cdot T_k(x)\]
\[T_{k+1}(x) = 2x \cdot T_k(x) + T_{k-1}(x)\]
\[\implies T_{k}(x) = 2x \cdot T_{k-1}(x) - T_{k-2}(x)\]
Using the first order Approximation i.e. K=1
dataset = Planetoid(root="tutorial1", name="Cora")
print(f"Cora Graphs, {len(dataset)}")
print(f"Classes: {dataset.num_classes}")
print(f"Features: {dataset.num_features}")Cora Graphs, 1
Classes: 7
Features: 1433
GCN (Using K=1)
This type of GCN uses the following propagation proposed by (Kipf & Welling, 2017).
\[H^{(l+1)} = \sigma \bigg( \tilde{D}^{-\frac{1}{2}} \tilde{A}\tilde{D}^{-\frac{1}{2}} H^{(l)} W^{(l)}\bigg)\]
We get this by:
\[g_{\theta} * x \approx \sum_{k=0}^{K} \theta_k \cdot T_k(\tilde{L})x\]
at \(k=0\), \(g_{\theta} * x \approx \theta_{0} \cdot T_{0}(\tilde{L}) = I\)
at \(k=1\), \(g_{\theta} * x \approx \theta_{1} \cdot T_{1}(\tilde{L}) = \tilde{L}\)
The equation becomes:
\[g_{\theta} * x \approx \theta_{0} \cdot I \cdot x + \theta_{1}\cdot\tilde{L}\cdot x \tag{i}\]
\(\lambda\) represents the Eigenvalues (frequencies) of the Laplacian: - \(\lambda_{min}\) = 0: (Lowest frequency, i.e. constant signal) - \(\lambda_{max} =\): Highest frequency (sharpest changes between neighbours)
The Chebyshev polynomials are only valid and stable if the input \(x\) is between -1 and 1. However, eigenvalues are between 0 to \(\lambda_{max}\). We fix it by \(\frac{1}{\lambda_{max}}\).
\(\tilde{L} = \frac{2}{\lambda_{max}} L - I\), It was shown by (Kipf & Welling, 2017) that \(\lambda_{max} \approx 2\).
substituting in \((\text{i})\)
\[g_{\theta} * x \approx \theta_{0} \cdot x + \theta_{1}\cdot (L - I) \cdot x\]
since, \(L = I - D^{-\frac{1}{2}} A D^{-\frac{1}{2}}\)
\[g_{\theta} * x \approx \theta_{0} \cdot x + \theta_{1}\cdot (I - D^{-\frac{1}{2}} A D^{-\frac{1}{2}} - I) \cdot x\]
\[g_{\theta} * x \approx \theta_{0} \cdot x - \theta_{1}\cdot ( D^{-\frac{1}{2}} A D^{-\frac{1}{2}}) \cdot x\]
sharing the weight parameters \(\theta = \theta_0 = -\theta_1\)
\[g_{\theta} * x \approx \theta(I + D^{-\frac{1}{2}} A D^{-\frac{1}{2}}) \cdot x\]
We add the Identity matrix before normalisation instead of adding identity marix \(I\) separately we add the self-loop as just another edge in the graph. \(\tilde{A}\): New adjacency matrix with self loop and \(\tilde{D}\): New diagonal matrix.
The final propagation becomes:
\[ H(l+1) = \sigma\bigg( \tilde{D}^{-\frac{1}{2}} \tilde{A} \tilde{D}^{-\frac{1}{2}} H^{l} W^{l}\bigg)\]
## Model using GCNN
class GCNN_Network(torch.nn.Module):
def __init__(self):
super(GCNN_Network, self).__init__()
self.hid = 32
self.conv1 = GCNConv(dataset.num_features, self.hid)
self.conv2 = GCNConv(self.hid, dataset.num_classes)
def forward(self, data):
x, edge_index = data.x, data.edge_index
x = F.dropout(x, p=0.6, training=self.training)
# print(f"Input shape Conv1: {x.shape}")
x = self.conv1(x, edge_index)
# print(f"Output shape Conv1, hidden->{self.hid} : output->{x.shape}")
x = F.elu(x)
x = F.dropout(x, p=0.6, training=self.training)
x = self.conv2(x, edge_index)
# print(f"Output shape Conv2, heads : {x.shape}")
return F.log_softmax(x, dim=1)model = GCNN_Network()
optimizer = torch.optim.Adam(model.parameters(), lr=0.005, weight_decay=5e-4)
data = dataset[0]for epoch in range(1000):
model.train()
optimizer.zero_grad()
out = model(data)
loss = F.nll_loss(out[data.train_mask], data.y[data.train_mask])
if epoch%200 == 0:
print(loss)
loss.backward()
optimizer.step()tensor(1.9614, grad_fn=<NllLossBackward0>)
tensor(0.0526, grad_fn=<NllLossBackward0>)
tensor(0.0276, grad_fn=<NllLossBackward0>)
tensor(0.0156, grad_fn=<NllLossBackward0>)
tensor(0.0157, grad_fn=<NllLossBackward0>)
Chebyshev Conv (using K>1)
It considers more neighbours, as K>1.
## Model using GCNN
class ChebConv_Network(torch.nn.Module):
def __init__(self):
super(ChebConv_Network, self).__init__()
self.hid = 32
self.conv1 = ChebConv(dataset.num_features, self.hid, K=2)
self.conv2 = ChebConv(self.hid, dataset.num_classes, K=2)
def forward(self, data):
x, edge_index = data.x, data.edge_index
x = F.dropout(x, p=0.6, training=self.training)
# print(f"Input shape Conv1: {x.shape}")
x = self.conv1(x, edge_index)
# print(f"Output shape Conv1, hidden->{self.hid} : output->{x.shape}")
x = F.elu(x)
x = F.dropout(x, p=0.6, training=self.training)
x = self.conv2(x, edge_index)
# print(f"Output shape Conv2, heads : {x.shape}")
return F.log_softmax(x, dim=1)model = ChebConv_Network()
optimizer = torch.optim.Adam(model.parameters(), lr=0.005, weight_decay=5e-4)
data = dataset[0]for epoch in range(1000):
model.train()
optimizer.zero_grad()
out = model(data)
loss = F.nll_loss(out[data.train_mask], data.y[data.train_mask])
if epoch%200 == 0:
print(loss)
loss.backward()
optimizer.step()tensor(2.1998, grad_fn=<NllLossBackward0>)
tensor(0.0090, grad_fn=<NllLossBackward0>)
tensor(0.0222, grad_fn=<NllLossBackward0>)
tensor(0.0101, grad_fn=<NllLossBackward0>)
tensor(0.0084, grad_fn=<NllLossBackward0>)