[人工智能]图卷积神经网络(GCN)综述与实现（PyTorch版)

图卷积神经网络(GCN)综述与实现（PyTorch版)

本文的实验环境为PyTorch = 1.11.0 + cu113，PyG = 2.0.4，相关依赖库和数据集的下载请见链接。

一、图卷积神经网络介绍

1.1 传统图像卷积

卷积神经网络中的卷积(Convolution)指的是在图像上进行的输入和卷积核之间离散内积运算，其本质上就是利用共享参数的滤波器，通过计算中心值以及相邻节点的值进行加权获得带有局部空间特征的特征提取器。

其具有三个重要的特征，分别为：

• 稀疏连接
  • 相较于全连接层，卷积层输入和输出间的连接是稀疏的，能够大大减少参数的数量，加快网络的训练速度。
• 参数共享
  • 卷积核的权重参数可以被多个函数或操作共享，这样只需要训练一个参数集，而不需要对每个位置都训练一个参数集。此外，由于卷积核的大小一般是小于输入大小，也能起到减少参数数量的作用
• 等变表示
  • 事实上，每个卷积层可以通过多个卷积核来进行特征提取，并且在卷积运算后，卷积神经网络对输入的图像具有平移不变性(有严格的数学论证)
    

一般在卷积层后，会通过一个池化层进行降维，进一步降低网络的复杂度和非线性程度。在那之后，可以将通过卷积池化层后的特征输入全连接层或反卷积层进行进一步的分类、分割、重建工作。当然，传统的卷积操作一般适用于结构化数据。

1.2 图结构

图作为一种典型的非结构化非线性数据(非欧几里得数据)，其可以表示一对一、一对多、多对多的关系，因而常被用于描述复杂的数据对象，譬如社交网络、知识图谱、城市路网、3D点云等。与结构化数据不同，图的局部输入维度可变，即每个节点的邻居节点数量不同；图具有无序性，即节点间并不存在先后关系，仅存在连接关系(点云是置换不变性，在无序性的基础上，交换两点或多点不会影响整体结果)。由于图结构的特殊性，传统CNN和RNN对其的表征能力并不理想。

对于图结构，我们可以将其抽象表示为：
$$G=(V,E)$$
在这里$V$表示图中节点的集合，而$E$为边的集合。对于图特征，我们一般有三个重要矩阵进行表示。

• 邻接矩阵$A$：adjacency matrix用来表示节点间的连接关系。
  

$$\left|\begin{array}{cccccc}& 0& 1& 0& 0& \\ & 0& 0& 1& 0& \\ & 1& 1& 0& 0& \\ & 0& 1& 1& 0& \end{array}\right|$$
? 对于带权的图，邻接矩阵将把1替换为对应的权重

• 度矩阵$D$：degree matrix用来表示节点的连接数，可以表征某个节点在图中的重要程度，是一个对角矩阵，例如针对上图的入度矩阵：
  $$\left|\begin{array}{cccccc}& 1& 0& 0& 0& \\ & 0& 3& 0& 0& \\ & 0& 0& 2& 0& \\ & 0& 0& 0& 0& \end{array}\right|$$

• 特征矩阵$X$：feature matrix用来表示节点的特征

1.3 图卷积神经网络

目前主流的图卷积基本上可以分为两类，一种是基于谱的图卷积，一种是基于空域的图卷积。

基于谱的图卷积通过傅里叶变换(FFT干的一件事情就是连接空域和频域)将节点映射到频域空间，通过频域空间上的卷积来实现时域上的卷积，最后将特征映射回空域。而基于空域的图卷积则是直接基于节点与邻居进行卷积提取特征，没有做域上的变换。

