最长公共子序列(LCS)

1.摘要：

继上篇最长上升子序列后，本篇主要讲述最长公共子序列 (LCS) 。

 

2.LCS定义：

       最长公共子序列，英文缩写为LCS（Longest Common Subsequence）。其定义是，一个序列 S ，如果分别是两个或多个已知序列的子序列，且是所有符合此条件序列中最长的，则 S 称为已知序列的最长公共子序列。

       如果觉得抽象不好理解，那么咱们还是采用学习LIS的时候的方式。首先，让我们先来看一下子串、子序列还有公共子序列的概念(在上篇LIS中也曾涉及过) ，我们以字符子串和字符子序列为例，更为形象，也能顺带着理解字符的子串和子序列：

     （1）字符子串：指的是字符串中连续的n个字符，如abcdefg中，ab，cde，fg等都属于它的字串。

     （2）字符子序列：指的是字符串中不一定连续但先后顺序一致的n个字符，即可以去掉字符串中的部分字符，但不可改变其前后顺序。如abcdefg中，acdg，bdf属于它的子序列，而bac，dbfg则不是，因为它们与字符串的字符顺序不一致。

       (3)  公共子序列：如果序列C既是序列A的子序列，同时也是序列B的子序列，则称它为序列A和序列B的公共子序列。如对序列 1 , 3 , 5 , 4 , 2 , 6 , 8 , 7 和序列 1 , 4 , 8 , 6 , 7 , 5  来说，序列1,8,7是它们的一个公共子序列。

       那么现在，我们再通俗的总结一下最长公共子序列（LCS）：就是A和B的公共子序列中长度最长的（包含元素最多的）
仍然用序列 1 , 3 , 5 , 4 , 2 , 6 , 8 , 7 和序列 1 , 4 , 8 , 6 , 7 , 5 为例，它们的最长公共子序列有 1 , 4 , 8 , 7 和 1 , 4 , 6 , 7 两种，但最长公共子序列的长度是4。由此可见，最长公共子序列（LCS）也不一定唯一。

      请大家用集合的观点来理解这些概念，子序列、公共子序列以及最长公共子序列都不唯一，所以我们通常特判取一个最长公共子序列，但很显然，对于固定的两个数组，虽然最LCS不一定唯一，但LCS的长度是一定的。查找最长公共子序列与查找最长公共子串的问题不同的地方在于：子序列不需要在原序列中占用连续的位置。最长公共子串（要求连续）和最长公共子序列是不同的。

那么该如何求出两个序列的最长公共子序列长度呢？请继续往下看～

 

3.LCS长度求法：

       你首先能想到的恐怕是暴力枚举？那我们先来看看：序列A有 2^n 个子序列，序列B有 2^m 个子序列，如果任意两个子序列一一比较，比较的子序列高达 2^(n+m) 对，这还没有算具体比较的复杂度。或许你说，只有长度相同的子序列才会真正进行比较。那么忽略空序列，我们来看看：对于A长度为1的子序列有C(n,1)个，长度为2的子序列有C(n,2)个，……长度为n的子序列有C(n,n)个。对于B也可以做类似分析，即使只对序列A和序列B长度相同的子序列做比较，那么总的比较次数高达：C(n,1)*C(m,1)*1 + C(n,2) * C(m,2) * 2+ …+C(n,p) * C(m,p)*p，其中p = min(m, n)。

       吓着了吧？怎么办？我们试试使用动态规划算法！

       我们用Ax表示序列A的连续前x项构成的子序列，即Ax= a1,a2,……ax, By= b1,b2,……by, 我们用LCS(x, y)表示它们的最长公共子序列长度，那原问题等价于求LCS(m,n)。为了方便我们用L(x, y)表示Ax和By的一个最长公共子序列。让我们来看看如何求LCS(x, y)。我们令x表示子序列考虑最后一项

