文本相似度计算/文本比较算法

参考： 

文本比较算法Ⅰ——LD算法

文本比较算法Ⅱ——Needleman/Wunsch算法

文本比较算法Ⅲ——计算文本的相似度

文本比较算法Ⅳ——Nakatsu算法

目录：

问题

LD算法

Needleman/Wunsch算法

Nakatsu算法

问题

字符串s1 和 字符串s2 的比较算法 ==> 相似度 or 差异性。

主流的算法有两大类：

• 基于编辑距离 ( Edit Distance)，例如：LD算法；
• 基于最长公共子串 ( Longest Common Subsequence），例如：Needleman/Wunsch算法等。

LD算法

LD算法（Levenshtein Distance）又称为编辑距离算法（Edit Distance）：以字符串A通过插入字符、删除字符、替换字符变成另一个字符串B，其中，这三类操作的总次数可表示两个字符串的差异。

　　例如：字符串A：kitten如何变成字符串B：sitting。

　　　　第一步：kitten——》sitten。k替换成s

　　　　第二步：sitten——》sittin。e替换成i

　　　　第三步：sittin——》sitting。在末尾插入g

　　故kitten和sitting的编辑距离为3。

定义说明：

　　LD(A,B)表示字符串A和字符串B的编辑距离。若LD(A, B)=0表示字符串A和字符串B完全相同；

　　Rev(A)表示反转字符串A；

　　Len(A)表示字符串A的长度；

　　A+B表示连接字符串A和字符串B；

　　A=a₁a₂……a_N，表示A是由a₁a₂……a_N这N个字符组成，Len(A)=N

　　B=b₁b₂……b_M，表示B是由b₁b₂……b_M这M个字符组成，Len(B)=M

　　定义LD(i, j) = LD(a₁a₂……a_i, b₁b₂……b_j)，其中0≤i≤N，0≤j≤M， ===> LD(N,M)=LD(A,B)

　　

　　对于1≤i≤N，1≤j≤M，有如下公式（动态规划的思想）：

　　若a_i=b_j，则LD(i,j) = LD(i-1,j-1) [无需替换 删除 插入]

　　若a_i≠b_j，则LD(i,j) = Min(LD(i-1,j-1)[可以在此基础上经过替换变成另一字符串],  LD(i-1,j),LD(i,j-1)[可以在此基础上经过插入or删除变成另一字符串])+1

　　

举例说明：

A=GGATCGA，B=GAATTCAGTTA，计算LD(A,B)

　　第一步：初始化LD矩阵　　

　　

LD算法矩阵

		G	A	A	T	T	C	A	G	T	T	A
	0	1	2	3	4	5	6	7	8	9	10	11
G	1
G	2
A	3
T	4
C	5
G	6
A	7

　　第二步：利用上述的公式，计算第一行

LD算法矩阵

		G	A	A	T	T	C	A	G	T	T	A
	0	1	2	3	4	5	6	7	8	9	10	11
G	1	0	1	2	3	4	5	6	7	8	9	10
G	2
A	3
T	4
C	5
G	6
A	7

　　第三步，利用上述的公式，计算其余各行 

LD算法矩阵

		G	A	A	T	T	C	A	G	T	T	A
	0	1	2	3	4	5	6	7	8	9	10	11
G	1	0	1	2	3	4	5	6	7	8	9	10
G	2	1	1	2	3	4	5	6	6	7	8	9
A	3	2	1	1	2	3	4	5	6	7	8	8
T	4	3	2	2	1	2	3	4	5	6	7	8
C	5	4	3	3	2	2	2	3	4	5	6	7
G	6	5	4	4	3	3	3	3	3	4	5	6
A	7	6	5	4	4	4	4	3	4	4	5	5

　　则LD(A,B) = LD(7,11) = 5 

代码实现（C++）：

class Solution {
public:
    int minDistance(string word1, string word2) {
        int n = word1.length();
        int m = word2.length();

        // 有一个字符串为空串
        if (n * m == 0) return n + m;

        // DP 数组
        int D[n + 1][m + 1];

        // 边界状态初始化
        for (int i = 0; i < n + 1; i++) {
            D[i][0] = i;
        }
        for (int j = 0; j < m + 1; j++) {
            D[0][j] = j;
        }

        // 计算所有 DP 值
        for (int i = 1; i < n + 1; i++) {
            for (int j = 1; j < m + 1; j++) {
                int left = D[i - 1][j] + 1;
                int down = D[i][j - 1] + 1;
                int left_down = D[i - 1][j - 1];
                if (word1[i - 1] != word2[j - 1]) left_down += 1;
                D[i][j] = min(left, min(down, left_down));

            }
        }
        return D[n][m];
    }
};

　　

　　我们往往不仅仅是计算出字符串A和字符串B的编辑距离，还要能得出它们的匹配结果。

　　以上面为例A=GGATCGA，B=GAATTCAGTTA，LD(A,B)=5

　　他们的匹配为：

　　　　A：GGA_TC_G__A

　　　　B：GAATTCAGTTA

　　如上面所示，蓝色表示完全匹配，黑色表示编辑操作，_表示插入字符或者是删除字符操作。如上面所示，黑色字符有5个，表示编辑距离为5。

　　利用上面的LD矩阵，通过回溯，能找到匹配字串。

　　第一步：定位在矩阵的右下角　　

LD算法矩阵

		G	A	A	T	T	C	A	G	T	T	A
	0	1	2	3	4	5	6	7	8	9	10	11
G	1	0	1	2	3	4	5	6	7	8	9	10
G	2	1	1	2	3	4	5	6	6	7	8	9
A	3	2	1	1	2	3	4	5	6	7	8	8
T	4	3	2	2	1	2	3	4	5	6	7	8
C	5	4	3	3	2	2	2	3	4	5	6	7
G	6	5	4	4	3	3	3	3	3	4	5	6
A	7	6	5	4	4	4	4	3	4	4	5	5

　　第二步：回溯单元格，至矩阵的左上角

　　　　若a_i=b_j，则回溯到左上角单元格

LD算法矩阵

		G	A	A	T	T	C	A	G	T	T	A
	0	1	2	3	4	5	6	7	8	9	10	11
G	1	0	1	2	3	4	5	6	7	8	9	10
G	2	1	1	2	3	4	5	6	6	7	8	9
A	3	2	1	1	2	3	4	5	6	7	8	8
T	4	3	2	2	1	2	3	4	5	6	7	8
C	5	4	3	3	2	2	2	3	4	5	6	7
G	6	5	4	4	3	3	3	3	3	4	5	6
A	7	6	5	4	4	4	4	3	4	4	5	5

　　　　若a_i≠b_j，回溯到 (左上角、上边、左边) 里值最小的单元格，若有相同最小值的单元格，优先级按照左上角、上边、左边的顺序

LD算法矩阵

		G	A	A	T	T	C	A	G	T	T	A
	0	1	2	3	4	5	6	7	8	9	10	11
G	1	0	1	2	3	4	5	6	7	8	9	10
G	2	1	1	2	3	4	5	6	6	7	8	9
A	3	2	1	1	2	3	4	5	6	7	8	8
T	4	3	2	2	1	2	3	4	5	6	7	8
C	5	4	3	3	2	2	2	3	4	5	6	7
G	6	5	4	4	3	3	3	3	3	4	5	6
A	7	6	5	4	4	4	4	3	4	4	5	5

LD算法矩阵

		G	A	A	T	T	C	A	G	T	T	A
	0	1	2	3	4	5	6	7	8	9	10	11
G	1	0	1	2	3	4	5	6	7	8	9	10
G	2	1	1	2	3	4	5	6	6	7	8	9
A	3	2	1	1	2	3	4	5	6	7	8	8
T	4	3	2	2	1	2	3	4	5	6	7	8
C	5	4	3	3	2	2	2	3	4	5	6	7
G	6	5	4	4	3	3	3	3	3	4	5	6
A	7	6	5	4	4	4	4	3	4	4	5	5

　　　　依照上面的回溯法则，回溯到矩阵的左上角

　　第三步：根据三个不同方向的回溯路径，写出匹配字串

　　　　若回溯到左上角单元格，将a_i添加到匹配字串A，将b_j添加到匹配字串B 

　　　　若回溯到上边单元格，将a_i添加到匹配字串A，将_添加到匹配字串B

　　　　若回溯到左边单元格，将_添加到匹配字串A，将b_j添加到匹配字串B

　　　　搜索晚整个匹配路径，匹配字串也就完成了

　　从上面可以看出，LD算法在不需要计算出匹配字串的话，时间复杂度为O(MN)，空间复杂度经优化后为O(M)

　　不过，如果要计算匹配字符串的话，时间复杂度为O(MN)，空间复杂度由于需要利用LD矩阵计算匹配路径，故空间复杂度仍然为O(MN)。这个在两个字符串都比较短小的情况下，能获得不错的性能。不过，如果字符串比较长的情况下，就需要极大的空间存放矩阵。例如：两个字符串都是20000字符，则LD矩阵的大小为20000*20000*2=800000000Byte=800MB。在比较长字符串的时候，还有其他性能更好的算法。

Needleman/Wunsch算法 

该算法是基于最长公共子串(不一定连续出现)的文本比较算法。

实例说明：

对于字符串A=kitten，字符串B=sitting

最长公共子串为ittn（注：最长公共子串不需要连续出现，但一定是出现的顺序一致），最长公共子串长度为4。

定义：

　　LCS(A,B)表示字符串A和字符串B的最长公共子串的长度，LSC(A,B)=0表示两个字符串没有公共部分。

　　Rev(A)表示反转字符串A

　　Len(A)表示字符串A的长度

　　A+B表示连接字符串A和字符串B

　　

　　A=a₁a₂……a_N，表示A是由a₁a₂……a_N这N个字符组成，Len(A)=N

　　B=b₁b₂……b_M，表示B是由b₁b₂……b_M这M个字符组成，Len(B)=M

　　定义LCS(i,j)=LCS(a₁a₂……a_i,b₁b₂……b_j)，其中0≤i≤N，0≤j≤M

　　对于1≤i≤N，1≤j≤M，有：

　　若a_i=b_j，则LCS(i,j)=LCS(i-1,j-1)+1

　　若a_i≠b_j，则LCS(i,j)=Max(LCS(i-1,j-1), LCS(i-1,j), LCS(i,j-1))    ==> 类似编辑距离

计算LCS(A,B)的算法有很多，下面介绍的Needleman/Wunsch算法是其中的一种。和LD算法类似，Needleman/Wunsch算法用的都是动态规划的思想。在Needleman/Wunsch算法中还设定了一个权值，用以区分三种操作（插入、删除、更改）的优先级。在下面的算法中，认为三种操作的优先级都一样。故权值默认为1。

举例说明：

A=GGATCGA，B=GAATTCAGTTA，计算LCS(A,B)

　　第一步：初始化LCS矩阵

Needleman/Wunsch算法矩阵

		G	A	A	T	T	C	A	G	T	T	A
	0	0	0	0	0	0	0	0	0	0	0	0
G	0
G	0
A	0
T	0
C	0
G	0
A	0

　　第二步：计算矩阵的第一行

Needleman/Wunsch算法矩阵

		G	A	A	T	T	C	A	G	T	T	A
	0	0	0	0	0	0	0	0	0	0	0	0
G	0	1	1	1	1	1	1	1	1	1	1	1
G	0
A	0
T	0
C	0
G	0
A	0

 　　第三步：计算矩阵的其余各行

Needleman/Wunsch算法矩阵

		G	A	A	T	T	C	A	G	T	T	A
	0	0	0	0	0	0	0	0	0	0	0	0
G	0	1	1	1	1	1	1	1	1	1	1	1
G	0	1	1	1	1	1	1	1	2	2	2	2
A	0	1	2	2	2	2	2	2	2	2	2	2
T	0	1	2	2	3	3	3	3	3	3	3	3
C	0	1	2	2	3	3	4	4	4	4	4	4
G	0	1	2	2	3	3	3	4	5	5	5	5
A	0	1	2	3	3	3	3	4	5	5	5	6

　　则，LCS(A,B)=LCS(7,11)=6

　　可以看出，Needleman/Wunsch算法实际上和LD算法是非常接近的。故他们的时间复杂度和空间复杂度也一样。时间复杂度为O(MN)，空间复杂度为O(MN)。空间复杂度经过优化，可以优化到O(M)，但是一旦优化就丧失了计算匹配字串的机会了。代码和LD算法几乎一样。

代码（c++):

