分类目录:《算法设计与分析》总目录
相关文章:
· 动态规划(一):基础知识
· 动态规划(二):钢条切割
· 动态规划(三):矩阵链乘法
· 动态规划(四):动态规划详解
· 动态规划(五):最长公共子序列
在生物应用中,经常需要比较两个(或多个)不同生物体的DNA。一个DNA串由一串称为碱基的分子组成,碱基有腺嘌呤、鸟嘌呤、胞嘧啶和胸腺嘧啶4种类型。我们用英文单词首字母表示4种碱基,这样就可以将一个DNA串表示为有限集 A , C , G , T {A,C,G,T} A,C,G,T上的一个字符串。例如,某种生物的DNA可能为 S 1 = A C C G G T C G A G T G C G C G G A A G C C G G C C G A A S_1= ACCGGTCGAGTGCGCGGAAGCCGGCCGAA S1=ACCGGTCGAGTGCGCGGAAGCCGGCCGAA,另一种生物的DNA可能为 S 2 = G T C G T T C G G A A T G C C G T T G C T C T G T A A A S_2= GTCGTTCGGAATGCCGTTGCTCTGTAAA S2=GTCGTTCGGAATGCCGTTGCTCTGTAAA。我们比较两个DNA串的一个原因是希望确定它们的“相似度”,作为度量两种生物相近程度的指标。我们可以用很多不同的方式来定义相似度,实际上也确实已经出现了很多相似度的定义。例如,如果一个DNA串是另一个DNA串的子串,那么可以说它们是相似的。但在我们的例子中, S 1 S_1 S1和 S 2 S_2 S2都不是对方的子串。我们还可以这样来定义相似性:如果将一个串转换为另一个串所需的操作很少,那么可以说两个串是相似的。另一种衡量 S 1 S_1 S1和 S 2 S_2 S2的相似度的方式是:寻找第三个串 S 3 S_3 S3,它的所有碱基也都出现在 S 1 S_1 S1和 S 2 S_2 S2中,且在三个串中出现的顺序都相同,但在 S 1 S_1 S1和 S 2 S_2 S2中不要求连续出现。可以找到的 S 3 S_3 S3越长,就可以认为 S 1 S_1 S1和 S 2 S_2 S2的相似度越高。在我们的例子中,最长的 S 3 S_3 S3为 G T C G T C G G A A G C C G G C C G A A GTCGTCGGAAGCCGGCCGAA GTCGTCGGAAGCCGGCCGAA。
我们将最后一种相似度的概念命名为最长公共子序列问题。一个给定序列的子序列,就是将给定序列中零个或多个元素去掉之后得到的结果。其形式化定义如下:给定一个序列 X = < x 1 , x 2 , ⋯ , x n > X=<x_1, x_2, \cdots, x_n> X=<x1,x2,⋯,xn>,另一个序列 Z = < z 1 , z 2 , ⋯ , z m > Z=<z_1, z_2, \cdots, z_m> Z=<z1,z2,⋯,zm>满足如下条件时称为 X X X的子序列,即存在一个严格递增的 X X X的下标序列 < i 1 , i 2 , ⋯ , i k > <i_1, i_2, \cdots, i_k> <i1,i2,⋯,ik>,对所有 j = 1 , 2 , ⋯ , k j=1, 2, \cdots, k j=1,2,⋯,k满足 X i k = Z j X_{i_k}=Z_j Xik=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 , 5 2, 3, 5 2,3,5。
给定两个序列 X X X和 Y Y Y,如果 Z Z Z既是 X X X的子序列,也是 Y Y Y的子序列,我们称它是 X X X和 Y Y Y的公共子序列。例如,如果 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,>就是 X X X和 Y Y Y的公共子序列。但它不是 X X X和 Y Y Y的最长公共子序列(LCS),因为它长度为2,而 < B , C , E > <B, C, E> <B,C,E>也是 X X X和 Y Y Y的公共子序列,其长度为3。 < B , C , D , E , G > <B, C, D, E, G> <B,C,D,E,G>是 X X X和 Y Y Y的最长公共子序列。最长公共子序列问题( longest-common-subsequence problem)说的是给定两个序列 X = < x 1 , x 2 , ⋯ , x m > X=<x_1, x_2, \cdots, x_m> X=<x1,x2,⋯,xm>和 Y = < x 1 , x 2 , ⋯ , x n > Y=<x_1, x_2, \cdots, x_n> Y=<x1,x2,⋯,xn>,求 X X X和 Y Y Y长度最长的公共子序列。下面我们讲如何用动态规划方法高效地求解LCS问题:
步骤1:刻画最长公共子序列的特征
如果用暴力搜索方法求解LCS问题,就要穷举 X X X的所有子序列,对每个子序列检查它是否也是 Y Y Y的子序列,记录找到的最长子序列。 X X X的每个子序列对应 X X X的下标集合 1 , 2 , ⋯ , m {1, 2, \cdots, m} 1,2,⋯,m的一个子集,所以 X X X有 2 m 2^m 2m个子序列,因此暴力方法的运行时间为指数阶,对较长的序列是不实用的。
根据定义,我们可以很直观地得到最长公共子序列的最优质子结构性质:
令 X = < x 1 , x 2 , ⋯ , x m > X=<x_1, x_2, \cdots, x_m> X=<x1,x2,⋯,xm>和 Y = < x 1 , x 2 , ⋯ , x n > Y=<x_1, x_2, \cdots, x_n> Y=<x1,x2,⋯,xn>为两个序列, Z = < z 1 , z 2 , ⋯ , z k > Z=<z_1, z_2, \cdots, z_k> Z=<z1,z2,⋯,zk>为 X X X和 Y Y Y的任意LCS。
- 如果 X m = Y n X_m=Y_n Xm=Yn,则 Z k = X m = Y n Z_k=X_m=Y_n Zk=Xm=Yn且 Z k − 1 Z_{k-1} Zk−1是 X m − 1 X_{m-1} Xm−1和 Y n − 1 Y_{n-1} Yn−1的一个LCS。
- 如果 X m ≠ Y n X_m≠Y_n Xm=Yn,且 Z k ≠ X m Z_k≠X_m Zk=Xm意味着 Z Z Z是 X m − 1 X_{m-1} Xm−1和 Y Y Y的一个LCS。
- 如果 X m ≠ Y n X_m≠Y_n Xm=Yn,且 Z k ≠ Y n Z_k≠Y_n Zk=Yn意味着 Z Z Z是 X X X和 Y n − 1 Y_{n-1} Yn−1的一个LCS。
如上文所示,LCS问题具有最优子结构性质。我们将看到,子问题的自然分类对应两个输入序列的“前缀”对。前缀的严谨定义如下:给定一个序列 X = < x 1 , x 2 , ⋯ , x m > X=<x_1, x_2, \cdots, x_m> X=<x1,x2,⋯,xm>,对 i = 1 , 2 , ⋯ , m i = 1, 2, \cdots, m i=1,2,⋯,m,定义 X X X的第i前缀为 X = < x 1 , x 2 , ⋯ , x i > 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>,则 X 3 = < A , B , C > X_3=<A, B, C> X3=<A,B,C>, X 0 X_0 X0为空串。所以,两个序列的LCS包含两个序列的前缀的LCS。因此,LCS问题具有最优子结构性质。
步骤2:一个递归解
有步骤1可以得出,在求 X = < x 1 , x 2 , ⋯ , x m > X=<x_1, x_2, \cdots, x_m> X=<x1,x2,⋯,xm>和 Y = < x 1 , x 2 , ⋯ , x n > Y=<x_1, x_2, \cdots, x_n> Y=<x1,x2,⋯,xn>的一个LCS时,我们需要求解一个或两个子问题。如果 X m = Y n X_m=Y_n Xm=Yn,我们应该求解 X m − 1 X_{m-1} Xm−1和 Y n − 1 Y_{n-1} Yn−1的一个LCS。将 X m = Y n X_m=Y_n Xm=Yn追加到这个LCS的末尾,就得到X和Y的一个LCS。如果 X m ≠ Y n X_m≠Y_n Xm=Yn,我们必须求解两个子问题:求 X m − 1 X_{m-1} Xm−1和 Y Y Y的一个LCS与 X X X和 Y n − 1 Y_{n-1} Yn−1的一个LCS。两个LCS较长者即为 X X X和 Y Y Y的一个LCS。由于这些情况覆盖了所有可能性,因此我们知道必然有一个子问题的最优解出现在 X X X和 Y Y Y的LCS中。
我们可以很容易看出LCS问题的重叠子问题性质。为了求 X X X和 Y Y Y的一个LCS,我们可能需要求 X m − 1 X_{m-1} Xm−1和 Y Y Y的一个LCS与 X X X和 Y n − 1 Y_{n-1} Yn−1的一个LCS。但是这几个子问题都包含求解 X m − 1 X_{m-1} Xm−1和 Y n − 1 Y_{n-1} Yn−1的ICS的子子问题。很多其他子问题也都共享子子问题。
与矩阵链乘法问题相似,设计LCS问题的递归算法首先要建立最优解的递归式。我们定义
c
[
i
,
j
]
c[i, j]
c[i,j]表示
X
i
X_i
Xi和
Y
j
Y_j
Yj的LCS的长度。如果
i
=
0
i=0
i=0或
j
=
0
j=0
j=0,即一个序列长度为0,那么LCS的长度为0。根据ICS问题的最优子结构性质,可得如下公式:
观察到在递归公式中,我们通过限制条件限定了需要求解哪些子问题。当
X
i
=
Y
j
X_i=Y_j
Xi=Yj时,我们可以而且应该求解子问题:
X
i
−
1
X_i-1
Xi−1和
Y
j
−
1
Y_{j-1}
Yj−1的一个LCS。否则,应该求解两个子问题:
X
i
X_i
Xi和
Y
j
−
1
Y_{j-1}
Yj−1的一个LCS及
X
i
−
1
X_i-1
Xi−1和
Y
j
Y_j
Yj的一个LCS。在之前讨论过的钢条切割问题和矩阵链乘法问题的动态规划算法中,根据问题的条件,我们没有排除任何子问题。
步骤3:计算LCS的长度
根据上面的分析,我们可以很容易地写出一个指数时间的递归算法来计算两个序列的LCS的长度。但是,由于LCS问题只有 m n mn mn个不同的子问题,我们可以用动态规划方法自底向上地计算。
过程lce_length(X, Y)
接受两个序列
X
=
[
x
1
,
x
2
,
⋯
,
x
m
]
X=[x_1, x_2, \cdots, x_m]
X=[x1,x2,⋯,xm]和
Y
=
[
x
1
,
x
2
,
⋯
,
x
n
]
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]时所选择的子问题最优解。过程返回表
b
b
b和表
c
c
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
(
m
n
)
O(mn)
O(mn)。
其中,b
中的
1
1
1表示
↖
↖
↖、
2
2
2表示
↑
↑
↑、
3
3
3表示
←
←
←。
步骤4:构造LCS
我们可以用print_lcs(X, b, m, n)
返回的表b
快速构造
X
=
<
x
1
,
x
2
,
⋯
,
x
m
>
X=<x_1, x_2, \cdots, x_m>
X=<x1,x2,⋯,xm>和
Y
=
<
x
1
,
x
2
,
⋯
,
x
n
>
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)