Alan Lee

理解最长公共子序列算法

2020/04/07 Share

最近看一篇 Google 的论文:《Encode, Tag, Realize: High-Precision Text Editing》,看源码的时候发现其前期预处理的时候用了最长公共子序列(Longest Common Subsequence,LCS)算法来生成词汇表。之前只是知道这个算法,但是具体解法并不是很清楚。如果你经常刷题,可能会很熟悉了,但是我对刷题一直保持排斥态度,所以不太懂。趁这个机会来看看这个 LCS。

定义

根据 Wikipedia

The longest common subsequence (LCS) problem is the problem of finding the longest subsequence common to all sequences in a set of sequences (often just two sequences).

也就是说寻找多个序列(通常是两个)之间的最长的公共子序列,输入是多个序列,输出是一个子序列。还有一个和这个很相似的问题:最长公共子串(Longest Common Substring),和 LCS 的区别就在于子序列不要求连续,而子串是要求连续的,即 subsequence 和 substring 的区别。例如 ABCBDABBDCABA,最长公共子序列是 BCBA,就是不连续的。

解法

其解法大致有两种:穷举法和动态规划法。

其实还有一种方法是递归法,动态规划法是对其的一种改进(也使用了递归),就暂且不说这种方法了。

以下均以寻找两个字符串之间的 LCS 为例。

穷举

穷举是最直接的方法,也是最耗时的。其思想是:先找到两个字符串中较短的字符串,找出其所有子序列,然后依次查看每个子序列是否在较长字符串中,最后得出一个公共子序列列表,进而得到最长的公共子序列。

主函数如下:

为了尽量保持文章简洁,其他辅助函数未列出,完整代码见文章底部 GitHub 地址。

1
2
3
4
5
6
7
8
9
10
11
12
13
def lcs_brute_force(s1, s2):
if len(s1) < len(s2):
short = s1
long_ = s2
else:
short = s2
long_ = s1
cs = []
for sub in get_all_subsequences(short):
if is_subsequence(sub, long_):
cs.append(sub)
max_len = len(max(cs, key=len))
return (s for s in cs if len(s) == max_len)

一个长度为 $n$ 的字符串,共有 $2^n$ 个子序列:

假设字符串 $X$ 的长度为 $m$,字符串 $Y$ 的长度为 $n$,且 $m > n$,那么寻找其 LCS 就需要先找到 $Y$ 的所有子序列,然后去 $X$ 中查。显然时间复杂度为 $O(2^n)$,指数级的复杂度,这很明显不是一个适合大规模应用的算法。

动态规划

动态规划法(Dynamic Programming,DP)要含蓄得多,并不直接大海捞针式地穷举所有可能,而是根据将这个大问题分解为很多小问题,逐个解决小问题从而得到大问题的解。那么在这里,「小问题」就指的是诸多子串之间的 LCS

总体上来说,动态规划法的思路分两步走:构建和重建。构建是要构建一个表 $T$,记录 $X_{:i}$ 和 $Y_{:j}$ 之间的 LCS 长度。重建是要根据表 $T$ 去重建恢复具体的 LCS,即具体的字符。

构建

构建部分是基于以下几点的:

  • 若 $X_i = Y_j = t$,则 $t$ 是 $X_{:i}$ 和 $Y_{:j}$ 之间的 LCS 的最后一个字符。
  • 若 $X_i \neq Y_j$,则寻找 $X_{:(i-1)}$ 和 $Y_{:j}$,或者 $X_{:i}$ 和 $Y_{:(j-1)}$ 之间的 LCS,这两者之间的最长者,即为 $X$ 和 $Y$ 的 LCS。

可以看到这其实是一个递归的思路,每次都寻找 $X$ 和 $Y$ 中相等元素加入 LCS,LCS 的长度 +1。这一步结束之后便会生成一个 $(m+1) \times (n+1)$ 的表 $T$,这个表的第一行和第一列是初始值,均为 $0$,其他值记录的是相应子串之间的 LCS 长度。

代码实现如下(来自 lasertagger):

1
2
3
4
5
6
7
8
9
10
11
12
def lcs_table(source, target):
"""构建表 T。"""
rows = len(source)
cols = len(target)
lcs_table = [[0] * (cols + 1) for _ in range(rows + 1)]
for i in range(1, rows + 1):
for j in range(1, cols + 1):
if source[i - 1] == target[j - 1]:
lcs_table[i][j] = lcs_table[i - 1][j - 1] + 1
else:
lcs_table[i][j] = max(lcs_table[i - 1][j], lcs_table[i][j - 1])
return lcs_table

为了更直观的展示构建过程,我制作了一个 GIF 来说明:

source 为 $X$,target 为 $Y$。

construct table

重建

也叫回溯,指从表 $T$ 的右下角开始向上回溯,重建整个 LCS。

记 $\text{LCS}_{X, Y}$ 为 $X$ 和 $Y$ 的之间的 LCS。

