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

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

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

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

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

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

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

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

令X=<x1,x2,⋯,xm>X=<x_1, x_2, \cdots, x_m>X=<x1,x2,⋯,xm>和Y=<x1,x2,⋯,xn>Y=<x_1, x_2, \cdots, x_n>Y=<x1,x2,⋯,xn>为两个序列，Z=<z1,z2,⋯,zk>Z=<z_1, z_2, \cdots, z_k>Z=<z1,z2,⋯,zk>为XXX和YYY的任意LCS。

如果Xm=YnX_m=Y_nXm=Yn，则Zk=Xm=YnZ_k=X_m=Y_nZk=Xm=Yn且Zk−1Z_{k-1}Zk−1是Xm−1X_{m-1}Xm−1和Yn−1Y_{n-1}Yn−1的一个LCS。

如果Xm≠YnX_m≠Y_nXm=Yn，且Zk≠XmZ_k≠X_mZk=Xm意味着ZZZ是Xm−1X_{m-1}Xm−1和YYY的一个LCS。

如果Xm≠YnX_m≠Y_nXm=Yn，且Zk≠YnZ_k≠Y_nZk=Yn意味着ZZZ是XXX和Yn−1Y_{n-1}Yn−1的一个LCS。

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

步骤2：一个递归解

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

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

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

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

步骤3：计算LCS的长度

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

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

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] + 1b[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)O(mn)O(mn)。

其中，b中的111表示↖↖↖、222表示↑↑↑、333表示←←←。

步骤4：构造LCS

我们可以用print_lcs(X, b, m, n)返回的表b快速构造X=<x1,x2,⋯,xm>X=<x_1, x_2, \cdots, x_m>X=<x1,x2,⋯,xm>和Y=<x1,x2,⋯,xn>Y=<x_1, x_2, \cdots, x_n>Y=<x1,x2,⋯,xn>的LCS，只需简单地从b[m,n]b[m, n]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)