LINE 图嵌入算法介绍与源码浅析

LINE 图嵌入算法介绍与源码浅析

前言 (可以略过~)

最近阅读了一些 Paper 和代码实现, 打算及时记录下来, 以便日后查阅和回忆; 写关于 Paper 的博客, 可以对文章进行翻译, 也可以只描述中心思想+自己的感悟, 两者均有好处: 前者内容详尽, 事无巨细, 查阅起来方便快捷, 不足之处是我不太擅长, 没有动力做~ 🤣; 后者方便了写作者, 但对阅读者不太友好, 如果时间一长, 其实对自己也不太友好, 因为文章的主要内容估计忘得差不多了, 又得返回去看原文 … 联想到我写的 PyGCN 源码阅读, 感觉还是应该多介绍一些理论内容, 再配合源码进行分析, 事半功倍. 后面写作的时候, 尽可能使用理论+代码的形式吧, 如何恰到好处, 得多写写文章来把握.

广而告之

可以在微信中搜索 “珍妮的算法之路” 或者 “world4458” 关注我的微信公众号；另外可以看看知乎专栏 PoorMemory-机器学习, 以后文章也会发在知乎专栏中；

LINE 图嵌入算法介绍

文章信息

• 论文标题: LINE: Large-scale Information Network Embedding
• 论文地址: https://arxiv.org/abs/1503.03578
• 代码地址: https://github.com/tangjianpku/LINE, 另外作者所在的组还发布了一个图框架 GraphVite, 支持大规模的 embedding learning, 内置了多款图算法.
• 发表时间: 2015
• 论文作者: Jian Tang, Meng Qu , Mingzhe Wang, Ming Zhang, Jun Yan, Qiaozhu Mei
• 作者单位: Microsoft Research Asia + Peking University + University of Michigan

主要内容

文章介绍了 LINE 算法, 引入了两个核心概念: first-order proximity (一阶邻近) 和 second-order proximity (二阶邻近). 其中一阶邻近描述的是直接相连的节点之间的关系, 而二阶邻近则指的是节点之间虽然不直接相连, 但是它们会有一些共同的邻居. 在具体实现上, 为了描述二阶邻近, 作者用了两个 Embedding Layer, 其中一个用来描述节点本身, 另一个用来描述节点作为其他节点的 context 时的情况. Loss 在设计时借鉴了 Word2Vec 中的 Negative Sampling. 此外, 通过 AliasMethod 来对边进行采样, AliasMethod 采样方法的特点是它的时间复杂度为 $O(1)$, 权重越大的边有更大的几率被采样到.

下图描述了一阶邻近和二阶邻近:

其中节点 6 和节点 7 是直接相连的 (一阶邻近), 它们在低维空间的距离应该是接近的; 另外节点 6 和节点 5 虽然不直接相连, 然后它们有共同的邻居, 因此它们的低维特征表达也应该是接近的.

为了形式化的表达一阶和二阶邻近, 作者分别它们俩进行建模.

一阶邻近

其中一阶邻近表述如下:

对于每条无向边 $(i,j)$, 定义顶点 (vertex) $v_i$ 和 $v_j$ 的联合概率为:

$$p_1\left(v_i,v_j\right)=\frac{1}{1+\exp \left(-{\vec{u}}_i^T\cdot {\vec{u}}_j\right)}$$

其中 ${\vec{u}}_i\in {\mathbb{R}}^d$ 为顶点 $v_i$ 的低维特征表达. 上式定义了顶点间的联合概率, 而顶点间的经验概率定义为

$${\hat{p}}_1\left(i,j\right)=\frac{w_{ij}}{W},W=\sum _{(i,j)\in E}{w}_{ij}$$

其中 $E$ 为 Graph 中边的集合, $w_{ij}$ 为顶点 $v_i$ 和 $v_j$ 之间边的权重. 为了保留顶点间的一阶邻近关系, 文中采用计算两个联合概率分布和经验概率分布之间的距离, 而分布之间的距离一般采用 KL 散度 (divergence) 来进行计算. KL 散度定义为:

