图神经网络（GNN）入门笔记（2）——从谱域理解图卷积，ChebNet和GCN实现

一、谱域图卷积（Spectral Domain Graph Convolution）

与谱域图卷积（Spectral Domain Graph Convolution）对应的是空间域（Spatial Domain）图卷积。本节学习的谱域图卷积指的是通过频率来理解卷积的方法。

二、ChebNet

上一节结尾，我们将谱域图卷积写作：
$g *_g f)= U W U^T f$

这其实已经可以做一些任务了，例如对于一个三维点云图，特征为每个点的坐标或者法向量，进行低通滤波可以把这个三维模型变得更加平滑。
而对于另一些任务，就不大合适。例如，我们有一张论文互相引用图，特征为每篇论文的种类，想要预测其中一些未标记的论文的种类。这时，这种图卷积就暴露了一些问题，例如，参数量和图规模相关以至于不适合处理不同的图，以及该处理方法过于全局，缺乏局部信息。
于是我们重新考虑拉普拉斯矩阵。乘一次拉普拉斯矩阵，相当于每个点聚合一次邻居的信息。那么乘两次，就是聚合2跳邻居的信息。乘k次，是k-hop信息。而：
$L^k = (U\Lambda U^T)^k = U\Lambda^kU^T$
因此，我们可以这样修改信号处理方法。我们学习参数 $\theta_0,...,\theta_{k-1}$ ，代表不同距离邻居的重要程度，令 $W=\theta_0 \Lambda^0 +...+\theta_{k-1} \Lambda^{k-1}$ ，就是一个对这个任务不错的滤波函数。
不过ChebNet之所以叫ChebNet，就是因为它出于种种复杂的原因使用了一个名为Chebyshev polynomial的技巧来拟合上述的 $W$ 。具体地：
$\approx \sum_{k=0}^{K-1} \theta_k T_k(\tilde \Lambda)$
其中 $T_0(x)=1,T_1(x)=x,T_k(x)=2xT_{k-1}(x)-T_{k-2}(x)$

为什么要使用切比雪夫多项式？我看网上有些人说是为了降低复杂度，但实际上计算 $T_k(L)$ 应该并不会比计算 $L^k$ 复杂度更低。实际上应该和切比雪夫多项式在信号处理中的性质有关，由于我相关知识不足，所以暂且略过。总之，目前需要学到的是ChebNet引入的把参数量从和图中点数有关的 $O (n)$ 减少到 $O (1)$ 级别的思想。

而 $\tilde \Lambda$ 是对原 $\Lambda$ 进行放缩的值，因为切比雪夫多项式要求参数取值在 $[- 1, 1]$ 之间，所以对其进行了一个 $\tilde \Lambda=\frac{2\Lambda}{\lambda_{max}}-I$ 这样的缩放。接着代回到原式，得：
$(g*_g f)=\sum_{k=0}^{K-1} \theta_kT_k(\tilde L)f$
其中 $\tilde L = \frac{2L}{\lambda_{max}}-I$
不过等等，这么说我们还是得花费高额复杂度求特征值吗？其实也不必，因为我们可以相信神经网络参数对规模缩放的自适应性（或许也可以使用一些估计方法？），取 $l_{max} \approx 2$ 即可。那么此时， $\tilde L=D^{-\frac{1}{2}}AD^{-\frac{1}{2}}$ 。
至此，我们已经解决了上一章讨论的所有问题：

$W$ 不再与图结构相关。
不需要计算特征向量。
因为不需要计算特征向量，不需求对称性保证正交基存在，可以扩展用于有向图）。
可以拟合局部信息。

三、ChebNet的实现

接下来我们尝试实现一个简单的ChebNet。
在此我使用了小规模论文类别-引用关系数据集Cora，可以使用torch_geometric.datasets来下载这个数据集。下载时如果出现超时问题，可以参考这篇博客。
另外，实际问题中的图一般都是稀疏图，拉普拉斯矩阵也是稀疏矩阵，可以用一些稀疏矩阵乘法方法加速过程。因为本篇笔记尚未讨论相关问题（~~我还没学会~~），此处给出使用邻接矩阵进行完整矩阵乘法的实现。

卷积核：

