最长公共子序列（LCS）问题

转自博客：http://blog.chinaunix.net/uid-26548237-id-3374211.html

一、什么是最长公共子序列

   

    什么是最长公共子序列呢？举个简单的例子吧，一个数列S，若分别是两个或多个已知序列的子序列，且是所有符合条件序列中最长的，则S称为已知序列的最长公共子序列。

  举例如下，如：有两个随机数列，1 2 3 4 5 6 和 3 4 5 8 9，则它们的最长公共子序列便是：3 4 5。

  一直不明白：最长公共子串和最长公共子序列的区别。

  

   上网查了下，最长公共子串（Longest Common Substirng）和最长公共子序列（Longest Common Subsequence，LCS）的区别为：子串是串的一个连续的部分，子序列则是从不改变序列的顺序，而从序列中去掉任意的元素而获得新的序列；也就是说，子串中字符的位置必须是连续的，子序列则可以不必连续。

二、蛮力法

   蛮力法是解决最长公共子序列问题最容易想到的方法，即对S的每一个子序列，检查是否为T的子序列，从而确定它是否为S和T的公共子序列，并且选出最长的公共子序列。

 

   S和T的所有子序列都检查过后即可求出S和T的最长公共子序列。S的一个子序列相应于下标序列1,2，...，n的一个子序列。因此，S共有2^n个子序列。当然，T也有2^m个子序列。

   因此，蛮力法的时间复杂度为O(2^n * 2^m)，这可是指数级别的啊。

三、动态规划方法

   1、序列str1和序列str2

 

  ·长度分别为m和n；

  ·创建1个二维数组L[m.n]；

    ·初始化L数组内容为0

    ·m和n分别从0开始，m++，n++循环：

       - 如果str1[m] == str2[n]，则L[m,n] = L[m - 1, n -1] + 1；

       - 如果str1[m] != str2[n]，则L[m,n] = max{L[m,n - 1]，L[m - 1, n]}

    ·最后从L[m,n]中的数字一定是最大的，且这个数字就是最长公共子序列的长度

    ·从数组L中找出一个最长的公共子序列

   2、从数组L中查找一个最长的公共子序列

   i和j分别从m，n开始，递减循环直到i = 0，j = 0。其中，m和n分别为两个串的长度。

  ·如果str1[i] == str2[j]，则将str[i]字符插入到子序列内，i--，j--；

  ·如果str1[i] != str[j]，则比较L[i,j-1]与L[i-1,j]，L[i,j-1]大，则j--，否则i--；（如果相等，则任选一个）

图1 效果演示图

   

   根据上图，我们可以得到其中公共子串：B C B A 和 B D A B。

   总感觉，上面这个过程说的不是很清楚，但是不知道怎么才能更加清楚的表述？？纠结啊。

四、代码实现

1.   //代码实现比较简单，有可能不符合规矩，如有哪位前辈看到后，可以指出，我会虚心学习。
2.   1 #include <iostream>
   
3.   2 #include <string>
   
4.   3 using namespace std;
   
5.   4 int main(int argc, char **argv)
   
6.   5 {
   
7.   6 string str1 = "ABCBDAB";
   
8.   7 string str2 = "BDCABA";
   
9.   8
   
10.   9 int x_len = str1.length();
    
11.  10 int y_len = str2.length();
    
12.  11
    
13.  12 int arr[50][50] = {{0,0}};
    
14.  13
    
15.  14 int i = 0;
    
16.  15 int j = 0;
    
17.  16
    
18.  17 for(i = 1; i <= x_len; i++)
    
19.  18 {
    
20.  19 for(j = 1; j <= y_len; j++)
    
21.  20 {
    
22.  21 if(str1[i - 1] == str2[j - 1])
    
23.  22 {
    
24.  23 arr[i][j] = arr[i - 1][j - 1] + 1;
    
25.  24 }
    
26.  25 else
    
27.  26 {
    
28.  27
    
29.  28 if(arr[i][j - 1] >= arr[i - 1][j])
    
30.  29 {
    
31.  30 arr[i][j] = arr[i][j - 1];
    
32.  31 }
    
33.  32 else
    
34.  33 {
    
35.  34 arr[i][j] = arr[i -1][j];
    
36.  35 }
    
37.  36 }
    
38.  37
    
39.  38 }
    
40.  39 }
41. 
    
42. 41     for(i = 0 ; i <= x_len; i++)
43. 42     {
44. 43         for( j = 0; j <= y_len; j++)
45. 44         {
46. 45             cout << arr[i][j] << "  ";
47. 46         }
48. 47         cout << endl;
49. 48     }
50. 49     for(i = x_len, j = y_len; i >= 1 && j >= 1;)
51. 50     {
52. 51             if(str1[i - 1] == str2[j - 1])
53. 52             {
54. 53                 cout << str1[i - 1] << " ";//倒序打印的
55. 54                 i--;
56. 55                 j--;
57. 56             }
58. 57             else
59. 58             {
60. 59             //  if(arr[i][j -1] >= arr[i - 1][j])//打印：B A D B
61. 60                 if(arr[i][j -1] > arr[i - 1][j]) //打印:A B C B
62. 61                 {
63. 62                     j--;
64. 63                 }
65. 64                 else
66. 65                 {
67. 66                     i--;
68. 67                 }
69. 68             }
70. 69     }
71. 70     cout << endl;
72. 71     return 0;
73. 72 }
74.  

    运行结果如下所示。

