图神经网络(GNN)入门笔记(2)——从谱域理解图卷积,ChebNet和GCN实现
一、谱域图卷积(Spectral Domain Graph Convolution)
与谱域图卷积(Spectral Domain Graph Convolution)对应的是空间域(Spatial Domain)图卷积。本节学习的谱域图卷积指的是通过频率来理解卷积的方法。
二、ChebNet
上一节结尾,我们将谱域图卷积写作:
(
g
∗
g
f
)
=
U
W
U
T
f
(g *_g f)= U W U^T f
(g∗gf)=UWUTf
这其实已经可以做一些任务了,例如对于一个三维点云图,特征为每个点的坐标或者法向量,进行低通滤波可以把这个三维模型变得更加平滑。
而对于另一些任务,就不大合适。例如,我们有一张论文互相引用图,特征为每篇论文的种类,想要预测其中一些未标记的论文的种类。这时,这种图卷积就暴露了一些问题,例如,参数量和图规模相关以至于不适合处理不同的图,以及该处理方法过于全局,缺乏局部信息。
于是我们重新考虑拉普拉斯矩阵。乘一次拉普拉斯矩阵,相当于每个点聚合一次邻居的信息。那么乘两次,就是聚合2跳邻居的信息。乘k次,是k-hop信息。而:
L
k
=
(
U
Λ
U
T
)
k
=
U
Λ
k
U
T
L^k = (U\Lambda U^T)^k = U\Lambda^kU^T
Lk=(UΛUT)k=UΛkUT
因此,我们可以这样修改信号处理方法。我们学习参数
θ
0
,
.
.
.
,
θ
k
−
1
\theta_0,...,\theta_{k-1}
θ0,...,θk−1,代表不同距离邻居的重要程度,令
W
=
θ
0
Λ
0
+
.
.
.
+
θ
k
−
1
Λ
k
−
1
W=\theta_0 \Lambda^0 +...+\theta_{k-1} \Lambda^{k-1}
W=θ0Λ0+...+θk−1Λk−1,就是一个对这个任务不错的滤波函数。
不过ChebNet之所以叫ChebNet,就是因为它出于种种复杂的原因使用了一个名为Chebyshev polynomial的技巧来拟合上述的
W
W
W。具体地:
W
≈
∑
k
=
0
K
−
1
θ
k
T
k
(
Λ
~
)
W \approx \sum_{k=0}^{K-1} \theta_k T_k(\tilde \Lambda)
W≈k=0∑K−1θkTk(Λ~)
其中
T
0
(
x
)
=
1
,
T
1
(
x
)
=
x
,
T
k
(
x
)
=
2
x
T
k
−
1
(
x
)
−
T
k
−
2
(
x
)
T_0(x)=1,T_1(x)=x,T_k(x)=2xT_{k-1}(x)-T_{k-2}(x)
T0(x)=1,T1(x)=x,Tk(x)=2xTk−1(x)−Tk−2(x)
为什么要使用切比雪夫多项式?我看网上有些人说是为了降低复杂度,但实际上计算 T k ( L ) T_k(L) Tk(L)应该并不会比计算 L k L^k Lk复杂度更低。实际上应该和切比雪夫多项式在信号处理中的性质有关,由于我相关知识不足,所以暂且略过。总之,目前需要学到的是ChebNet引入的把参数量从和图中点数有关的 O ( n ) O(n) O(n)减少到 O ( 1 ) O(1) O(1)级别的思想。
而
Λ
~
\tilde \Lambda
Λ~是对原
Λ
\Lambda
Λ进行放缩的值,因为切比雪夫多项式要求参数取值在
[
−
1
,
1
]
[-1,1]
[−1,1]之间,所以对其进行了一个
Λ
~
=
2
Λ
λ
m
a
x
−
I
\tilde \Lambda=\frac{2\Lambda}{\lambda_{max}}-I
Λ~=λmax2Λ−I这样的缩放。接着代回到原式,得:
(
g
∗
g
f
)
=
∑
k
=
0
K
−
1
θ
k
T
k
(
L
~
)
f
(g*_g f)=\sum_{k=0}^{K-1} \theta_kT_k(\tilde L)f
(g∗gf)=k=0∑K−1θkTk(L~)f
其中
L
~
=
2
L
λ
m
a
x
−
I
\tilde L = \frac{2L}{\lambda_{max}}-I
L~=λmax2L−I
不过等等,这么说我们还是得花费高额复杂度求特征值吗?其实也不必,因为我们可以相信神经网络参数对规模缩放的自适应性(或许也可以使用一些估计方法?),取
l
m
a
x
≈
2
l_{max} \approx 2
lmax≈2即可。那么此时,
L
~
=
D
−
1
2
A
D
−
1
2
\tilde L=D^{-\frac{1}{2}}AD^{-\frac{1}{2}}
L~=D−21AD−21。
至此,我们已经解决了上一章讨论的所有问题:
- W W 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_channels
self.out_channles = out_channles
self.K = K
self.use_bias = use_bias
self.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 terms
terms[1] = laplacian
for 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.bias
return 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 = droput
def 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 adj
def 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 - adj
L -= I
return 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 l
def 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_acc
epoch = 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
K=2
K=2的情况下,卷积的定义为:
g
θ
∗
f
=
θ
0
f
−
θ
1
D
−
1
2
A
D
−
1
2
g_\theta * f=\theta_0 f - \theta_1 D^{-\frac{1}{2}}AD^{-\frac{1}{2}}
gθ∗f=θ0f−θ1D−21AD−21
如果我们进一步减少参数量,取
θ
1
=
−
θ
0
\theta_1=-\theta_0
θ1=−θ0,则此时有:
g
θ
∗
f
=
θ
(
I
+
D
−
1
2
A
D
−
1
2
)
f
g_\theta * f=\theta (I+D^{-\frac{1}{2}}AD^{-\frac{1}{2}})f
gθ∗f=θ(I+D−21AD−21)f
到这一步已经离GCN很近了,只缺少最后一点:Renormalization Trick。
在上式中,
I
+
D
−
1
2
A
D
−
1
2
I+D^{-\frac{1}{2}}AD^{-\frac{1}{2}}
I+D−21AD−21的特征值范围在
[
0
,
2
]
[0,2]
[0,2]之间(可以使用盖尔圆理论进行估计),大于1的数反复相乘有可能引起数值不稳定和梯度爆炸的问题。而Renormalization Trick就是通过对图中所有点添加自环,然后再统一归一化,来估计一个特征值在
[
−
1
,
1
]
[-1,1]
[−1,1]之间的
L
~
\tilde L
L~。具体地:
A
~
=
A
+
I
,
D
~
i
i
=
∑
j
A
~
i
j
,
L
~
=
D
~
−
1
2
A
~
D
~
−
1
2
\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+I,D~ii=j∑A~ij,L~=D~−21A~D~−21
盖尔圆理论
若 A = ( a i j ) ∈ C A=(a_{ij}) \in \mathbb{C} A=(aij)∈C,则第 i i i的盖尔圆为以 a i i a_{ii} aii为圆心,以 R i = ∑ j ≠ i ∣ a i j ∣ R_i=\sum_{j \not=i} |a_{ij}| Ri=∑j=i∣aij∣为半径的圆。 A A A的所有特征值落在盖尔圆的并集之内。
I + D − 1 2 A D − 1 2 I+D^{-\frac{1}{2}}AD^{-\frac{1}{2}} I+D−21AD−21的每一个盖尔圆都是以1为圆心,以1为半径,特征值取值范围: [ 0 , 2 ] [0,2] [0,2]。而 L ~ \tilde L L~的盖尔圆是以 1 D i i \frac{1}{D_{ii}} Dii1为圆心,以 1 − 1 D i i 1-\frac{1}{D_{ii}} 1−Dii1为半径的,特征值取值范围: [ − 1 , 1 ] [-1,1] [−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_channels
self.out_channles = out_channles
self.use_bias = use_bias
self.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.bias
return 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 = droput
def 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 + I
D = 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 = droput
def 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