文本相似度量方法(2): LCS 和编辑距离

对给定序列 A 和 B, 满足以下条件的一个序列 C 被称为 A 和 B 的公共子序列:

1. C 中每一个元素都对应 A 和 B 中一个元素
2. 从 C 中挑选两个元素 \(C_{i}\) 和 \(C_{j}\) ，其中 \(i\) 和 \(j\) 表示这两个元素在 C 中的序号(从左至右)，假设这两个元素分别对应 \(A_{m}\) 和 \(A_{n}\) ，那么有 \((j - i)\cdot(n - m)\gt 0\)，在 B 中对应的两个元素同理

比如说给定 \(A =\) "打南边来了个喇嘛，手里提拉着五斤鳎目", \(B =\) "打北边来了个哑巴，腰里别着个喇叭" ，以下都是 A 和 B 的子序列:

• "打边"
• "打边来了个"
• "边个着"

在所有可能的 C 中，长度最大的即所谓 "最长公共子序列"。上述例子中 A 和 B 的最长公共子序列是: "打边来了个，里着"。如下图所示:

直观感受上，我们可以认为，如果 A 和 B 的最长公共子序列越长，A 和 B 就越相似。这也是用最长公共子序列来度量文本相似程度的思想。

LCS 是一个比较经典的算法问题，了解动态规划的读者应该对它不会陌生。

要求解 LCS，朴素的想法是将 A 的所有子序列都枚举一遍，看是否是 B 的子序列，然后从中挑选出最长的。对给定字符串 A ，假定其长度为 \(L\)，其所有可能的子序列数量为 \(\sum_{i=1}^{L}\binom{L}{i}\) 也就是 \(2^{L} - 1\) ，所以暴力求解方法的复杂度至少是指数级的，这显然是不可取的。

通常我们用动态规划方法来求解 LCS 问题。

不难发现 LCS 的求解可以按照以下方法来得到:

1. 对比 A 和 B 的第一个字符，如果相等，则转入步骤 2，否则转入步骤 3
2. 将 A 和 B 的第一个字符记录为 LCS 的第一个字符，求 A 和 B 剩下部分的 LCS (转回步骤 1，下同)
3. 将 A 去掉第一个字符，用剩下的部分和 B 一起计算 LCS，转入步骤 4
4. 将 B 去掉第一个字符，用剩下的部分和 A 一起计算 LCS，转入步骤 5
5. 比较第 3 步和第 4 步得到的 LCS，其中较长的就是 A 和 B 的 LCS.

上述过程很明显是一个递归过程，所以可以简单地实现出来。

# coding: utf-8

def lcs(a, b):
    if not a or not b:
        return []

    res = []
    if a[0] == b[0]:
        res.append(a[0])
        res.extend(lcs(a[1:], b[1:]))
    else:
        lcs_one = lcs(a[1:], b)
        lcs_two = lcs(a, b[1:])
        res = lcs_one if len(lcs_one) > len(lcs_two) else lcs_two

    return res

res = lcs('hello world', 'hero word')
print ''.join(res)

上述代码可以解决问题，看起来行数也不多，但其实还是有问题的，那就是递归树的存在导致了大量的重复计算、以及可能存在的栈溢出风险。举个栗子，计算 "hello" 和 "hero" 的过程如下图所示:

可以看到 "lcs('lo', 'o')" 被计算了两次。

解决上面提到的问题的方法是将递归方法改为循环方法。用 \(m\) 表示 A 的长度，用 \(n\) 表示 B 的长度，我们可以构建一个 \(m\times n\) 的矩阵 \(D\) ，用来保存上面递归过程中计算到的 A 和 B 的公共子序列。

具体方法是， \(D_{i,j}\) 表示去除 A 的前 i 个元素后的子串 \(A^{\prime}\) 和 B 去掉前 j 个元素后的子串 \(B^{\prime}\) 的 LCS。这样我们首先计算 \(D_{m,n}\) ，然后计算 \(D_{m-1,n}\) 和 \(D_{m,n-1}\) ，再计算 \(D_{m-1,n-1}\)，依次类推。为什么可以这么做呢？仔细看看之前提到的递归过程中的 2、3、4 步，以及上面那张计算 "hello" 和 "hero" 的 LCS 的递归树图，可以发现计算最后都可以归结为计算 A 和 B 尾部片段的 LCS 。

我们发现上面的形式话表述不太直观，如果能先计算 \(D_{0,0}\) 是不是更好一些呢？是的，只要把之前提到的递归规则中的 "比较 A 和 B 的第一个字符" 改为 "比较 A 和 B 的最后一个字符" 并对步骤 2、3、4 做出相应的改变即可。

这样，LCS 的算法可以改写为:

# coding: utf-8

import numpy as np


