分类目录:《算法设计与分析》总目录
相关文章:
· 动态规划(一):基础知识
· 动态规划(二):钢条切割
· 动态规划(三):矩阵链乘法
· 动态规划(四):动态规划详解
· 动态规划(五):最长公共子序列
在生物应用中,经常需要比较两个(或多个)不同生物体的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 npdef 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)