（1） Ax ＝ By

         那么它们L(Ax, By)的最后一项一定是这个元素！

       为什么呢？为了方便，我们令t = Ax = By, 我们用反证法：假设L(x,y)最后一项不是t，则要么L(x,y)为空序列（别忘了这个），要么L(x,y)的最后一项是Aa＝Bb ≠ t, 且显然有a < x, b < y。无论是哪种情况我们都可以把t接到这个L(x,y)后面,从而得到一个更长的公共子序列。矛盾！
       如果我们从序列Ax中删掉最后一项ax得到Ax-1,从序列By中也删掉最后一项by得到By-1，(多说一句角标为0时，认为子序列是空序列)，则我们从L(x,y)也删掉最后一项t得到的序列是L(x – 1, y - 1)。为什么呢？和上面的道理相同，如果得到的序列不是L(x - 1, y - 1)，则它一定比L(x - 1, y - 1)短（注意L（，）是个集合！），那么它后面接上元素t得到的子序列L(x,y)也比L(x - 1, y - 1)接上元素t得到的子序列短，这与L(x, y)是最长公共子序列矛盾。因此L(x, y) = L(x - 1, y - 1) 最后接上元素t，LCS(Ax, By) = LCS(x - 1, y - 1) + 1。

（2）  Ax ≠ By

        仍然设t = L(Ax, By), 或者L(Ax, By)是空序列（这时t是未定义值不等于任何值）。则t  ≠ Ax和t  ≠ By至少有一个成立，因为t不能同时等于两个不同的值嘛！

（2.1）如果t  ≠ Ax，则有L(x, y)= L(x - 1, y)，因为根本没Ax的事嘛。

            LCS(x,y) = LCS(x – 1, y)
（2.2）如果t  ≠ By,l类似L(x, y)= L(x , y - 1)

            LCS(x,y) = LCS(x, y – 1)
       可是，我们事先并不知道t，由定义，我们取最大的一个，因此这种情况下,有LCS(x,y) = max(LCS(x – 1, y) , LCS(x, y – 1))。看看目前我们已经得到了什么结论：

LCS(x,y) = 
(1) LCS(x - 1,y - 1) + 1      （Ax ＝ By）
(2) max(LCS(x – 1, y) , LCS(x, y – 1))    （Ax ≠ By）

这时一个显然的递推式，光有递推可不行，初值是什么呢？显然，一个空序列和任何序列的最长公共子序列都是空序列！所以我们有:

LCS(x,y) = 
(1) LCS(x - 1,y - 1) + 1 如果Ax ＝ By
(2) max(LCS(x – 1, y) , LCS(x, y – 1)) 如果Ax ≠ By
(3) 0 如果x = 0或者y = 0

到此我们求出了计算最长公共子序列长度的递推公式。我们实际上计算了一个(n + 1)行(m + 1)列的表格（行是0..n，列是0..m)，也就这个二维度数组LCS(,)。

n²算法c++模版：

#include <iostream>
using namespace std;
const int MAXN = 1000 + 10;
int n, data[MAXN];
int dp[MAXN]; 
int from[MAXN]; 
void output(int x)
{
    if(!x)return;
    output(from[x]);
    cout<<data[x]<<" ";
    //迭代输出 
}
int main()
{
    cin>>n;
    for(int i=1;i<=n;i++)cin>>data[i];
    
    // DP
    for(int i=1;i<=n;i++)
    {
        dp[i]=1;
        from[i]=0;
        for(int j=1;j<i;j++)
        if(data[j]<data[i] && dp[i]<dp[j]+1)
        {
            dp[i]=dp[j]+1;
            from[i]=j;//逐个记录前驱 
        }
    }
    
    int ans=dp[1], pos=1;
    for(int i=1;i<=n;i++)
        if(ans<dp[i])
        {
            ans=dp[i];
            pos=i;//由于需要递归输出
    //所以要记录最长上升子序列的最后一
    //个元素，来不断回溯出路径来 
        }
    cout<<ans<<endl;
    output(pos);
    
    return 0;
}

nlog(n)算法c++模版（以上讲的）：