图卷积算子可表示为：
$$h_i^{l+1}=\sigma (\sum _{j\in N_i}\frac{1}{C_{ij}}{h}_j^l{w}_{R_j}^l)$$
其中，设中心结点为$i$；$h_i^l$为结点$i$在$l$层的特征表达；$\sigma$是非线性激活函数；$C_{ij}$则是归一化因子，譬如结点度的倒数、反距离权重、高斯衰减权重等；$N_i$是结点$i$的邻接节点(包括自身)；$R_i$表示节点$i$的类型；$W_{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\times 3$网络结构，有部分是带有标识的，而有部分则是无标识的。

GCN通过考虑节点本身以及邻居节点的特征信息来提取潜在的关系。比如我们在中心节点拼接了邻居的特征后，使用平均池化的方式对这些特征进行聚合，再通过浅层网络进行学习训练得到新的数据。

每个GCN层做的事情：

• 获取当前节点和邻接节点特征
• 通过聚合函数获取局部特征(带有拓扑关系)
• 浅层学习训练，获取高维特征
  

数学推导

对于一个如上图所示的无向图，我们要怎么样才能获取到某个节点以及其邻居节点的特征呢？一种非常直观的想法💡是，在当前节点的特征之后，根据权重拼接与之相邻节点的特征。此时，空域上的距离往往成为了权重的影响因素。

那么，邻接表就成了我们考虑节点间拓扑关系最重要的结构。我们写成$A\times X\times W$，$A$是邻接矩阵，$X$是特征矩阵，$W$是权重。右乘相当于控制行，左乘相当于控制列，$A$左乘$X$，相当于在对应节点处，使用哪些节点特征构建新的特征：
$$\left|\begin{array}{cccccc}& 0& 1& 0& 0& \\ & 0& 0& 1& 1& \\ & 0& 1& 0& 1& \\ & 1& 1& 0& 0& \end{array}\right|\times \left|\begin{array}{cccccc}& 1& 1& 1& 1& \\ & 2& 2& 2& 2& \\ & 3& 3& 3& 3& \\ & 40& 40& 40& 40& \end{array}\right|=\left|\begin{array}{cccccc}& 2& 2& 2& 2& \\ & 43& 43& 43& 43& \\ & 42& 42& 42& 42& \\ & 3& 3& 3& 3& \end{array}\right|$$
诶，糟糕，这矩阵数数数着数着咋把自己给忘了呀！而且似乎…有些小朋友发育的过于良好，这样会导致梯度消化不良的！(梯度爆炸或消失问题)

为了处理这种情况，我们引入了新的邻接矩阵$\tilde{A}$用于承载当前节点信息~
$$\tilde{A}=A+\lambda I_N$$
当权重因子$\lambda$取$1$时，$\tilde{A}=A+I_N$，此时意味着节点中心和邻居一样重要啦。🎉

值得注意的是，$\lambda$本身也是一个可以由训练得到的参数。

针对梯度问题呢，我们可以利用对角矩阵实现标准化，譬如利用度矩阵的逆${\tilde{D}}^{?1}$进行平均化的操作：
$$\left|\begin{array}{cccccc}& 2& 0& 0& 0& \\ & 0& 3& 0& 0& \\ & 0& 0& 4& 0& \\ & 0& 0& 0& 1& \end{array}\right|\stackrel{inverse}{}\left|\begin{array}{cccccc}& 1/2& 0& 0& 0& \\ & 0& 1/3& 0& 0& \\ & 0& 0& 1/4& 0& \\ & 0& 0& 0& 1& \end{array}\right|$$
于是乎，新的公式诞生啦：
$${\tilde{D}}^{?1}(\tilde{A}X)W$$
但我们在这里仅仅对行进行了标准化，在列方向上并没有被标准化，所以，修正后的公式为：
$$\begin{gathered} ({\tilde{D}}^{?1}\tilde{A}{\tilde{D}}^{?1})XW \\ Let‘scallA=({\tilde{D}}^{?1}\tilde{A}{\tilde{D}}^{?1}) \end{gathered}$$
当然了，行做了一次标准化，列做了一次标准化，那么对于元素$X_{ij}$来说，那就是做了两次标准化诶！我们再开个方吧:happy:
$$\begin{gathered} ({\tilde{D}}^{?1/2}\tilde{A}{\tilde{D}}^{?1/2})XW \\ Let‘scallA=({\tilde{D}}^{?1/2}\tilde{A}{\tilde{D}}^{?1/2}) \end{gathered}$$
于是乎，最终可以得到以下的标准格式：
$$Z=f(X,A)=softmax(AReLU(AX{W}^{(0)})W^{(1)})$$
$W$是可训练的权重，$A$是标准化后的邻接矩阵。

此时我们在回过头来看这个图卷积算子，是不是有些恍然大悟了呢🤔
$$h_i^{l+1}=\sigma (\sum _{j\in N_i}\frac{1}{C_{ij}}{h}_j^l{w}_{R_j}^l)$$

Thomas经过实验比对，浅层的GCN(一般为2层)取得的效果往往要比没有加入残差块的深层GCN来的好。GCN的网络层数代表着节点特征所能到达的最远距离，或者说节点信息的传播距离。深层的网络可能会令节点信息传播到整个网络，反而效果不那么好。

举个简单的栗子~🐈

假设有邻接矩阵$A$：
$$A=\left|\begin{array}{cccccc}& 0& 1& 0& 0& \\ & 1& 0& 1& 1& \\ & 0& 1& 0& 1& \\ & 0& 1& 1& 1& \end{array}\right|$$
对于节点$A_1$来说，跟他接轨的只有节点$A_2$，那么在第一层计算后，节点$A_1$聚合了来自节点$A_2$的特征。而对于节点$A_2$，与之接轨的则有节点$A_1,A_3,A_4$，经过第一层后，该节点聚合了所有节点的信息。那么在第二层，由于图结构不是动态更新的，$A_1$又会聚合$A_2$的特征信息，但此时的$A_2$已经聚合了$A_3,A_4$的特征信息！换言之，$A_1$感受到了来自远方的召唤！($A_3,A_4$的信息经过两层网络后传递到了$A_1$处，传递过程类似于并查集啦)

最后对于这些分类任务，采用交叉熵函数作为损失函数，计算获取最可能的情况即可。

二、基于PyG的图卷积神经网络半监督分类

2.1 模型准备

模块导入阶段，在本实验中，networkx主要用于图结构可视化，argparse用于参数的读写

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import torch
import torch.nn.functional as F
from torch_geometric.nn import GCNConv
from torch_geometric.datasets import Planetoid
from torch_geometric.utils import to_networkx
import networkx as nx
from sklearn.metrics import accuracy_score
from sklearn.manifold import TSNE
from sklearn.svm import SVC
from sklearn.semi_supervised import _label_propagation
import 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，用于后面的归纳学习设置。

通过以下语句下载或读取，若当前路径下没有找到文件，则会自动下载。

# 导入Cora数据集
dataset=Planetoid(root=r"./Cora",name="Cora") # root: 指定路径 name: 数据集名称
# 查看数据的基本情况
print("网络数据包含的类数量:",dataset.num_classes)
print("网络数据边的特征数量:",dataset.num_edge_features)
print("网络数据边的数量:",dataset.data.edge_index.shape[1]/2) # 除以2是OOC的组织形式
print("网络数据节点的特征数量:",dataset.num_node_features)
print("网络数据节点的数量:",dataset.data.x.shape[0])

输出的结果为：

网络数据包含的类数量: 7
网络数据边的特征数量: 0
网络数据边的数量: 5278.0
网络数据节点的特征数量: 1433
网络数据节点的数量: 2708
网络的边的数量:  5278
网络的节点的数量:  2708

我们可以看到，Cora数据集囊括了七个类共2708个节点，每个节点都带有1433个特征，一共有5278条边，边没有包含特征。

重点来看dataset.data：

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库中的有向图数据，方便进一步的分析和可视化。

# 通过networkx库进行可视化
CoraNet=to_networkx(dataset.data) # type : networkx.classes.graph.Graph
CoraNet=CoraNet.to_undirected() # 转化为无向图

# 查看网络情况
print("网络的边的数量: ",CoraNet.number_of_edges())
print("网络的节点的数量: ",CoraNet.number_of_nodes())
Node_class=dataset.data.y.data.numpy()
print("节点分类:",Node_class)

结果为：

网络的边的数量:  5278
网络的节点的数量:  2708
节点分类: [3 4 4 ... 3 3 3]

此时数据已经由Tensor格式转化为Graph格式啦，并且我们将标签y单独提取出来，方便接下来的操作。

在Graph类中，我们可以使用Graph.degree方法计算度，度越大的节点，说明跟其他节点的连通性越好，在整个网络中具有更大的贡献度，即越重要。下面我们对前三十的节点进行可视化：

# 计算节点的度
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)用于绘制图结构的边

