最长公共子序列（LCS）

0、前言

    程序员编程艺术系列重新开始创作了（前十章，请参考程序员编程艺术第一~十章集锦与总结）。回顾之前的前十章，有些代码是值得商榷的，因当时的代码只顾阐述算法的原理或思想，所以，很多的与代码规范相关的问题都未能做到完美。日后，会着力修缮之。

    搜遍网上，讲解这个LCS问题的文章不计其数，但大多给读者一种并不友好的感觉，稍感晦涩，且代码也不够清晰。本文力图避免此些情况。力保通俗，阐述详尽。同时，经典算法研究系列的第三章（三、dynamic
programming）也论述了此LCS问题。有任何问题，欢迎不吝赐教。

第一节、问题描述

    什么是最长公共子序列呢?好比一个数列 S，如果分别是两个或多个已知数列的子序列，且是所有符合此条件序列中最长的，则S 称为已知序列的最长公共子序列。

    举个例子，如：有两条随机序列，如 1 3 4 5 5 ，and 2 4 5 5 7 6，则它们的最长公共子序列便是：4 5 5。

    注意最长公共子串（Longest CommonSubstring）和最长公共子序列（LongestCommon Subsequence, LCS）的区别：子串（Substring）是串的一个连续的部分，子序列（Subsequence）则是从不改变序列的顺序，而从序列中去掉任意的元素而获得的新序列；更简略地说，前者（子串）的字符的位置必须连续，后者（子序列LCS）则不必。比如字符串acdfg同akdfc的最长公共子串为df，而他们的最长公共子序列是adf。LCS可以使用动态规划法解决。下文具体描述。

第二节、LCS问题的解决思路

穷举法   

    解最长公共子序列问题时最容易想到的算法是穷举搜索法，即对X的每一个子序列，检查它是否也是Y的子序列，从而确定它是否为X和Y的公共子序列，并且在检查过程中选出最长的公共子序列。X和Y的所有子序列都检查过后即可求出X和Y的最长公共子序列。X的一个子序列相应于下标序列{1, 2, …, m}的一个子序列，因此，X共有2m个不同子序列（Y亦如此，如为2^n），从而穷举搜索法需要指数时间（2^m * 2^n）。
动态规划算法

    事实上，最长公共子序列问题也有最优子结构性质。

记:

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-1与Yn-1的LCS（LCS(X , Y)的长度等于LCS(Xm-1 , Yn-1)的长度加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-1与Y的LCS及X与Yn-1的LCS。LCS(X , Y)的长度为：max{LCS(Xm-1 , Y)的长度, LCS(X , Yn-1)的长度}。

    由于上述当xm≠yn的情况中，求LCS(Xm-1 , Y)的长度与LCS(X , Yn-1)的长度，这两个问题不是相互独立的：两者都需要求LCS(Xm-1，Yn-1)的长度。另外两个序列的LCS中包含了两个序列的前缀的LCS，故问题具有最优子结构性质考虑用动态规划法。

    也就是说，解决这个LCS问题，你要求三个方面的东西：1、LCS（Xm-1，Yn-1）+1；2、LCS（Xm-1，Y），LCS（X，Yn-1）；3、max{LCS（Xm-1，Y），LCS（X，Yn-1）}。

    行文至此，其实对这个LCS的动态规划解法已叙述殆尽，不过，为了成书的某种必要性，下面，我试着再多加详细阐述这个问题。

第三节、动态规划算法解LCS问题

3.1、最长公共子序列的结构

    最长公共子序列的结构有如下表示：

    设序列X=<x1, x2, …, xm>和Y=<y1, y2, …, yn>的一个最长公共子序列Z=<z1, z2, …, zk>，则：
若xm=yn，则zk=xm=yn且Zk-1是Xm-1和Yn-1的最长公共子序列；
若xm≠yn且zk≠xm ，则Z是Xm-1和Y的最长公共子序列；
若xm≠yn且zk≠yn ，则Z是X和Yn-1的最长公共子序列。

    其中Xm-1=<x1, x2, …, xm-1>，Yn-1=<y1, y2, …, yn-1>，Zk-1=<z1, z2, …, zk-1>。