#include<iostream>
#include<cstdio>
using namespace std;
int a[100001],b[100001],map[100001],f[100001];
int main()
{
    int n;
    cin>>n;
    for(int i=1;i<=n;i++){scanf("%d",&a[i]);map[a[i]]=i;}
    for(int i=1;i<=n;i++){scanf("%d",&b[i]);f[i]=0x7fffffff;}
    int len=0;
    f[0]=0;
    for(int i=1;i<=n;i++)
    {
        int l=0,r=len,mid;
        if(map[b[i]]>f[len])f[++len]=map[b[i]];
        else 
        {
        while(l<r)
        {    
            mid=(l+r)/2;
            if(f[mid]>map[b[i]])r=mid;
            else l=mid+1; 
        }
        f[l]=min(map[b[i]],f[l]);
         }
    }
    cout<<len;
    return 0
}

 学完以上知识后，你就可以刷洛谷 P1349，以下一篇通俗易懂的题解，共还没理解的OIER继续阅读。

 

一、$DP$的意义以及线性动规简介

动态规划自古以来是$DALAO$凌虐萌新的分水岭，但有些OIer认为并没有这么重要——会打暴力，大不了记忆化。但是其实，动态规划学得好不好，可以彰显出一个$OIer$的基本素养——能否富有逻辑地思考一些问题，以及更重要的——能否将数学、算筹学（决策学）、数据结构合并成一个整体并且将其合理运用$qwq$。

而我们首先要了解的，便是综合难度在所有动规题里最为简单的线性动规了。线性动规既是一切动规的基础，同时也可以广泛解决生活中的各项问题——比如在我们所在的三维世界里，四维的时间就是不可逆式线性，比如我们需要决策在相同的时间内做价值尽量大的事情，该如何决策，最优解是什么——这就引出了动态规划的真正含义：

在一个困难的嵌套决策链中，决策出最优解。

二、动态规划性质浅谈

首先，动态规划和递推有些相似（尤其是线性动规），但是不同于递推的是：

递推求出的是数据，所以只是针对数据进行操作；而动态规划求出的是最优状态，所以必然也是针对状态的操作，而状态自然可以出现在最优解中，也可以不出现——这便是决策的特性（布尔性）。

其次，由于每个状态均可以由之前的状态演变形成，所以动态规划有可推导性，但同时，动态规划也有无后效性，即每个当前状态会且仅会决策出下一状态，而不直接对未来的所有状态负责。

现在决定未来，未来与过去无关。

三、扯正题——子序列问题

（一）一个序列中的最长上升子序列（$LIS$）

例：由6个数，分别是： 1 7 6 2 3 4，求最长上升子序列。

评析：首先，我们要理解什么叫做最长上升子序列：1、最长上升子序列的元素不一定相邻 2、最长上升子序列一定是原序列的子集。所以这个例子中的$LIS$就是：1 2 3 4，共4个

1、$n^2$做法

首先我们要知道，对于每一个元素来说，最长上升子序列就是其本身。那我们便可以维护一个$dp$数组，使得**$dp[i]$表示以第$i$元素为结尾的最长上升子序列长度**，那么对于每一个$dp[i]$而言，初始值即为$1$；

那么dp数组怎么求呢？我们可以对于每一个$i$，枚举在$i$之前的每一个元素$j$，然后对于每一个$dp[j]$,如果元素$i$大于元素$j$，那么就可以考虑继承，而最优解的得出则是依靠对于每一个继承而来的$dp$值，取$max$.

	for(int i=1;i<=n;i++)
	{
		dp[i]=1;//初始化 
		for(int j=1;j<i;j++)//枚举i之前的每一个j 
		if(data[j]<data[i] && dp[i]<dp[j]+1)
		//用if判断是否可以拼凑成上升子序列，
		//并且判断当前状态是否优于之前枚举
		//过的所有状态,如果是，则↓ 
		dp[i]=dp[j]+1;//更新最优状态 
		
	}