# nx图结构可视化
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()

结果如下图所示：

我们可以再对训练集单独绘制一下：

# 可视化训练集的节点分布
nodelabel=np.arange(0,140) # 训练集节点位置
Node_class=dataset.data.y.data.numpy()[0:140]
# nx可视化
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进行搭建：

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
            # 这里我们用了两个图卷积层
            # 1433->32
            self.conv1=GCNConv(input_feature,32)
            # 32->num_classes
            self.conv2=GCNConv(32,num_classes)

        def forward(self,data):
            x,edge_index=data.x,data.edge_index
            # 需要输入的是节点的数据(index,feature)
            # 以及边索引COO，用于构建邻接矩阵和入度矩阵
            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库进行操作。

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()

我们可以定义一个保存路径，用于存放最优的模型参数和最优的优化器参数，譬如：

# 环境
    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

那首先进行网络的初始化

 # 初始化网络
    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

接着咱定义优化器

# 优化器
    if args.optimizer == 'Adam':
        optimizer = torch.optim.Adam(
            model.parameters(),
            lr=args.learning_rate,
            # betas=(0.9, 0.999),
            # eps=1e-08,
            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)

此时的输出为：

'''
GCNNet(
  (conv1): GCNConv(1433, 32)
  (conv2): GCNConv(32, 7)
)
'''