int needleman(char s1[], char s2[]){
    int len1 = s1.size();
    int len2 = s2.size();
    int i,j;
    vector<vector<int>> count(len1+1,vector<int>(len2+1, 0) );

    for(i=1;i<len1+1;i++){
        for(j=1;j<len2+1;j++){
            if(s1[i] == s2[j]){
                count[i][j] = count[i-1][j-1]+1;
            } else {
                count[i][j] = max(max(count[i-1][j-1],count[i][j-1]),count[i-1][j]);
            }
        }
    }
    return count[len1][len2];
}

　　LD算法和Needleman/Wunsch算法的回溯路径是一样的。这样找到的匹配字串也是一样的。

　　不过，Needleman/Wunsch算法和LD算法一样，若要找出匹配字串，空间的复杂度就一定是O(MN)，在文本比较长的时候，是极为耗用存储空间的。故若要计算出匹配字串，还得用其他的算法。

　　

 计算文本相似度

　　在给定的字符串A和字符串B，LD(A,B)表示编辑距离，LCS(A,B)表示最长公共子串的长度。

　　如何来度量它们之间的相似度呢？

　　不妨设S(A,B)来表示字符串A和字符串B的相似度。那么，比较合理的相似度应该满足下列性质。

　　性质一：0≤S(A,B)≤100%，0表示完全不相似，100%表示完全相等