图2 运行效果

   

   最后输出为A B C B，则最大子串为B C B A。

   其实，应该将结果保存起来，然后，正序打印呢。

 

    备注：

   我想还应该有更加优化的方法：求取最大公共子序列。关门了，明天看看。

   如果哪位前辈知道哪里有更加详细的介绍动态规划求最大公共子序列的资料，可以告诉我下，不胜感激。

   如果看到这些文字，对于文字中，有哪些不合适合理的地方，请定要指出，也不胜感激，这样我才可以更好的进步。

        梦醒潇湘

                  2012/10/15  22:42

*******************************更新于10月25日***********************************

   下面的文字来自于网络，最后给出引用位置。

   

一、动态规划算法

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

   记：

   Xi = <x1,x2,x3,....xi> 即X序列的前i个字符（1<= i <= m）(前缀)

   Yj = <y1,y2,y3,....yi> 即Y序列的前j个字符（1<= j <= m）(前缀)

   

   假定Z = <z1,z2,z3,...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问题

2.1 最长公共子序列的结构

 

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

   

   设序列X=<x1, x2, …, xm>和Y=<y1, y2, …, yn>的一个最长公共子序列Z=<z1, z2, …, zk>，则：

   1> 若 xm=yn，则 zk=xm=yn，且Zk-1是Xm-1和Yn-1的最长公共子序列； 

   2> 若 xm≠yn且 zk≠xm ，则 Z是 Xm-1和 Y的最长公共子序列；

   3> 若 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>。

2.2 子问题的递归结构

   由最长公共子序列问题的最优子结构性质可知，要找出 Xm=<x1, x2, …, xm>和 Yn=<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。其他情况下，可得递归关系如下所示：

   

2.3 计算最优值

 

   直接利用上节节末的递归式，我们将很容易就能写出一个计算c[i,j]的递归算法，但其计算时间是随输入长度指数增长的。由于在所考虑的子问题空间中，总共只有O(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]中。

   伪代码如下所示。

1. Procedure LCS_LENGTH(X,Y);LCS_LENGTH(X,Y);
2. begin 
3.    m:=length[X];
4.    n:=length[Y];
5. for i:=1 to m do c[i,0]:=0;
6. for j:=1 to n do c[0,j]:=0;
7. for i:=1 to m do 
8.    for j:=1 to n do 
9.       if x[i]=y[j] then 
10.         begin
11.            c[i,j]:= c[i-1,j -1]+ 1; 
12.            b[i,j]:="↖"; 
13.         end 
14.       else if c[i -1,j]≥ c[i,j -1] then
15.         begin
16.            c[i,j]:= c[i-1,j];
17.            b[i,j]:= "↑" ; 
18.         end 
19.       else 
20.         begin 
21.            c[i,j]:= c[j-1];
22.            b[i,j]:="←" 
23.         end;
24.    return(c,b);
25. 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的每一个元素的。由于每个数组单元的计算花费O(1)时间，算法LCS_Length耗时O(mn)。

2.4 构造最长公共子序列

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

   伪代码如下所示。

1. Procedure LCS(b,X,i,j);
   
2. begin
   
3.   if i=0 or j=0 then return;
   
4.   if b[i,j]=”↖” then
   
5.     begin
   
6.       LCS(b,X,i-1,j-1);
   
7.       print(x[i]); {打印x[i]}
   
8.     end
   
9.   else if b[i,j]=”↑” then LCS(b,X,i-1,j)
   
10.                       else LCS(b,X,i,j-1);
    
11. 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”。

   

   上述文字转自：http://blog.csdn.net/v_july_v/article/details/6695482

   让我对动态规划求LCS问题，有了更加深刻的理解。感谢作者！  

   

  

   根据以上思路，写出代码如下所示。

1. //只能打印一个最长公共子序列
2. #include <iostream>
   
3. using namespace std;
   
4.  
   
5.  const int X = 100, Y = 100;        //串的最大长度
   