$$\sum _i{p}_i\left(\log p_i-\log q_i\right)$$

将 $p_1\left(v_i,v_j\right)$ 和 ${\hat{p}}_1\left(i,j\right)$ 代入上式, 并省去一些常数得到一阶邻近的优化目标为:

$$O_1=-\sum _{(i,j)\in E}{w}_{ij}\log p_1\left(v_i,v_j\right)$$

注意一阶邻近只适用于无向图, 而不能用于有向图. 具体原因见下面我的想法.

插入我的一点想法:

• 首先, 使用 Sigmoid 定义两个节点之间的联合概率, 输入是两个节点 embedding 的内积, 注意到 Sigmoid 是递增函数, 那么如果两个向量的内积越大, 即两个节点比较相似, 那么 p 的概率就越大. 另外作者后面定义了经验概率 $\hat{p}$, 如果两个节点之间的边权越大, 表示两个节点的联系越紧密.
• 对于一阶邻近只适用于无向图的问题, 是因为在计算 $p_1\left(v_i,v_j\right)$ 时的公式是满足对称性的, 即交换 $u_i$ 和 $u_j$ 的位置, 结果是不变的; 但是对于 ${\hat{p}}_1\left(i,j\right)$ 的结果, 对有向图来说结果会发生变化, 因为 $w_{ij}$ 不一定等于 $w_{ji}$; 这样联合概率的分布和经验概率的分布是有差异的, 再去计算它们之间的距离可能学习不到正确的表达.
• 另外需要注意的是, 使用 Sigmoid 计算出来的值不算分布, 毕竟经验分布 ${\hat{p}}_1\left(i,j\right)$ 是可以归一化的, 而 $p_1\left(v_i,v_j\right)$ 的结果不行.

二阶邻近

二阶邻近可以作用于有向/无向/加权/无权图. 二阶邻近认为, 如果节点 A 和节点 B 共享部分邻居, 那么它们之间也是相似的. 在这种情况下, 每个节点有两种角色, 其中一个即节点自身, 而另一个角色是作为其他节点的 “context”, 即作为其他节点的上下文, 如果节点 A 和 B 的上下文分布相似, 那么 A 和 B 也应该是相似的.

因此作者对节点 $v_i$ 引入了两个特征表达, 使用 ${\vec{u}}_i$ 表示它作为节点时的特征表达, 另外使用 ${\vec{u}}_i^{\prime }$ 表示它作为上下文 context 时的特征表达. 对于有向边 $(i,j)$, 节点 $v_j$ 是作为节点 $u_i$ 的上下文的, 定义条件概率:

$$p_2\left(v_j\mid v_i\right)=\frac{\exp \left({\vec{u}}_j^{\prime T}\cdot {\vec{u}}_i\right)}{{\sum }_{k=1}^{|V|}\exp \left({\vec{u}}_k^{\prime T}\cdot {\vec{u}}_i\right)}$$

其中 $|V|$ 为上下文节点的数量. 对每一个源节点 $v_i$ 来说, 上式定义了它在上下文节点上的条件概率分布 $p_2\left(\cdot \mid v_i\right)$. 此外, 经验概率分布定义为:

$${\hat{p}}_2\left(v_j\mid v_i\right)=\frac{w_{ij}}{d_i},d_i=\sum _{k\in \mathcal{N}(i)}{w}_{ik}$$

其中 $\mathcal{N}(i)$ 表示节点 $v_i$ 指向的邻居集合, $d_i$ 表示该节点的出度 (out-degree). 为了刻画二阶邻近, 计算条件概率分布和经验概率分布之间的 KL 散度, 去掉一些常数项后得到优化目标为:

$$O_2=-\sum _{(i,j)\in E}{w}_{ij}\log p_2\left(v_j\mid v_i\right)$$

通过对 { u ⃗ i } i = 1 , … , ∣ V ∣ \left\{\vec{u}_{i}\right\}_{i = 1, \ldots, |V|} {u i​}i=1,…,∣V∣​ 和 { u ⃗ i ′ } i = 1 , … , ∣ V ∣ \left\{\vec{u}^\prime_{i}\right\}_{i = 1, \ldots, |V|} {u i′​}i=1,…,∣V∣​ 的学习, 最小化目标函数, 这样就可以得到每个节点 $v_i$ 的表达 ${\vec{u}}_i$.

