最长公共子序列问题(LCS) (生物信息学中常用算法)
子序列的概念:
设 X=< x1, x2,┅, xm>,若有 1≤i1< i2< ┅ <ik≤m,使得
Z=< z1, z2,┅, zk> = < xi1, xi2,┅, xik>,则称 Z 是 X 的子序列,
记为 Z<X。
则有 Z<X。
e.g. X=<A,B,C,B,D,A,B>, Z=<B,C,B,A>,
公共子序列的概念:
设 X,Y 是两个序列,且有 Z<X 和 Z<Y,
则称 Z 是 X 和 Y 的公共子序列。
最长公共子序列的概念:
若 Z<X,Z<Y,且不存在比 Z 更长的 X 和 Y 的公共子序列,
则称 Z 是 X 和 Y 的最长公共子序列,记为 Z∈LCS(X , Y)。
最长公共子序列往往不止一个。
Y=<B,D,C,A,B,A>, 则
e.g. X=<A,B,C,B,D,A,B>,
Z=<B,C,B,A>, Z’=<B,C,A,B>, Z’’=<B,D,A,B>
均属于 LCS(X , Y) ,即 X,Y 有 3 个 LCS。
如何找出 X 和 Y 的一个最长公共子序列?
Brute-force 法:列出 X 的所有长度不超过 n(即Y)的子序列,
从长到短逐一进行检查,看其是否为 Y 的子序列,
直到找到第一个最长公共子序列。
由于 X 共有 2m 个子序列,故此方法对较大的 m 没有实用价值。
是否能使用动态规划法?如何用?
分析:
记 Xi=<x1,...,xi>即 X 序列的前 i 个字符 (1≤i≤m)(前缀)
Yj=<y1,...,yj>即 Y 序列的前 j 个字符 (1≤j≤n)(前缀)
假定 Z=<z1,...,zk>∈LCS(X , Y)。
若 xm=yn(最后一个字符相同) ,则不难用反证法证明:
该字符必是 X 与 Y 的任一最长公共子序列 Z(设长度为 k)的
最后一个字符,即有 zk = xm = yn 且显然有 Zk-1∈LCS(Xm-1 , Yn-1)
即 Z 的前缀 Zk-1 是 Xm-1 与 Yn-1 的最长公共子序列。
若 xm≠yn,则亦不难用反证法证明:
要么 Z∈LCS(Xm-1, Y),要么 Z∈LCS(X , Yn-1)。
由于 zk≠xm 与 zk≠yn 其中至少有一个必成立,
若 zk≠xm 则有 Z∈LCS(Xm-1 , Y),类似的,
若 zk≠yn 则有 Z∈LCS(X , Yn-1)
∴若 xm=yn,则问题化归成求 Xm-1 与 Yn-1 的 LCS
(LCS(X , Y)的长度等于 LCS(Xm-1 , Yn-1)的长度加 1)
若 xm≠yn 则问题化归成求 Xm-1 与 Y 的 LCS 及 X 与 Yn-1 的 LCS
LCS(X , Y)的长度为:
max{LCS(Xm-1 , Y)的长度, LCS(X , Yn-1)的长度}
求 LCS(Xm-1 , Y)的长度与 LCS(X , Yn-1)的长度
这两个问题不是相互独立的:
∵两者都需要求 LCS(Xm-1,Yn-1)的长度,
因而具有重叠性。
另外两个序列的 LCS 中包含了两个序列的前缀的 LCS,
故问题具有最优子结构性质 ⇒ 考虑用动态规划法。
引进一个二维数组 C,用 C[i,j]记录 Xi 与 Yj 的 LCS 的长度
如果我们是自底向上进行递推计算,那么在计算 C[i,j]之前,
C[i-1,j-1], C[i-1,j]与 C[i,j-1]均已计算出来。此时我们
根据 X[i]=Y[j]还是 X[i]≠Y[j],就可以计算出 C[i,j]:
若 X[i]=Y[j],则执行 C[i,j]←C[i-1,j-1]+1;若 X[i]≠Y[j],则根据:
C[i-1,j]≥C[i,j-1],则 C[i,j]取 C[i-1,j];否则 C[i,j]取 C[i,j-1]。即有
0 若i = 0或j = 0
C[i − 1, j − 1] + 1 若i, j > 0且 xi = y j
C[i,j]=
y j
max{C[i − 1, j ], C[i, j − 1]} 若i, j > 0且 xi ≠
e.g. 如下图:
0 1 2 345 6
yj B DCABA
0 xi 0 0 00 00 0
↑ ↑ ↑↖ ↖
0 ←0 ←0←0
1 A 1←1 1
↖ ↑↖
1 ←1←1 ←1 2 ←2
2 B 0
↑ ↑↖ ↑↑
1 ←1 2 ←2←2 ←2
3 C 0
↖ ↑ ↑ ↑↖
1 ←1 2 ←2 3 ←3
4 B 0
↑↖ ↑ ↑↑ ↑
2←2 ←2 3 ←3
5 D 0 1
↑ ↑ ↑↖ ↑↖
6 A 0 1 2←2 3←3 4
↖ ↑ ↑ ↑↖ ↑
3 4 ←4
7 B 0 1 2←2
为了构造出 LCS,使用一个 m×n 的二维数组 b,
b[i,j]记录 C[i,j]是通过哪一个子问题的值求得的,
以决定搜索的方向:
若 C[i-1,j]≥C[i,j-1],则 b[i,j]中记入“↑” ;
若 C[i-1,j] < C[i,j-1],则 b[i,j]中记入“←” ;
算法 LCS_L(X,Y,m,n,C)
for i=0 to m do C[i,0]←0
for j=1 to n do C[0,j]←0
for i=1 to m do {
for j=1 to n do{
if x[i]=y[j] then {C[i,j]←C[i-1,j-1]+1;
b[i,j]←“↖” ;}
}else if C[i-1,j]≥C[i,j-1] then {C[i,j]←C[i-1,j];
b[i,j]←“↑” ;}
}else{C[i,j]←C[i,j-1];
b[i,j]←“←” ;}
算法的时间复杂度显然是Θ(m×n)的。
输出一个 LCS(X,Y)的递归算法:
LCS_Output(b,X,i,j)
If i=0 or j=0 then return;
If b[i,j]=“↖”then /*X[i]=Y[j]*/
{LCS_Output(b,X,i-1,j-1);
输出 X[i];}
else if b[i,j]=“↑”then /*C[i-1,j]≥C[i,j-1]*/
LCS_Output(b,X,i-1,j)
else if b[i,j]=“←” then /*C[i-1,j]<C[i,j-1]*/
LCS_Output(b,X,i,j-1)
e.g. 对上述例子调用 LCS_Output(b,X,7,6)。
算法分析:
由于每次调用至少向上或向左(或向上向左同时)移动一步,
故最多调用(m+n)次就会遇到 i=0 或 j=0 的情况,
此时开始返回。返回时与递归调用时方向相反,步数相同,
故算法时间复杂度为Θ(m+n)。
注意:仅用“↑” , “←” , “↖”是搜索不到所有的 LCS 的;
主要是由于在上述 LCS_L 算法里,
对于语句 if C[i-1] ≥ C[i,j-1] then ....,
没有区分 C[i-1,j]>C[i,j-1] 和 C[i-1,j]=C[i,j-1]
这两种不同的情况。因此,要找出所有的 LCS,
当满足 X[i]≠Y[j]且 C[i-1,j]=C[i,j-1]时,要执行 b[i,j]←“←↑” ,
即记录两个搜索方向,才能够据此找出所有的 LCS。
为节省空间,数组 b 亦可不用,直接
根据 X[i]=Y[j]还是 X[i]≠Y[j]以及 C[i,j-1], C[i-1,j]来找出搜索方向:
,
X[i]=Y[j]等价于 b[i,j]=“↖”
X[i]≠Y[j]且 C[i,j-1] > C[i-1,j] 等价于 b[i,j]=“←” ,
X[i]≠Y[j]且 C[i,j-1] < C[i-1,j] 等价于 b[i,j]=“↑” ,
X[i]≠Y[j]且 C[i,j-1] = C[i-1,j] 等价于 b[i,j]=“←↑” ,
不用 b 数组可节省Θ(m×n)的空间,但算法写起来稍烦一些。
执行时间的量级两者相同。
如果只需计算 LCS 的长度而无需求出具体的 LCS,
则二维数组 C 也可以不用,
只用一个长为 min{m,n}的数组即可。
作业:设计一个算法求出全部的 LCS,分析最坏情况。
用”会计方法”证明,利用 b[i,j]来求所有 LCS 的算法
m
在最坏情况下时间复杂度为 O( C )。m+ n
子序列的概念:
设 X=< x1, x2,┅, xm>,若有 1≤i1< i2< ┅ <ik≤m,使得
Z=< z1, z2,┅, zk> = < xi1, xi2,┅, xik>,则称 Z 是 X 的子序列,
记为 Z<X。
则有 Z<X。
e.g. X=<A,B,C,B,D,A,B>, Z=<B,C,B,A>,
公共子序列的概念:
设 X,Y 是两个序列,且有 Z<X 和 Z<Y,
则称 Z 是 X 和 Y 的公共子序列。
最长公共子序列的概念:
若 Z<X,Z<Y,且不存在比 Z 更长的 X 和 Y 的公共子序列,
则称 Z 是 X 和 Y 的最长公共子序列,记为 Z∈LCS(X , Y)。
最长公共子序列往往不止一个。
Y=<B,D,C,A,B,A>, 则
e.g. X=<A,B,C,B,D,A,B>,
Z=<B,C,B,A>, Z’=<B,C,A,B>, Z’’=<B,D,A,B>
均属于 LCS(X , Y) ,即 X,Y 有 3 个 LCS。
如何找出 X 和 Y 的一个最长公共子序列?
Brute-force 法:列出 X 的所有长度不超过 n(即Y)的子序列,
从长到短逐一进行检查,看其是否为 Y 的子序列,
直到找到第一个最长公共子序列。
由于 X 共有 2m 个子序列,故此方法对较大的 m 没有实用价值。
是否能使用动态规划法?如何用?
分析:
记 Xi=<x1,...,xi>即 X 序列的前 i 个字符 (1≤i≤m)(前缀)
Yj=<y1,...,yj>即 Y 序列的前 j 个字符 (1≤j≤n)(前缀)
假定 Z=<z1,...,zk>∈LCS(X , Y)。
若 xm=yn(最后一个字符相同) ,则不难用反证法证明:
该字符必是 X 与 Y 的任一最长公共子序列 Z(设长度为 k)的
最后一个字符,即有 zk = xm = yn 且显然有 Zk-1∈LCS(Xm-1 , Yn-1)
即 Z 的前缀 Zk-1 是 Xm-1 与 Yn-1 的最长公共子序列。
若 xm≠yn,则亦不难用反证法证明:
要么 Z∈LCS(Xm-1, Y),要么 Z∈LCS(X , Yn-1)。
由于 zk≠xm 与 zk≠yn 其中至少有一个必成立,
若 zk≠xm 则有 Z∈LCS(Xm-1 , Y),类似的,
若 zk≠yn 则有 Z∈LCS(X , Yn-1)
∴若 xm=yn,则问题化归成求 Xm-1 与 Yn-1 的 LCS
(LCS(X , Y)的长度等于 LCS(Xm-1 , Yn-1)的长度加 1)
若 xm≠yn 则问题化归成求 Xm-1 与 Y 的 LCS 及 X 与 Yn-1 的 LCS
LCS(X , Y)的长度为:
max{LCS(Xm-1 , Y)的长度, LCS(X , Yn-1)的长度}
求 LCS(Xm-1 , Y)的长度与 LCS(X , Yn-1)的长度
这两个问题不是相互独立的:
∵两者都需要求 LCS(Xm-1,Yn-1)的长度,
因而具有重叠性。
另外两个序列的 LCS 中包含了两个序列的前缀的 LCS,
故问题具有最优子结构性质 ⇒ 考虑用动态规划法。
引进一个二维数组 C,用 C[i,j]记录 Xi 与 Yj 的 LCS 的长度
如果我们是自底向上进行递推计算,那么在计算 C[i,j]之前,
C[i-1,j-1], C[i-1,j]与 C[i,j-1]均已计算出来。此时我们
根据 X[i]=Y[j]还是 X[i]≠Y[j],就可以计算出 C[i,j]:
若 X[i]=Y[j],则执行 C[i,j]←C[i-1,j-1]+1;若 X[i]≠Y[j],则根据:
C[i-1,j]≥C[i,j-1],则 C[i,j]取 C[i-1,j];否则 C[i,j]取 C[i,j-1]。即有
0 若i = 0或j = 0
C[i − 1, j − 1] + 1 若i, j > 0且 xi = y j
C[i,j]=
y j
max{C[i − 1, j ], C[i, j − 1]} 若i, j > 0且 xi ≠
e.g. 如下图:
0 1 2 345 6
yj B DCABA
0 xi 0 0 00 00 0
↑ ↑ ↑↖ ↖
0 ←0 ←0←0
1 A 1←1 1
↖ ↑↖
1 ←1←1 ←1 2 ←2
2 B 0
↑ ↑↖ ↑↑
1 ←1 2 ←2←2 ←2
3 C 0
↖ ↑ ↑ ↑↖
1 ←1 2 ←2 3 ←3
4 B 0
↑↖ ↑ ↑↑ ↑
2←2 ←2 3 ←3
5 D 0 1
↑ ↑ ↑↖ ↑↖
6 A 0 1 2←2 3←3 4
↖ ↑ ↑ ↑↖ ↑
3 4 ←4
7 B 0 1 2←2
为了构造出 LCS,使用一个 m×n 的二维数组 b,
b[i,j]记录 C[i,j]是通过哪一个子问题的值求得的,
以决定搜索的方向:
若 C[i-1,j]≥C[i,j-1],则 b[i,j]中记入“↑” ;
若 C[i-1,j] < C[i,j-1],则 b[i,j]中记入“←” ;
算法 LCS_L(X,Y,m,n,C)
for i=0 to m do C[i,0]←0
for j=1 to n do C[0,j]←0
for i=1 to m do {
for j=1 to n do{
if x[i]=y[j] then {C[i,j]←C[i-1,j-1]+1;
b[i,j]←“↖” ;}
}else if C[i-1,j]≥C[i,j-1] then {C[i,j]←C[i-1,j];
b[i,j]←“↑” ;}
}else{C[i,j]←C[i,j-1];
b[i,j]←“←” ;}
算法的时间复杂度显然是Θ(m×n)的。
输出一个 LCS(X,Y)的递归算法:
LCS_Output(b,X,i,j)
If i=0 or j=0 then return;
If b[i,j]=“↖”then /*X[i]=Y[j]*/
{LCS_Output(b,X,i-1,j-1);
输出 X[i];}
else if b[i,j]=“↑”then /*C[i-1,j]≥C[i,j-1]*/
LCS_Output(b,X,i-1,j)
else if b[i,j]=“←” then /*C[i-1,j]<C[i,j-1]*/
LCS_Output(b,X,i,j-1)
e.g. 对上述例子调用 LCS_Output(b,X,7,6)。
算法分析:
由于每次调用至少向上或向左(或向上向左同时)移动一步,
故最多调用(m+n)次就会遇到 i=0 或 j=0 的情况,
此时开始返回。返回时与递归调用时方向相反,步数相同,
故算法时间复杂度为Θ(m+n)。
注意:仅用“↑” , “←” , “↖”是搜索不到所有的 LCS 的;
主要是由于在上述 LCS_L 算法里,
对于语句 if C[i-1] ≥ C[i,j-1] then ....,
没有区分 C[i-1,j]>C[i,j-1] 和 C[i-1,j]=C[i,j-1]
这两种不同的情况。因此,要找出所有的 LCS,
当满足 X[i]≠Y[j]且 C[i-1,j]=C[i,j-1]时,要执行 b[i,j]←“←↑” ,
即记录两个搜索方向,才能够据此找出所有的 LCS。
为节省空间,数组 b 亦可不用,直接
根据 X[i]=Y[j]还是 X[i]≠Y[j]以及 C[i,j-1], C[i-1,j]来找出搜索方向:
,
X[i]=Y[j]等价于 b[i,j]=“↖”
X[i]≠Y[j]且 C[i,j-1] > C[i-1,j] 等价于 b[i,j]=“←” ,
X[i]≠Y[j]且 C[i,j-1] < C[i-1,j] 等价于 b[i,j]=“↑” ,
X[i]≠Y[j]且 C[i,j-1] = C[i-1,j] 等价于 b[i,j]=“←↑” ,
不用 b 数组可节省Θ(m×n)的空间,但算法写起来稍烦一些。
执行时间的量级两者相同。
如果只需计算 LCS 的长度而无需求出具体的 LCS,
则二维数组 C 也可以不用,
只用一个长为 min{m,n}的数组即可。
作业:设计一个算法求出全部的 LCS,分析最坏情况。
用”会计方法”证明,利用 b[i,j]来求所有 LCS 的算法
m
在最坏情况下时间复杂度为 O( C )。m+ n