　　性质二：S(A,B)=S(B,A)

　　目前，网上介绍的各种相似度的计算，都有各自的不尽合理的地方。

　　计算公式一：S(A,B)=1/(LD(A,B)+1)

　　　　能完美的满足性质二。

　　　　当LD(A,B)=0时，S(A,B)=100%，不过无论LD(A,B)取任何值，S(A,B)＞0，不能满足性质一。

　　计算公式二：S(A,B)=1－LD(A,B)/Len(A)

　　　　当Len(B)＞Len(A)时，S(A,B)＜0。不满足性质一。

　　　　有人会说，当S(A,B)<0时，强制指定S(A,B)=0就解决问题了。

　　　　问题是，S(A,B)=1－LD(A,B)/Len(A)，而S(B,A)=1－LD(B,A)/Len(B)。当Len(A)≠Len(B)时，S(A,B)≠S(B,A)。不满足性质二

　　　　还有一个例子可以说明问题

　　　　A="BC"，B="CD"，C="EF"

　　　　S(A,B)=1－LD(A,B)/Len(A)=1－2/2=0

　　　　S(A,C)=1－LD(A,C)/Len(A)=1－2/2=0

　　　　A和B的相似度与A和C的相似度是一样的。不过很明显的是B比C更接近A

　　计算公式三：S(A,B)=LCS(A,B)/Len(A)