模型优化以及边采样

文章提到, 优化 $O_2$ 的代价比较大, 因为对每个节点都需要计算 $p_2\left(\cdot \mid v_i\right)$. 为此, 作者借鉴 Word2Vec 中的 Negative Sampling, 对负样本进行采样. 因此最终对每一条边 $(i,j)$, 优化如下目标函数:

$$\log \sigma \left({\vec{u}}_j^{\prime T}\cdot {\vec{u}}_i\right)+\sum _{n=1}^K{E}_{v_n\sim P_n(v)}\left[\log \sigma \left(-{\vec{u}}_n^{\prime T}\cdot {\vec{u}}_i\right)\right]$$

其中 $\sigma (x)=1/(1+\exp (-x))$ 为 sigmoid 函数, 第一项对正样本进行建模, 第二项描述负样本, $K$ 表示负样本的数量, 负样本的采样分布 $P_n(v)\propto d_v^{3/4}$, 其中 $d_v$ 表示节点 $v$ 的出度.

文章采用 ASGD (Asynchronous Stochastic Gradient Algorithm) 来优化上式. 注意一阶邻近和二阶邻近均采用上式作为目标函数.

由于目标函数 $O_1$ 和 $O_2$ 中均含有边权 $w_{ij}$, 如果边权过大, 学习速率就需要进行调小, 以避免梯度爆炸, 然而学习速率过小对于那些边权较小的边是不利的. 为了处理这个问题, 作者采用 AliasMethod 算法进行边的采样, 该方法的时间复杂度为 $O(1)$, 对于那些边权较大的边, 有更大的概率被采样到.

好吧, 论文介绍就到这吧, 下面看源码. 理论介绍是一回事, 代码实现又是另一回事 🤣🤣🤣

LINE 算法源码分析

代码编译

LINE 的 C++ 源码地址在 https://github.com/tangjianpku/LINE, 我看的是 Linux 版本, 但我是在 MacOS 上编译的. 编译时需要 GSL (GNU Scientific Library) 科学计算库, 以产生随机数. 在 MacOS 使用 brew install gsl 就能方便的安装. 之后运行 train_youtube.sh 即可.

LINE 算法的核心实现主要在 line.cpp 中, 因此抓大放小, 下面只介绍 line.cpp 中的代码. 另外如果有 Word2Vec 源码阅读经验的话, 看 LINE 的源码会有帮助, 毕竟 LINE 最后采用的是 Negative Sampling 来进行模型的优化.

Hash Table

这个不多介绍, 图文件读入进来, 节点使用字符串的表示, 需要转换为节点 id 的形式. 对于大规模的图网络, 肯定是需要用到 Hash 的. 代码中实现的 InsertHashTable 以及 SearchHashTable 考虑了哈希冲突的问题. 如图:

代码如下:

// 节点使用 ClassVertex 表示, 包含名字和边权 degree
struct ClassVertex {
	double degree;
	char *name;
};

void InsertHashTable(char *key, int value)
{
	int addr = Hash(key);
	while (vertex_hash_table[addr] != -1) addr = (addr + 1) % hash_table_size;
	vertex_hash_table[addr] = value;
}

int SearchHashTable(char *key)
{
	int addr = Hash(key);  // 获取节点 id
	while (1)
	{
		if (vertex_hash_table[addr] == -1) return -1;
		// 在 addr 位置上查找节点, 还需要保证名字相同, 否则就查找下一个位置
		if (!strcmp(key, vertex[vertex_hash_table[addr]].name)) return vertex_hash_table[addr];
		addr = (addr + 1) % hash_table_size;
	}
	return -1;
}

文件读取 - 获取边的关系

在 ReadData() 中实现对网络文件的读取, 达到对节点个数以及边的个数的统计, 使用 edge_source_id 保存所有源节点, edge_target_id 保存所有目标节点, edge_weight 保存所有边的权重.

