首页 > 算法 数据结构 > 最长公共子序列(Longest-Common-Subsequence,LCS)

最长公共子序列(Longest-Common-Subsequence,LCS)

2012年6月12日 发表评论 阅读评论
文章作者:Yx.Ac   文章来源:勇幸|Thinking (http://www.ahathinking.com)   转载请注明,谢谢合作。

---

一个字符串S,去掉零个或者多个元素所剩下的子串称为S的子序列。最长公共子序列就是寻找两个给定序列的子序列,该子序列在两个序列中以相同的顺序出现,但是不必要是连续的。

例如序列X=ABCBDAB,Y=BDCABA。序列BCA是X和Y的一个公共子序列,但是不是X和Y的最长公共子序列,子序列BCBA是X和Y的一个LCS,序列BDAB也是。

寻找LCS的一种方法是枚举X所有的子序列,然后注意检查是否是Y的子序列,并随时记录发现的最长子序列。假设X有m个元素,则X有2^m个子序列,指数级的时间,对长序列不实际。

使用动态规划求解这个问题,先寻找最优子结构。设X=<x1,x2,…,xm>和Y=<y1,y2,…,yn>为两个序列,LCS(X,Y)表示X和Y的一个最长公共子序列,可以看出

  1. 如果xm=yn,则LCS ( X,Y ) = xm + LCS ( Xm-1,Yn-1 )。
  2. 如果xm!=yn,则LCS( X,Y )= max{ LCS ( Xm-1, Y ), LCS ( X, Yn-1 ) }

LCS问题也具有重叠子问题性质:为找出X和Y的一个LCS,可能需要找X和Yn-1的一个LCS以及Xm-1和Y的一个LCS。但这两个子问题都包含着找Xm-1和Yn-1的一个LCS,等等.

DP最终处理的还是数值(极值做最优解),找到了最优值,就找到了最优方案;为了找到最长的LCS,我们定义dp[i][j]记录序列LCS的长度,合法状态的初始值为当序列X的长度为0或Y的长度为0,公共子序列LCS长度为0,即dp[i][j]=0,所以i和j分别表示序列X的长度和序列Y的长度,状态转移方程为

  1. dp[i][j] = 0  如果i=0或j=0
  2. dp[i][j] = dp[i-1][j-1] + 1  如果X[i-1] = Y[i-1]
  3. dp[i][j] = max{ dp[i-1][j], dp[i][j-1] }  如果X[i-1] != Y[i-1]

求出了最长公共子序列的长度后,输出LCS就是输出dp的最优方案了,这在01背包中已经讲过,既可以用一个额外的矩阵存储路径,也可以直接根据状态转移矩阵倒推最优方案。代码如下:

#include <iostream>
using namespace std;

/* LCS
 * 设序列长度都不超过20
*/

int dp[21][21]; /* 存储LCS长度, 下标i,j表示序列X,Y长度 */
char X[21];
char Y[21];
int i, j;

void main()
{
	cin.getline(X,20);
	cin.getline(Y,20);

	int xlen = strlen(X);
	int ylen = strlen(Y);

	/* dp[0-xlen][0] & dp[0][0-ylen] 都已初始化0 */
	for(i = 1; i <= xlen; ++i)
	{
		for(j = 1; j <= ylen; ++j)
		{
			if(X[i-1] == Y[j-1])
			{
				dp[i][j] = dp[i-1][j-1] + 1;
			}else if(dp[i][j-1] > dp[i-1][j])
			{
				dp[i][j] = dp[i][j-1];
			}else
			{
				dp[i][j] = dp[i-1][j];
			}
		}
	}
	printf("len of LCS is: %d\n", dp[xlen][ylen]);

	/* 输出LCS 本来是逆序打印的,可以写一递归函数完成正序打印
	   这里采用的方法是将Y作为临时存储LCS的数组,最后输出Y
	*/
	i = xlen;
	j = ylen;
	int k = dp[i][j];
	Y[k] = '\0';
	while(i && j)
	{
		if(dp[i][j] == dp[i-1][j-1] + 1)
		{
			Y[--k] = X[i-1];
			--i; --j;
		}else if(dp[i-1][j] > dp[i][j-1])
		{
			--i;
		}else
		{
			--j;
		}
	}
	printf("%s\n",Y);
}

在LCS问题中,如果仅仅要求求出LCS的长度,而不要求输出序列,那么由于每步迭代都只用到了前面的状态,之前的信息便无用了,我们就可以使用滚动数组了,代码如下:

#include <iostream>
using namespace std;

/* 滚动数组 */

int dp[2][21];  /* 存储LCS长度 */
char X[21];
char Y[21];
int i, j, k;

void main()
{
	cin.getline(X,20);
	cin.getline(Y,20);

	int xlen = strlen(X);
	int ylen = strlen(Y);

	for(i = 1; i <= xlen; ++i)
	{
		k = i & 1;
		for(j = 1; j <= ylen; ++j)
		{
			if(X[i-1] == Y[j-1])
			{
				dp[k][j] = dp[k^1][j-1] + 1;
			}else if(dp[k][j-1] > dp[k^1][j])
			{
				dp[k][j] = dp[k][j-1];
			}else
			{
				dp[k][j] = dp[k^1][j];
			}
		}
	}
	printf("len of LCS is: %d\n", dp[k][ylen]);
}

本文相关代码可以到这里下载。

(全文完)

参考资料:算法导论 15.4节

勘误:(感谢@jiawei)

LCS在输出子序列时,忽略了测试条件,勘误如下

	i = xlen;
	j = ylen;
	int k = dp[i][j];
	char lcs[21] = {'\0'};
	while(i && j)
	{
		if(X[i-1] == Y[j-1] && dp[i][j] == dp[i-1][j-1] + 1)
		{
			lcs[--k] = X[i-1];
			--i; --j;
		}else if(X[i-1] != Y[j-1] && dp[i-1][j] > dp[i][j-1])
		{
			--i;
		}else
		{
			--j;
		}
	}
	printf("%s\n",lcs);
  • yongmi

    第三个等式同理

  • yongmi

    dp[i][j] = dp[i-1][j-1] + 1 如果X[i-1] = Y[i-1]
    应该是下面的等式。
    dp[i][j] = dp[i-1][j-1] + 1 如果X[i] = Y[j]

    • Yx.Ac

      这里是这样,我们“用i和j分别表示序列X的长度和序列Y的长度”,又因为数组下标是从0开始的,所以是这样,并没有错。您实现一下代码就明白为啥这样写了,欢迎继续讨论 :-)

      • SCaffrey

        他的意思是x[i - 1] == y[j - 1]不是==y[i - 1]?

  • Yx.Ac

    再次感谢@jiawei,又是实现的错误,忽视了条件,当时测试了一条通过了就没再去做过多测试,请见勘误@jiawei

  • jiawei

    运行第一段代码时,输入:
    abcd
    bd

    结果:len of LCS is: 2
    cd
    有误!