　　　　这个公式能完美的满足的性质一

　　　　不过当Len(A)≠Len(B)时，S(A,B)≠S(B,A)。不满足性质二

　　　　用一个例子说明问题

　　　　A="BC"，B="BCD"，C="BCEF"

　　　　S(A,B)=LCS(A,B)/Len(A)=2/2=100%

　　　　S(A,C)=LCS(A,C)/Len(A)=2/2=100%

　　　　A和B的相似度与A和C的相似度是一样的。不过很明显的是B比C更接近A

　　上面是网上能找到的三个计算公式，从上面的分析来看都有各自的局限性。

　　我们看例子：

　　A=GGATCGA，B=GAATTCAGTTA，LD(A,B)=5，LCS(A,B)=6

　　他们的匹配为：

　　　　A：GGA_TC_G__A

　　　　B：GAATTCAGTTA

　　给出一个新的公式

　　S(A,B)=LCS(A,B)/(LD(A,B)+LCS(A,B))

　　这个公式能解决上述三个公式的种种不足。

　　而LD(A,B)+LCS(A,B)表示两个字符串A、B的最佳匹配字串的长度。这个是唯一的。

　　还有注意的是LD(A,B)+LCS(A,B)和Max(Len(A),Len(B))这两个并不完全相等。

　　

Nakatsu算法

LD算法和LCS算法都是基于动态规划的。它们的时间复杂度O(MN)、空间复杂度O(MN)（在基于计算匹配字符串情况下，是不可优化的。如果只是计算LD和LCS，空间占用可以优化到O(M)）。

Nakatsu算法在计算匹配字符串的情况下，有着良好的时间复杂度O(N(M-P))和空间复杂度O(N²)，而且在采取适当的优化手段时，可以将空间复杂度优化到O(N)。

定义说明：

1. 设M=Len(A)，N=Len(B)，不失一般性，假设M≤N。

2. A=a₁a₂……a_M，表示A是由a₁a₂……a_M这M个字符组成

B=b₁b₂……b_N，表示B是由b₁b₂……b_N这N个字符组成

LCS(i,j)=LCS(a₁a₂……a_i,b₁b₂……b_j)，其中1≤i≤M，1≤j≤N

3. L(k,i)=Min｛j｝ Where LCS(i,j)=k 表示，所有与字符串a₁a₂……a_i有长度为k的LCS的字符串b₁b₂……b_j中j的最小值（这个可以看后面的例子好好理解）

为了推导L的计算，有下面几个定理。

　　定理一：任意的i，1≤i≤M。有L(1,i)＜L(2,i)＜L(3,i)……

　　定理二：任意的i，1≤i≤M-1。任意的k，1≤k≤M。有L(k,i+1)≤L(k,i)