接下来，我们再进行网络的训练：

 # 网络训练
    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))

可以发现，到200个epoch后，模型已经收敛。当然有很大程度上是因为我们并没有随机取出数据集(用的是固定顺序)，最终的精度结果为：

Accuracy: 0.8060

为了直观体现网络的特征提取能力，我们对隐藏层获得的32维特征在空间中的分布情况进行可视化,与原始数据1433维进行比较，统一使用TSNE算法降到二维。首先我们先定义一个用于TSNE绘制的函数：

 # 对原始数据分布进行展示
    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()

我们来绘制原始数据：

x_tsne=TSNE(n_components=2).fit_transform(dataset.data.x.data.numpy())
    draw(x_tsne,"Original feature TSNE")

接着呢，利用hook函数获取中间层的输入，板子如下：

 # 对中间层32维数据进行可视化
    activation={} # 保存不同层的输出
    def get_activation(name):
        def hook(model,input,output):
            activation[name]=output.detach()
        return hook
    
    model.conv1.register_forward_hook(get_activation("conv1"))

# 绘制32维特征
    _=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)算法:

 # 使用SVM进行计算
    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))

最终的结果为：

SVM的预测精度： 0.56

 # 使用LabelPropagation模型进行计算
 # 对于不是有监督的训练数据的样本标签使用-1表示
    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]))

结果为：

LP模型的预测精度： 0.45

GCN在精度上远高于上面两个模型，具体原因在前面已经介绍过了。

2.7 完整代码

# 完整代码
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import torch
import torch.nn.functional as F
from torch_geometric.nn import GCNConv
from torch_geometric.datasets import Planetoid
from torch_geometric.utils import to_networkx
import networkx as nx
from sklearn.metrics import accuracy_score
from sklearn.manifold import TSNE
from sklearn.svm import SVC
from sklearn.semi_supervised import _label_propagation
import argparse
import os

dataset=Planetoid(root=r"./Cora",name="Cora") # root: 指定路径 name: 数据集名称

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
            # 这里我们用了两个图卷积层
            # 1433->32
            self.conv1=GCNConv(input_feature,32)
            # 32->num_classes
            self.conv2=GCNConv(32,num_classes)

        def forward(self,data):
            x,edge_index=data.x,data.edge_index
            # 需要输入的是节点的数据(index,feature)
            # 以及边索引COO，用于构建邻接矩阵和入度矩阵
            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,
            # betas=(0.9, 0.999),
            # eps=1e-08,
            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))

    # 为了直观体现网络的特征提取能力，我们对隐藏层获得的32维特征在空间中的分布情况进行可视化,
    # 与原始数据1433维进行比较，统一使用TSNE算法降维

    # 对原始数据分布进行展示
    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")

    # 对中间层32维数据进行可视化
    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")

    # 使用SVM进行计算
    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))

    # 使用LabelPropagation模型进行计算

    # 对于不是有监督的训练数据的样本标签使用-1表示
    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结构进行分析：

