动态规划-最长公共子序列问题

最长公共子序列问题
    最长公共子序列问题(longest-common-subsequence problem)给定两个序列X=<x1,x2,...,xm>和Y=<y1,y2,...,yn>,求X和Y长度最小的公共子序列(可以不连续)。接下来将展示如何用动态规划方法高效地求解LCS问题。
步骤1：刻画最长公共子序列的特征
     如果用暴力搜索方法求解LCS问题，就要穷举X的所有子序列，对每个子序列检查它是否也是Y的子序列，记录找到的最长子序列。X的每个子序列对应X的下标集合<1,2,...,m>，所以X有2^m个子序列(例序列<1,2,3>的子序列分别是1,2,3,12,13,23,123,所以应该是(2^m)-1个子序列)，因此暴力方法的运行时间为指数阶，对较长的序列是不实用的。
     但是LCS问题具有最优子结构性质。我们将看到，子问题的自然分类对应两个输入序列的“前缀”对。前缀的严谨定义如下：给定一个序列X=<X1,X2，...Xm>，对i=0,1,...,m，定义X的第i前缀为Xi=<x1,x2,...,xi>。例如，若X=<A,B,C,B,D,A,B>，则X4=<A,B,C,B>，X0为空串。
     存在以下定理(LCS的最优子结构) 令X=<x1,x2,...,xm>和Y=<y1,y2,...,yn>为两个序列，Z=<z1,z2,...,zk>为X和Y的任意LCS。
   1.如果x(m)=y(n)，则z(k)=x(m)=y(n)且Z(k-1)是X(m-1)和Y(n-1)的一个LCS。
   2.如果x(m)≠y(n)，那么z(k)≠x(m)意味着Z是X(m-1)和Y的一个LCS。
   3.如果x(m)≠y(n)，那么z(k)≠y(n)意味着Z是X和Y(n-1)的一个LCS。
步骤2：一个递归解
   上述定理意味着，在求序列X和Y的一个LCS时，我们需要求解一个或两个子问题。如果x(m)=y(n)，我们应该求解X(m-1)和Y(n-1)的一个LCS。将x(m)=y(n)追加到这个LCS的末尾，就得到X和Y的一个LCS。如果x(m)≠y(n)，我们必须求解两个子问题：求X(m-1)和Y的一个LCS与X和Y(n-1)的一个LCS。两个LCS较长者即为X和Y的一个LCS。由于这些情况覆盖了所有可能性，因此我们知道必然有一个子问题的最优解出现在X和Y的LCS中。
   我们可以很容易看出LCS问题的重叠子问题性质。为了求X和Y的一个LCS，我们可能要求 X和Y(n-1)的一个LCS以及X(m-1)和Y的一个LCS。但这几个子问题都包含求解X(m-1)和Y(n-1)的LCS的子子问题。很多其他子问题也都享有子子问题。
   与矩阵链乘法问题相似，设计LCS问题的递归算法首先要建立最优解的递归式。我们定义c[i,j]表示X(i)和Y(j)的LCS的长度。如果i=0或j=0，即一个序列的长度为0，那么LCS的长度为0。根据LCS问题的最优子结构性质，可得到下面公式：
   ①若i=0或j=0，c[i,j]=0；
   ②若i,j>0且x(i)=y(j)，c[i,j]=c[i-1,j-1]+1。
   ③若i,j>0且x(i)≠y(j)，c[i,j]=max(c[i,j-1],c[i-1,j])
   观察上述公式，我们通过限制条件限定了需要求解那些子问题。当x(i)=y(j)我们可以而且应该求解子问题：X(i-1)和Y(j-1)的一个LCS。否则，应该求解两个子问题：X(i)和Y(j-1)的一个LCS及X(i-1)和Y(j)和一个LCS。在之前的钢条切割和矩阵链乘法问题的动态规划算法中，根据问题的条件，我们没有排除任何子问题。不过，LCS问题并非唯一根据条件排除子问题的动态规划算法。
步骤3：计算LCS的长度
    我们可以很容易地写出一个指数时间的递归算法来计算两个序列的LCS的长度。但是由于LCS问题只有Θ(mn)个不同的子问题，我们可以用动态规划方法自底向上地计算。

过程LCS_LENGTH接受两个序列X=<x1,x2,...xm>和Y=<y1,y2,...y2>为输入。它将c[i,j]的值保存在表c中，并按行主次序(row-major order)计算表项(即首先由左至右计算c的第一行，然后计算第二行，依次类推)。过程还维护一个表b，帮助构造最优解。b[i,j]指向表项对应计算c[i,j]时选择的子问题最优解。过程返回表b和c，c[m,n]保存了X和Y的LCS的长度。