edge_source_id = (int *)malloc(num_edges*sizeof(int));
edge_target_id = (int *)malloc(num_edges*sizeof(int));
edge_weight = (double *)malloc(num_edges*sizeof(double));

网络文件的格式是 name_v1 name_v2 weight, 读取文件的每一行, 并将节点名字映射为 hash id, 插入到 哈希表以及 edge_source_id 数组中:

vid = SearchHashTable(name_v1); // 现在哈希表中查找节点 name_v1, 找寻其 id
if (vid == -1) vid = AddVertex(name_v1); // 如果没有找到 name_v1, 那么就插入到 hash_table 中
vertex[vid].degree += weight;
edge_source_id[k] = vid;

对目标节点也做同样的处理.

初始化 AliasTable

初始化 AliasTable, 以便后面进行边的采样. 关于 AliasMethod, 强推 【数学】时间复杂度O(1)的离散采样算法—— Alias method/别名采样方法 这篇文章. 原来大致是, 对于输入的边权分布 $[p{w}_1,p{w}_2,\dots ,p{w}_N]$, 乘上 N 之后, 其中总有大于 1 的值, 也有小于 1 的值, 这些值的求和肯定等于 N. 如下图原始边权分布, 有很多矩阵块, 这些矩阵块的面积总和为 N. 现在要设计一种算法, 将大于面积大于 1 的矩阵块的部分面积, 填补到面积小于 1 的矩阵块上, 比如下图中将矩阵 4 的面积分别填补到了矩阵块 1 和 3 上, 将矩阵 2 的部分面积填充到矩阵块 5 上. 这样就得到了 AliasTable.

那么之后如何采样呢? 需要掷两次骰子 🎲, 第一次投掷, 选择 AliasTable 的第 i 个位置; 之后第二次投掷, 确定到底使用哪一个矩阵块 (即哪条边), 因为第二次投掷时, AliasTable 的第 i 个位置可能会有两种颜色的矩阵块 (比如第 0 个位置上, 有红色 4 和黄色 1 的矩阵块), 需要确定使用哪种颜色. 具体代码就不贴了, 理解了 AliasTable 的原理, 代码实现不难.

采样边

因为要投掷两次骰子 🎲, 因此 SampleAnEdge 函数需要传入两个随机数.

long long SampleAnEdge(double rand_value1, double rand_value2)
{
	long long k = (long long)num_edges * rand_value1;
	return rand_value2 < prob[k] ? k : alias[k];
}

Embedding Layer 初始化

为了实现二阶邻近, 每个节点都会有两个特征表示, 其中一个为节点自身的表达 emb_vertex, 另一个是节点作为其他节点的上下文时的表达 emb_context.

Negative Table 初始化

负样本的采样需要使用到 Negative Table, 采样的概率分布服从 $P_n(v)\propto d_v^{3/4}$, 这就是说, 对于节点 u 来说, 其出度为 $d_u^{3/4}$, 那么其出现的概率为:

$$P_n(u)\propto \frac{d_u^{3/4}}{\sum _{i=1}^{|V|}{d}_i^{3/4}}$$

假设 Negative Table 的大小为 $N$, 那么节点 $u$ 出现的次数应该是:

$$N\times \frac{d_u^{3/4}}{\sum _{i=1}^{|V|}{d}_i^{3/4}}$$

下面的代码就是为了在 Negative Table 中填充节点 id 值, 对于出度较大的节点, 在 Negative Table 中出现的次数越多.

/* Sample negative vertex samples according to vertex degrees */
// NEG_SAMPLING_POWER 定义为 0.75
void InitNegTable()
{
	double sum = 0, cur_sum = 0, por = 0;
	int vid = 0;
	neg_table = (int *)malloc(neg_table_size * sizeof(int));
	for (int k = 0; k != num_vertices; k++) sum += pow(vertex[k].degree, NEG_SAMPLING_POWER);
	for (int k = 0; k != neg_table_size; k++)
	{
		if ((double)(k + 1) / neg_table_size > por)
		{
			cur_sum += pow(vertex[vid].degree, NEG_SAMPLING_POWER);
			por = cur_sum / sum;
			vid++;
		}
		neg_table[k] = vid - 1;
	}
}