回溯时:

  • 和构建时相同,若 $X_{i-1} = Y_{j-1} = t$,则将 $t$ 加入当前 LSC 列表,然后继续从 $T[i-1, j-1]$ 回溯,即向左上方走。
  • 若 $X_i \neq Y_j$:
    • 若 $T[i, j-1] > T[i-1, j]$,即 $\text{LCS}_{X_{i}, Y_{j-1}} > \text{LCS}_{X_{i-1}, Y_{j}}$,说明 $\text{LCS}_{X, Y}$ 肯定是包括 $\text{LCS}_{X_{i}, Y_{j-1}}$ 的,所以从 $T[i, j-1]$ 处继续回溯。
    • 否则,从 $T[i-1, j]$ 处开始回溯。

代码实现如下(改自 lasertagger):

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
def _backtrack(table, source, target, i, j):
"""回溯 LCS 表来重建 LCS。Backtracks the Longest Common Subsequence table to reconstruct the LCS.

Args:
table: LCS 表,注意该表的大小为 (len(source) + 1) × (len(target) + 1)。
source: 源字符串,list 形式。
target: 目标字符串,list 形式。
i: 当前行的索引。
j: 当前列的索引。

Returns:
LCS 列表。
"""
if i == 0 or j == 0:
# 如果i或者j为0,则说明有一方已经是空字符串,此时直接返回空
return []
if source[i - 1] == target[j - 1]: # 为什么不是 source[i] 和 target[j] 呢?实际上就是,因为 table 表的大小是要比实际 source 和 target 长度大 1 的,所以要减 1,不然就 index out of range 了。
# Append the aligned token to output.
# 如果 source[i - 1] 和 target[j - 1] 相等,则表示该值便是当前的 source 和 target 的 LCS 的最后一个 token,
# 然后 source 和 target 分别往前推一个 token 再去寻找,也是一个递归的过程
return _backtrack(table, source, target, i - 1, j - 1) + [target[j - 1]]
if table[i][j - 1] > table[i - 1][j]:
# 如果不相等的话,接下来就是寻找相等的元素,这些元素才可能属于 LCS。
# 如果不等且 source[:i] 和 target[:j-1] 的 LCS 长度,大于 source[:i-1] 和 target[:j] 的 LCS 长度,
# 那么就从 source[:i] 和 target[:j-1] 寻找相等的元素
return _backtrack(table, source, target, i, j - 1)
else:
# 和上面同理
return _backtrack(table, source, target, i - 1, j)

同样为了直观理解,我制作了一个 GIF:

source 为 $X$,target 为 $Y$。

reconstruct

LCS 并不一定唯一

以上图中的两个字符串为例,ABCBDABBDCABA 之间的 LCS 并不唯一,共有 ['BDAB', 'BCAB', 'BCBA'] 三个,穷举法可以都找到。然而上述的动态规划法却只能找到其中一个。

为什么?

两个数的大小存在三种情况,然而在判断 $T[i, j-1]$ 和 $T[i-1, j]$ 的大小时,上述算法只是比较前者是否大于后者,而相等和小于的情况被同等对待。这就造成了不能找到全部 LCS 的情况。

若 $T[i, j-1] = T[i-1, j]$,说明 $\text{LCS}_{X_{i}, Y_{j-1}}$ 和 $\text{LCS}_{X_{i-1}, Y_{j}}$ 都有可能是 $\text{LCS}_{X, Y}$ 的一部分,因此需要分两条路走,然后合并两条路得到的结果,这样就可以得到全部的 LCS。

构建部分不用改动,重建部分我用一个简单的树结构来存储 LCS 字符,对上述结构进行了改造:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
class Tree:
"""一个简单的树结构。"""
def __init__(self, data=None, left=None, right=None):
self.data = data
self.left = left
self.right = right


def backtrack(table, source, target, i, j):
"""改造的回溯函数。"""
tree = Tree()
if i == 0 or j == 0:
return Tree()
if source[i - 1] == target[j - 1]:
tree.data = target[j - 1]
tree.left = backtrack(table, source, target, i - 1, j - 1)
return Tree(
data=target[j - 1], left=backtrack(table, source, target, i - 1, j - 1)
)
if table[i][j - 1] > table[i - 1][j]:
return Tree(left=backtrack(table, source, target, i, j - 1))
elif table[i][j - 1] == table[i - 1][j]:
left = backtrack(table, source, target, i, j - 1)
right = backtrack(table, source, target, i - 1, j)
return Tree(left=left, right=right)
else:
return Tree(right=backtrack(table, source, target, i - 1, j))

重建完这个「树」之后,对齐进行类似 NLR(前序)的遍历即可得到全部 LCS:

由于节点默认的 dataNone,因此遍历结束后需要去除这些 None,详细可见 GitHub 完整代码。

1
2
3
4
5
6
7
8
def traverse_tree(root):
if root is None:
return []
if root.left == None and root.right == None:
return [str(root.data)]
return [
str(root.data) + p for p in traverse_tree(root.left) + traverse_tree(root.right)
]

完整代码

见 GitHub

Reference

END

CATALOG
  1. 1. 定义
  2. 2. 解法
    1. 2.1. 穷举
    2. 2.2. 动态规划
      1. 2.2.1. 构建
      2. 2.2.2. 重建
  3. 3. LCS 并不一定唯一
  4. 4. 完整代码
  5. 5. Reference
  6. 6. END