图卷积神经网络(GCN)综述与实现(PyTorch版)
本文的实验环境为PyTorch = 1.11.0 + cu113,PyG = 2.0.4,相关依赖库和数据集的下载请见链接 。
一、图卷积神经网络介绍
1.1 传统图像卷积
卷积神经网络中的卷积(Convolution)指的是在图像上进行的输入和卷积核之间 离散内积运算 ,其本质上就是利用共享参数 的滤波器,通过计算中心值以及相邻节点的值进行加权 获得带有局部空间特征 的特征提取器。
其具有三个重要的特征,分别为:
稀疏连接
相较于全连接层,卷积层输入和输出间的连接是稀疏的,能够大大减少参数的数量,加快网络的训练速度。
参数共享
卷积核的权重参数可以被多个函数或操作共享,这样只需要训练一个参数集,而不需要对每个位置都训练一个参数集。此外,由于卷积核的大小一般是小于输入大小,也能起到减少参数数量的作用
等变表示
事实上,每个卷积层可以通过多个卷积核来进行特征提取,并且在卷积运算后,卷积神经网络对输入的图像具有平移不变性(有严格的数学论证)
一般在卷积层后,会通过一个池化层 进行降维,进一步降低网络的复杂度和非线性程度。在那之后,可以将通过卷积池化层后的特征输入全连接层或反卷积层进行进一步的分类、分割、重建工作。当然,传统的卷积操作一般适用于结构化数据。
1.2 图结构
图作为一种典型的非结构化非线性数据(非欧几里得数据),其可以表示一对一、一对多、多对多的关系,因而常被用于描述复杂的数据对象,譬如社交网络、知识图谱、城市路网、3D点云等。与结构化数据不同,图的局部输入维度可变,即每个节点的邻居节点数量不同;图具有无序性,即节点间并不存在先后关系,仅存在连接关系(点云是置换不变性,在无序性的基础上,交换两点或多点不会影响整体结果)。由于图结构的特殊性,传统CNN和RNN对其的表征能力并不理想。
对于图结构,我们可以将其抽象表示为:
G = ( V , E ) G=(V,E)
G = ( V , E )
在这里V V V 表示图中节点的集合,而E E E 为边的集合。对于图特征,我们一般有三个重要矩阵进行表示。
邻接矩阵A A A :adjacency matrix用来表示节点间的连接关系。
$$
\begin{vmatrix}
&0 & 1 & 0 & 0&\\
&0 & 0 & 1 & 0&\\
&1 & 1 & 0 & 0&\\
&0 & 1 & 1 & 0&\\
\end{vmatrix}
$$
对于带权的图,邻接矩阵将把1替换为对应的权重
度矩阵D D D :degree matrix用来表示节点的连接数,可以表征某个节点在图中的重要程度,是一个对角矩阵,例如针对上图的入度矩阵:
∣ 1 0 0 0 0 3 0 0 0 0 2 0 0 0 0 0 ∣ \begin{vmatrix}
&1 & 0 & 0 & 0&\\
&0 & 3 & 0 & 0&\\
&0 & 0 & 2 & 0&\\
&0 & 0 & 0 & 0&\\
\end{vmatrix}
∣ ∣ 1 0 0 0 0 3 0 0 0 0 2 0 0 0 0 0 ∣ ∣
特征矩阵X X X :feature matrix用来表示节点的特征
1.3 图卷积神经网络
目前主流的图卷积基本上可以分为两类,一种是基于谱的图卷积,一种是基于空域的图卷积。
基于谱的图卷积通过傅里叶变换(FFT干的一件事情就是连接空域和频域)将节点映射到频域空间,通过频域空间上的卷积来实现时域上的卷积,最后将特征映射回空域。而基于空域的图卷积则是直接基于节点与邻居进行卷积提取特征,没有做域上的变换。
图卷积算子可表示为:
h i l + 1 = σ ( ∑ j ∈ N i 1 C i j h j l w R j l ) h_i^{l+1}=\sigma(\sum_{j\in N_i}\frac{1}{C_{ij}}h_j^lw_{R_j}^l)
h i l + 1 = σ ( j ∈ N i ∑ C ij 1 h j l w R j l )
其中,设中心结点为i i i ;h i l h_i^l h i l 为结点i i i 在l l l 层的特征表达;σ \sigma σ 是非线性激活函数;C i j C_{ij} C ij 则是归一化因子,譬如结点度的倒数、反距离权重、高斯衰减权重等;N i N_i N i 是结点i i i 的邻接节点(包括自身);R i R_i R i 表示节点i i i 的类型;W R j W_{R_j} W R j 表示R j R_j R j 类型的节点变换权重参数。
下面围绕Semi-supervised Classification with Graph Convolutional Networks 一文中提出的GCN结构进行分析。
该篇文章由阿姆斯特丹(the University of Amsterdam)大学机器学习专业的Thomas Kipf 博士于2016年提出,并于2017年被深度学习的顶会ICLR(International Conference on Learning Representations)接收!这位大佬的研究方向是学习结构化数据和结构化表示/计算,包括推理、(多智能体)强化学习和结构化深度生成模型。
核心思想
该篇文章提出了一种新的网络结构,用于处理非结构化的图数据,并解决了在一个图中,只有少部分节点的标签是已知情况下的节点分类问题(半监督学习)。
对于带有特征的图结构,例如下图中的n × 3 n\times3 n × 3 网络结构,有部分是带有标识的,而有部分则是无标识的。
GCN通过考虑节点本身以及邻居节点的特征信息来提取潜在的关系。比如我们在中心节点拼接了邻居的特征后,使用平均池化的方式对这些特征进行聚合,再通过浅层网络进行学习训练得到新的数据。
每个GCN层做的事情:
获取当前节点和邻接节点特征
通过聚合函数获取局部特征(带有拓扑关系)
浅层学习训练,获取高维特征
数学推导
对于一个如上图所示的无向图,我们要怎么样才能获取到某个节点以及其邻居节点的特征呢?一种非常直观的想法💡是,在当前节点的特征之后,根据权重拼接 与之相邻节点 的特征。此时,空域上的距离往往成为了权重的影响因素。
那么,邻接表就成了我们考虑节点间拓扑关系最重要的结构。我们写成A × X × W A\times X\times W A × X × W ,A A A 是邻接矩阵,X X X 是特征矩阵,W W W 是权重。右乘相当于控制行,左乘相当于控制列,A A A 左乘X X X ,相当于在对应节点处,使用哪些节点特征构建新的特征:
∣ 0 1 0 0 0 0 1 1 0 1 0 1 1 1 0 0 ∣ × ∣ 1 1 1 1 2 2 2 2 3 3 3 3 40 40 40 40 ∣ = ∣ 2 2 2 2 43 43 43 43 42 42 42 42 3 3 3 3 ∣ \begin{vmatrix}
&0 & 1 & 0&0&\\
&0 & 0 & 1&1&\\
&0 & 1 & 0&1&\\
&1 & 1 & 0&0&\\
\end{vmatrix}\times
\begin{vmatrix}
&1 & 1 & 1&1&\\
&2 & 2 & 2&2&\\
&3 & 3 & 3&3&\\
&40 & 40 & 40&40&\\
\end{vmatrix}=
\begin{vmatrix}
&2 & 2 & 2&2&\\
&43 & 43 & 43&43&\\
&42 & 42 & 42&42&\\
&3 & 3 & 3&3&\\
\end{vmatrix}
∣ ∣ 0 0 0 1 1 0 1 1 0 1 0 0 0 1 1 0 ∣ ∣ × ∣ ∣ 1 2 3 40 1 2 3 40 1 2 3 40 1 2 3 40 ∣ ∣ = ∣ ∣ 2 43 42 3 2 43 42 3 2 43 42 3 2 43 42 3 ∣ ∣
诶,糟糕,这矩阵数数数着数着咋把自己给忘了呀!而且似乎…有些小朋友发育的过于良好,这样会导致梯度消化不良的!(梯度爆炸或消失问题)
为了处理这种情况,我们引入了新的邻接矩阵A ~ \tilde A A ~ 用于承载当前节点信息~
A ~ = A + λ I N \tilde A =A+\lambda I_N
A ~ = A + λ I N
当权重因子λ \lambda λ 取1 1 1 时,A ~ = A + I N \tilde A =A+I_N A ~ = A + I N ,此时意味着节点中心和邻居一样重要啦。🎉
值得注意的是,λ \lambda λ 本身也是一个可以由训练得到的参数。
针对梯度问题呢,我们可以利用对角矩阵实现标准化,譬如利用度矩阵的逆D ~ − 1 \tilde D^{-1} D ~ − 1 进行平均化的操作:
∣ 2 0 0 0 0 3 0 0 0 0 4 0 0 0 0 1 ∣ → i n v e r s e ∣ 1 / 2 0 0 0 0 1 / 3 0 0 0 0 1 / 4 0 0 0 0 1 ∣ \begin{vmatrix}
&2 & 0 & 0&0&\\
&0 & 3 & 0&0&\\
&0 & 0 & 4&0&\\
&0 &0 & 0&1&\\
\end{vmatrix}\xrightarrow{inverse}
\begin{vmatrix}
&1/2 & 0 & 0&0&\\
&0 & 1/3 & 0&0&\\
&0 & 0 & 1/4&0&\\
&0 &0 &0&1&\\
\end{vmatrix}
∣ ∣ 2 0 0 0 0 3 0 0 0 0 4 0 0 0 0 1 ∣ ∣ in v erse ∣ ∣ 1/2 0 0 0 0 1/3 0 0 0 0 1/4 0 0 0 0 1 ∣ ∣
于是乎,新的公式诞生啦:
D ~ − 1 ( A ~ X ) W \tilde D^{-1}(\tilde AX)W
D ~ − 1 ( A ~ X ) W
但我们在这里仅仅对行进行了标准化,在列方向上并没有被标准化,所以,修正后的公式为:
( D ~ − 1 A ~ D ~ − 1 ) X W L e t ‘ s c a l l A ^ = ( D ~ − 1 A ~ D ~ − 1 ) (\tilde D^{-1}\tilde A\tilde D^{-1})XW
\\
Let‘s\quad call\quad \widehat{A}=(\tilde D^{-1}\tilde A\tilde D^{-1})
( D ~ − 1 A ~ D ~ − 1 ) X W L e t ‘ s c a ll A = ( D ~ − 1 A ~ D ~ − 1 )
当然了,行做了一次标准化,列做了一次标准化,那么对于元素X i j X_{ij} X ij 来说,那就是做了两次标准化诶!我们再开个方吧:happy:
( D ~ − 1 / 2 A ~ D ~ − 1 / 2 ) X W L e t ‘ s c a l l A ^ = ( D ~ − 1 / 2 A ~ D ~ − 1 / 2 ) (\tilde D^{-1/2}\tilde A\tilde D^{-1/2})XW\\
Let‘s\quad call\quad \widehat{A}=(\tilde D^{-1/2}\tilde A\tilde D^{-1/2})
( D ~ − 1/2 A ~ D ~ − 1/2 ) X W L e t ‘ s c a ll A = ( D ~ − 1/2 A ~ D ~ − 1/2 )
于是乎,最终可以得到以下的标准格式:
Z = f ( X , A ) = s o f t m a x ( A ^ R e L U ( A ^ X W ( 0 ) ) W ( 1 ) ) Z=f(X,A)=softmax(\widehat A ReLU(\widehat AXW^{(0)})W^{(1)})
Z = f ( X , A ) = so f t ma x ( A R e LU ( A X W ( 0 ) ) W ( 1 ) )
W W W 是可训练的权重,A ^ \widehat A A 是标准化后的邻接矩阵。
此时我们在回过头来看这个图卷积算子,是不是有些恍然大悟了呢🤔
h i l + 1 = σ ( ∑ j ∈ N i 1 C i j h j l w R j l ) h_i^{l+1}=\sigma(\sum_{j\in N_i}\frac{1}{C_{ij}}h_j^lw_{R_j}^l)
h i l + 1 = σ ( j ∈ N i ∑ C ij 1 h j l w R j l )
Thomas经过实验比对,浅层的GCN(一般为2层)取得的效果往往要比没有加入残差块的深层GCN来的好。GCN的网络层数代表着节点特征所能到达的最远距离,或者说节点信息的传播距离。深层的网络可能会令节点信息传播到整个网络,反而效果不那么好。
举个简单的栗子~🐈
假设有邻接矩阵A A A :
A = ∣ 0 1 0 0 1 0 1 1 0 1 0 1 0 1 1 1 ∣ A=\begin{vmatrix}
&0&1&0&0&\\
&1&0&1&1&\\
&0&1&0&1&\\
&0&1&1&1&\\
\end{vmatrix}
A = ∣ ∣ 0 1 0 0 1 0 1 1 0 1 0 1 0 1 1 1 ∣ ∣
对于节点A 1 A_1 A 1 来说,跟他接轨的只有节点A 2 A_2 A 2 ,那么在第一层计算后,节点A 1 A_1 A 1 聚合了来自节点A 2 A_2 A 2 的特征。而对于节点A 2 A_2 A 2 ,与之接轨的则有节点A 1 , A 3 , A 4 A_1,A_3,A_4 A 1 , A 3 , A 4 ,经过第一层后,该节点聚合了所有节点的信息。那么在第二层,由于图结构不是动态更新的,A 1 A_1 A 1 又会聚合A 2 A_2 A 2 的特征信息,但此时的A 2 A_2 A 2 已经聚合了A 3 , A 4 A_3,A_4 A 3 , A 4 的特征信息!换言之,A 1 A_1 A 1 感受到了来自远方的召唤!(A 3 , A 4 A_3,A_4 A 3 , A 4 的信息经过两层网络后传递到了A 1 A_1 A 1 处,传递过程类似于并查集啦)
最后对于这些分类任务,采用交叉熵函数作为损失函数,计算获取最可能的情况即可。
二、基于PyG的图卷积神经网络半监督分类
2.1 模型准备
模块导入阶段,在本实验中,networkx主要用于图结构可视化,argparse用于参数的读写
1 2 3 4 5 6 7 8 9 10 11 12 13 14 import numpy as npimport pandas as pdimport matplotlib.pyplot as pltimport torchimport torch.nn.functional as Ffrom torch_geometric.nn import GCNConvfrom torch_geometric.datasets import Planetoidfrom torch_geometric.utils import to_networkximport networkx as nxfrom sklearn.metrics import accuracy_scorefrom sklearn.manifold import TSNEfrom sklearn.svm import SVCfrom sklearn.semi_supervised import _label_propagationimport argparse
2.2 数据探查
我们导入Cora数据集。
Cora数据集由2078篇机器学习领域的论文构成,每个样本点都是一篇论文,这些论文主要分为了7个类别,分别为基于案例、遗传算法、神经网络、概率方法、强化学习、规则学习与理论。在该数据集中,每篇论文都至少引用了该数据集中的另一篇论文,对每个节点所代表的论文,都由一个1433维的词向量表示,即该图上每个节点都具有1433个特征,词向量的每个元素都对应一个词,且该元素仅有0或1两个取值,取0表示该元素对应的词不在论文中,取1表示在。
Cora参数
ind.cora.x : 训练集节点特征向量,保存对象为:scipy.sparse.csr.csr_matrix,实际展开后大小为: (140, 1433)
ind.cora.tx : 测试集节点特征向量,保存对象为:scipy.sparse.csr.csr_matrix,实际展开后大小为: (1000, 1433)
ind.cora.allx : 包含有标签和无标签的训练节点特征向量,保存对象为:scipy.sparse.csr.csr_matrix,实际展开后大小为:(1708, 1433),可以理解为除测试集以外的其他节点特征集合,训练集是它的子集
ind.cora.y : one-hot表示的训练节点的标签,保存对象为:numpy.ndarray
ind.cora.ty : one-hot表示的测试节点的标签,保存对象为:numpy.ndarray
ind.cora.ally : one-hot表示的ind.cora.allx对应的标签,保存对象为:numpy.ndarray
ind.cora.graph : 保存节点之间边的信息,保存格式为:{ index : [ index_of_neighbor_nodes ] }
ind.cora.test.index : 保存测试集节点的索引,保存对象为:List,用于后面的归纳学习设置。
通过以下语句下载或读取,若当前路径下没有找到文件,则会自动下载。
1 2 3 4 5 6 7 8 dataset=Planetoid(root=r"./Cora" ,name="Cora" ) print ("网络数据包含的类数量:" ,dataset.num_classes)print ("网络数据边的特征数量:" ,dataset.num_edge_features)print ("网络数据边的数量:" ,dataset.data.edge_index.shape[1 ]/2 ) print ("网络数据节点的特征数量:" ,dataset.num_node_features)print ("网络数据节点的数量:" ,dataset.data.x.shape[0 ])
输出的结果为:
1 2 3 4 5 6 7 网络数据包含的类数量: 7 网络数据边的特征数量: 0 网络数据边的数量: 5278.0 网络数据节点的特征数量: 1433 网络数据节点的数量: 2708 网络的边的数量: 5278 网络的节点的数量: 2708
我们可以看到,Cora数据集囊括了七个类共2708个节点,每个节点都带有1433个特征,一共有5278条边,边没有包含特征。
重点来看dataset.data:
1 Data(x=[2708 , 1433 ], edge_index=[2 , 10556 ], y=[2708 ], train_mask=[2708 ], val_mask=[2708 ], test_mask=[2708 ])
该数据中的x是节点和对应的特征矩阵,y是标签值,edge_index则是COO的网络连接格式(第一行表示节点顺序,第二行表示连接对象,因而每条边会出现两次),mask对象则是一个逻辑向量(与节点数量相同),表示该节点的用途。
2.3 结构可视化
为了更直观地分析数据,我们可以考虑将图结构可视化,但使用张量格式的数据十分不方便,PyG提供了to_networkx()函数,可用于将torch_geometric.data.Data对象转化为networkx库中的有向图数据,方便进一步的分析和可视化。
1 2 3 4 5 6 7 8 9 CoraNet=to_networkx(dataset.data) CoraNet=CoraNet.to_undirected() print ("网络的边的数量: " ,CoraNet.number_of_edges())print ("网络的节点的数量: " ,CoraNet.number_of_nodes())Node_class=dataset.data.y.data.numpy() print ("节点分类:" ,Node_class)
结果为:
1 2 3 网络的边的数量: 5278 网络的节点的数量: 2708 节点分类: [3 4 4 ... 3 3 3 ]
此时数据已经由Tensor格式转化为Graph格式啦,并且我们将标签y单独提取出来,方便接下来的操作。
在Graph类中,我们可以使用Graph.degree方法计算度,度越大的节点,说明跟其他节点的连通性越好,在整个网络中具有更大的贡献度,即越重要。下面我们对前三十的节点进行可视化:
1 2 3 4 5 6 7 8 9 10 Node_degree=pd.DataFrame(data=CoraNet.degree,columns=["Node" ,"Degree" ]) Node_degree=Node_degree.sort_values(by=["Degree" ],ascending=False ) Node_degree=Node_degree.reset_index(drop=True ) Node_degree.iloc[0 :30 ,:].plot(x="Node" ,y="Degree" ,kind="bar" ,figsize=(10 ,7 )) plt.xlabel("Node" ,size=12 ) plt.ylabel("Degree" ,size=12 ) plt.show()
接着借助networkx库对网络结构进行可视化,nx.spring_layout(Graph)会自动计算网络节点中的布局方式,我们需要做的事情有:
为不同类的节点指定不同的颜色
draw_networkx_nodes(Graph,pos,nodelist,nodesize,node_color,alpha)可用于以特定颜色绘制一张图里面的指定索引的节点
为网络添加边
draw_networkx_edges(Graph,pos,width,edge_color)用于绘制图结构的边
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 pos=nx.spring_layout(CoraNet) nodecolor=["red" ,"blue" ,"green" ,"yellow" ,"peru" ,"violet" ,"cyan" ] nodelabel=np.array(list (CoraNet.nodes)) plt.figure(figsize=(16 ,12 )) for ii in np.arange(len (np.unique(Node_class))): nodelist=nodelabel[Node_class==ii] nx.draw_networkx_nodes(CoraNet,pos,nodelist=list (nodelist),node_size=20 , node_color=nodecolor[ii],alpha=0.8 ) nx.draw_networkx_edges(CoraNet,pos,width=1 ,edge_color="gray" ) plt.show()
结果如下图所示:
我们可以再对训练集单独绘制一下:
1 2 3 4 5 6 7 8 9 10 11 nodelabel=np.arange(0 ,140 ) Node_class=dataset.data.y.data.numpy()[0 :140 ] plt.figure(figsize=(8 ,6 )) for ii in np.arange(len (np.unique(Node_class))): nodelist=nodelabel[Node_class==ii] nx.draw_networkx_nodes(CoraNet,pos,nodelist=list (nodelist),node_size=20 , node_color=nodecolor[ii],alpha=0.8 ) plt.show()
简单来说,训练集的分布大致体现了原始数据的分布
2.4 网络搭建
我们这边直接使用PyG的GCNConv进行搭建:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 class GCNNet (torch.nn.Module): def __init__ (self,input_feature,num_classes ): super (GCNNet, self).__init__() self.input_feature=input_feature self.num_classes=num_classes self.conv1=GCNConv(input_feature,32 ) self.conv2=GCNConv(32 ,num_classes) def forward (self,data ): x,edge_index=data.x,data.edge_index x=F.relu(self.conv1(x,edge_index)) x=F.relu(self.conv2(x,edge_index)) return F.softmax(x,dim=1 )
2.5 网络训练
在进行训练前,需要先设定好全局变量,这里我们使用argparse库进行操作。
1 2 3 4 5 6 7 8 9 10 11 def parse_args (): '''PARAMETERS''' parser=argparse.ArgumentParser("GCN" ) parser.add_argument('-e' ,'--epochs' ,type =int ,default=500 ,help ='number of epoch in training [default: 200]' ) parser.add_argument('-lr' ,'--learning_rate' ,type =float ,default=0.01 , help ='learning rate in training [default: 0.01]' ) parser.add_argument('-op' ,'--optimizer' ,type =str , default='Adam' , help ='optimizer for training [default: Adam]' ) parser.add_argument('-g' ,'--gpu' ,type =str , default='0' , help ='specify gpu device [default: 0]' ) parser.add_argument('-p' ,'--path' ,type =str , default='./save' , help ='the path of file saving [default: ./save]' ) parser.add_argument('-dr' ,'--decay_rate' , type =float , default=5e-4 , help ='decay rate [default: 5e-4]' ) return parser.parse_args()
我们可以定义一个保存路径,用于存放最优的模型参数和最优的优化器参数,譬如:
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 args = args savepath = args.path if not os.path.exists(savepath): os.makedirs(savepath) if 当前精度>最优精度: sp=savepath+"/best_model.pth" state={ "epoch" :epoch, "accuracy" :acc, "model_state_dict" :model.state_dict(), "optimizer_state_dict" :optimizer.state_dict() } torch.save(state,sp) best_acc = acc try : checkpoint=torch.load(savepath+"/best_model.pth" ) mygcn.load_state_dict(checkpoint['model_state_dict' ]) start_epoch = checkpoint['epoch' ] optimizer.load_state_dict(checkpoint['optimizer_state_dict' ]) except : mygcn=GCNNet(input_feature,num_classes) start_epoch=0
那首先进行网络的初始化
1 2 3 4 5 6 7 8 9 10 11 12 13 args = args savepath = args.path if not os.path.exists(savepath): os.makedirs(savepath) input_feature = dataset.num_node_features num_classes = dataset.num_classes mygcn = GCNNet(input_feature, num_classes) model=mygcn data = dataset[0 ] os.environ["CUDA_VISIBLE_DEVICES" ] = args.gpu
接着咱定义优化器
1 2 3 4 5 6 7 8 9 10 11 12 if args.optimizer == 'Adam' : optimizer = torch.optim.Adam( model.parameters(), lr=args.learning_rate, weight_decay=args.decay_rate ) else : optimizer = torch.optim.SGD(model.parameters(), lr=0.01 , momentum=0.9 )
再来读取模型
1 2 3 4 5 6 7 8 9 10 11 try : checkpoint=torch.load(savepath+"/best_model.pth" ) mygcn.load_state_dict(checkpoint['model_state_dict' ]) start_epoch = checkpoint['epoch' ] optimizer.load_state_dict(checkpoint['optimizer_state_dict' ]) except : mygcn=GCNNet(input_feature,num_classes) start_epoch=0 print (mygcn)
此时的输出为:
1 2 3 4 5 6 ''' GCNNet( (conv1): GCNConv(1433, 32) (conv2): GCNConv(32, 7) ) '''
接下来,我们再进行网络的训练:
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 train_loss_all=[] val_loss_all=[] best_acc=0 train_acc_all=[] model.train() for epoch in range (start_epoch,args.epochs): optimizer.zero_grad() out=model(data) loss=F.cross_entropy(out[data.train_mask],data.y[data.train_mask]) loss.backward() optimizer.step() train_loss_all.append(loss.data.cpu().numpy()) loss=F.cross_entropy(out[data.val_mask],data.y[data.val_mask]) val_loss_all.append(loss.data.cpu().numpy()) _,pred=out.max (dim=1 ) acc=float (pred[data.train_mask].eq(data.y[data.train_mask]).sum ().item())/data.train_mask.sum ().item() train_acc_all.append(acc) if acc>best_acc: sp=savepath+"/best_model.pth" state={ "epoch" :epoch, "accuracy" :acc, "model_state_dict" :model.state_dict(), "optimizer_state_dict" :optimizer.state_dict() } torch.save(state,sp) best_acc = acc if epoch%20 ==0 : print ("Epoch:" ,epoch,";Train Loss:" ,train_loss_all[-1 ],";Val Loss:" ,val_loss_all[-1 ],"Train acc:" ,train_acc_all[-1 ])
最后自然是结果的可视化啦:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 plt.figure(figsize=(10 ,6 )) plt.plot(train_loss_all,"ro-" ,label="Train Loss" ) plt.plot(val_loss_all,"bs-" ,label="Val Loss" ) plt.legend() plt.grid() plt.xlabel("epoch" ,size=13 ) plt.ylabel("loss" ,size=13 ) plt.title("Graph Convolutional Networks" ,size=14 ) plt.show() model.eval () _,pred=model(data).max (dim=1 ) correct=float (pred[data.test_mask].eq(data.y[data.test_mask]).sum ().item()) acc=correct/data.test_mask.sum ().item() print ("Accuracy: {:.4f}" .format (acc))
可以发现,到200个epoch后,模型已经收敛。当然有很大程度上是因为我们并没有随机取出数据集(用的是固定顺序),最终的精度结果为:
为了直观体现网络的特征提取能力,我们对隐藏层获得的32维特征在空间中的分布情况进行可视化,与原始数据1433维进行比较,统一使用TSNE算法降到二维。首先我们先定义一个用于TSNE绘制的函数:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 def draw (x_tsne,text ): plt.figure(figsize=(12 , 8 )) axl = plt.subplot(1 , 1 , 1 ) X = x_tsne[:, 0 ] Y = x_tsne[:, 1 ] axl.set_xlim([min (X), max (X)]) axl.set_ylim([min (Y), max (Y)]) for ii in range (x_tsne.shape[0 ]): text = data.y.data.numpy()[ii] axl.text(X[ii], Y[ii], str (text), fontsize=5 , bbox=dict (boxstyle="round" , facecolor=plt.cm.Set1(text), alpha=0.7 )) axl.set_xlabel("TSNE Feature 1" , size=13 ) axl.set_xlabel("TSNE Feature 2" , size=13 ) axl.set_title(t, size=15 ) plt.show()
我们来绘制原始数据:
1 2 x_tsne=TSNE(n_components=2 ).fit_transform(dataset.data.x.data.numpy()) draw(x_tsne,"Original feature TSNE" )
接着呢,利用hook函数获取中间层的输入,板子如下:
1 2 3 4 5 6 7 8 activation={} def get_activation (name ): def hook (model,input ,output ): activation[name]=output.detach() return hook model.conv1.register_forward_hook(get_activation("conv1" ))
1 2 3 4 5 6 7 _=model(data) conv1=activation["conv1" ].data.cpu().numpy() print ("conv1.shape: " ,conv1.shape) x_tsne=TSNE(n_components=2 ).fit_transform(conv1) draw(x_tsne,"GCNConv Feature TSNE" )
可以发现,在高维特征中挖掘出了每个类的潜在关系。
2.6 模型比较
我们再用GCN对比下传统的SVM和LP(Label Propagation)算法:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 X=dataset.data.x.data.numpy() Y=dataset.data.y.data.numpy() train_mask=dataset.data.train_mask.data.numpy() test_mask=dataset.data.test_mask.data.numpy() train_x=X[:140 ,:] train_y=Y[train_mask] test_x=X[1708 :2709 ,:] test_y=Y[test_mask] svmmodel=SVC() svmmodel.fit(train_x,train_y) prelab=svmmodel.predict(test_x) print ("SVM的预测精度:" ,accuracy_score(test_y,prelab))
最终的结果为:
1 2 3 4 5 6 7 8 9 10 train_y=Y.copy() train_y[test_mask==True ]=-1 test_y=Y[test_mask] lp_model=_label_propagation.LabelPropagation(kernel="knn" ,n_neighbors=3 ) lp_model.fit(X,train_y) prelab=lp_model.transduction_ print ("LP模型的预测精度:" ,accuracy_score(Y[test_mask],prelab[test_mask]))
结果为:
GCN在精度上远高于上面两个模型,具体原因在前面已经介绍过了。
2.7 完整代码
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 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 import numpy as npimport pandas as pdimport matplotlib.pyplot as pltimport torchimport torch.nn.functional as Ffrom torch_geometric.nn import GCNConvfrom torch_geometric.datasets import Planetoidfrom torch_geometric.utils import to_networkximport networkx as nxfrom sklearn.metrics import accuracy_scorefrom sklearn.manifold import TSNEfrom sklearn.svm import SVCfrom sklearn.semi_supervised import _label_propagationimport argparseimport osdataset=Planetoid(root=r"./Cora" ,name="Cora" ) def parse_args (): '''PARAMETERS''' parser=argparse.ArgumentParser("GCN" ) parser.add_argument('-e' ,'--epochs' ,type =int ,default=500 ,help ='number of epoch in training [default: 200]' ) parser.add_argument('-lr' ,'--learning_rate' ,type =float ,default=0.01 , help ='learning rate in training [default: 0.01]' ) parser.add_argument('-op' ,'--optimizer' ,type =str , default='Adam' , help ='optimizer for training [default: Adam]' ) parser.add_argument('-g' ,'--gpu' ,type =str , default='0' , help ='specify gpu device [default: 0]' ) parser.add_argument('-p' ,'--path' ,type =str , default='./save' , help ='the path of file saving [default: ./save]' ) parser.add_argument('-dr' ,'--decay_rate' , type =float , default=5e-4 , help ='decay rate [default: 5e-4]' ) return parser.parse_args() class GCNNet (torch.nn.Module): def __init__ (self,input_feature,num_classes ): super (GCNNet, self).__init__() self.input_feature=input_feature self.num_classes=num_classes self.conv1=GCNConv(input_feature,32 ) self.conv2=GCNConv(32 ,num_classes) def forward (self,data ): x,edge_index=data.x,data.edge_index x=F.relu(self.conv1(x,edge_index)) x=F.relu(self.conv2(x,edge_index)) return F.softmax(x,dim=1 ) def main (dataset,args ): args = args savepath = args.path if not os.path.exists(savepath): os.makedirs(savepath) input_feature = dataset.num_node_features num_classes = dataset.num_classes data = dataset[0 ] os.environ["CUDA_VISIBLE_DEVICES" ] = args.gpu if args.optimizer == 'Adam' : optimizer = torch.optim.Adam( model.parameters(), lr=args.learning_rate, weight_decay=args.decay_rate ) else : optimizer = torch.optim.SGD(model.parameters(), lr=0.01 , momentum=0.9 ) try : checkpoint=torch.load(savepath+"/best_model.pth" ) mygcn.load_state_dict(checkpoint['model_state_dict' ]) start_epoch = checkpoint['epoch' ] optimizer.load_state_dict(checkpoint['optimizer_state_dict' ]) except : mygcn=GCNNet(input_feature,num_classes) start_epoch=0 print (mygcn) model=mygcn train_loss_all=[] val_loss_all=[] best_acc=0 train_acc_all=[] model.train() for epoch in range (start_epoch,args.epochs): optimizer.zero_grad() out=model(data) loss=F.cross_entropy(out[data.train_mask],data.y[data.train_mask]) loss.backward() optimizer.step() train_loss_all.append(loss.data.cpu().numpy()) loss=F.cross_entropy(out[data.val_mask],data.y[data.val_mask]) val_loss_all.append(loss.data.cpu().numpy()) _,pred=out.max (dim=1 ) acc=float (pred[data.train_mask].eq(data.y[data.train_mask]).sum ().item())/data.train_mask.sum ().item() train_acc_all.append(acc) if acc>best_acc: sp=savepath+"/best_model.pth" state={ "epoch" :epoch, "accuracy" :acc, "model_state_dict" :model.state_dict(), "optimizer_state_dict" :optimizer.state_dict() } torch.save(state,sp) best_acc = acc if epoch%20 ==0 : print ("Epoch:" ,epoch,";Train Loss:" ,train_loss_all[-1 ],";Val Loss:" ,val_loss_all[-1 ],"Train acc:" ,train_acc_all[-1 ]) plt.figure(figsize=(10 ,6 )) plt.plot(train_loss_all,"ro-" ,label="Train Loss" ) plt.plot(val_loss_all,"bs-" ,label="Val Loss" ) plt.legend() plt.grid() plt.xlabel("epoch" ,size=13 ) plt.ylabel("loss" ,size=13 ) plt.title("Graph Convolutional Networks" ,size=14 ) plt.show() model.eval () _,pred=model(data).max (dim=1 ) correct=float (pred[data.test_mask].eq(data.y[data.test_mask]).sum ().item()) acc=correct/data.test_mask.sum ().item() print ("Accuracy: {:.4f}" .format (acc)) def draw (x_tsne,t ): plt.figure(figsize=(12 , 8 )) axl = plt.subplot(1 , 1 , 1 ) X = x_tsne[:, 0 ] Y = x_tsne[:, 1 ] axl.set_xlim([min (X), max (X)]) axl.set_ylim([min (Y), max (Y)]) for ii in range (x_tsne.shape[0 ]): text = data.y.data.numpy()[ii] axl.text(X[ii], Y[ii], str (text), fontsize=5 , bbox=dict (boxstyle="round" , facecolor=plt.cm.Set1(text), alpha=0.7 )) axl.set_xlabel("TSNE Feature 1" , size=13 ) axl.set_xlabel("TSNE Feature 2" , size=13 ) axl.set_title(t, size=15 ) plt.show() x_tsne=TSNE(n_components=2 ).fit_transform(dataset.data.x.data.numpy()) draw(x_tsne,"Original feature TSNE" ) activation={} def get_activation (name ): def hook (model,input ,output ): activation[name]=output.detach() return hook model.conv1.register_forward_hook(get_activation("conv1" )) _=model(data) conv1=activation["conv1" ].data.cpu().numpy() print ("conv1.shape: " ,conv1.shape) x_tsne=TSNE(n_components=2 ).fit_transform(conv1) draw(x_tsne,"GCNConv Feature TSNE" ) X=dataset.data.x.data.numpy() Y=dataset.data.y.data.numpy() train_mask=dataset.data.train_mask.data.numpy() test_mask=dataset.data.test_mask.data.numpy() train_x=X[:140 ,:] train_y=Y[train_mask] test_x=X[1708 :2709 ,:] test_y=Y[test_mask] svmmodel=SVC() svmmodel.fit(train_x,train_y) prelab=svmmodel.predict(test_x) print ("SVM的预测精度:" ,accuracy_score(test_y,prelab)) train_y=Y.copy() train_y[test_mask==True ]=-1 test_y=Y[test_mask] lp_model=_label_propagation.LabelPropagation(kernel="knn" ,n_neighbors=3 ) lp_model.fit(X,train_y) prelab=lp_model.transduction_ print ("LP模型的预测精度:" ,accuracy_score(Y[test_mask],prelab[test_mask])) if __name__ == '__main__' : args=parse_args() main(dataset,args)
三、图卷积层底层原理
3.1 基于Python的Graph Convolution结构
我们对PyG的Graph Convolution结构进行分析:
1 2 3 4 5 def glorot (tensor ): if tensor is not None : stdv = math.sqrt(6.0 / (tensor.size(-2 ) + tensor.size(-1 ))) tensor.data.uniform_(-stdv, stdv)
1 2 3 4 def zeros (tensor ): if tensor is not None : tensor.data.fill_(0 )
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 105 class GCNConv (MessagePassing ): def __init__ (self, in_channels, out_channels, improved=False , cached=False , bias=True , normalize=True , **kwargs ): super (GCNConv, self).__init__(aggr='add' , **kwargs) self.in_channels = in_channels self.out_channels = out_channels self.improved = improved self.cached = cached self.normalize = normalize self.weight = Parameter(torch.Tensor(in_channels, out_channels)) if bias: self.bias = Parameter(torch.Tensor(out_channels)) else : self.register_parameter('bias' , None ) self.reset_parameters() def reset_parameters (self ): glorot(self.weight) zeros(self.bias) self.cached_result = None self.cached_num_edges = None @staticmethod def norm (edge_index, num_nodes, edge_weight=None , improved=False , dtype=None ): if edge_weight is None : edge_weight = torch.ones((edge_index.size(1 ), ), dtype=dtype, device=edge_index.device) fill_value = 1 if not improved else 2 edge_index, edge_weight = add_remaining_self_loops( edge_index, edge_weight, fill_value, num_nodes) row, col = edge_index deg = scatter_add(edge_weight, row, dim=0 , dim_size=num_nodes) deg_inv_sqrt = deg.pow (-0.5 ) deg_inv_sqrt[deg_inv_sqrt == float ('inf' )] = 0 return edge_index, deg_inv_sqrt[row] * edge_weight * deg_inv_sqrt[col] def forward (self, x, edge_index, edge_weight=None ): """""" x = torch.matmul(x, self.weight) if self.cached and self.cached_result is not None : if edge_index.size(1 ) != self.cached_num_edges: raise RuntimeError( 'Cached {} number of edges, but found {}. Please ' 'disable the caching behavior of this layer by removing ' 'the `cached=True` argument in its constructor.' .format ( self.cached_num_edges, edge_index.size(1 ))) if not self.cached or self.cached_result is None : self.cached_num_edges = edge_index.size(1 ) if self.normalize: edge_index, norm = self.norm(edge_index, x.size(self.node_dim), edge_weight, self.improved, x.dtype) else : norm = edge_weight self.cached_result = edge_index, norm edge_index, norm = self.cached_result return self.propagate(edge_index, x=x, norm=norm) def message (self, x_j, norm ): return norm.view(-1 , 1 ) * x_j if norm is not None else x_j def update (self, aggr_out ): if self.bias is not None : aggr_out = aggr_out + self.bias return aggr_out def __repr__ (self ): return '{}({}, {})' .format (self.__class__.__name__, self.in_channels, self.out_channels)
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 105 106 107 108 109 def propagate (self, edge_index: Adj, size: Size = None , **kwargs ): decomposed_layers = 1 if self._explain else self.decomposed_layers for hook in self._propagate_forward_pre_hooks.values(): res = hook(self, (edge_index, size, kwargs)) if res is not None : edge_index, size, kwargs = res size = self.__check_input__(edge_index, size) if (isinstance (edge_index, SparseTensor) and self.fuse and not self._explain): coll_dict = self.__collect__(self.__fused_user_args__, edge_index, size, kwargs) msg_aggr_kwargs = self.inspector.distribute( 'message_and_aggregate' , coll_dict) for hook in self._message_and_aggregate_forward_pre_hooks.values(): res = hook(self, (edge_index, msg_aggr_kwargs)) if res is not None : edge_index, msg_aggr_kwargs = res out = self.message_and_aggregate(edge_index, **msg_aggr_kwargs) for hook in self._message_and_aggregate_forward_hooks.values(): res = hook(self, (edge_index, msg_aggr_kwargs), out) if res is not None : out = res update_kwargs = self.inspector.distribute('update' , coll_dict) out = self.update(out, **update_kwargs) elif isinstance (edge_index, Tensor) or not self.fuse: if decomposed_layers > 1 : user_args = self.__user_args__ decomp_args = {a[:-2 ] for a in user_args if a[-2 :] == '_j' } decomp_kwargs = { a: kwargs[a].chunk(decomposed_layers, -1 ) for a in decomp_args } decomp_out = [] for i in range (decomposed_layers): if decomposed_layers > 1 : for arg in decomp_args: kwargs[arg] = decomp_kwargs[arg][i] coll_dict = self.__collect__(self.__user_args__, edge_index, size, kwargs) msg_kwargs = self.inspector.distribute('message' , coll_dict) for hook in self._message_forward_pre_hooks.values(): res = hook(self, (msg_kwargs, )) if res is not None : msg_kwargs = res[0 ] if isinstance (res, tuple ) else res out = self.message(**msg_kwargs) for hook in self._message_forward_hooks.values(): res = hook(self, (msg_kwargs, ), out) if res is not None : out = res if self._explain: edge_mask = self._edge_mask if self._apply_sigmoid: edge_mask = edge_mask.sigmoid() if out.size(self.node_dim) != edge_mask.size(0 ): edge_mask = edge_mask[self._loop_mask] loop = edge_mask.new_ones(size[0 ]) edge_mask = torch.cat([edge_mask, loop], dim=0 ) assert out.size(self.node_dim) == edge_mask.size(0 ) out = out * edge_mask.view([-1 ] + [1 ] * (out.dim() - 1 )) aggr_kwargs = self.inspector.distribute('aggregate' , coll_dict) for hook in self._aggregate_forward_pre_hooks.values(): res = hook(self, (aggr_kwargs, )) if res is not None : aggr_kwargs = res[0 ] if isinstance (res, tuple ) else res out = self.aggregate(out, **aggr_kwargs) for hook in self._aggregate_forward_hooks.values(): res = hook(self, (aggr_kwargs, ), out) if res is not None : out = res update_kwargs = self.inspector.distribute('update' , coll_dict) out = self.update(out, **update_kwargs) if decomposed_layers > 1 : decomp_out.append(out) if decomposed_layers > 1 : out = torch.cat(decomp_out, dim=-1 ) for hook in self._propagate_forward_hooks.values(): res = hook(self, (edge_index, size, kwargs), out) if res is not None : out = res return out
好复杂呀!我们只需要其中的精华就行啦,那么有一种比较简单的实现方式:
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 import torch.nn as nnimport torchclass GCNConv (nn.Module): def __init__ (self, input_dim, output_dim, use_bias=True ): super (GCNConv, self).__init__() self.input_dim = input_dim self.output_dim = output_dim self.use_bias = use_bias self.weight = nn.Parameter(torch.Tensor(input_dim, output_dim)) if self.use_bias: self.bias = nn.Parameter(torch.Tensor(output_dim)) else : self.register_parameter('bias' , None ) self.reset_parameters() def reset_parameters (self ): nn.init.kaiming_uniform_(self.weight) if self.use_bias: nn.init.zeros_(self.bias) def forward (self, adjacency, input_feature,l=1 ): ''' :param adjacency: 邻接矩阵 :param input_feature: 输入特征 :param l: lambda 影响自环权重值 :return: ''' size = adjacency.shape[0 ] support = torch.mm(input_feature, self.weight) A=adjacency+l*torch.eye(size) SUM=A.sum (dim=1 ) D=torch.diag_embed(SUM) D=D.__pow__(-0.5 ) D[D==float ("inf" )]=0 adjacency=torch.sparse.mm(D,adjacency) adjacency=torch.sparse.mm(adjacency,D) output = torch.sparse.mm(adjacency, support) if self.use_bias: output += self.bias return output def __repr__ (self ): return self.__class__.__name__ + ' (' + str (self.input_dim) + ' -> ' + str (self.output_dim) + ')'
四、参考
本文中的部分图片参考自博客:https://blog.csdn.net/qq_43787862/article/details/113830925,非常感谢大佬的翻译!