算法设计与分析——动态规划（五）：最长公共子序列

分类目录：《算法设计与分析》总目录
相关文章：
· 动态规划（一）：基础知识
· 动态规划（二）：钢条切割
· 动态规划（三）：矩阵链乘法
· 动态规划（四）：动态规划详解
· 动态规划（五）：最长公共子序列

---

在生物应用中，经常需要比较两个（或多个）不同生物体的DNA。一个DNA串由一串称为碱基的分子组成，碱基有腺嘌呤、鸟嘌呤、胞嘧啶和胸腺嘧啶4种类型。我们用英文单词首字母表示4种碱基，这样就可以将一个DNA串表示为有限集$A,C,G,T$上的一个字符串。例如，某种生物的DNA可能为$S_1=ACCGGTCGAGTGCGCGGAAGCCGGCCGAA$，另一种生物的DNA可能为$S_2=GTCGTTCGGAATGCCGTTGCTCTGTAAA$。我们比较两个DNA串的一个原因是希望确定它们的“相似度”，作为度量两种生物相近程度的指标。我们可以用很多不同的方式来定义相似度，实际上也确实已经出现了很多相似度的定义。例如，如果一个DNA串是另一个DNA串的子串，那么可以说它们是相似的。但在我们的例子中，$S_1$和$S_2$都不是对方的子串。我们还可以这样来定义相似性：如果将一个串转换为另一个串所需的操作很少，那么可以说两个串是相似的。另一种衡量$S_1$和$S_2$的相似度的方式是：寻找第三个串$S_3$，它的所有碱基也都出现在$S_1$和$S_2$中，且在三个串中出现的顺序都相同，但在$S_1$和$S_2$中不要求连续出现。可以找到的$S_3$越长，就可以认为$S_1$和$S_2$的相似度越高。在我们的例子中，最长的$S_3$为 $GTCGTCGGAAGCCGGCCGAA$。

我们将最后一种相似度的概念命名为最长公共子序列问题。一个给定序列的子序列，就是将给定序列中零个或多个元素去掉之后得到的结果。其形式化定义如下：给定一个序列$X=<x_1,x_2,\cdots ,x_n>$，另一个序列$Z=<z_1,z_2,\cdots ,z_m>$满足如下条件时称为$X$的子序列，即存在一个严格递增的$X$的下标序列$<i_1,i_2,\cdots ,i_k>$，对所有$j=1,2,\cdots ,k$满足$X_{i_k}=Z_j$。例如，$Z=<B,C,E>$是$X=<A,B,C,D,E,F,G>$的子序列，对应的下标序列为$2,3,5$。

给定两个序列$X$和$Y$，如果$Z$既是$X$的子序列，也是$Y$的子序列，我们称它是$X$和$Y$的公共子序列。例如，如果$X=<A,B,C,D,E,F,G>$，$X=<B,B,D,C,D,E,G>$，那么序列$Z=<B,C,>$就是$X$和$Y$的公共子序列。但它不是$X$和$Y$的最长公共子序列（LCS），因为它长度为2，而$<B,C,E>$也是$X$和$Y$的公共子序列，其长度为3。$<B,C,D,E,G>$是$X$和$Y$的最长公共子序列。最长公共子序列问题（ longest-common-subsequence problem）说的是给定两个序列$X=<x_1,x_2,\cdots ,x_m>$和$Y=<x_1,x_2,\cdots ,x_n>$，求$X$和$Y$长度最长的公共子序列。下面我们讲如何用动态规划方法高效地求解LCS问题：

步骤1：刻画最长公共子序列的特征

如果用暴力搜索方法求解LCS问题，就要穷举$X$的所有子序列，对每个子序列检查它是否也是$Y$的子序列，记录找到的最长子序列。$X$的每个子序列对应$X$的下标集合$1,2,\cdots ,m$的一个子集，所以$X$有$2^m$个子序列，因此暴力方法的运行时间为指数阶，对较长的序列是不实用的。

根据定义，我们可以很直观地得到最长公共子序列的最优质子结构性质：

令$X=<x_1,x_2,\cdots ,x_m>$和$Y=<x_1,x_2,\cdots ,x_n>$为两个序列，$Z=<z_1,z_2,\cdots ,z_k>$为$X$和$Y$的任意LCS。

1. 如果$X_m=Y_n$，则$Z_k=X_m=Y_n$且$Z_{k-1}$是$X_{m-1}$和$Y_{n-1}$的一个LCS。
2. 如果$X_m\ne Y_n$，且$Z_k\ne X_m$意味着$Z$是$X_{m-1}$和$Y$的一个LCS。
3. 如果$X_m\ne Y_n$，且$Z_k\ne Y_n$意味着$Z$是$X$和$Y_{n-1}$的一个LCS。

