LCS公共子串问题求解

本文转自csdn

0、前言

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

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

第一节、问题描述

??? 什么是最长公共子序列呢?好比一个数列?

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]中。

1. Procedure?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[i,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的每一个元素的。由于每个数组单元的计算耗费Ο(1)时间，算法LCS_LENGTH耗时Ο(mn)。

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

??? 下面的算法LCS(b,X,i,j)实现根据b的内容打印出Xi与Yj的最长公共子序列。通过算法的调用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”。

3、5.算法的改进

??? 对于一个具体问题，按照一般的算法设计策略设计出的算法，往往在算法的时间和空间需求上还可以改进。这种改进，通常是利用具体问题的一些特殊性。

??? 例如，在算法LCS_LENGTH和LCS中，可进一步将数组b省去。事实上，数组元素c[i,j]的值仅由c[i-1,j-1]，c[i-1,j]和c[i,j-1]三个值之一确定，而数组元素b[i,j]也只是用来指示c[i,j]究竟由哪个值确定。因此，在算法LCS中，我们可以不借助于数组b而借助于数组c本身临时判断c[i,j]的值是由c[i-1,j-1]，c[i-1,j]和c[i,j-1]中哪一个数值元素所确定，代价是Ο(1)时间。既然b对于算法LCS不是必要的，那么算法LCS_LENGTH便不必保存它。这一来，可节省θ(mn)的空间，而LCS_LENGTH和LCS所需要的时间分别仍然是Ο(mn)和Ο(m+n)。不过，由于数组c仍需要Ο(mn)的空间，因此这里所作的改进，只是在空间复杂性的常数因子上的改进。

??? 另外，如果只需要计算最长公共子序列的长度，则算法的空间需求还可大大减少。事实上，在计算c[i,j]时，只用到数组c的第i行和第i-1行。因此，只要用2行的数组空间就可以计算出最长公共子序列的长度。更进一步的分析还可将空间需求减至min(m, n)。

第四节、编码实现LCS问题

??? 动态规划的一个计算最长公共子序列的方法如下，以两个序列?

1. ???
2. import?java.util.Random;??
3. ???
4. public?class?LCS{??
5. ????public?static?void?main(String[]?args){??
6. ???
7. ????????//设置字符串长度??
8. ????????int?substringLength1?=?20;??
9. ????????int?substringLength2?=?20;??//具体大小可自行设置??
10. ???
11. ????????//?随机生成字符串??
12. ????????String?x?=?GetRandomStrings(substringLength1);??
13. ????????String?y?=?GetRandomStrings(substringLength2);??
14. ???
15. ????????Long?startTime?=?System.nanoTime();??
16. ????????//?构造二维数组记录子问题x[i]和y[i]的LCS的长度??
17. ????????int[][]?opt?=?new?int[substringLength1?+?1][substringLength2?+?1];??
18. ???
19. ????????//?动态规划计算所有子问题??
20. ????????for?(int?i?=?substringLength1?-?1;?i?>=?0;?i--){??
21. ????????????for?(int?j?=?substringLength2?-?1;?j?>=?0;?j--){??
22. ????????????????if?(x.charAt(i)?==?y.charAt(j))??
23. ????????????????????opt[i][j]?=?opt[i?+?1][j?+?1]?+?1;?????????????????????????????????//参考上文我给的公式。??
24. ????????????????else??
25. ????????????????????opt[i][j]?=?Math.max(opt[i?+?1][j],?opt[i][j?+?1]);????????//参考上文我给的公式。??
26. ????????????}??
27. ????????}??
28. ???
29. ????????-------------------------------------------------??
30. ???
31. ????????理解上段，参考上文我给的公式：??
32. ???
33. ????????根据上述结论，可得到以下公式，??
34. ???
35. ????????如果我们记字符串Xi和Yj的LCS的长度为c[i,j]，我们可以递归地求c[i,j]：??
36. ???
37. ??????????????????/??????0???????????????????????????????if?i<0?or?j<0??
38. ????????c[i,j]=??????????c[i-1,j-1]+1????????????????????if?i,j>=0?and?xi=xj??
39. ?????????????????/???????max(c[i,j-1],c[i-1,j]???????????if?i,j>=0?and?xi≠xj??
40. ???
41. ????????-------------------------------------------------??
42. ???
43. ????????System.out.println("substring1:"+x);??
44. ????????System.out.println("substring2:"+y);??
45. ????????System.out.print("LCS:");??
46. ???
47. ????????int?i?=?0,?j?=?0;??
48. ????????while?(i?<?substringLength1?&&?j?<?substringLength2){??
49. ????????????if?(x.charAt(i)?==?y.charAt(j)){??
50. ????????????????System.out.print(x.charAt(i));??
51. ????????????????i++;??
52. ????????????????j++;??
53. ????????????}?else?if?(opt[i?+?1][j]?>=?opt[i][j?+?1])??
54. ????????????????i++;??
55. ????????????else??
56. ????????????????j++;??
57. ????????}??
58. ????????Long?endTime?=?System.nanoTime();??
59. ????????System.out.println("?Totle?time?is?"?+?(endTime?-?startTime)?+?"?ns");??
60. ????}??
61. ???
62. ????//取得定长随机字符串??
63. ????public?static?String?GetRandomStrings(int?length){??
64. ????????StringBuffer?buffer?=?new?StringBuffer("abcdefghijklmnopqrstuvwxyz");??
65. ????????StringBuffer?sb?=?new?StringBuffer();??
66. ????????Random?r?=?new?Random();??
67. ????????int?range?=?buffer.length();??
68. ????????for?(int?i?=?0;?i?<?length;?i++){??
69. ????????????sb.append(buffer.charAt(r.nextInt(range)));??
70. ????????}??
71. ????????return?sb.toString();??
72. ????}??
73. }??

第五节、改进的算法

??? 下面咱们来了解一种不同于动态规划法的一种新的求解最长公共子序列问题的方法,该算法主要是把求解公共字符串问题转化为求解矩阵L(p,m)的问题，在利用定理求解矩阵的元素过程中（1）while(i<k),L(k,i)=null，
????????????????? （2）while(L(k,i)=k),L(k,i+1)=L(k,i+2)=…L(k,m)=k；

??? 求出每列元素，一直到发现第p+1 行都为null 时退出循环，得出矩阵L(k,m)后，B[L(1,m-p+1)]B[L(2,m-p+2)]…B[L(p,m)]即为A 和B 的LCS，其中p 为LCS 的长度。

6.1 主要定义及定理

定义 1 子序列(Subsequence)：给定字符串A=A[1]A[2]…A[m]，(A[i]是A 的第i 个字母，A[i]∈字符集Σ，l<= i<m = A ， A 表示字符串A 的长度)，字符串B 是A 的子序列是指B=A[ 1 i ]A[ 2 i ]…A[ k i ],其中1 i < 2 i <…< k i 且k<=m.定义2 公共子序列(Common Subsequence)：给定字符串A、B、C，C 称为A 和B 的公共子序列是指C 既是A 的子序列，又是B 的子序列。定义3 最长公共子序列(Longest Common Subsequence 简称LCS)：给定字符串A、B、C，C 称为A 和B 的最长公共子序列是指C 是A 和B 的公共子序列，且对于A 和B 的任意公共子序列D，都有D <= C 。给定字符串A 和B，A =m，B =n，不妨设m<=n，LCS 问题就是要求出A 和B 的LCS。定义4 给定字符串A=A[1]A[2]…A[m]和字符串B=B[1]B[2]…[n]，A( 1:i)表示A 的连续子序列A[1]A[2]…A[i]，同样B(1:j)表示B 的连续子序列B[1]B[2]…[j]。Li(k)表示所有与字符串A(1:i) 有长度为k 的LCS 的字符串B(l:j) 中j 的最小值。用公式表示就是Li(k)=Minj(LCS(A(1:i)，B(l:j))=k) [3]。

定理1 ? i∈[1，m]，有Li(l)<Li(2)<Li(3)<…<Li(m) .
定理2 ?i∈[l，m-1]，?k∈[l，m]，有i 1 L + (k)<= i L (k).
定理3 ? i∈[l，m-1]， ? k∈[l，m-l]，有i L (k)< i 1 L + (k+l).
以上三个定理都不考虑Li(k)无定义的情况。
定理4[3] i 1 L + (k)如果存在，那么它的取值必为: i 1 L + (k)=Min(j, i L (k))。这里j 是满足以下条件的最小整数:A[i+l]=B[j]且j> i L (k-1)。

??? 矩阵中元素L(k，i)=Li(k)，这里(1<i<=m，1<k<=m)，null 表示L(k,i)不存在。当i<k 时，显然L(k，i)不存在。
??? 设p=Maxk(L(k ， m) ≠ null) ， 可以证明L 矩阵中L(p,m) 所在的对角线,L(1,m-p+1),L(2,m-p+2)…L(p-1,m-1),L(p,m) 所对应的子序列B[L(1,m-p+1)]B[L(2,m-p+2)]…B[L(p,m)]即为A 和B 的LCS，p 为该LCS 的长度。这样，LCS 问题的求解就转化为对m m L × 矩阵的求解。

6.2 算法思想
??? 根据定理,第一步求出第一行元素,L(1,1),L(1,2),…L(1,m),第二步求第二行,一直到发现第p+1 行都为null 为止。在计算过程中遇到i<k 时,L(k,i)=null, 及L(k,i)=k时,L(k,i+1)=L(k,i+2)=…L(k,m)=k。这样,计算每行的时间复杂度为O(n),则整个时间复杂度为O(pn)。在求L 矩阵的过程中不用存储整个矩阵,只需存储当前行和上一行即可。空间复杂度为O(m+n)。

??? 下面给出一个例子来说明:给定字符串A 和B，A=acdabbc，B=cddbacaba，(m= A =7，n= B =9)。按照定理给出的递推公式，求出A 和B 的L 矩阵如图2，其中的$表示NULL。

??? 则A 和B 的LCS 为B[1]B[2]B[4]B[6]=cdbc,LCS 的长度为4。

6.3 算法伪代码
算法 L(A,B,L)
输入 长度分别为m,n 的字符串A,B
输出 A,B 的最长公共子序列LCS

1. L(A,B,L){//字符串A，B，所求矩阵L??
2. ??for(k=1;k<=m;k++){?//m?为A?的长度??
3. ????for(i=1;i<=m;i++){??
4. ??????if(i<k)?L[k][i]=N;//i<k?时,L(k,i)=null，N?代表无穷大??
5. ??????if(L[k][i]==k)//L(k,i)=k?时,L(k,i+1)=L(k,i+2)=…L(k,m)=k??
6. ??????for(l=i+1;l<=m;l++)??
7. ???????{?L[k][l]=k;??
8. ?????????Break;}??
9. ??????for(j=1;j<=n;j++){//定理4?的实现??
10. ???????if(A[i+1]==B[j]&&j>L[k-1][i]){??
11. ????????L[k][i+1]=(j<L[k][i]?j:L[k][i]);??
12. ????????break;??
13. ??????}??
14. ??????if(L[k][i+1]==0)??
15. ????????L[k][i]=N;??
16. ?????}??
17. ?????if(L[k][m]==N)??
18. ??????{p=k-1;break;}??
19. ??}??
20. ??p=k-1;??
21. }??

6.4 结语
??? 本节主要描述区别于动态规划法的一种新的求解最长公共子序列问题的方法，在不影响精确度的前提下，提高序列匹配的速度，根据定理i 1 L + (k)=Min(j, i L (k))得出矩阵，在求解矩阵的过程中对最耗时的L(p,m)进行条件约束优化。我们在Intel(R) Core(TM)2 Quad 双核处理器、1G 内存，软件环境：windows xp 下试验结果证明，本文算法与其他经典的比对算法相比,不但能够取得准确的结果,而且速度有了较大的提高（本节参考了刘佳梅女士的论文）。

??? 若有任何问题，恳请不吝指正。谢谢各位。完。