3、2.子问题的递归结构

    由最长公共子序列问题的最优子结构性质可知，要找出X=<x1, x2, …, xm>和Y=<y1, y2, …, yn>的最长公共子序列，可按以下方式递归地进行：当xm=yn时，找出Xm-1和Yn-1的最长公共子序列，然后在其尾部加上xm(=yn)即可得X和Y的一个最长公共子序列。当xm≠yn时，必须解两个子问题，即找出Xm-1和Y的一个最长公共子序列及X和Yn-1的一个最长公共子序列。这两个公共子序列中较长者即为X和Y的一个最长公共子序列。

    由此递归结构容易看到最长公共子序列问题具有子问题重叠性质。例如，在计算X和Y的最长公共子序列时，可能要计算出X和Yn-1及Xm-1和Y的最长公共子序列。而这两个子问题都包含一个公共子问题，即计算Xm-1和Yn-1的最长公共子序列。

    与矩阵连乘积最优计算次序问题类似，我们来建立子问题的最优值的递归关系。用c[i,j]记录序列Xi和Yj的最长公共子序列的长度。其中Xi=<x1, x2, …, xi>，Yj=<y1, y2, …, yj>。当i=0或j=0时，空序列是Xi和Yj的最长公共子序列，故c[i,j]=0。其他情况下，由定理可建立递归关系如下：

3、3.计算最优值

    直接利用上节节末的递归式，我们将很容易就能写出一个计算c[i,j]的递归算法，但其计算时间是随输入长度指数增长的。由于在所考虑的子问题空间中，总共只有θ(m*n)个不同的子问题，因此，用动态规划算法自底向上地计算最优值能提高算法的效率。

    计算最长公共子序列长度的动态规划算法LCS_LENGTH(X,Y)以序列X=<x1, x2, …, xm>和Y=<y1, y2, …, yn>作为输入。输出两个数组c[0..m ,0..n]和b[1..m ,1..n]。其中c[i,j]存储Xi与Yj的最长公共子序列的长度，b[i,j]记录指示c[i,j]的值是由哪一个子问题的解达到的，这在构造最长公共子序列时要用到。最后，X和Y的最长公共子序列的长度记录于c[m,n]中。

Procedure LCS_LENGTH(X,Y);  

begin  

  m:=length[X];  

  n:=length[Y];  

  for i:=1 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  

        begin  

          c[i,j]:=c[i-1,j-1]+1;  

          b[i,j]:="↖";  

        end  

      else if c[i-1,j]≥c[i,j-1] then  

        begin  

          c[i,j]:=c[i-1,j];  

          b[i,j]:="↑";  

        end  

      else  

        begin  

          c[i,j]:=c[i,j-1];  

          b[i,j]:="←"  

        end;  

  return(c,b);  

end;   

    由算法LCS_LENGTH计算得到的数组b可用于快速构造序列X=<x1, x2, …, xm>和Y=<y1, y2, …, yn>的最长公共子序列。首先从b[m,n]开始，沿着其中的箭头所指的方向在数组b中搜索。
当b[i,j]中遇到"↖"时（意味着xi=yi是LCS的一个元素），表示Xi与Yj的最长公共子序列是由Xi-1与Yj-1的最长公共子序列在尾部加上xi得到的子序列；
当b[i,j]中遇到"↑"时，表示Xi与Yj的最长公共子序列和Xi-1与Yj的最长公共子序列相同；
当b[i,j]中遇到"←"时，表示Xi与Yj的最长公共子序列和Xi与Yj-1的最长公共子序列相同。

    这种方法是按照反序来找LCS的每一个元素的。由于每个数组单元的计算耗费Ο(1)时间，算法LCS_LENGTH耗时Ο(mn)。

3、4.构造最长公共子序列

    下面的算法LCS(b,X,i,j)实现根据b的内容打印出Xi与Yj的最长公共子序列。通过算法的调用LCS(b,X,length[X],length[Y])，便可打印出序列X和Y的最长公共子序列。

Procedure LCS(b,X,i,j);  

begin  

  if i=0 or j=0 then return;  

  if b[i,j]="↖" then  

    begin  

      LCS(b,X,i-1,j-1);  

      print(x[i]); {打印x[i]}  

    end  

  else if b[i,j]="↑" then LCS(b,X,i-1,j)   

                      else LCS(b,X,i,j-1);  