def lcs(a, b):
    m, n = len(a), len(b)

    # 为便于计算，为 D 多增加一行一列
    # 其中第一行和第一列中的元素保持为空字符串
    D = np.zeros((m+1, n+1), dtype=object)
    D[:] = ''                   # 初始化为空字符串

    for idx_a, ch_a in enumerate(a, 1):
        for idx_b, ch_b in enumerate(b, 1):
            # 若 D 不增加一行一列，下标 idx_a-1/idx_b-1 要判断是否非负
            if ch_a == ch_b:
                D[idx_a, idx_b] = D[idx_a-1, idx_b-1] + ch_a
            else:
                lcs_one = D[idx_a, idx_b-1]
                lcs_two = D[idx_a-1, idx_b]
                if len(lcs_one) > len(lcs_two):
                    D[idx_a, idx_b] = lcs_one
                else:
                    D[idx_a, idx_b] = lcs_two

    return D[m, n]


print lcs('hello', 'hero')

当然实际上我们不会去用一个二维数组来保存计算过程中用到的(非最长)公共子序列，这样虽然很直观，但是在内存使用上有点丑陋。标准的做法是只记录这些公共子序列的长度，计算完整个长度矩阵后，再从最后的位置回溯取得 LCS 。

先观察一下计算 "lcs('hello', 'hero')" 时得到的公共子序列矩阵:

矩阵中出现过的公共子序列有: 'h', 'he' 和 'heo'。从中我们 似乎可以发现这么一个规律: 从上往下逐行查看，这三个公共子序列 第一次 出现的时候，恰好就是 'hello' 和 'hero' 中有字符相等的时候，换成记录长度后，也就对应某个特定的长度值第一次出现的时候。

就这个例子而言，我们 似乎 可以这样来根据长度矩阵得到 LCS (行/列序号从 0 开始，后同):

• 在第 1 行第 1 列找到 1 ，对应的字符是 'h'
• 在第 2 行第 2 列找到 2 ，对应的字符是 'e'
• 在第 5 行第 4 列找到 3 ，对应的字符是 'o'

这种方法 直观上感觉是对的，但实际上是有问题的 ，下面是计算 'GCGGACTG' 和 'GCCCTAGCG' 时得到的长度矩阵。

如果按照刚才的做法来反推 LCS ，会得到下面的结果:

1. 在第 1 行第 1 列找到 1 ，对应的字符是 'G'
2. 在第 2 行第 2 列找到 2 ，对应的字符是 'C'
3. 在第 3 行第 7 列找到 3 ，对应的字符是 'G'
4. 在第 4 行第 9 列找到 4 ，对应的字符是 'G'
5. 第 4 步找到的 'G' 是 'GCCCTAGCG' 的最后一个字符(表中最后一列)，因此停止

实际上真正求得的 LCS 长度应该是 5 ，为 'GCGCG'、'GCACG' 或 'GCCTG'，而不是 'GCGG'。问题出在哪呢？ LCS 可以看成是一个最优解，但到达最优解的路径可能有不止一条，而且局部的最优解并不一定是最优解的组成部分，所以前面提到的贪心方法在有些情况下可以得到正确的结果，但有的情况下就会出错。

正确的做法是从长度矩阵右下角，根据长度矩阵的计算规则往前反推，这样就能保证得到的结果是最长的公共子序列。

先把子序列长度矩阵的计算方法实现:

import numpy as np


def lcs_matrix(a, b):
    m, n = len(a), len(b)
    matrix = np.zeros((m+1, n+1), dtype=int)

    for idx_a, ch_a in enumerate(a, 1):
        for idx_b, ch_b in enumerate(b, 1):
            if ch_a == ch_b:
                matrix[idx_a, idx_b] = matrix[idx_a-1, idx_b-1] + 1
            else:
                matrix[idx_a, idx_b] = max(
                    matrix[idx_a, idx_b-1],
                    matrix[idx_a-1, idx_b]
                )

    return matrix


print lcs_matrix('GCGGACTG', 'GCCCTAGCG')

根据长度矩阵的计算规则，可以按照以下步骤来反推出 LCS:

1. 首先定位到长度矩阵右下角位置
2. 如果当前位置的值为 0 ，则停止；否则转到步骤 3
3. 如果当前位置对应的 A 和 B 的元素相等，则向当前位置的左上角后退一步(行号和列号各减 1)，并回到步骤 2，否则转到步骤 4
4. 检查矩阵当前位置左边的值和上边的值，跳转到其中值更大的那个位置(如果相等，则在往上和往左中选择一个方向)，回到步骤 2

用代码实现出来大概是这样:

import numpy as np


def lcs_backtrace(a, b, matrix):
    idx_a, idx_b = len(a) - 1, len(b) - 1

    lcs_list = []
    while matrix[idx_a+1, idx_b+1] > 0:
        if a[idx_a] == b[idx_b]:
            lcs_list.append(a[idx_a])
            idx_a -= 1
            idx_b -= 1
        else:
            upper_value = matrix[idx_a, idx_b+1]
            left_value = matrix[idx_a+1, idx_b]
            if upper_value > left_value:
                idx_a -= 1
            else:
                idx_b -= 1

    lcs_list.reverse()
    return lcs_list


