在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 xVert_1 &= sum_{i=1}^n |x_i|\ Vert xVert_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=sumlimits_{i=1}^n p_i=sumlimits_{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_iequiv frac{1}{n}, q_jequiv 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_iequiv frac{1}{n}, q_jequiv frac{1}{n'}或许太粗糙，原论文提出的方法是令每个词的权重p_i=d_i=frac{c_i}{sumlimits_{j=1}^nc_j}，即这个词的权重为该词在整个文档中出现的频率，此时
d与d'表示两个文档，但实际上大部分情况我们不会划分两个文档，而是只用一个文档
这便是Word Mover's Distance(WMD)（推词机距离？？），大概可以理解为将一个句子变为另一个句子的最短路径，某种意义上也可以理解为编辑距离的光滑版。实际使用的时候，通常会去掉停用词再计算WMD
参考实现
参考实现如下：
代码语言：javascript复制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} leftVert boldsymbol{w}_i - boldsymbol{w}'_jrightVert_2 =& sum_{i,j} leftVert gamma_{i,j}(boldsymbol{w}_i - boldsymbol{w}'_j)rightVert_2\ geq& leftVert sum_{i,j}gamma_{i,j}(boldsymbol{w}_i - boldsymbol{w}'_j)rightVert_2\ =& leftVert sum_ileft(sum_jgamma_{i,j}right)boldsymbol{w}_i - sum_jleft(sum_igamma_{i,j}right)boldsymbol{w}'_jrightVert_2\ =& {color{red} {leftVert frac{1}{n}sum_iboldsymbol{w}_i - frac{1}{n'}sum_jboldsymbol{w}'_jrightVert_2}}\ end{aligned}tag{8}end{equation} $$
也就是说，WMD大于两个句子的平均向量的欧式距离，所以欧式距离大的两个句子，WMD一定大，因此我们要检索WMD比较小的句子时，可以先用欧式距离过滤掉距离比较大的句子，剩下的再采用WMD进行比较
Word Rotator's Distance
WMD其实已经听不错了，但非要鸡蛋里挑骨头的话，还是能挑出一些缺点来：
- 它使用的是欧式距离作为语义差距度量，但从Word2Vec的经验我们知道，用cos往往比欧式距离要好
- WMD理论上是一个无上界的量，这意味着我们不太容易直观感知相似程度，从而不能很好调整相似与否的阈值
为了解决这两个问题，一个 比较朴素的想法是将所有向量除以各自的模长归一化后再算WMD，但这样就完全失去模长信息了。最近的论文《Word Rotator's Distance: Decomposing Vectors Gives Better Representations》则巧妙的提出，在归一化的同时可以把模长融入到约束条件p,q里去，这就形成了WRD
基本形式
首先，WRD提出了"词向量的模长正相关于这个词的重要程度"的观点，并通过一些实验结果验证了这个观点。在WMD中，
$$ begin{equation}begin{aligned}&p_i = frac{leftVert boldsymbol{w}_irightVert _2}{Z}, &Z=sum_{i=1}^n leftVertboldsymbol{w}_irightVert_2\ &q_j = frac{leftVert boldsymbol{w}'_jrightVert_2}{Z'}, &Z'=sum_{j=1}^{n'}leftVertboldsymbol{w}'_jrightVert_2 end{aligned}tag{9}end{equation} $$
然后d_{i,j}就用余弦距离：
得到
注意Vert boldsymbol{w}_iVert_2和Vert boldsymbol{w}_j'Vert_2都是数，并且Vert boldsymbol{w}_iVert_2times Vert boldsymbol{w}_j'Vert_2得到的结果也是数，不是向量 boldsymbol{w}_i^{top}cdot boldsymbol{w_j'}是向量内积，得到的结果也是一个数
这就是Word Rotator's Distance（WRD）了。由于使用的度量是余弦距离，所以两个向量之间的变换更像是一种旋转（rotate）而不是移动（move），所以有了这个命名；同样由于使用了余弦距离，所以它的结果在[0,2]内，相对来说更容易去感知其相似程度
参考实现
参考实现如下：
代码语言：javascript复制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的一个下界公式：
$$ begin{equation}begin{aligned} 2sum_{i,j} gamma_{i,j} left(1 - frac{boldsymbol{w}_i^{top}cdot boldsymbol{w}'_j}{leftVertboldsymbol{w}_irightVert_2times leftVertboldsymbol{w}'_jrightVert_2}right)=&,sum_{i,j} gamma_{i,j} leftVert frac{boldsymbol{w}_i}{leftVert boldsymbol{w}_irightVert_2} - frac{boldsymbol{w}'_j}{leftVert boldsymbol{w}'_jrightVert_2}rightVert^2 \ geq&, sum_{i,j} leftVert gamma_{i,j}left(frac{boldsymbol{w}_i}{leftVert boldsymbol{w}_irightVert_2} - frac{boldsymbol{w}'_j}{leftVert boldsymbol{w}'_jrightVert_2}right)rightVert^2\ geq&, leftVert sum_{i,j}gamma_{i,j}left(frac{boldsymbol{w}_i}{leftVert boldsymbol{w}_irightVert_2} - frac{boldsymbol{w}'_j}{leftVert boldsymbol{w}'_jrightVert_2}right)rightVert^2\ =& ,leftVert sum_ileft(sum_jgamma_{i,j}right)frac{boldsymbol{w}_i}{leftVert boldsymbol{w}_irightVert_2} - sum_jleft(sum_igamma_{i,j}right)frac{boldsymbol{w}'_j}{leftVert boldsymbol{w}'_jrightVert_2}rightVert^2\ =&, leftVert frac{1}{Z}sum_iboldsymbol{w}_i - frac{1}{Z'}sum_jboldsymbol{w}'_jrightVert^2\ end{aligned}tag{11}end{equation} $$
其中第一个等号基于一个简单的数学常识(x-y)^2=x^2-2xy y^2。我们从右往左推，即 $$ begin{aligned} leftVertfrac{boldsymbol{w}_i}{leftVert boldsymbol{w}_irightVert_2} - frac{boldsymbol{w}'_j}{leftVert boldsymbol{w}'_jrightVert_2}rightVert^2 =leftVert frac{boldsymbol{w}_i}{leftVert boldsymbol{w}_irightVert_2}rightVert^2 - 2cdotfrac{boldsymbol{w}_i^{top}cdot boldsymbol{w}_j'}{leftVertboldsymbol{w}_i rightVert_2times leftVert boldsymbol{w}_j'rightVert_2} leftVert frac{boldsymbol{w}_j'}{leftVert boldsymbol{w}_j'rightVert_2}rightVert^2 end{aligned} $$ 且（下面式子是定理，大家可以随便带个例子进去算一下是成立的）
所以 $$ begin{aligned} leftVertfrac{boldsymbol{w}_i^{top}}{leftVert boldsymbol{w}_irightVert_2} - frac{boldsymbol{w}'_j}{leftVert boldsymbol{w}'_jrightVert_2}rightVert^2 =&, 2 - 2cdotfrac{boldsymbol{w}_i^{top}cdot boldsymbol{w}_j'}{leftVertboldsymbol{w}_i rightVert_2times leftVert boldsymbol{w}_j'rightVert_2}\ =&, 2left(1 - frac{boldsymbol{w}_i^{top}cdot boldsymbol{w}'_j}{leftVertboldsymbol{w}_irightVert_2times leftVertboldsymbol{w}'_jrightVert_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}{leftVert boldsymbol{w}_irightVert_2} - frac{boldsymbol{w}'_j}{leftVert boldsymbol{w}'_jrightVert_2}right)这一项无法保证一定大于0，所以先对这一项进行norm再求和，肯定是大于等于先求和再norm的
参考实现
下面对下界公式给出一个代码实现：
代码语言：javascript复制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 disReferences
- 从EMD、WMD到WRD：文本向量序列的相似度计算
- Word Rotator‘s Distance——WRD算法应用