void LCS_LENGTH(string X,string Y,int **c, int **b)
{
	int m=X.length();
	int n=Y.length();
	for(int i=1;i<=m;i++)   c[i][0]=0;
	for(int j=0;j<=n;j++)  c[0][j]=0; 
	for(int i=1;i<=m;i++)
	{
		for(int j=1;j<=n;j++)
		 {
			if(X[i]==Y[j])
			{
				c[i][j]=c[i-1][j-1]+1;
				b[i][j]='↖';  /* 这个符号是int类似 */
			}
			 else if(c[i-1][j]>=c[i][j-1])
			{
				c[i][j]=c[i-1][j];
				b[i][j]='↑';
			 }
			else 
			{
				c[i][j]=c[i][j-1];
				b[i][j]='←';
			 }
		  }
	}
}

      下图显示了LCS_LENGTH对输入序列X=<A,B,C,B,D,A,B>和Y=<B,D,C,A,B,A>生成的过程。过程运行的时间为Θ(mn)，因为每个表项的计算时间为Θ(1)。

      其中表项c[7,6]中的4即为X和Y的一个LCS <B,C,B,A>的长度。对所有i,j>0，表项c[i,j]仅依赖于x(i)=y(j)以及c[i-1][j]、c[i][j-1]和c[i-1][j-1]的值买这些值都会在c[i][j]之前计算出来。为了构造LCS中的元素，从右下角开始沿着b[i,j]的箭头前进即可，如图中阴影方格序列。阴影序列中每个“↖”对应的表项表示x(i)=y(j)是LCS的一个元素。
步骤4：构造LCS

   我们可以用LCS_LENGTH返回的表b快速构造X=<x1,x2,...,xm>和Y=<y1,y2,...,yn>的LCS，只需简单地从b[m,n]开始，并按箭头方向追踪下去即可。当在表项b[i,j]中遇到一个“↖”时，意味着x(i)=y(j)是LCS的一个元素。按照这种方法，我们可以按逆序依次构造出LCS的所有元素。下面的递归过程会按正确的顺序打印出X和Y的一个LCS。对它的起始调用为PRINT_LCS

void PRINT_LCS(int **b,string X,int i,int j)

{
	if(i==0||j==0) return ;
	if(b[i][j]=='↖') 
	{
		PRINT_LCS(b,X,i-1,j-1);
		cout<<X[i];
	}
	else if(b[i][j]=='↑')
	{
		PRINT_LCS(b,X,i-1,j);
	}
	else  PRINT_LCS(b,X,i,j-1);
}    

   对于上图中表b，此过程会打印出BCBA。过程运行时间为O(m+n)，每次递归调用i和j至少有一个会减少1。

下面给出完整的代码实现，假定给定序列X=<A,B,C,B,D,A,B>和Y=<B,D,C,A,B,A>。

#include<iostream>
#include<string>
using namespace std;
void LCS_LENGTH(string X,string Y,int **c, int **b)
{
	int m=X.length();
	int n=Y.length();
	for(int i=1;i<=m;i++)   c[i][0]=0;
	for(int j=0;j<=n;j++)  c[0][j]=0; 
	for(int i=1;i<=m;i++)
	{
		for(int j=1;j<=n;j++)
		 {
			if(X[i]==Y[j])
			{
				c[i][j]=c[i-1][j-1]+1;
				b[i][j]='↖';  /* 这个符号是int类似 */
			}
			 else if(c[i-1][j]>=c[i][j-1])
			{
				c[i][j]=c[i-1][j];
				b[i][j]='↑';
			 }
			else 
			{
				c[i][j]=c[i][j-1];
				b[i][j]='←';
			 }
		  }
	}
}	
void PRINT_LCS(int **b,string X,int i,int j)
{
	if(i==0||j==0) return ;
	if(b[i][j]=='↖') 
	{
		PRINT_LCS(b,X,i-1,j-1);
		cout<<X[i];
	}
	else if(b[i][j]=='↑')
	{
		PRINT_LCS(b,X,i-1,j);
	}
	else  PRINT_LCS(b,X,i,j-1);
}  
int main()
{
	string X="0ABCBDAB";  /* 因为对序列的操作是从1开始，所以X[0]=0 */
	string Y="0BDCABA";
	int m=X.length();
	int n=Y.length();
	/* 动态申请 */
	int **c=new int *[m+1];
	for(int i=0;i<=m;i++) c[i]=new int[n+1];
	int **b=new int *[m+1];
	for(int i=0;i<=m;i++) b[i]=new int[n+1];
	/* 计算表b和c */
	LCS_LENGTH(X,Y,c,b);
	/* 打印 */
	PRINT_LCS(b,X,m,n);
	/* 回收内存 */
	for (int i = 0; i < m+1; i++)     
       {  
          delete b[i];     
          b[i] = NULL;
	  delete c[i];
	  c[i] = NULL;
       }  
       delete []b;     
        b = NULL; 
	delete []c;
	c = NULL;
	return 0;
}