def glorot(tensor):
    # 参数标准化(方差)
    if tensor is not None:
        stdv = math.sqrt(6.0 / (tensor.size(-2) + tensor.size(-1)))
        tensor.data.uniform_(-stdv, stdv)

def zeros(tensor):
    # 归零
    if tensor is not None:
        tensor.data.fill_(0)

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 # 影响自环权重，也就是hat A = A+\lambda I 中的\lambda
        self.cached = cached # 缓存计算值
        self.normalize = normalize # 是否添加自环并对称标准化

        self.weight = Parameter(torch.Tensor(in_channels, out_channels)) # 权重参数，size为(in,out)

        if bias:
            self.bias = Parameter(torch.Tensor(out_channels)) # 偏置参数,size为(out)
        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)
		# 这玩意就是lambda , 影响自环权重
        fill_value = 1 if not improved else 2
        
        # 加入自环，做的就是：A' = A + \lambda I
        edge_index, edge_weight = add_remaining_self_loops(
            edge_index, edge_weight, fill_value, num_nodes)

        row, col = edge_index
        # 构建对角矩阵(degree) D
        deg = scatter_add(edge_weight, row, dim=0, dim_size=num_nodes)
        # D=(D^-0.5)
        deg_inv_sqrt = deg.pow(-0.5)
        # 把某些离谱的家伙修正为0
        # 会出现inf也是因为浮点的计算可能出现问题
        deg_inv_sqrt[deg_inv_sqrt == float('inf')] = 0
		# return : COO边,(D^-0.5) W (D^-0.5)
        return edge_index, deg_inv_sqrt[row] * edge_weight * deg_inv_sqrt[col]

    def forward(self, x, edge_index, edge_weight=None):
        """"""
        
        # 公式： (D^-0.5) A' (D^-0.5) X W
        # A' = A +\lambda I
        
        # 这里做的就是 X W
        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:
                # 是否进行标准化(利用D矩阵，相当于权重)
                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
		
        # Message.propagate 主要是利用邻接矩阵完成一些信息的堆叠
        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)

    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)

        # Run "fused" message and aggregation (if applicable).
        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)

        # Otherwise, run both functions in separation.
        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

                # For `GNNExplainer`, we require a separate message and
                # aggregate procedure since this allows us to inject the
                # `edge_mask` into the message passing computation scheme.
                if self._explain:
                    edge_mask = self._edge_mask
                    if self._apply_sigmoid:
                        edge_mask = edge_mask.sigmoid()
                    # Some ops add self-loops to `edge_index`. We need to do
                    # the same for `edge_mask` (but do not train those).
                    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

好复杂呀！我们只需要其中的精华就行啦，那么有一种比较简单的实现方式：

import torch.nn as nn
import torch
class 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:
            # 偏置先全给0
            nn.init.zeros_(self.bias)

    def forward(self, adjacency, input_feature,l=1):
        '''

        :param adjacency: 邻接矩阵
        :param input_feature: 输入特征
        :param l: lambda 影响自环权重值
        :return:
        '''
        # 公式: (D^-0.5) A' (D^-0.5) X W
        size = adjacency.shape[0]
        # X W
        support = torch.mm(input_feature, self.weight)

        # A' = A + \lambda I
        A=adjacency+l*torch.eye(size)

        # D: degree
        SUM=A.sum(dim=1)
        D=torch.diag_embed(SUM)
        # D'=D^(-0.5)
        D=D.__pow__(-0.5)
        # 让inf值变成0
        D[D==float("inf")]=0

        # (D^-0.5) A' (D^-0.5)
        adjacency=torch.sparse.mm(D,adjacency)
        adjacency=torch.sparse.mm(adjacency,D)
        
        # (D^-0.5) A' (D^-0.5) X W
        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，非常感谢大佬的翻译！