end;   

在算法LCS中，每一次的递归调用使i或j减1，因此算法的计算时间为O(m+n)。

例如，设所给的两个序列为X=<A，B，C，B，D，A，B>和Y=<B，D，C，A，B，A>。由算法LCS_LENGTH和LCS计算出的结果如下图所示：



    我来说明下此图（参考算法导论）。在序列X={A，B，C，B，D，A，B}和 Y={B，D，C，A，B，A}上，由LCS_LENGTH计算出的表c和b。第i行和第j列中的方块包含了c[i，j]的值以及指向b[i，j]的箭头。在c[7,6]的项4，表的右下角为X和Y的一个LCS<B，C，B，A>的长度。对于i，j>0，项c[i，j]仅依赖于是否有xi=yi，及项c[i-1，j]和c[i，j-1]的值，这几个项都在c[i，j]之前计算。为了重构一个LCS的元素，从右下角开始跟踪b[i，j]的箭头即可，这条路径标示为阴影，这条路径上的每一个“↖”对应于一个使xi=yi为一个LCS的成员的项（高亮标示）。

    所以根据上述图所示的结果，程序将最终输出：“B C B A”。

//
// main.cpp
// LCS
//
// Created by Bryan on 14-6-26.
// Copyright (c) 2014年 Bryan. All rights reserved.
//

#include <iostream>

using namespace std;

void Print_LCS(int **LCS_Dir,char *s,int row,int col)
{
if(s == 0)
return;
if(strlen(s) == 0 || row < 0 || col < 0)
return;
if(LCS_Dir[row][col] == 1)
{
if(row > 0 && col > 0)
Print_LCS(LCS_Dir, s, row - 1, col -1);
cout<<s[row]<<" ";
}
else if(LCS_Dir[row][col] == 2)
{
Print_LCS(LCS_Dir, s, row - 1, col);
}
else
Print_LCS(LCS_Dir, s, row, col - 1);

}

int LCS(char *s1,char *s2)
{
int len1 = strlen(s1);
int len2 = strlen(s2);
if(len1 == 0 || len2 == 0)
{
return 0;
}

//===============Init begin====================
int **lcs_Len = new int*[len1];
int **lcs_Dir = new int*[len1];

for(int i = 0;i<len1;i++)
{
lcs_Len[i] = new int[len2];
lcs_Dir[i] = new int[len2];
}

for(int i = 0;i<len1;i++)
{
for(int j = 0;j<len2;j++)
lcs_Dir[i][j] = 0;
}
for(int i = 0;i<len1;i++)
lcs_Len[i][0] = 0;
for(int j = 0;j<len2;j++)
lcs_Len[0][j] = 0;
//==============Init end===============

//========Find LCS_length start=======
for(int i = 0;i<len1;i++)
{
for(int j = 0;j<len2;j++)
{
if(i == 0 || j==0)
{
if (s1[i] == s2[j])
{
lcs_Len[i][j] = 1;
lcs_Dir[i][j] = 1;
}
}
else if(s1[i] == s2[j])
{
lcs_Len[i][j] = lcs_Len[i - 1][j - 1] + 1;
lcs_Dir[i][j] = 1;
}
else if(lcs_Len[i - 1][j] > lcs_Len[i][j - 1])
{
lcs_Len[i][j] = lcs_Len[i - 1][j];
lcs_Dir[i][j] = 2;
}
else
{
lcs_Len[i][j] = lcs_Len[i][j - 1];
lcs_Dir[i][j] = 3;
}

}
}

//========Find LCS_length end=======

Print_LCS(lcs_Dir,s1,len1 -1,len2 - 1);

int lcs_length = lcs_Len[len1 - 1][len2 - 1];
for(int i = 0;i<len1;i++)
{
delete [] lcs_Dir[i];
delete [] lcs_Len[i];
}

delete [] lcs_Dir;
delete [] lcs_Len;

return lcs_length;

}

int main(int argc, const char * argv[])
{

cout<<LCS("ABCBDAB", "BDCABA")<<endl;
// insert code here...
std::cout << "Hello, World!\n";
return 0;
}