可以看下图对上面代码进行理解:

Sigmoid Table 的使用

由于 Sigmoid 中含有指数函数, 计算起来相对耗时, 可以使用 Sigmoid Table 提前算好一些值, 后面只需查询即可. 假设下面代码中, 使用 $B$ 来表示 SIGMOID_BOUND, 使用 $S$ 来表示 sigmoid_table_size, 那么计算 $[-B,B]$ 范围内的值的 Sigmoid 结果, 需要先对这个区间进行离散化, 长度为 $\frac{2\times B}{S}$, 那么给定索引 $k$, 得到对应的输入值 $x$ 为:

$$x=\frac{2\times B}{S}\times k-B$$

那么给定实数 $x$, 计算它在 Sigmoid Table 中的索引为:

$$k=\frac{(x+B)\times S}{2\times B}$$

/* Fastly compute sigmoid function */
void InitSigmoidTable()
{
	real x;
	sigmoid_table = (real *)malloc((sigmoid_table_size + 1) * sizeof(real));
	for (int k = 0; k != sigmoid_table_size; k++)
	{
		x = 2.0 * SIGMOID_BOUND * k / sigmoid_table_size - SIGMOID_BOUND;
		sigmoid_table[k] = 1 / (1 + exp(-x));
	}
}

real FastSigmoid(real x)
{
	if (x > SIGMOID_BOUND) return 1;
	else if (x < -SIGMOID_BOUND) return 0;
	int k = (x + SIGMOID_BOUND) * sigmoid_table_size / SIGMOID_BOUND / 2;
	return sigmoid_table[k];
}

负样本采用随机数

由于要在 Negative Table 中采样负样本, 需要随机产生索引:

/* Fastly generate a random integer */
int Rand(unsigned long long &seed)
{
	seed = seed * 25214903917 + 11;
	return (seed >> 16) % neg_table_size;
}

LINE 模型训练以及梯度更新

由于 LINE 使用 ASGD 对梯度进行更新, 代码中使用多线程来进行实现, 因此核心代码在 TrainLINEThread 中.

下面代码中使用 count 对样本进行计数, 第二个 if 用来对学习速率进行调整, 采用 Learning Rate Decay 的思路. 其中 rho 表示学习速率.

//judge for exit
		if (count > total_samples / num_threads + 2) break;

		if (count - last_count > 10000)
		{
			current_sample_count += count - last_count;
			last_count = count;
			printf("%cRho: %f  Progress: %.3lf%%", 13, rho, (real)current_sample_count / (real)(total_samples + 1) * 100);
			fflush(stdout);
			rho = init_rho * (1 - current_sample_count / (real)(total_samples + 1));
			if (rho < init_rho * 0.0001) rho = init_rho * 0.0001;
		}

核心算法:

		// 使用 AliasMethod 进行边的采样, 得到源节点 u, 正样本 v 
		curedge = SampleAnEdge(gsl_rng_uniform(gsl_r), gsl_rng_uniform(gsl_r));
		u = edge_source_id[curedge];
		v = edge_target_id[curedge];

		// dim 为 embedding 的大小, lu 表示 u 的 embedding 的起始位置, 因为矩阵使用数组来
		// 表示, embedding 矩阵大小为 [N, emb], 表示成一维数组的大小为 N*emb, 要找到 u 所
		// 代表的 emb, 需要用索引值 u 乘上 dim
		lu = u * dim;
		for (int c = 0; c != dim; c++) vec_error[c] = 0;

		// NEGATIVE SAMPLING
		for (int d = 0; d != num_negative + 1; d++)
		{
			if (d == 0)
			{
				target = v;
				label = 1;
			}
			else
			{
				target = neg_table[Rand(seed)]; // 从 Negative Table 采样负样本
				label = 0;
			}
			lv = target * dim;
			// 如果使用一阶邻近, 那么 Update 的输入均为 emb_vertex
			if (order == 1) Update(&emb_vertex[lu], &emb_vertex[lv], vec_error, label);
			// 如果使用二阶邻近, 那么 Update 的输入分别为 emb_vertex 以及 emb_context
			// emb_vertex 是节点作为自身所代表的 embedding,
			// emb_context 是节点作为其他节点的 context 所代表的 embedding
			if (order == 2) Update(&emb_vertex[lu], &emb_context[lv], vec_error, label);
		}
		// embedding 更新, 这个要根据对 loss 求导来得到, 具体见下面的讲解
		for (int c = 0; c != dim; c++) emb_vertex[c + lu] += vec_error[c];