最后，因为我们对于$dp$数组的定义是到i为止的最长上升子序列长度，所以我们最后对于整个序列，只需要输出$dp[n]$($n$为元素个数)即可。

从这个题我们也不难看出，状态转移方程可以如此定义：

下一状态最优值=最优比较函数（已经记录的最优值，可以由先前状态得出的最优值）

——即动态规划具有 判断性继承思想

2、$nlogn$ 做法

我们其实不难看出，对于$n^{}$做法而言，其实就是暴力枚举：将每个状态都分别比较一遍。但其实有些没有必要的状态的枚举，导致浪费许多时间，当元素个数到了$10^4-10^5$以上时，就已经超时了。而此时，我们可以通过另一种动态规划的方式来降低时间复杂度：

将原来的dp数组的存储由数值换成该序列中，上升子序列长度为i的上升子序列，的最小末尾数值

这其实就是一种几近贪心的思想：我们当前的上升子序列长度如果已经确定，那么如果这种长度的子序列的结尾元素越小，后面的元素就可以更方便地加入到这条我们臆测的、可作为结果、的上升子序列中。

qwq一定要好好看注释啊！

int n;
	cin>>n;
	for(int i=1;i<=n;i++)
	{
		cin>>a[i];
		f[i]=0x7fffffff;
		//初始值要设为INF
		/*原因很简单，每遇到一个新的元素时，就跟已经记录的f数组当前所记录的最长
		上升子序列的末尾元素相比较：如果小于此元素，那么就不断向前找，直到找到
		一个刚好比它大的元素，替换；反之如果大于，么填到末尾元素的下一个q，INF
                就是为了方便向后替换啊！*/ 
	}
	f[1]=a[1];
	int len=1;//通过记录f数组的有效位数，求得个数 
	/*因为上文中所提到我们有可能要不断向前寻找，
	所以可以采用二分查找的策略，这便是将时间复杂
    度降成nlogn级别的关键因素。*/ 
	for(int i=2;i<=n;i++)
	{
		int l=0,r=len,mid;
		if(a[i]>f[len])f[++len]=a[i];
		//如果刚好大于末尾，暂时向后顺次填充 
		else 
		{
		while(l<r)
		{	
		    mid=(l+r)/2;
		    if(f[mid]>a[i])r=mid;
	//如果仍然小于之前所记录的最小末尾，那么不断
	//向前寻找(因为是最长上升子序列，所以f数组必
	//然满足单调) 
			else l=mid+1; 
		}
		f[l]=min(a[i],f[l]);//更新最小末尾 
     	}
    }
    cout<<len;

---

$AnotherSituation$

但是事实上，$nlogn$做法偷了个懒，没有记录以每一个元素结尾的最长上升子序列长度。那么我们对于$n^2$的统计方案数，有很好想的如下代码（再对第一次的$dp$数组$dp$一次）：

for(i = 1; i <= N; i ++){
	if(dp[i] == 1) f[i] = 1 ;
	for(j = 1; j <= N: j ++)
		if(base[i] > base[j] && dp[j] == dp[i] - 1) f[i] += f[j] ;
		else if(base[i] == base[j] && dp[j] == dp[i]) f[i] = 0 ;
	if(f[i] == ans) res ++ ;
	}

但是$nlogn$呢？虽然好像也可以做，但是想的话会比较麻烦，在这里就暂时不讨论了$qwq$，但笔者说这件事的目的是为了再次论证一个观点:时间复杂度越高的算法越全能

---

$3$、输出路径

只要记录前驱，然后递归输出即可（也可以用栈的）

下面贴出$n^2$的完整代码qwq