class ChebConv(nn.Module):def __init__(self, in_channels, out_channles, K=2,use_bias=True):super(ChebConv, self).__init__()self.in_channels = in_channelsself.out_channles = out_channlesself.K = Kself.use_bias = use_biasself.weights = nn.Parameter(torch.Tensor(K, 1, in_channels, out_channles))nn.init.xavier_normal_(self.weights)if(use_bias):self.bias = nn.Parameter(torch.FloatTensor(out_channles))else:self.register_parameter('bias', None)# 切比雪夫多项式，其实可以提前实现而不必在卷积核里每次都重算一遍# 放在此处只是为了演示清晰def cheb_polynomial(self, laplacian):N = laplacian.size(0)terms = torch.zeros([self.K, N, N], device=laplacian.device, dtype=torch.float32)terms[0] = torch.eye(N, device=device, dtype=torch.float32)if(self.K == 1):return termsterms[1] = laplacianfor k in range(2, self.K):terms[k] = 2.0 * torch.mm(laplacian, terms[k-1]) - terms[k-2]return terms # [K, N, N]def forward(self, inputs, laplacian):# inputs: [B, N, in_channels]# outputs: [B, N, out_channels]terms = self.cheb_polynomial(laplacian).unsqueeze(1) # [K, 1, N, N]outputs = torch.matmul(terms, inputs) # [K, B, N, in_channels]outputs = torch.matmul(outputs, self.weights) # [K, B, N, out_channels]outputs = torch.sum(outputs, dim=0) #[B, N, out_channels]if(self.use_bias):outputs += self.biasreturn outputs

网络架构：

class ChebNet(nn.Module):def __init__(self, in_channels, hidden_channels, out_channles, K=2, droput=0.5):super(ChebNet, self).__init__()self.conv1 = ChebConv(in_channels, hidden_channels, K)self.conv2 = ChebConv(hidden_channels, out_channles, K)self.dropout = droputdef forward(self, x, laplacian):outputs1 = self.conv1(x, laplacian)outputs1 = outputs1.relu()outputs1 = F.dropout(outputs1, p=self.dropout, training=self.training) outputs2 = self.conv2(outputs1, laplacian)outputs2 = outputs2.relu()return outputs2

数据集处理及训练和测试：

# 获取数据集，Cora只有一组数据，且不分训练和测试
dataset = Planetoid(root='./Cora', name='Cora', transform=NormalizeFeatures())
data = dataset[0]device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')# 定义模型
model = ChebNet(in_channels=dataset.num_node_features,hidden_channels=16,out_channles=dataset.num_classes,K=3, droput=0).to(device)# 邻接矩阵
def edge_index_to_adj(edge_index, num_nodes):  # 构建一个大小为 (num_nodes, num_nodes) 的零矩阵  adj = torch.zeros(num_nodes, num_nodes, dtype=torch.float32)# 使用索引广播机制，一次性将边索引映射到邻接矩阵的相应位置上  adj[edge_index[0], edge_index[1]] = 1.adj[edge_index[1], edge_index[0]] = 1.return adjdef get_laplacian(adj, use_normalize=True):I = torch.eye(adj.size(0), device=adj.device, dtype=adj.dtype)if(use_normalize):D = torch.diag(torch.sum(adj, dim=-1) ** (-1 / 2))L = I - torch.mm(torch.mm(D, adj), D)else:D = torch.diag(torch.sum(adj, dim=-1))L = D - adjL -= Ireturn L# 提前计算拉普拉斯矩阵
adj = edge_index_to_adj(data.edge_index, data.num_nodes)
laplacian = get_laplacian(adj).to(device)# 定义损失函数和优化器
loss = torch.nn.CrossEntropyLoss()
optimizer = torch.optim.Adam(params=model.parameters(), lr = 0.001, weight_decay=5e-4)def train():model.train()optimizer.zero_grad()  # 梯度清理y_hat = model(data.x.unsqueeze(0).to(device), laplacian).squeeze(0)l = loss(y_hat[data.train_mask].to(device), data.y[data.train_mask].to(device))l.backward() # 误差反向传播optimizer.step()return ldef test():model.eval()pred = model(data.x.unsqueeze(0).to(device), laplacian).squeeze(0)pred = pred.argmax(dim=1)test_correct = pred[data.test_mask] == data.y[data.test_mask].to(device)test_acc = int(test_correct.sum()) / int(data.test_mask.sum())return test_accepoch = 501
for epoch in range(1, epoch):train_loss = train()test_acc = test()if epoch % 10 == 0:print(f"epoch:{epoch}  loss:{train_loss}  test_acc:{test_acc}")

四、从ChebNet的first-order近似到GCN：从谱域理解GCN

