从EMD、WMD、WRD：文本向量序列的相似度计算

在NLP中，我们经常要比较两个句子的相似度，其标准方法是将句子编码为固定大小的向量，然后用某种几何距离（欧氏距离、cos距离等）作为相似度。这种方案相对来说比较简单，而且检索起来比较快速，一定程度上能满足工程需求

此外，还可以直接比较两个变长序列的差异性，比如编辑距离，它通过动态规划找出两个字符串之间的最优映射，然后算不匹配程度；现在我们还有Word2Vec、BERT等工具，可以将文本序列转换为对应的向量序列，所以也可以直接比较这两个向量序列的差异，而不是先将向量序列弄成单个向量。

后一种方案速度相对慢一点，但可以比较得更精细一些，并且理论比较优雅，所以也有一定的应用场景。本文就来简单介绍一下属于后者的两个相似度指标，分别简称为WMD、WRD

Earth Mover's Distance

假设现在有两个概率分布p({x}),q({x})，那么Wasserstein距离的定义为

我们来逐项看看

成本函数

首先d(x,y)，它不一定是距离，其准确含义应该是一个成本函数，代表着从x运输到y的成本。常用的d是基于l范数衍生出来的，例如

都是常见的选择，其中

$$ \begin{aligned} \Vert x\Vert_1 &= \sum_{i=1}^n |x_i|\\ \Vert x\Vert_2 &= \sqrt{\sum_{i=1}^n x_i^2} \end{aligned} \tag{3} $$

特别地，其实哪种距离并不是特别重要，因为很多范数都是相互等价的，范数的等价性表明其实最终定义出来的\mathcal{W}距离都差不多

成本最小化

然后来看\gamma，条件\gamma \in \Pi[p,q]意味着：

$$ \begin{aligned} \int \gamma (x,y)dy &= p(x)\\ \int \gamma (x,y)dx &= p(y) \end{aligned} \tag{4} $$

也就是说\gamma是联合分布，它的边缘分布就是原来的p和q。事实上\gamma就描述了一种运输方案，不失一般性，设p是原始分布，设q是目标分布，p(x)的意思是原来在位置x处有p(x)量的货物，q(x)是指最终x处要存放的货物量，如果p(x)>q(x)x处的一部分运到别处；反之如果p(x)<q(x)x处。而\gamma (x,y)的意思是指，要从x处搬\gamma (x,y)dx那么多的东西到y处

最后是\inf，这表示下确界，简单来说就是取最小，也就是说，要从所有的运输方案中，找出总运输成本\iint \gamma (x,y)d(x,y)dxdy最小的方案，这个方案的成本，就是我们要算的\mathcal{W}[p,q]。如果将上述比喻中的“货物”换成“沙土”，那么Wasserstein距离就是在求最省力的“搬土”方案了，所以Wasserstein距离也被称为“推土机距离”（Earth Mover's Distance）

左边p(x)每处的沙土被分为若干份，然后运输到右端q(x)同色的位置（或者不动）

最优传输

假设在位置i=1,2,...,n处我们分布有p_1,p_2,...,p_n那么多的土，简单起见我们设土的总数量为1，即p_1+p_2+\cdots +p_n=1，现在要将土推到位置j=1,2,...,n'上，每处的量为q_1,q_2,...,q_{n'}，而从i推到j的成本为d_{ij}，求成本最低的方案以及对应的最低成本

这其实就是一个经典的最优传输问题。我们将最优方案表示为\gamma_{i,j}，表示这个方案中要从i把\gamma_{i,j}数量的土推到j处，很明显我们有约束

所以我们的优化问题是

参考实现

看上去很复杂，但认真观察下就能发现上式其实就是一个线性规划问题——在线性约束下求线性函数的极值。而scipy本身自带了线性规划求解函数linprog，因此我们可以利用它实现求Wasserstein距离的函数

import numpy as np
from scipy.optimize import linprog
def wasserstein_distance(p, q, D):

'''

通过线性规划求Wasserstein距离

p.shape=[m], q.shape=[n], D.shape=[m,n]

p.sum()=1, q.sum()=1, p∈[0,1], q∈[0,1]

'''

