最长公共子序列问题

      最长公共子序列问题(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