6.  char result[X+1];                    //用于保存结果
   
7.  int count=0;                        //用于保存公共最长公共子串的个数
   
8. 
   
9.  /*功能：计算最优值
   
10.   *参数：
    
11.   *        x:字符串x
    
12.   *        y:字符串y
    
13.   *        b:标志数组
    
14.   *        xlen:字符串x的长度
    
15.   *        ylen:字符串y的长度
    
16.   *返回值：最长公共子序列的长度
    
17.   *
    
18.   */
    
19.  int Lcs_Length(string x, string y, int b[][Y+1],int xlen,int ylen)
    
20.  {
    
21.      int i = 0;
    
22.      int j = 0;
    
23. 
    
24.      int c[X+1][Y+1];
    
25.      for (i = 0; i<=xlen; i++)
    
26.      {
    
27.          c[i][0]=0;
    
28.      }
    
29.      for (i = 0; i <= ylen; i++ )
    
30.      {
    
31.          c[0][i]=0;
    
32.      }
    
33.      for (i = 1; i <= xlen; i++)
    
34.      {
    
35.         
    
36.          for (j = 1; j <= ylen; j++)
    
37.          {
    
38.              if (x[i - 1] == y[j - 1])
    
39.              {
    
40.                  c[i][j] = c[i-1][j-1]+1;
    
41.                  b[i][j] = 1;
    
42.              }
    
43.              else
    
44.                  if (c[i-1][j] > c[i][j-1])
    
45.                  {
    
46.                      c[i][j] = c[i-1][j];
    
47.                      b[i][j] = 2;
    
48.                  }
    
49.                  else
    
50.                      if(c[i-1][j] <= c[i][j-1])
    
51.                      {
    
52.                          c[i][j] = c[i][j-1];
    
53.                          b[i][j] = 3;
    
54.                      }
    
55.                      
    
56.          }
    
57.      }
    
58.      cout << "计算最优值效果图如下所示：" << endl;
    
59.      for(i = 1; i <= xlen; i++)
    
60.      {
    
61.          for(j = 1; j < ylen; j++)
    
62.          {
    
63.              cout << c[i][j] << " ";
    
64.          }
    
65.          cout << endl;
    
66.      }
    
67.      return c[xlen][ylen];
    
68.  }
    
69.  
    
70. void Display_Lcs(int i, int j, string x, int b[][Y+1],int current_Len)
    
71.  {
    
72.      if (i ==0 || j==0)
    
73.      {
    
74.          return;
    
75.      }
    
76.      if(b[i][j]== 1)
    
77.      {
    
78.          current_Len--;
    
79.          result[current_Len]=x[i- 1];
    
80.          Display_Lcs(i-1, j-1, x, b, current_Len);
    
81.      }
    
82.      else
    
83.      {
    
84.          if(b[i][j] == 2)
    
85.          {
    
86.              Display_Lcs(i-1, j, x, b, current_Len);
    
87.          }
    
88.          else
    
89.          {
    
90.              if(b[i][j]==3)
    
91.              {
    
92.                  Display_Lcs(i, j-1, x, b, current_Len);
    
93.              }
    
94.              else
    
95.              {
    
96.                  Display_Lcs(i-1,j,x,b, current_Len);
    
97.              }
    
98.          }
    
99.      }
    
100.  }
     
101.  
     
102.  int main(int argc, char* argv[])
     
103.  {
     
104.      string x = "ABCBDAB";
     
105.      string y = "BDCABA";
     
106.      int xlen = x.length();
     
107.      int ylen = y.length();
     
108. 
     
109.      int b[X + 1][Y + 1];
     
110. 
     
111.      int lcs_max_len = Lcs_Length( x, y, b, xlen,ylen );
     
112.      cout << lcs_max_len << endl;
     
113. 
     
114.      Display_Lcs( xlen, ylen, x, b, lcs_max_len );
     
115.      
     
116.      //打印结果如下所示
     
117.     for(int i = 0; i < lcs_max_len; i++)
     
118.     {
     
119.         cout << result[i];
     
120.     }
     
121.     cout << endl;
     
122.      return 0;
     
123.  }

    运行结果如下所示。

   由于有时并不是只有一个最长公共子序列，所以，对上面的代码进行改进，增加一个数组保存结果等....代码如下所示。

1. //求取所有的最长公共子序列
2. #include <iostream>
   
3. using namespace std;
   
4.  
   
5.  const int X = 100, Y = 100;        //串的最大长度
   
6.  char result[X+1];                    //用于保存结果
   
7.  int count=0;                        //用于保存公共最长公共子串的个数
   
8. 
   