A_eq = []

for i in range(len(p)):

A = np.zeros_like(D)

A[i, :] = 1

A_eq.append(A.reshape(-1))

for i in range(len(q)):

A = np.zeros_like(D)

A[:, i] = 1

A_eq.append(A.reshape(-1))

A_eq = np.array(A_eq)

b_eq = np.concatenate([p, q])

D = np.array(D)

D = D.reshape(-1)

result = linprog(D, A_eq=A_eq[:-1], b_eq=b_eq[:-1])

return result.fun

读者会发现，在传入约束的时候用的是A_eq=A_eq[:-1], b_eq=b_eq[:-1]，也就是去掉了最后一个约束。这是因为1=\sum\limits_{i=1}^n p_i=\sum\limits_{j=1}^{n'}q_j，所以(1)中的等式约束本身存在冗余，而实际计算中有时候存在浮点误差，导致冗余的约束之间相互矛盾，从而使得线性规划的求解失败，所以干脆去掉最后一个冗余的约束，减少出错的可能性

Word Mover's Distance

很明显，Wasserstein距离适合于用来计算两个长度不同的序列的差异性，而我们要做语义相似度的时候，两个句子的长度通常也是不一样的，刚好对应这个特性，因此很自然地就会联想到Wasserstein距离也许可以用来比较句子相似度，首次进行这个尝试的是论文《From Word Embeddings To Document Distances》

基本形式

设有两个句子s=(t_1,t_2,...,t_n), s'=(t_1',t_2',...,t_{n'}')，经过某种映射（比如Word2Vec或BERT）后，它们变成了对应的向量序列(\boldsymbol{w}_1,\boldsymbol{w}_2,\dots,\boldsymbol{w}_n), (\boldsymbol{w}'_1, \boldsymbol{w}'_2, \dots, \boldsymbol{w}'_{n'})，现在我们就想办法用Wasserstein距离来比较这两个序列的相似度

根据前一节的介绍，Wasserstein距离需要知道p_i,q_j,d_{i,j}三个量，我们逐一把它们都定义好即可。由于没有什么先验知识，所以我们可以很朴素地令p_i\equiv \frac{1}{n}, q_j\equiv \frac{1}{n'}，因此现在还剩d_{i,j}。显然，d_{i,j}代表着第一个序列的向量\boldsymbol{w}_1与第二个序列的向量\boldsymbol{w}_j'的某种差异性，简单起见我们可以用欧式距离\Vert \boldsymbol{w}_i - \boldsymbol{w}_j'\Vert_2，所以两个句子的差异程度可以表示为

直接朴素地令p_i\equiv \frac{1}{n}, q_j\equiv \frac{1}{n'}或许太粗糙，原论文提出的方法是令每个词的权重p_i=d_i=\frac{c_i}{\sum\limits_{j=1}^nc_j}，即这个词的权重为该词在整个文档中出现的频率，此时

d与d'表示两个文档，但实际上大部分情况我们不会划分两个文档，而是只用一个文档

这便是Word Mover's Distance(WMD)（推词机距离？？），大概可以理解为将一个句子变为另一个句子的最短路径，某种意义上也可以理解为编辑距离的光滑版。实际使用的时候，通常会去掉停用词再计算WMD

参考实现

参考实现如下：

def word_mover_distance(x, y):
    '''WMD (Word Mover's Distance) 参考实现
    x.shape=[m, d], y.shape=[n, d]
    '''
    x, y = np.array(x), np.array(y)
    p = np.ones(x.shape[0]) / x.shape[0]
    q = np.ones(y.shape[0]) / y.shape[0]
    D = np.sqrt(np.square(x[:, None] - y[None, :]).mean(axis=2))
    return wasserstein_distance(p, q, D)

下界公式

如果是检索场景，要将输入句子跟数据库里的所有句子一一算WMD并排序的话，那计算成本是相当大的，所以我们要尽量减少算WMD的次数，比如通过一些更简单高效的指标来过滤掉一些样本，然后再对剩下的样本算WMD

幸运的是，我们确实可以推导出WMD的一个下界公式，原论文称之为Word Centroid Distance（WCD）：

$$ \begin{equation}\begin{aligned} \sum_{i,j} \gamma_{i,j} \left\Vert \boldsymbol{w}_i - \boldsymbol{w}'_j\right\Vert_2 =& \sum_{i,j} \left\Vert \gamma_{i,j}(\boldsymbol{w}_i - \boldsymbol{w}'_j)\right\Vert_2\\ \geq& \left\Vert \sum_{i,j}\gamma_{i,j}(\boldsymbol{w}_i - \boldsymbol{w}'_j)\right\Vert_2\\ =& \left\Vert \sum_i\left(\sum_j\gamma_{i,j}\right)\boldsymbol{w}_i - \sum_j\left(\sum_i\gamma_{i,j}\right)\boldsymbol{w}'_j\right\Vert_2\\ =& {\color{red} {\left\Vert \frac{1}{n}\sum_i\boldsymbol{w}_i - \frac{1}{n'}\sum_j\boldsymbol{w}'_j\right\Vert_2}}\\ \end{aligned}\tag{8}\end{equation} $$

也就是说，WMD大于两个句子的平均向量的欧式距离，所以欧式距离大的两个句子，WMD一定大，因此我们要检索WMD比较小的句子时，可以先用欧式距离过滤掉距离比较大的句子，剩下的再采用WMD进行比较

Word Rotator's Distance

WMD其实已经听不错了，但非要鸡蛋里挑骨头的话，还是能挑出一些缺点来：

1. 它使用的是欧式距离作为语义差距度量，但从Word2Vec的经验我们知道，用cos往往比欧式距离要好
2. WMD理论上是一个无上界的量，这意味着我们不太容易直观感知相似程度，从而不能很好调整相似与否的阈值

为了解决这两个问题，一个 比较朴素的想法是将所有向量除以各自的模长归一化后再算WMD，但这样就完全失去模长信息了。最近的论文《Word Rotator's Distance: Decomposing Vectors Gives Better Representations》则巧妙的提出，在归一化的同时可以把模长融入到约束条件p,q里去，这就形成了WRD

基本形式

首先，WRD提出了"词向量的模长正相关于这个词的重要程度"的观点，并通过一些实验结果验证了这个观点。在WMD中，

$$ \begin{equation}\begin{aligned}&p_i = \frac{\left\Vert \boldsymbol{w}_i\right\Vert _2}{Z}, &Z=\sum_{i=1}^n \left\Vert\boldsymbol{w}_i\right\Vert_2\\ &q_j = \frac{\left\Vert \boldsymbol{w}'_j\right\Vert_2}{Z'}, &Z'=\sum_{j=1}^{n'}\left\Vert\boldsymbol{w}'_j\right\Vert_2 \end{aligned}\tag{9}\end{equation} $$

然后d_{i,j}就用余弦距离：

得到

注意\Vert \boldsymbol{w}_i\Vert_2和\Vert \boldsymbol{w}_j'\Vert_2都是数，并且\Vert \boldsymbol{w}_i\Vert_2\times \Vert \boldsymbol{w}_j'\Vert_2得到的结果也是数，不是向量 \boldsymbol{w}_i^{\top}\cdot \boldsymbol{w_j'}是向量内积，得到的结果也是一个数

这就是Word Rotator's Distance（WRD）了。由于使用的度量是余弦距离，所以两个向量之间的变换更像是一种旋转（rotate）而不是移动（move），所以有了这个命名；同样由于使用了余弦距离，所以它的结果在[0,2]内，相对来说更容易去感知其相似程度

参考实现

参考实现如下：

def word_rotator_distance(x, y):
    """WRD（Word Rotator's Distance）的参考实现
    x.shape=[m,d], y.shape=[n,d]
    """
    x, y = np.array(x), np.array(y)
    x_norm = (x**2).sum(axis=1, keepdims=True)**0.5
    y_norm = (y**2).sum(axis=1, keepdims=True)**0.5
    p = x_norm[:, 0] / x_norm.sum()
    q = y_norm[:, 0] / y_norm.sum()
    D = 1 - np.dot(x / x_norm, (y / y_norm).T)
    return wasserstein_distance(p, q, D)
def word_rotator_similarity(x, y):

"""1 - WRD

x.shape=[m,d], y.shape=[n,d]

"""

return 1 - word_rotator_distance(x, y)

下界公式

同WMD一样，我们也可以推导出WRD的一个下界公式：

其中第一个等号基于一个简单的数学常识(x-y)^2=x^2-2xy+y^2。我们从右往左推，即

\begin{aligned} \left\Vert\frac{\boldsymbol{w}_i}{\left\Vert \boldsymbol{w}_i\right\Vert_2} - \frac{\boldsymbol{w}&#39;_j}{\left\Vert \boldsymbol{w}&#39;_j\right\Vert_2}\right\Vert^2 =\left\Vert \frac{\boldsymbol{w}_i}{\left\Vert \boldsymbol{w}_i\right\Vert_2}\right\Vert^2 - 2\cdot\frac{\boldsymbol{w}_i^{\top}\cdot \boldsymbol{w}_j&#39;}{\left\Vert\boldsymbol{w}_i \right\Vert_2\times \left\Vert \boldsymbol{w}_j&#39;\right\Vert_2}+\left\Vert \frac{\boldsymbol{w}_j&#39;}{\left\Vert \boldsymbol{w}_j&#39;\right\Vert_2}\right\Vert^2 \end{aligned}

且（下面式子是定理，大家可以随便带个例子进去算一下是成立的）

所以 $$ \begin{aligned} \left\Vert\frac{\boldsymbol{w}_i^{\top}}{\left\Vert \boldsymbol{w}_i\right\Vert_2} - \frac{\boldsymbol{w}'_j}{\left\Vert \boldsymbol{w}'_j\right\Vert_2}\right\Vert^2 =&\, 2 - 2\cdot\frac{\boldsymbol{w}_i^{\top}\cdot \boldsymbol{w}_j'}{\left\Vert\boldsymbol{w}_i \right\Vert_2\times \left\Vert \boldsymbol{w}_j'\right\Vert_2}\\ =&\, 2\left(1 - \frac{\boldsymbol{w}_i^{\top}\cdot \boldsymbol{w}'_j}{\left\Vert\boldsymbol{w}_i\right\Vert_2\times \left\Vert\boldsymbol{w}'_j\right\Vert_2}\right) \end{aligned} $$

第一个不等号基于Jensen不等式，实际上Jensen不等式的特殊形式大家应该都见过：给定下凸函数f，则对于给定区间内任意两点x_1,x_2，有不等式\frac{f(x_1)+f(x_2)}{2}\ge f(\frac{x_1+x_2}{2})

第二个不等号基于三角不等式，或者简单的判断一下，因为\gamma_{i,j}\left(\frac{\boldsymbol{w}_i}{\left\Vert \boldsymbol{w}_i\right\Vert_2} - \frac{\boldsymbol{w}'_j}{\left\Vert \boldsymbol{w}'_j\right\Vert_2}\right)这一项无法保证一定大于0，所以先对这一项进行norm再求和，肯定是大于等于先求和再norm的

参考实现

下面对下界公式给出一个代码实现：

def dis_lower_boundary(x,y):
    """
    WRD的一个下界距离
    """
    x, y = np.array(x), np.array(y)
    x_norm = (x ** 2).sum(axis=1, keepdims=True) ** 0.5
    y_norm = (y ** 2).sum(axis=1, keepdims=True) ** 0.5
    z_x = x.sum(axis=0)/x_norm.sum()
    z_y = y.sum(axis=0)/y_norm.sum()
    dis = ((z_x-z_y) ** 2).sum()**0.5 * 0.5 # 别忘了最后要乘以1/2
    return dis

References

• 从EMD、WMD到WRD：文本向量序列的相似度计算
• Word Rotator‘s Distance——WRD算法应用