对于ChebNet，在 $K = 2$ 的情况下，卷积的定义为：
$g_\theta * f=\theta_0 f - \theta_1 D^{-\frac{1}{2}}AD^{-\frac{1}{2}}$
如果我们进一步减少参数量，取 $\theta_1=-\theta_0$ ，则此时有：
$g_\theta * f=\theta (I+D^{-\frac{1}{2}}AD^{-\frac{1}{2}})f$

到这一步已经离GCN很近了，只缺少最后一点：Renormalization Trick。
在上式中， $I+D^{-\frac{1}{2}}AD^{-\frac{1}{2}}$ 的特征值范围在 $[0, 2]$ 之间（可以使用盖尔圆理论进行估计），大于1的数反复相乘有可能引起数值不稳定和梯度爆炸的问题。而Renormalization Trick就是通过对图中所有点添加自环，然后再统一归一化，来估计一个特征值在 $[- 1, 1]$ 之间的 $\tilde L$ 。具体地：
$\tilde A = A+I, \tilde D_{ii}=\sum_j \tilde A_{ij}, \tilde L =\tilde D^{-\frac{1}{2}}\tilde A\tilde D^{-\frac{1}{2}}$

盖尔圆理论
若 $A=(a_{ij}) \in \mathbb{C}$ ，则第 $i$ 的盖尔圆为以 $a_{ii}$ 为圆心，以 $R_i=\sum_{j \not=i} |a_{ij}|$ 为半径的圆。 $A$ 的所有特征值落在盖尔圆的并集之内。
$I+D^{-\frac{1}{2}}AD^{-\frac{1}{2}}$ 的每一个盖尔圆都是以1为圆心，以1为半径，特征值取值范围： $[0, 2]$ 。而 $\tilde L$ 的盖尔圆是以 $\frac{1}{D_{ii}}$ 为圆心，以 $1-\frac{1}{D_{ii}}$ 为半径的，特征值取值范围： $[- 1, 1]$ 。

好的，恭喜你也学会GCN了，让我们来实现吧！

五、GCN的实现

卷积核：

class GCNConv(nn.Module):def __init__(self, in_channels, out_channles, use_bias=True):super(GCNConv, self).__init__()self.in_channels = in_channelsself.out_channles = out_channlesself.use_bias = use_biasself.weights = nn.Parameter(torch.Tensor(in_channels, out_channles))nn.init.xavier_normal_(self.weights)if(use_bias):self.bias = nn.Parameter(torch.FloatTensor(out_channles))else:self.register_parameter('bias', None)def forward(self, inputs, laplacian):# inputs: [B, N, in_channels]# outputs: [B, N, out_channels]outputs = torch.matmul(inputs, self.weights) # [B, N, out_channels]outputs = torch.matmul(laplacian, outputs) # [B, N, out_channels]if(self.use_bias):outputs += self.biasreturn outputs

网络：

class GCN(nn.Module):def __init__(self, in_channels, hidden_channels, out_channles, droput=0.5):super(GCN, self).__init__()self.conv1 = GCNConv(in_channels, hidden_channels)self.conv2 = GCNConv(hidden_channels, out_channles)self.dropout = droputdef forward(self, x, laplacian):outputs1 = self.conv1(x, laplacian)outputs1 = outputs1.relu()outputs1 = F.dropout(outputs1, p=self.dropout, training=self.training) outputs2 = self.conv2(outputs1, laplacian)outputs2 = F.softmax(outputs2, dim=2)return outputs2

renormalization trick后的拉普拉斯矩阵计算：

def get_laplacian(adj):I = torch.eye(adj.size(0), device=adj.device, dtype=adj.dtype)adj = adj + ID = torch.diag(torch.sum(adj, dim=-1) ** (-1 / 2))L = torch.mm(torch.mm(D, adj), D)return L

除了自己手写GCN卷积核以外，还可以使用torch_geometric中的GCNConv实现，此时传入的不再是拉普拉斯矩阵，而是所有的边集（edge index），包中会用一些针对拉普拉斯矩阵的稀疏性质的方法加速计算：

class GCN(nn.Module):def __init__(self, in_channels, hidden_channels, out_channles, droput=0.5):super(GCN, self).__init__()self.conv1 = GCNConv(in_channels, hidden_channels)self.conv2 = GCNConv(hidden_channels, out_channles)self.dropout = droputdef forward(self, x, edge_index):outputs1 = self.conv1(x, edge_index)outputs1 = outputs1.relu()outputs1 = F.dropout(outputs1, p=self.dropout, training=self.training) outputs2 = self.conv2(outputs1, edge_index)outputs2 = F.softmax(outputs2, dim=2)return outputs2