LCS最长子序列匹配算法讲解

首先将要看到如何运用动态编程查找两个DNA 序列的最长公共子序列（longest common subsequence，LCS）。

发现了新的基因序列的生物学家通常想知道该基因序列与其他哪个序列最相似。

查找LCS 是计算两个序列相似程度的一种方法：LCS 越长，两个序列越相似。

子序列中的字符与子字符串中的字符不同，它们不需要是连续的。

例如，ACE是ABCDE的子序列，但不是它的子字符串。

请看下面两个DNA 序列：∙S1 = DE>GCCCTAGCGDE>∙S2 = DE>GCGCAATGDE>这两个序列的LCS 是GCCAG。

（请注意，这仅是一个LCS，而不是唯一的LCS，因为可能存在其他长度相同的公共子序列。

这种最优化问题和其他最优化问题的解可能不止一个。

）LCS 算法首先，考虑如何递归地计算LCS。

令：∙C1是S1最右侧的字符∙C2是S2最右侧的字符∙S1'是S1中“切掉” C1的部分∙S2'是S2中“切掉” C2的部分有三个递归子问题：∙L1 = LCS(S1', S2)∙L2 = LCS(S1, S2')∙L3 = LCS(S1', S2')结果表明（而且很容易使人相信）原始问题的解就是下面三个子序列中最长的一个：∙L1∙L2∙如果C1等于C2，则为L3后端加上C1，如果C1不等于C2，则为L3。

（基线条件（base case）是S1或S2为长度为零的字符串的情形。

在这种情况下，S1和S2的LCS 显然是长度为零的字符串。

）但是，就像计算斐波纳契数的递归过程一样，这个递归解需要多次计算相同的子问题。

可以证明，这种递归解法需要耗费指数级的时间。

相比之下，这一问题的动态编程解法的运行时间是Θ(mn)，其中m和n分别是两个序列的长度。

为了用动态编程有效地计算LCS，首先需要构建一个表格，用它保存部分结果。

沿着顶部列出一个序列，再沿着左侧从上到下列出另一个序列，如图2 所示：图2. 初始LCS 表格这种方法的思路是：将从上向下、从左到右填充表格，每个单元格包含一个数字，代表该行和该列之前的两个字符串的LCS 的长度。

也就是说，每个单元格包含原始问题的一个字问题的解。

例如，请看第6 行第7 列的单元格：它在GCGCAATG 序列第二个C 的右侧，在GCCCTAGCG 的T 的下面。

这个单元格最终包含的数字就是GCGC 和GCCCT 的LCS 的长度。

首先看一下表格的第二行中应该是什么条目。

这一行的单元格保存的LCS 长度对应的是序列GCGCAATA 的零长前端和序列GCCCTAGCG 的LCS。

显然，这些LCS 的值都是0。

类似的，沿着第二列向下的值也都是0，这与递归解的基线条件对应。

现在表格如图3 所示：图3.填充了基线条件的LCS 表格接下来，要实现与递归算法中递归子问题对应的情景，但这时使用的是表格中已经填充的值。

在图4 中，我已经填充了一半左右的单元格：图4.填充了一半的LCS 表格在填充单元格时，需要考虑以下条件：∙它左侧的单元格∙它上面的单元格∙它左上侧的单元格下面三个值分别对应着我在前面列出的三个递归子问题返回的值。

∙V1 = 左侧单元格的值∙V2 = 上面单元格的值∙V3 = 左上侧单元格的值在空单元格中填充下面3 个数字中的最大值：∙V1∙V2∙如果C1等于C2则为V3 + 1，如果C1不等于C2，则为V3，其中C1是当前单元格上面的字符，C2是当前单元格左侧的字符请注意，我在图中还添加了箭头，指向当前单元格值的来源。

后面的“回溯” 一节将用这些箭头建立实际的LCS（与仅仅发现LCS 长度相反）。

现在填充图4 中接下来的空单元格—在GCCCTAGCG 中第三个C 下面和GCGCAATG 第二个C 的右侧的单元格。

它上面的值是2，左侧的值是3，左上侧的值是2。

这个单元格上面的字符和左侧的字符相等（都是C），所以必须选择2、3 和3（左上侧单元格中的2 + 1）的最大值。

所以，这个单元格的值为3。

绘制一个箭头，从该单元格指向其中的值的源单元格。

在这个示例中，新的数值可能来自不止一个单元格，所以可以任选一个：例如左上侧单元格。

作为练习，您可以尝试填充表格的余下部分。

如果在关联过程中，一直按照左上侧-上侧-左侧的顺序选择单元格，那么会得到如图5 所示的表格。

（当然，如果在关联过程中做了不同的选择，那么箭头就会不同，但是数字是相同的。