#include <iostream>
using namespace std;
const int MAXN = 1000 + 10;
int n, data[MAXN];
int dp[MAXN]; 
int from[MAXN]; 
void output(int x)
{
	if(!x)return;
	output(from[x]);
	cout<<data[x]<<" ";
	//迭代输出 
}
int main()
{
	cin>>n;
	for(int i=1;i<=n;i++)cin>>data[i];
	
	// DP
	for(int i=1;i<=n;i++)
	{
		dp[i]=1;
		from[i]=0;
		for(int j=1;j<i;j++)
		if(data[j]<data[i] && dp[i]<dp[j]+1)
		{
			dp[i]=dp[j]+1;
			from[i]=j;//逐个记录前驱 
		}
	}
	
	int ans=dp[1], pos=1;
	for(int i=1;i<=n;i++)
		if(ans<dp[i])
		{
			ans=dp[i];
			pos=i;//由于需要递归输出
	//所以要记录最长上升子序列的最后一
	//个元素，来不断回溯出路径来 
		}
	cout<<ans<<endl;
	output(pos);
	
	return 0;
}

（二）两个序列中的最长公共子序列（$LCS$）

1、譬如给定2个序列：

1 2 3 4 5

3 2 1 4 5

试求出最长的公共子序列。

$qwq$显然长度是$3$，包含$345$ 三个元素（不唯一）

解析：我们可以用$dp[i][j]$来表示第一个串的前$i$位，第二个串的前j位的$LCS$的长度，那么我们是很容易想到状态转移方程的：

如果当前的$A1[i]$和$A2[j]$相同（即是有新的公共元素） 那么

$dp[i][j]=max(dp[i][j],dp[i-1][j-1]+1);$

如果不相同，即无法更新公共元素，考虑继承：

$dp[i][j]=max(dp[i-1][j],dp[i][j-1]$

那么代码:

#include<iostream>
using namespace std;
int dp[1001][1001],a1[2001],a2[2001],n,m;
int main()
{
   //dp[i][j]表示两个串从头开始，直到第一个串的第i位 
   //和第二个串的第j位最多有多少个公共子元素 
   cin>>n>>m;
   for(int i=1;i<=n;i++)scanf("%d",&a1[i]);
   for(int i=1;i<=m;i++)scanf("%d",&a2[i]);
   for(int i=1;i<=n;i++)
    for(int j=1;j<=m;j++)
     {
     	dp[i][j]=max(dp[i-1][j],dp[i][j-1]);
     	if(a1[i]==a2[j])
     	dp[i][j]=max(dp[i][j],dp[i-1][j-1]+1);
     	//因为更新，所以++； 
     }
   cout<<dp[n][m];
}

$2$、而对于洛谷$P1439$而言，不仅是卡上面的朴素算法，也考察到了全排列的性质：

对于这个题而言，朴素算法是$n^2$的，会被$10^5$卡死，所以我们可以考虑$nlogn$的做法：

因为两个序列都是$1n$的全排列，那么两个序列元素互异且相同，也就是说只是位置不同罢了，那么我们通过一个$map$数组将$A$序列的数字在$B$序列中的位置表示出来——

因为最长公共子序列是按位向后比对的，所以a序列每个元素在b序列中的位置如果递增，就说明b中的这个数在a中的这个数整体位置偏后，可以考虑纳入$LCS$——那么就可以转变成$nlogn$求用来记录新的位置的map数组中的**$LIS$**。

最后贴$AC$代码：

#include<iostream>
#include<cstdio>
using namespace std;
int a[100001],b[100001],map[100001],f[100001];
int main()
{
	int n;
	cin>>n;
	for(int i=1;i<=n;i++){scanf("%d",&a[i]);map[a[i]]=i;}
	for(int i=1;i<=n;i++){scanf("%d",&b[i]);f[i]=0x7fffffff;}
	int len=0;
	f[0]=0;
	for(int i=1;i<=n;i++)
	{
		int l=0,r=len,mid;
		if(map[b[i]]>f[len])f[++len]=map[b[i]];
		else 
		{
		while(l<r)
		{	
		    mid=(l+r)/2;
		    if(f[mid]>map[b[i]])r=mid;
			else l=mid+1; 
		}
		f[l]=min(map[b[i]],f[l]);
     	}
    }
    cout<<len;
    return 0
}