GCN学习笔记：第一部分，手把手用Numpy实现GCN

第一部分：图卷积网络到底怎么做，这是一份极简的Numpy实现

由于图结构非常复杂且信息量很大，因此对于图的机器学习是一项艰巨的任务。本文介绍了如何使用图卷积网络（GCN）对图进行深度学习，GCN 是一种可直接作用于图并利用其结构信息的强大神经网络。

本文将介绍 GCN，并使用代码示例说明信息是如何通过 GCN 的隐藏层传播的。读者将看到 GCN 如何聚合来自前一层的信息，以及这种机制如何生成图中节点的有用特征表征

何为图卷积网络？

GCN 是一类非常强大的用于图数据的神经网络架构。事实上，它非常强大，即使是随机初始化的两层 GCN 也可以生成图网络中节点的有用特征表征。下图展示了这种两层 GCN 生成的每个节点的二维表征。请注意，即使没有经过任何训练，这些二维表征也能够保存图中节点的相对邻近性。

更形式化地说，图卷积网络（GCN）是一个对图数据进行操作的神经网络。给定图 G = (V, E)，GCN 的输入为：

• 1，一个输入维度为 N × F⁰ 的特征矩阵 X，其中 N 是图网络中的节点数而 F⁰ 是每个节点的输入特征数。

• 2，一个图结构的维度为 N × N 的矩阵表征，例如图 G 的邻接矩阵 A。[1]

因此，GCN 中的隐藏层可以写作 Hⁱ = f(Hⁱ⁻¹, A))。其中，H⁰ = X，f 是一种传播规则 [1]。每一个隐藏层 Hⁱ 都对应一个维度为 N × Fⁱ 的特征矩阵，该矩阵中的每一行都是某个节点的特征表征。在每一层中，GCN 会使用传播规则 f 将这些信息聚合起来，从而形成下一层的特征。这样一来，在每个连续的层中特征就会变得越来越抽象。在该框架下，GCN 的各种变体只不过是在传播规则 f 的选择上有所不同 [1]。

传播规则的简单示例

下面，本文将给出一个最简单的传播规则示例 [1]：

f(Hⁱ, A) = σ(AHⁱWⁱ)

其中，Wⁱ 是第 i 层的权重矩阵，σ 是非线性激活函数（如 ReLU 函数）。权重矩阵的维度为 Fⁱ × Fⁱ⁺¹，即权重矩阵第二个维度的大小决定了下一层的特征数。如果你对卷积神经网络很熟悉，那么你会发现由于这些权重在图中的节点间共享，该操作与卷积核滤波操作类似。

简化

接下来我们在最简单的层次上研究传播规则。令：

• 1，i = 1，（约束条件 f 是作用于输入特征矩阵的函数）

• 2，σ 为恒等函数

• 3，选择权重（约束条件： AH⁰W⁰ =AXW⁰ = AX）

换言之，f(X, A) = AX。该传播规则可能过于简单，本文后面会补充缺失的部分。此外，AX 等价于多层感知机的输入层。

简单的图示例

我们将使用下面的图作为简单的示例：

                                          

                                                                                 一个简单的有向图

使用 numpy 编写的上述有向图的邻接矩阵表征如下：

A = np.matrix([
[0, 1, 0, 0],
[0, 0, 1, 1],
[0, 1, 0, 0],
[1, 0, 1, 0]],
dtype=float
)

接下来，需要抽取出特征！我们基于每个节点的索引为其生成两个整数特征，这简化了本文后面手动验证矩阵运算的过程。

In [3]: X = np.matrix([
[i, -i]
for i in range(A.shape[0])
], dtype=float)
X
Out[3]: matrix([
[ 0.,  0.],
[ 1., -1.],
[ 2., -2.],
[ 3., -3.]
])

应用传播规则

现在已经建立了一个图，其邻接矩阵为 A，输入特征的集合为 X。下面让我们来看看，当我们对其应用传播规则后会发生什么：