如上文所示，LCS问题具有最优子结构性质。我们将看到，子问题的自然分类对应两个输入序列的“前缀”对。前缀的严谨定义如下：给定一个序列$X=<x_1,x_2,\cdots ,x_m>$，对$i=1,2,\cdots ,m$，定义$X$的第i前缀为$X=<x_1,x_2,\cdots ,x_i>$。例如，若$X=<A,B,C,D,E,F,G>$，则$X_3=<A,B,C>$，$X_0$为空串。所以，两个序列的LCS包含两个序列的前缀的LCS。因此，LCS问题具有最优子结构性质。

步骤2：一个递归解

有步骤1可以得出，在求$X=<x_1,x_2,\cdots ,x_m>$和$Y=<x_1,x_2,\cdots ,x_n>$的一个LCS时，我们需要求解一个或两个子问题。如果$X_m=Y_n$，我们应该求解$X_{m-1}$和$Y_{n-1}$的一个LCS。将$X_m=Y_n$追加到这个LCS的末尾，就得到X和Y的一个LCS。如果$X_m\ne Y_n$，我们必须求解两个子问题：求$X_{m-1}$和$Y$的一个LCS与$X$和$Y_{n-1}$的一个LCS。两个LCS较长者即为$X$和$Y$的一个LCS。由于这些情况覆盖了所有可能性，因此我们知道必然有一个子问题的最优解出现在$X$和$Y$的LCS中。

我们可以很容易看出LCS问题的重叠子问题性质。为了求$X$和$Y$的一个LCS，我们可能需要求$X_{m-1}$和$Y$的一个LCS与$X$和$Y_{n-1}$的一个LCS。但是这几个子问题都包含求解$X_{m-1}$和$Y_{n-1}$的ICS的子子问题。很多其他子问题也都共享子子问题。

与矩阵链乘法问题相似，设计LCS问题的递归算法首先要建立最优解的递归式。我们定义$c[i,j]$表示$X_i$和$Y_j$的LCS的长度。如果$i=0$或$j=0$，即一个序列长度为0，那么LCS的长度为0。根据ICS问题的最优子结构性质，可得如下公式：

观察到在递归公式中，我们通过限制条件限定了需要求解哪些子问题。当$X_i=Y_j$时，我们可以而且应该求解子问题：$X_i-1$和$Y_{j-1}$的一个LCS。否则，应该求解两个子问题：$X_i$和$Y_{j-1}$的一个LCS及$X_i-1$和$Y_j$的一个LCS。在之前讨论过的钢条切割问题和矩阵链乘法问题的动态规划算法中，根据问题的条件，我们没有排除任何子问题。

步骤3：计算LCS的长度

根据上面的分析，我们可以很容易地写出一个指数时间的递归算法来计算两个序列的LCS的长度。但是，由于LCS问题只有$mn$个不同的子问题，我们可以用动态规划方法自底向上地计算。

过程lce_length(X, Y)接受两个序列$X=[x_1,x_2,\cdots ,x_m]$和$Y=[x_1,x_2,\cdots ,x_n]$为输入。它将$c[i,j]$的值保存在表$c[0\cdots m,0\cdots n]$中，并按行主次序计算表项（即首先由左至右计算c的第一行，然后计算第二行，依此类推）。过程还维护一个表$b[i,j]$，帮助构造最优解。$b[i,j]$指向的表项对应计算$c[i,j]$时所选择的子问题最优解。过程返回表$b$和表$c$。

import numpy as np
def lcs_length(X, Y):
    m = len(X)
    n = len(Y)
    b = np.zeros([m + 1, n + 1]) 
    c = np.zeros([m + 1, n + 1])
    for i in range(1, m + 1):
        for j in range(1, n + 1):
            if X[i - 1] == Y[j - 1]:
                c[i, j] = c[i - 1, j - 1] + 1
                b[i, j] = '1'
            elif c[i - 1, j] >= c[i, j -1]:
                c[i, j] = c[i - 1, j]
                b[i, j] = '2'
            else:
                c[i, j] = c[i, j -1]
                b[i, j] = '3'
    return c, b

下图显示了ce_length(X, Y)对输入X = 'ABCBDAB'和Y = 'BDCABA'，过程的运行时间为$O(mn)$。

其中，b中的$1$表示$\nwarrow$、$2$表示$\uparrow$、$3$表示$\leftarrow$。

步骤4：构造LCS

我们可以用print_lcs(X, b, m, n)返回的表b快速构造$X=<x_1,x_2,\cdots ,x_m>$和$Y=<x_1,x_2,\cdots ,x_n>$的LCS，只需简单地从$b[m,n]$开始，并按箭头方向追踪下去即可。

def print_lcs(X, b, m, n):
    if m == 0 or n == 0:
        return ''
    if b[m, n] == 1:
        print_lcs(X, b, m - 1, n - 1)
        print(X[m - 1])
    elif b[m, n] == 2:
        print_lcs(X, b, m - 1, n)
    else:
        print_lcs(X, b, m, n - 1)