）图5.填充好的LCS 表格回想一下，任何单元格中的数字都是该单元格所在行之上和列之前的字符串的LCS 长度。

所以，表格右下角的数字就是字符串S1和S2 （在本例中是GCCCTAGCG 和GCGCAATG）的LCS 的长度。

所以，这两个序列的LCS 长度是5。

这是在所有动态编程算法中都要牢记的关键点。

表格中的每个单元格都包含单元格所在行上面和所在列左侧序列前端问题的解。

使用回溯方法寻找实际的LCS接下来要做的就是寻找实际的LCS。

使用单元格箭头进行回溯可以完成。

在构建表格的时候，请记住，如果箭头指向左上侧的单元格，那么当前单元格中的值要比左上侧单元格的值大1，这意味着左侧单元格和上面单元格中的字符相等。

构建LCS 时，这会将相应的字符添加到LCS 中。

所以，构建LCS 的途径就是从右下角的单元格开始，沿着箭头一路返回。

每当沿着对角箭头回到左上角的单元格而且该单元格的值比当前单元格的值小1 时，就要将对应的公共字符添加到正在构建的LCS 的前端。

请注意，之所以将字符放在LCS 前端，是因为我们是从LCS 末端开始的。

（在图5的示例中，右下角的5 与要添加的第5 个字符对应。

）依此类推，继续构建LCS。

从右下侧的单元格开始，看到单元格指针指向左上侧的单元格，而且当前单元格的值（5）比其左上侧单元格的值（4）大1。

所以将字符G 添加到最初的零长度的字符串之前。

下一个箭头，从包含4 的单元格开始，也指向左上侧，但是值没有变。

接着这个箭头也是如此。

下一个单元格的箭头还是指向左上侧，但是这次值从3 变为4。

这意味着需要将这一行和这一列中的公共字符A 添加到LCS 中。

所以现在的LCS 是AG。

接下来，沿着指针向左（对应着跳过上面的T）到达另一个3。

然后有一个对角指针指向2。

因此，又添加了在当前行和当前列中的公共字符C，生成的LCS 为CAG。

继续使用这种方式，直到最后到达0。

图6 显示了整个回溯过程：图6.在填满的LCS 表格上进行回溯通过这种回溯方法，得到的LCS 为GCCAG图5.填充好的LCS 表格回想一下，任何单元格中的数字都是该单元格所在行之上和列之前的字符串的LCS 长度。

所以，表格右下角的数字就是字符串S1和S2 （在本例中是GCCCTAGCG 和GCGCAATG）的LCS 的长度。

所以，这两个序列的LCS 长度是5。

这是在所有动态编程算法中都要牢记的关键点。

表格中的每个单元格都包含单元格所在行上面和所在列左侧序列前端问题的解。

题目这样：求两个字符串的最大公共子序列。

一、什么是最长公共子序列什么是最长公共子序列呢？举个简单的例子吧，一个数列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。

四、代码实现用C++：#include <iostream>#include <string>using namespace std;//求字串x和y的最长子序列void LCSLength(int m,int n,string x,string y,int **c)//m:数组b的行数n:数组b的列数x:第一个字符串y:第二个字符串b:标示数组c:子序列长度数组s:要返回的最大子序列字符串{int i,j;for(i=0;i<=m;i++)for(j=0;j<=n;j++){c[i][j]=0;}for(i=1;i<=m;i++){for(j=1;j<=n;j++){if(x[i-1]==y[j-1]){c[i][j]=c[i-1][j-1]+1;}else if(c[i-1][j]>=c[i][j-1]){c[i][j]=c[i-1][j];}else {c[i][j]=c[i][j-1];}}}}string LCS(int i,int j,string x,int **c,string s){if(i==0||j==0)return "";if(c[i][j]==(c[i-1][j-1]+1)){s=LCS(i-1,j-1,x,c,s);s=s+x[i-1];}else if(c[i][j]==c[i-1][j])s=LCS(i-1,j,x,c,s);else s=LCS(i,j-1,x,c,s);return s;}//定义必要变量，并调用LCS()输出最长子序列string Enter(string s1,string s2,string s)//s1:第一个字符串s2:第二个字符串s:要返回的最大子序列字符串{int l1=s1.length();int l2=s2.length();int **c=new int*[l1+1];for(int i=0;i<=l1;i++)c[i]=new int[l2+1];LCSLength(l1,l2,s1,s2,c);s=LCS(l1,l2,s1,c,s);return s;}void main(){string s1;string s2;string s="\0";cout<<"输入第一个字符串: ";cin>>s1;cout<<"输入第二个字符串: ";cin>>s2;s=Enter(s1,s2,s);cout<<endl<<"最长子序列长度:"<<s.length()<<endl<<"最长子序列："<<s<<endl;cout<<"程序结束！"<<endl;getchar();getchar();}。