我们现在已经建立了一个图，其邻接矩阵为 A，输入特征的集合为 X。下面让我们来看看，当我们对其应用传播规则后会发生什么：

In [6]: A * X
Out[6]: matrix([
[ 1., -1.],
[ 5., -5.],
[ 1., -1.],
[ 2., -2.]]

每个节点的表征（每一行）现在是其相邻节点特征的和！换句话说，图卷积层将每个节点表示为其相邻节点的聚合。大家可以自己动手验证这个计算过程。请注意，在这种情况下，如果存在从 v 到 n 的边，则节点 n 是节点 v 的邻居。

问题

你可能已经发现了其中的问题：

• 1，节点的聚合表征不包含它自己的特征！该表征是相邻节点的特征聚合，因此只有具有自环（self-loop）的节点才会在该聚合中包含自己的特征 [1]。

• 2，度大的节点在其特征表征中将具有较大的值，度小的节点将具有较小的值。这可能会导致梯度消失或梯度爆炸 [1, 2]，也会影响随机梯度下降算法（随机梯度下降算法通常被用于训练这类网络，且对每个输入特征的规模（或值的范围）都很敏感）。

接下来，本文将分别对这些问题展开讨论。

增加自环

为了解决第一个问题，我们可以直接为每个节点添加一个自环 [1, 2]。具体而言，这可以通过在应用传播规则之前将邻接矩阵 A 与单位矩阵 I 相加来实现。

In [4]: I = np.matrix(np.eye(A.shape[0]))
I
Out[4]: matrix([
[1., 0., 0., 0.],
[0., 1., 0., 0.],
[0., 0., 1., 0.],
[0., 0., 0., 1.]
])
In [8]: A_hat = A + I
A_hat * X
Out[8]: matrix([
[ 1., -1.],
[ 6., -6.],
[ 3., -3.],
[ 5., -5.]])

现在，由于每个节点都是自己的邻居，每个节点在对相邻节点的特征求和过程中也会囊括自己的特征！

对特征表征进行归一化处理

通过将邻接矩阵 A 与度矩阵 D 的逆相乘，对其进行变换，从而通过节点的度对特征表征进行归一化。因此，我们简化后的传播规则如下：

f(X, A) = D⁻¹AX

让我们看看发生了什么。我们首先计算出节点的度矩阵。注意：此处计算节点的度是用节点的入度，也可以根据自身的任务特点用出度，在本文中，这个选择是任意的。一般来说，您应该考虑节点之间的关系是如何与您的具体任务相关。例如，您可以使用in-degree来进行度矩阵计算，前提是只有关于节点的in-neighbors的信息与预测其具体任务中的标签相关。相反，如果只有关于外部邻居的信息是相关的，则可以使用out-degree。最后，如果节点的out-和in-邻居都与您的预测相关，那么您可以基于in-和out-度的组合来计算度矩阵。

正如我将在下一篇文章中讨论的那样，您还可以通过其他方法对表示进行归一化，而不是使用逆矩阵乘法。

计算度矩阵：

In [9]: D = np.array(np.sum(A, axis=0))[0]
D = np.matrix(np.diag(D))
D
Out[9]: matrix([
[1., 0., 0., 0.],
[0., 2., 0., 0.],
[0., 0., 2., 0.],
[0., 0., 0., 1.]
])

在应用传播规则之前，不妨看看我们对邻接矩阵进行变换后发生了什么。

变换之前：

A = np.matrix([
[0, 1, 0, 0],
[0, 0, 1, 1],
[0, 1, 0, 0],
[1, 0, 1, 0]],
dtype=float
)

变换之后：

In [10]: D**-1 * A
Out[10]: matrix([
[0. , 1. , 0. , 0. ],
[0. , 0. , 0.5, 0.5],
[0. , 0.5, 0. , 0. ],
[1. , 0. , 1. , 0. ]
])

可以观察到，邻接矩阵中每一行的权重（值）都除以该行对应节点的度。我们接下来对变换后的邻接矩阵应用传播规则：

得到与相邻节点的特征均值对应的节点表征。这是因为（变换后）邻接矩阵的权重对应于相邻节点特征加权和的权重。大家可以自己动手验证这个结果。

整合

现在，我们将把自环和归一化技巧结合起来。此外，我们还将重新介绍之前为了简化讨论而省略的有关权重和激活函数的操作。

添加权重

首先要做的是应用权重。请注意，这里的 D_hat 是 A_hat = A + I 对应的度矩阵，即具有强制自环的矩阵 A 的度矩阵。

In [45]: W = np.matrix([
[1, -1],
[-1, 1]
])
D_hat**-1 * A_hat * X * W
Out[45]: matrix([
[ 1., -1.],
[ 4., -4.],
[ 2., -2.],
[ 5., -5.]
])

如果我们想要减小输出特征表征的维度，我们可以减小权重矩阵 W 的规模：

In [46]: W = np.matrix([
[1],
[-1]
])
D_hat**-1 * A_hat * X * W
Out[46]: matrix([[1.],
[4.],
[2.],
[5.]]
)

添加激活函数:

本文选择保持特征表征的维度，并应用 ReLU 激活函数。Relu函数的公式是

，代码为：

def relu(x):
return (abs(x) + x) / 2

一个带有邻接矩阵、输入特征、权重和激活函数的完整隐藏层如下：

In [51]: W = np.matrix([
[1, -1],
[-1, 1]
])
relu(D_hat**-1 * A_hat * X * W)
Out[51]: matrix([[1., 0.],
[4., 0.],
[2., 0.],
[5., 0.]])

---

在真实场景下的应用

Zachary 空手道俱乐部

Zachary 空手道俱乐部是一个被广泛使用的社交网络，其中的节点代表空手道俱乐部的成员，边代表成员之间的相互关系。当年，Zachary 在研究空手道俱乐部的时候，管理员和教员发生了冲突，导致俱乐部一分为二。下图显示了该网络的图表征，其中的节点标注是根据节点属于俱乐部的哪个部分而得到的，「0」表示属于Mr. Hi部分的中心节点，[32」表示属于Officer阵营的中心节点，参考https://networkx.github.io/documentation/stable/_modules/networkx/generators/social.html#karate_club_graph。

def plot_graph(G):
'''
G: a networkx G
'''
%matplotlib notebook
import matplotlib.pyplot as plt
plt.figure()
pos = nx.spring_layout(G)
edges = G.edges()
nodelist1 = []
nodelist2 = []
for i in range (34):
if zkc.nodes[i]['club'] == 'Mr. Hi':
nodelist1.append(i)
else:
nodelist2.append(i)
nx.draw_networkx(G, pos, edges=edges);
nx.draw_networkx_nodes(G, pos, nodelist=nodelist1, node_size=300, node_color='r',alpha = 0.8)
nx.draw_networkx_nodes(G, pos, nodelist=nodelist2, node_size=300, node_color='b',alpha = 0.8)
#     nx.draw_networkx_edges(G, pos, edgelist=edges,alpha =0.4)
plot_graph(zkc)

构建 GCN：

接下来，我们将构建一个图卷积网络。我们并不会真正训练该网络，但是会对其进行简单的随机初始化，从而生成我们在本文开头看到的特征表征。我们将使用 networkx，它有一个可以很容易实现的 Zachary 空手道俱乐部的图表征。然后，我们将计算 A_hat 和 D_hat 矩阵。

import networx as nx
from networkx import to_numpy_matrix
zkc = nx.karate_club_graph()
order = sorted(list(zkc.nodes()))
A = to_numpy_matrix(zkc, nodelist=order)
I = np.eye(zkc.number_of_nodes())
A_hat = A + I
D_hat = np.array(np.sum(A_hat, axis=0))[0]
D_hat = np.matrix(np.diag(D_hat))

接下来，我们将随机初始化权重。

W_1 = np.random.normal(
loc=0, scale=1, size=(zkc.number_of_nodes(), 4))
W_1

array([[-1.60850738, -1.81286031, -1.21488112, -0.67512344],
[ 0.4728306 ,  0.48590465,  0.46768218, -1.28452497],
[ 0.48404091,  0.36123352, -0.28535282,  1.3567935 ],
[ 0.5454595 ,  1.47414699, -0.2436099 ,  1.8448905 ],
[ 0.35601084,  0.17007017, -1.08012856,  0.12243646],
[ 1.08218031, -0.86061025,  0.2198768 , -0.88965361],
[ 0.77173117,  1.0753156 ,  0.80771322, -0.44305522],
[ 2.32292136, -1.11497939, -0.08832334, -0.48259347],
[-1.14989919, -1.07299688,  0.35833233, -1.02271864],
[-1.26604098,  1.11403405, -1.27368651,  0.68473466],
[ 0.03711537, -0.75404804,  2.41704567, -0.26348625],
[-0.82774307, -1.27096915,  0.01414155, -0.26087749],
[-0.47360135,  2.2111597 , -0.2785656 , -1.19886492],
[-0.39920476,  0.37976039,  0.37631679, -0.85434706],
[-1.06777751,  1.37249251, -1.5528223 , -0.7082742 ],
[ 1.21119404, -1.00086211,  0.79906034,  1.51296967],
[-1.93495147, -1.41661229, -1.46847639, -0.88607286],
[-0.57851013,  0.91926198, -0.4840015 ,  0.3120885 ],
[-0.1431892 ,  0.66630315,  0.37068425,  0.65605851],
[ 1.73521822, -1.78790492, -0.08257555,  1.30993567],
[-1.45908455,  0.12757605,  0.57191581,  0.02318818],
[ 0.7415613 ,  1.26621066, -0.51246359,  0.8104835 ],
[ 0.91050774,  0.05019478, -0.95328867, -0.17870357],
[ 1.208966  ,  0.04708878,  0.2219424 , -1.04159587],
[ 0.1054095 , -0.41581914, -1.09173746, -0.55462029],
[-1.21917691,  0.9778717 , -0.54145143, -2.0725276 ],
[-0.61375621,  0.98126815,  0.07500444, -0.76939277],
[-0.94490006,  1.98768759,  0.1917666 , -0.64337003],
[-1.47937171,  0.09036739,  0.57247029,  1.32377866],
[ 0.24979308,  0.65461977,  0.84397059, -2.47825603],
[-1.65513804,  0.49083851, -1.52508566, -1.12353088],
[-0.42206369, -0.54915933, -0.03701411, -1.96757786],
[ 0.75976155,  1.28285357,  0.44180894,  1.92307764],
[-1.0557052 , -0.7521497 , -1.97767775,  0.01836153]])

W_2 = np.random.normal(
loc=0, size=(W_1.shape[1], 2))
W_2

array([[-0.96457918,  0.45900237],
[-0.55060714, -1.08875363],
[-1.07701631, -0.71475231],
[ 0.45233694, -0.26491149]])

接着，我们会堆叠 GCN 层。这里，我们只使用单位矩阵作为特征表征，即每个节点被表示为一个 one-hot 编码的类别变量。

def gcn_layer(A_hat, D_hat, X, W):
return relu(D_hat**-1 * A_hat * X * W)

H_1 = gcn_layer(A_hat, D_hat, I, W_1)
H_2 = gcn_layer(A_hat, D_hat, H_1, W_2)
output = H_2
output

matrix([[0.17394349, 0.12613323],
[0.21514969, 0.1248979 ],
[0.2200705 , 0.07164891],
[0.20434201, 0.16982548],
[0.17972395, 0.26591497],
[0.13314916, 0.17843399],
[0.08887072, 0.11894365],
[0.24585268, 0.25294519],
[0.17443997, 0.14387933],
[0.27246705, 0.07764017],
[0.18066635, 0.25648615],
[0.2912415 , 0.46735709],
[0.23301799, 0.20337298],
[0.20719444, 0.21210805],
[0.10766761, 0.02108594],
[0.1048165 , 0.02436026],
[0.01142164, 0.0006152 ],
[0.23678846, 0.06878243],
[0.14423298, 0.        ],
[0.16002963, 0.10097575],
[0.13788295, 0.11606215],
[0.34037411, 0.27910508],
[0.14169454, 0.09565058],
[0.16444515, 0.        ],
[0.13790554, 0.        ],
[0.13951285, 0.        ],
[0.08406175, 0.04714129],
[0.18639611, 0.04326001],
[0.20328647, 0.26678944],
[0.11548443, 0.04857866],
[0.11449148, 0.        ],
[0.17931359, 0.1911088 ],
[0.14008594, 0.        ],
[0.14183489, 0.        ]])

经过多次随机生成W_1和W_2权重矩阵，得到上图H_2，但是我发现经过激活函数relu之后，x轴与y轴有很多零值，导致可视化效果很差，可视化效果如下图，初步分析，可能的原因是权重矩阵是随机生成的，没有用后面的具体任务去更新权重矩阵，，画图代码及图片如下：

import matplotlib.pyplot as plt
%matplotlib notebook
for i in range (34):
if zkc.nodes[i]['club'] == 'Mr. Hi':
plt.scatter(np.array(output)[i,0],np.array(output)[i,1] ,label=str(i),color = 'b',alpha=0.5,s = 250)
plt.text(np.array(output)[i,0],np.array(output)[i,1] ,i, horizontalalignment='center',verticalalignment='center', fontdict={'color':'black'})
# 为每个点添加标签，一些形如（x轴，y轴，标签）的元组，水平及垂直位置，背景颜色
else:
plt.scatter(np.array(output)[i,0],np.array(output)[i,1] ,label = 'i',color = 'r',alpha=0.5,s = 250)
plt.text(np.array(output)[i,0],np.array(output)[i,1] ,i, horizontalalignment='center',verticalalignment='center', fontdict={'color':'black'})
# plt.scatter(np.array(output)[:,0],np.array(output)[:,1],label = 0:33)

我尝试去掉激活函数relu，重新运行一遍，发现效果反而更好

def gcn_layer(A_hat, D_hat, X, W):
return D_hat**-1 * A_hat * X * W
H_1 = gcn_layer(A_hat, D_hat, I, W_1)
H_2 = gcn_layer(A_hat, D_hat, H_1, W_2)
output = H_2
output

matrix([[ 0.98839777,  1.13734122],
[ 0.51031539,  0.47345377],
[ 0.77387213,  0.49289261],
[ 0.61797444,  1.11746982],
[ 2.12295151,  2.43676755],
[ 1.56395358,  1.84093717],
[ 1.57723571,  1.94485012],
[ 0.62368931,  1.0618257 ],
[ 0.75132512,  0.99456766],
[ 0.55593252,  0.29953574],
[ 1.93535421,  2.47420126],
[ 0.70997674,  1.38977267],
[ 0.68421621,  1.82005162],
[ 0.53495206,  1.01818959],
[ 0.72965886, -0.50122342],
[ 0.48213244,  0.19867772],
[ 1.41885305,  1.59985406],
[ 0.27414389,  0.76435937],
[ 0.60733412, -0.38820873],
[ 0.27476357,  0.48294864],
[ 0.51044012,  0.21351918],
[ 0.45876093,  0.63466151],
[ 0.64684601,  0.00845191],
[-0.20156276, -1.09758021],
[-0.17949636, -1.40743413],
[-0.41499531, -1.76745965],
[-0.48696632, -0.69577392],
[ 0.06628659, -0.50773281],
[ 0.91267734,  0.53837441],
[-0.30924486, -0.83232449],
[ 0.59774103,  0.78497008],
[ 0.29614438, -0.24454598],
[ 0.72875561, -0.55889401],
[ 0.60941032, -0.69855984]])

feature_representations = {
node: np.array(output)[node]
for node in zkc.nodes()}
feature_representations

{0: array([0.98839777, 1.13734122]),
1: array([0.51031539, 0.47345377]),
2: array([0.77387213, 0.49289261]),
3: array([0.61797444, 1.11746982]),
4: array([2.12295151, 2.43676755]),
5: array([1.56395358, 1.84093717]),
6: array([1.57723571, 1.94485012]),
7: array([0.62368931, 1.0618257 ]),
8: array([0.75132512, 0.99456766]),
9: array([0.55593252, 0.29953574]),
10: array([1.93535421, 2.47420126]),
11: array([0.70997674, 1.38977267]),
12: array([0.68421621, 1.82005162]),
13: array([0.53495206, 1.01818959]),
14: array([ 0.72965886, -0.50122342]),
15: array([0.48213244, 0.19867772]),
16: array([1.41885305, 1.59985406]),
17: array([0.27414389, 0.76435937]),
18: array([ 0.60733412, -0.38820873]),
19: array([0.27476357, 0.48294864]),
20: array([0.51044012, 0.21351918]),
21: array([0.45876093, 0.63466151]),
22: array([0.64684601, 0.00845191]),
23: array([-0.20156276, -1.09758021]),
24: array([-0.17949636, -1.40743413]),
25: array([-0.41499531, -1.76745965]),
26: array([-0.48696632, -0.69577392]),
27: array([ 0.06628659, -0.50773281]),
28: array([0.91267734, 0.53837441]),
29: array([-0.30924486, -0.83232449]),
30: array([0.59774103, 0.78497008]),
31: array([ 0.29614438, -0.24454598]),
32: array([ 0.72875561, -0.55889401]),
33: array([ 0.60941032, -0.69855984])}

import matplotlib.pyplot as plt
%matplotlib notebook
for i in range (34):
if zkc.nodes[i]['club'] == 'Mr. Hi':
plt.scatter(np.array(output)[i,0],np.array(output)[i,1] ,label=str(i),color = 'b',alpha=0.5,s = 250)
plt.text(np.array(output)[i,0],np.array(output)[i,1] ,i, horizontalalignment='center',verticalalignment='center', fontdict={'color':'black'})
# 为每个点添加标签，一些形如（x轴，y轴，标签）的元组，水平及垂直位置，背景颜色
else:
plt.scatter(np.array(output)[i,0],np.array(output)[i,1] ,label = 'i',color = 'r',alpha=0.5,s = 250)
plt.text(np.array(output)[i,0],np.array(output)[i,1] ,i, horizontalalignment='center',verticalalignment='center', fontdict={'color':'black'})
# plt.scatter(np.array(output)[:,0],np.array(output)[:,1],label = 0:33)

                                                                     Zachary 空手道俱乐部图网络中节点的特征表征。

你看，这样的特征表征可以很好地将 Zachary 空手道俱乐部的两个社区划分开来。至此，我们甚至都没有开始训练模型！我们应该注意到，在该示例中由于 ReLU 函数的作用，在 x 轴或 y 轴上随机初始化的权重很可能为 0。

结语

本文中对图卷积网络进行了高水平的的介绍，并说明了 GCN 中每一层节点的特征表征是如何基于其相邻节点的聚合构建的。读者可以从中了解到如何使用 numpy 构建这些网络，以及它们的强大：即使是随机初始化的 GCN 也可以将 Zachary 空手道俱乐部网络中的社区分离开来。 

在下一篇文章中，我将更详细地介绍技术细节，并展示如何使用半监督学习实现和训练最近发布的GCN。你可以在本人csdn找到下一篇文章。

参考文献

[1] Blog post on graph convolutional networks by Thomas Kipf.

[2] Paper called Semi-Supervised Classification with Graph Convolutional Networks by Thomas Kipf and Max Welling.

[3]https://blog.csdn.net/qq_36793545/article/details/84844867

[4]https://towardsdatascience.com/how-to-do-deep-learning-on-graphs-with-graph-convolutional-networks-7d2250723780