代码中的思路主要是, 先使用 AliasMethod 对边进行采样, 得到 $(u_i,v_j)$ 边, 从而获得源节点 $u_i$ 和正样本 $v_j$, 之后采样 num_negative 个负样本, 从 Negative Table 中进行采样. 得到节点的 id 后, 然后获取它们分别对应的 embedding, 再根据 Negative Sampling 的目标函数进行 loss 的计算 (loss 没有显式计算出来) 和梯度的更新, 最后再更新 embedding.

SGD 的计算实现在 Update 函数中:

/* Update embeddings */
void Update(real *vec_u, real *vec_v, real *vec_error, int label)
{
	real x = 0, g;
	for (int c = 0; c != dim; c++) x += vec_u[c] * vec_v[c];
	g = (label - FastSigmoid(x)) * rho;
	for (int c = 0; c != dim; c++) vec_error[c] += g * vec_v[c];
	for (int c = 0; c != dim; c++) vec_v[c] += g * vec_u[c];
}

vec_u 为源节点的 embeddng, 而 vec_v 为目标节点的 embedding. 代码用公式来表示, 其实就简单多了. 首先, Negative Sampling 的目标函数为:

$$loss=\log \sigma \left(v_j^T\cdot {\vec{u}}_i\right)+\sum _{n=1}^K{E}_{v_n\sim P_n(v)}\left[\log \sigma \left(-{\vec{v}}_n^{\prime T}\cdot {\vec{u}}_i\right)\right]$$

先来看对 $loss$ 关于 $u_i$ 和 $v_k$ 的求导 (之所以用 $k$, 考虑到正负样本的问题).

$$\frac{\partial loss}{\partial u_i}=\left(1-\sigma \left(v_j^T\cdot {\vec{u}}_i\right)\right)\cdot v_j+\sum _{n=1}^K\left(-\sigma \left(v_n^T\cdot {\vec{u}}_i\right)\right)\cdot v_n$$

代码中的

$$\begin{array}{rl}x& =v_k^T\cdot u_i\\ g& =\left[\text{label}-\sigma (x)\right]\cdot \rho \end{array}$$

为了更新 $u_i$, 需要对 $v_k$ 的结果进行累加, 因此需要使用 vec_error 对误差进行累计.
而要更新 $v_k$, 有:

$$\begin{array}{rl}\frac{\partial loss}{\partial v_j}& =\left[1-\sigma \left(v_j^T\cdot {\vec{u}}_i\right)\right]\cdot u_i\\ \frac{\partial loss}{\partial v_n}& =\left[0-\sigma \left(v_n^T\cdot {\vec{u}}_i\right)\right]\cdot u_i\end{array}$$

由于无需累加, 因此在 Update 函数中使用 vec_v[c] += g * vec_u[c]; 进行了更新.

代码总结

Good, 有对 Word2Vec 源码的参考, 以后我去看 Word2Vec 的 C 源码会方便很多. 写文章的时间真的太久了, 但为了以后查阅方便, 以及对自己的总结提高, 还是有必要花这些时间的. OK, 晚上吃火锅! 🤣🤣🤣

资料总结

• Word2Vec 论文: Word2Vec
• AliasMethod 的详细介绍: 【数学】时间复杂度O(1)的离散采样算法—— Alias method/别名采样方法
• LINE:Large-scale Information Network Embedding阅读笔记: 评论区非常热闹和精彩, 是关于 LINE 算法更为深入的解读