9.  /*功能：计算最优值
   
10.   *参数：
    
11.   *        x:字符串x
    
12.   *        y:字符串y
    
13.   *        b:标志数组
    
14.   *        xlen:字符串x的长度
    
15.   *        ylen:字符串y的长度
    
16.   *返回值：最长公共子序列的长度
    
17.   *
    
18.   */
    
19.  int Lcs_Length(string x, string y, int b[][Y+1],int xlen,int ylen)
    
20.  {
    
21.      int i = 0;
    
22.      int j = 0;
    
23. 
    
24.      int c[X+1][Y+1];
    
25.      for (i = 0; i<=xlen; i++)
    
26.      {
    
27.          c[i][0]=0;
    
28.      }
    
29.      for (i = 0; i <= ylen; i++ )
    
30.      {
    
31.          c[0][i]=0;
    
32.      }
    
33.      for (i = 1; i <= xlen; i++)
    
34.      {
    
35.         
    
36.          for (j = 1; j <= ylen; j++)
    
37.          {
    
38.              if (x[i - 1] == y[j - 1])
    
39.              {
    
40.                  c[i][j] = c[i-1][j-1]+1;
    
41.                  b[i][j] = 1;
    
42.              }
    
43.              else
    
44.                  if (c[i-1][j] > c[i][j-1])
    
45.                  {
    
46.                      c[i][j] = c[i-1][j];
    
47.                      b[i][j] = 2;
    
48.                  }
    
49.                  else
    
50.                      if(c[i-1][j] < c[i][j-1])
    
51.                      {
    
52.                          c[i][j] = c[i][j-1];
    
53.                          b[i][j] = 3;
    
54.                      }
    
55.                      else
    
56.                      {
    
57.                          c[i][j] = c[i][j-1]; //或者c[i][j]=c[i-1][j];
    
58.                          b[i][j] = 4;
    
59.                      }
    
60.          }
    
61.      }
    
62.      cout << "计算最优值效果图如下所示：" << endl;
    
63.      for(i = 1; i <= xlen; i++)
    
64.      {
    
65.          for(j = 1; j < ylen; j++)
    
66.          {
    
67.              cout << c[i][j] << " ";
    
68.          }
    
69.          cout << endl;
    
70.      }
    
71.      return c[xlen][ylen];
    
72.  }
    
73.  
    
74.  /*功能：计算最长公共子序列
    
75.   *参数：
    
76.   *        xlen:字符串x的长度
    
77.   *        ylen:字符串y的长度
    
78.   *        x    :字符串x
    
79.   *        b:标志数组
    
80.   *        current_len:当前长度
    
81.   *        lcs_max_len:最长公共子序列长度
    
82.   *
    
83.   */
    
84.  void Display_Lcs(int i, int j, string x, int b[][Y+1],int current_len,int lcs_max_len)
    
85.  {
    
86.      if (i ==0 || j==0)
    
87.      {
    
88.          for(int s=0; s < lcs_max_len; s++)
    
89.          {
    
90.              cout << result[s];
    
91.          }
    
92.          cout<<endl;
    
93.          count++;
    
94.          return;
    
95.      }
    
96.      if(b[i][j]== 1)
    
97.      {
    
98.          current_len--;
    
99.          result[current_len]=x[i- 1];
    
100.          Display_Lcs(i-1, j-1, x, b,current_len,lcs_max_len);
     
101.      }
     
102.      else
     
103.      {
     
104.          if(b[i][j] == 2)
     
105.          {
     
106.              Display_Lcs(i-1, j, x, b,current_len,lcs_max_len);
     
107.          }
     
108.          else
     
109.          {
     
110.              if(b[i][j]==3)
     
111.              {
     
112.                  Display_Lcs(i, j-1, x, b,current_len,lcs_max_len);
     
113.              }
     
114.              else
     
115.              {
     
116.                  Display_Lcs(i,j-1,x,b,current_len,lcs_max_len);
     
117.                  Display_Lcs(i-1,j,x,b,current_len,lcs_max_len);
     
118.              }
     
119.          }
     
120.      }
     
121.  }
     
122.  
     
123.  int main(int argc, char* argv[])
     
124.  {
     
125.      string x = "ABCBDAB";
     
126.      string y = "BDCABA";
     
127.      int xlen = x.length();
     
128.      int ylen = y.length();
     
129. 
     
130.      int b[X + 1][Y + 1];
     
131. 
     
132.      int lcs_max_len = Lcs_Length( x, y, b, xlen,ylen );
     
133.      cout << lcs_max_len << endl;
     
134. 
     
135.      Display_Lcs( xlen, ylen, x, b, lcs_max_len, lcs_max_len );
     
136.      cout << "共有：" << count << "种";
     
137. 
     
138.   
     
139.      return 0;
     
140.  }

   运行结果如下图所示。