a = 'GCGGACTG'
b = 'GCCCTAGCG'
matrix = np.array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
                   [0, 1, 1, 1, 1, 1, 1, 1, 1, 1],
                   [0, 1, 2, 2, 2, 2, 2, 2, 2, 2],
                   [0, 1, 2, 2, 2, 2, 2, 3, 3, 3],
                   [0, 1, 2, 2, 2, 2, 2, 3, 3, 4],
                   [0, 1, 2, 2, 2, 2, 3, 3, 3, 4],
                   [0, 1, 2, 3, 3, 3, 3, 3, 4, 4],
                   [0, 1, 2, 3, 3, 4, 4, 4, 4, 4],
                   [0, 1, 2, 3, 3, 4, 4, 5, 5, 5]], dtype=int)
print lcs_backtrace(a, b, matrix)

需要注意的是，LCS 结果可能并不唯一。

上述求解方法可以用数学语言表示如下:

\[\begin{eqnarray}LCS(A_{i}, B_{j}) = \begin{cases} 0 & if\quad i=0\quad or\quad j=0\\ LCS(A_{i-1}, B_{j-1})+a_{i} & if\quad a_{i}=b_{j}\\ longest(LCS(A_{i}, B_{j-1}), LCS(A_{i-1}, B_{j})) & if\quad a_{i}\neq b_{j} \end{cases}\end{eqnarray}\]

其中 \(A_{i}\) 表示 A 的前 \(i\) 个字符组成的字符串，\(B_{j}\) 同理。

直观上，我们可以认为 A 和 B 的 LCS 越长，那么 A 和 B 就越相似。为了使所有用于比较的 (A, B) 对得到的相似度量能进行横向比较，定义 LCS 相似度为:

\[S(A, B) = \frac{2\cdot |LCS(A, B)|}{|A| + |B|}\]

这样得到的相似度的值就被变换到 [0, 1] 区间中了。

编辑距离最早由俄罗斯科学家 Levenshtein 提出，故又称 "Levenshtein 距离"。其定义为: 给定两个字符串 A 和 B，将 A 通过删除、插入、替换操作转换为 B 所需要的最少操作次数。

比如将 "kitten" 转换为 "sitting" 需要进行如下操作:

1. 替换操作: kitten -> sitten (k -> s)
2. 替换操作: sitten -> sittin (e -> i)
3. 插入操作: sittin -> sitting (SPC -> g)

我们就说 "kitten" 相对 "sitting" 的编辑距离是 3。

编辑距离衡量的是两个字符串之间的差异程度，所以差异程度越小，相似程度就越大了。

编辑距离的核心过程 LCS 在我看来是一样的，可以看一下它的数学描述:

\[\begin{eqnarray}LEV(A_{i}, B_{j}) = \begin{cases} max(i, j) & if\quad min(i, j)=0\\ min(LEV(A_{i-1}, B_{j})+1, LEV(A_{i}, B_{j-1})+1, LEV(A_{i-1}, B_{j-1})) & if\quad a_{i}=b_{j}\\ min(LEV(A_{i-1}, B_{j})+1, LEV(A_{i}, B_{j-1})+1, LEV(A_{i-1}, B_{j-1})+1) & if\quad a_{i}\neq b_{j}\\ \end{cases}\end{eqnarray}\]

按照其数学描述，可以实现如下:

# coding: utf-8
import numpy as np


def edit_distance(a, b):
    m, n = len(a), len(b)
    dis_matrix = np.zeros((m+1, n+1), dtype=int)

    # 初始化距离矩阵的第 0 行和第 0 列
    dis_matrix[0, :] = np.arange(n+1)
    dis_matrix[:, 0] = np.arange(m+1)

    # 开始计算
    for idx_a, ch_a in enumerate(a, 1):
        for idx_b, ch_b in enumerate(b, 1):
            cur_dis = None

            dis_left = dis_matrix[idx_a, idx_b-1]
            dis_upper = dis_matrix[idx_a-1, idx_b]
            dis_upper_left = dis_matrix[idx_a-1, idx_b-1]
            if ch_a == ch_b:
                cur_dis = min(dis_left+1, dis_upper+1, dis_upper_left)
            else:
                cur_dis = min(dis_left+1, dis_upper+1, dis_upper_left + 1)

            dis_matrix[idx_a, idx_b] = cur_dis

    return dis_matrix[m, n]


print edit_distance('kitten', 'sitting')

与 LCS 不同的是，编辑距离衡量的是差异性，所以用编辑距离来表示相似程度，按照如下式子进行转换:

\[S(A, B)=1-\frac{LEV(A, B)}{max(|A|, |B|)}\]