　　定理三：任意的i，1≤i≤M-1。任意的k，1≤k≤M－1。有L(k,i)＜L(k+1,i+1)

　　定理四：如果L(k,i+1)存在，则L(k,i+1)的计算公式为

　　　　　　L(k,i+1)=Min｛Min｛j｝，L(k,i)｝ Where ｛a_i+1=b_j And j>L(k-1,i)｝

　　上面四个定理证明从略。可以从上面四个定理推导出L的计算。

　　故，L的计算公式为

　　　　L(1,1)=Min｛j｝ Where ｛a₁=b_j｝

　　　　L(1,i)=Min｛Min｛j｝ Where ｛a_i=b_j｝，L(1,i-1)｝ 　　此时，i＞1

　　　　L(k,i)=Min｛Min｛j｝ Where ｛a_i=b_j  And j＞L(k-1,i-1)｝，L(k,i-1)｝ 　　此时，i＞1，k＞1

　　　　注：以上公式中，若找不到满足Where后面条件的j，则j=MaxValue

　　　　　　当i＜k时，则L(k,i)=MaxValue

　　　　　　MaxValue是一个常量，表示“不存在”

　　举例说明：A=GGATCGA，B=GAATTCAGTTA，计算LCS(A,B)

　　第一步：初始化L矩阵，表格中V=MaxValue。

Nakatsu算法L矩阵

	i=1	i=2	i=3	i=4	i=5	i=6	i=7
k=1
k=2	V
k=3	V	V
k=4	V	V	V
k=5	V	V	V	V
k=6	V	V	V	V	V
k=7	V	V	V	V	V	V

　　第二步：依据上面的计算公式，计算表格的其余单元格

Nakatsu算法L矩阵

	i=1	i=2	i=3	i=4	i=5	i=6	i=7
k=1	1	1	1	1	1	1	1
k=2	V	8	2	2	2	2	2
k=3	V	V	11	4	4	4	3
k=4	V	V	V	V	6	6	6
k=5	V	V	V	V	V	8	7
k=6	V	V	V	V	V	V	11
k=7	V	V	V	V	V	V	V

　　第三步：在矩阵中找寻对角线

 　　　　1、先找如下的对角线，对角线中有四个单元格的值是V(MaxValue)。不是本算法的合适答案

Nakatsu算法L矩阵

	i=1	i=2	i=3	i=4	i=5	i=6	i=7
k=1	1	1	1	1	1	1	1
k=2	V	8	2	2	2	2	2
k=3	V	V	11	4	4	4	3
k=4	V	V	V	V	6	6	6
k=5	V	V	V	V	V	8	7
k=6	V	V	V	V	V	V	11
k=7	V	V	V	V	V	V	V

 　　　　2、再找右边的一条对角线。

Nakatsu算法L矩阵

	i=1	i=2	i=3	i=4	i=5	i=6	i=7
k=1	1	1	1	1	1	1	1
k=2	V	8	2	2	2	2	2
k=3	V	V	11	4	4	4	3
k=4	V	V	V	V	6	6	6
k=5	V	V	V	V	V	8	7
k=6	V	V	V	V	V	V	11
k=7	V	V	V	V	V	V	V

　　　　　　对角线上的所有单元格的值都不是V(MaxValue)。故本对角线就是算法的求解。

　　　　　　LCS(A,B)就是对角线的长度。故LCS(A,B)=6。

　　　　　　本算法的精妙之处就在于这六个单元格的值所对应的字符串B的字符就是最长公共子串。

　　　　　　最长公共子串：b₁b₂b₄b₆b₈b₁₁=GATCGA

　　　　　　再将最长公共子串在两个字符串中搜索一遍，能得出字符串的匹配字串。

　　　　　　　　A：GGA_TC_G__A

　　　　　　　　B：GAATTCAGTTA

　　　　　　　　注：原本以为能很容易得出匹配字符串。不过现在看来还需费一番周折，也是考虑不周。不过已经有大概的解决方案，留待后文介绍。

　　　　　　

　　

　　Nakatsu算法关键就是找寻满足条件对角线（对角线的值没有MaxValue）,故计算的过程可以沿着对角线进行，先计算第一条对角线，看是否满足对角线条件，满足则退出，不满足则继续计算下一条对角线，直到计算出满足条件的对角线。

　　假设LCS(A,B)=P，则一共需要计算M-P+1条对角线，每条对角线的比较次数为N，则Nakatsu算法的时间复杂度为O((M-P+1)N)，空间复杂度为O(M²)，但由于计算顺序的优化，可以将空间复杂度降为O(M)，这应该是令人满意的了。