10525人阅读
(转自：&）
　　数据挖掘方法的提出，让人们有能力最终认识数据的真正价值，即蕴藏在数据中的信息和知识。数据挖掘&，指的是从大型数据库或数据仓库中提取人们感兴趣的知识，这些知识是隐含的、事先未知的潜在有用信息，数据挖掘是目前国际上，数据库和信息决策领域的最前沿研究方向之一。因此分享一下很久以前做的一个小研究成果。也算是一个简单的数据挖掘处理的例子。
1.数据挖掘与聚类分析概述&
数据挖掘一般由以下几个步骤：&
(l)分析问题源数据数据库必须经过评估确认其是否符合数据挖掘标准。以决定预期结果，也就选择了这项工作的最优算法。&
(2)提取、清洗和校验数据提取的数据放在一个结构上与数据模型兼容的数据库中。以统一的格式清洗那些不一致、不兼容的数据。一旦提取和清理数据后，浏览所创建的模型，以确保所有的数据都已经存在并且完整。
(3)创建和调试模型将算法应用于模型后产生一个结构。浏览所产生的结构中数据，确认它对于源数据中“事实”的准确代表性，这是很重要的一点。虽然可能无法对每一个细节做到这一点，但是通过查看生成的模型，就可能发现重要的特征。&
(4)查询数据挖掘模型的数据一旦建立模型，该数据就可用于决策支持了。&
(5)维护数据挖掘模型数据模型建立好后，初始数据的特征，如有效性，可能发生改变。一些信息的改变会对精度产生很大的影响，因为它的变化影响作为基础的原始模型的性质。因而，维护数据挖掘模型是非常重要的环节。&
　　聚类分析是数据挖掘采用的核心技术，成为该研究领域中一个非常活跃的研究课题。聚类分析基于”物以类聚”的朴素思想，根据事物的特征，对其进行聚类或分类。作为数据挖掘的一个重要研究方向，聚类分析越来越得到人们的关注。聚类的输入是一组没有类别标注的数据，事先可以知道这些数据聚成几簇爪也可以不知道聚成几簇。通过分析这些数据，根据一定的聚类准则，合理划分记录集合，从而使相似的记录被划分到同一个簇中，不相似的数据划分到不同的簇中。
2.特征选择与聚类分析算法
Relief为一系列算法，它包括最早提出的以及后来拓展的和，其中算法是针对目标属性为连续值的回归问题提出的，下面仅介绍一下针对分类问题的和算法。
2.1&Relief算法&
Relief算法最早由提出，最初局限于两类数据的分类问题。算法是一种特征权重算法，根据各个特征和类别的相关性赋予特征不同的权重，权重小于某个阈值的特征将被移除。算法中特征和类别的相关性是基于特征对近距离样本的区分能力。算法从训练集中随机选择一个样本，然后从和同类的样本中寻找最近邻样本，称为，从和不同类的样本中寻找最近邻样本，称为，然后根据以下规则更新每个特征的权重：如果和在某个特征上的距离小于和上的距离，则说明该特征对区分同类和不同类的最近邻是有益的，则增加该特征的权重；反之，如果和在某个特征的距离大于和上的距离，说明该特征对区分同类和不同类的最近邻起负面作用，则降低该特征的权重。以上过程重复次，最后得到各特征的平均权重。特征的权重越大，表示该特征的分类能力越强，反之，表示该特征分类能力越弱。算法的运行时间随着样本的抽样次数和原始特征个数的增加线性增加，因而运行效率非常高。具体算法如下所示：
2.2 ReliefF算法
由于算法比较简单，但运行效率高，并且结果也比较令人满意，因此得到广泛应用，但是其局限性在于只能处理两类别数据，因此年对其进行了扩展，得到了作算法，可以处理多类别问题。该算法用于处理目标属性为连续值的回归问题。算法在处理多类问题时，每次从训练样本集中随机取出一个样本，然后从和同类的样本集中找出的个近邻样本，从每个的不同类的样本集中均找出个近邻样本，然后更新每个特征的权重，如下式所示：
Relief系列算法运行效率高，对数据类型没有限制，属于一种特征权重算法，算法会赋予所有和类别相关性高的特征较高的权重，所以算法的局限性在于不能有效的去除冗余特征。&
2.3 K-means聚类算法
由于聚类算法是给予数据自然上的相似划法，要求得到的聚类是每个聚类内部数据尽可能的相似而聚类之间要尽可能的大差异。所以定义一种尺度来衡量相似度就显得非常重要了。一般来说，有两种定义相似度的方法。第一种方法是定义数据之间的距离，描述的是数据的差异。第二种方法是直接定义数据之间的相似度。下面是几种常见的定义距离的方法：&
1.Euclidean距离，这是一种传统的距离概念，适合于、维空间。&
2.Minkowski距离，是距离的扩展，可以理解为维空间的距离。&
聚类算法有很多种，在需要时可以根据所涉及的数据类型、聚类的目的以及具的应用要求来选择合适的聚类算法。下面介绍&K-means聚类算法:&
K-means算法是一种常用的基于划分的聚类算法。算法是以为参数，把个对象分成个簇，使簇内具有较高的相似度，而簇间的相似度较低。的处理过程为：首先随机选择个对象作为初始的个簇的质心；然后将余对象根据其与各个簇的质心的距离分配到最近的簇；最后重新计算各个簇的质心。不断重复此过程，直到目标函数最小为止。簇的质心由公式下列式子求得：
在具体实现时，为了防止步骤中的条件不成立而出现无限循环，往往定义一个最大迭代次数。尝试找出使平方误差函数值最小的个划分。当数据分布较均匀，且簇与簇之间区别明显时，它的效果较好。面对大规模数据集，该算法是相对可扩展的，并且具有较高的效率。其中，为数据集中对象的数目，为期望得到的簇的数目，为迭代的次数。通常情况下，算法会终止于局部最优解。但用，例如涉及有非数值属性的数据。其次，这种算法要求事先给出要生成的簇的数目，显然这对用户提出了过高的要求，并且由于算法的初始聚类中心是随机选择的，而不同的初始中心对聚类结果有很大的影响。另外，算法不适用于发现非凸面形状的簇，或者大小差别很大的簇，而且它对于噪音和孤立点数据是敏感的。
3.一个医学数据分析实例
3.1&数据说明&
本文实验数据来自著名的UCI机器学习数据库，该数据库有大量的人工智能数据挖掘数据，网址为:http://archive.ics.uci.edu/ml/。该数据库是不断更新的，也接受数据的捐赠。数据库种类涉及生活、工程、科学各个领域，记录数也是从少到多，最多达几十万条。截止2010年底，数据库共有199个数据集，每个数据集合中有不同类型、时间的相关数据。可以根据实际情况进行选用。&
本文选用的数据来类型为：Breast&Cancer&Wisconsin&(Original)&Data&Set，中文名称为：威斯康星州乳腺癌数据集。这些数据来源美国威斯康星大学医院的临床病例报告，每条数据具有11个属性。下载下来的数据文件格式为“.data”,通过使用Excel和Matlab工具将其转换为Matlab默认的数据集保存，方便程序进行调用。
&下表是该数据集的11个属性名称及说明：
对上述数据进行转换后，以及数据说明可知，可以用于特征提取的有9个指标，样品编号和分类只是用于确定分类。本文的数据处理思路是先采用ReliefF特征提取算法计算各个属性的权重，剔除相关性最小的属性，然后采用K-means聚类算法对剩下的属性进行聚类分析。
3.2 数据预处理与程序&
本文在转换数据后，首先进行了预处理，由于本文的数据范围都是，因此不需要归一化，但是数据样本中存在一些不完整，会影响实际的程序运行，经过程序处理，将这一部分数据删除。这些不完整的数据都是由于实际中一些原因没有登记或者遗失的，以“”的形式代表。&
本文采用软件进行编程计算。根据第三章提到的算法过程，先编写函数程序，用来计算特征属性，再编写主程序，在主程序中调用该函数进行计算，并对结果进行分析，绘图，得到有用的结论。
程序统一在最后贴出。
3.3 乳腺癌数据集特征提取&
本文采用节中的算法来计算各个特征的权重，权重小于某个阈值的特征将被移除，针对本文的实际情况，将对权重最小的种剔除。由于算法在运行过程中，会选择随机样本，随机数的不同将导致结果权重有一定的出入，因此本文采取平均的方法，将主程序运行次，然后将结果汇总求出每种权重的平均值。如下所示，列为属性编号，行为每一次的计算结果：&
下面是特征提取算法计算的特征权重趋势图，计算次的结果趋势相同：
上述结果是否运行主程序所得的计算结果，看起来不直观，下面将其按照顺序绘图，可以直观显示各个属性权重的大小分布，如下图所示：
按照从小到大顺序排列，可知，各个属性的权重关系如下：
&&属性属性属性属性属性属性属性属性属性
我们选定权重阀值为，则属性、属性和属性剔除。
从上面的特征权重可以看出，属性裸核大小是最主要的影响因素，说明乳腺癌患者的症状最先表现了裸核大小上，将直接导致裸核大小的变化，其次是属性1和属性8等，后几个属性权重大小接近，但是从多次计算规律来看，还是能够说明其中不同的重要程度，下面是着重对几个重要的属性进行分析。下面是20次测试中，裸核大小（属性6）的权重变化：
从上图中可以看到该属性权重大部分在0.22-0.26左右，是权重最大的一个属性。下面看看属性1的权重分布：
块厚度属性的特征权重在0.19-25左右变动，也是权重极高的一个，说明该特征属性在乳腺癌患者检测指标中是相当重要的一个判断依据。进一步分析显示，在单独对属性6，和属性1进行聚类分析，其成功率就可以达到91.8%。本文将在下节中的Kmeans算法中详细介绍。
3.4 乳腺癌数据集聚类分析&
上一节中通过算法对数据集的分析，可以得到属性权重的重要程度，这些可以对临床诊断有一些参考价值，可以用来对实际案例进行分析，可以尽量的避免错误诊断，并提高诊断的速度和正确率。下面将通过聚类分析算法对数据进行分析。本小节将分为几个步骤来进行对比，确定聚类分析算法的结果以及与算法结合的结果等。
1.K-means算法单独分析数据集&
下面将采用算法单独对数据集进行分析。中已经包括了一些常规数据挖掘的算法，例如本文所用到的算法。该函数名为kmeans，可以对数据集进行聚类分析。首先本文对乳腺癌数据集的所有属性列(除去身份信息和分类列)直接进行分类，由于数据集结果只有2种类型，所以首先进行分2类的测试，结果如下：总体将683条数据分成了2类，总体的正确率为94.44%，其中第一类的正确率为93.56%，第二类的正确率为96.31%。下面是分类后对按照不同属性的绘制的属性值分布图：&
限于篇幅，只选择了上述3个特征属性进行图像绘制，从结果来看，&可以很直观的观察到K-means算法分类后的情况，第一类与第一类的分类界限比较清晰。但是不容易观察到正确和错误的情况。下表是分类结果中各个属性的聚类中心：
从K-means算法的效果来看，能够很准确的将数据集进行分类。一方面是由于该数据集，可能是该案例特征比较明显，另一方面是由于K-menas算法对这种2类的作用较大。&
2.K-means结合分析数据集&
单从分类正确率和结果方面来看，算法已经完全可以对乳腺癌数据集做出非常准确的判断。但是考虑算法对属性权重的影响，本小节将结合算法和算法来对该数据集进行分析，一方面得到处理该问题一些简单的结论，另外一方面可以得到一些对医学处理数据的方法研究方法。&
首先，本小节首先根据节中的一些结论，根据不同属性的权重来对分类数据进行预处理，以得到更精确的结论和对该数据更深度的特征规律。&
从节中，得知属性属性属性属性属性属性属性属性属性，根据算法原理本文可以认为，对于这种属性和属性重要的特征属性，应该对分类起到更加到的作用。所以下面将单独对各个属性的数据进行分类测试，详细结果如下表：
总的分类正确率中，属性最低，属性最高，这与算法测试的结果大致相似，但是由于算法中间部分权重接近，所以也区分不明显。说明特征属性权重的判断对分类是有影响的。上述单独分类中，只将需要分类的列数据取出来，输入到算法中即可。由于输入数据的变化，分类时结果肯定是有差距的，所以单独从一个属性判断其类型是不可靠的。下面选择了单个分类时最高和最低的情况，绘制其分类属性值分布图，如下图所示：
下面将对特征权重按照从大到小的顺序，选择相应的数据，进行聚类分析，结论如下：
1.直接选择全部种属性，分类成功率为：；
2.选择属性，属性，分类成功率为：；
3.选择属性，，，，分类成功率为：；
4.选择属性，，，，，，分类成功率为：；
5.选择属性，，，，，，，，分类成功率为：；
从上面的测试可以看出，选择特征权重最大的个属性，其正确率就达到选择所有属性的情况，因此我们可以认为特征权重最小的几个属性在乳腺癌诊断过程的作用实际可能比较小，实际有可能造成反作用，也就是这几个属性值与乳腺癌没有必然的联系。这一点可以给诊断参考，或者引起注意，进行进一步的研究，确认。&
3.&K-means分成类的情况&
虽然从上述小节的实验中可以得到该数据集的大部分结果和结论。但是为了将相同类型的数据更加准确的分出，下面将尝试分为类的情况。一方面，可以分析在乳腺癌良性和恶性情况下的显著特征属性；另一方面也可以根据此结果找到更加合理的解决方法。&
还是采用中的函数，将分类数改为，由于分为类后数据类型增多，判断较复杂，所以手动对数据进行分析，将所有特征属性加入进去。运行结果如下，测试数据中总共条，其中良性共条，恶性共条：&
1.分为第一类的记录中，良性占；&
2.分为第二类的记录中，恶性占&；&
3.分为第三类的记录中，恶性占&；&
根据上述结果可以认为第一类为良性的分类，第二类为恶性分类，第三类为混合类。对于混合类，说明里面的数据较其他数据更加接近于偏离病例的典型数据，所以进一步分析在第一类中和第二类中的分类正确率：&
1.第一类为良性，共条数据，分类正确率为；
2.第二类为恶性，共条数据，分类正确率为&；&
3.第三类为混合类，共条数据&
因此单独从分类后的正确率来看，效果有提高，说明对典型的病例数据分类更准确，但是对于第三类数据，而无法区分，因此这种情况下，其意义不在于分类的整体正确率，而在于在一些特殊情况下，可以根据一些重要的特征属性值就可以为患者确诊，从而提高效率和准确率，减少误诊断的几率。&
上面是将所有属性进行变换，下面将结合算法，先去掉一部分特征权重较小的特征属性后，再进行处理。根据节中的结论，下面提取权重最大的个属性进行测试，分别是：属性，属性&，属性&，属性&，属性，属性&。&
1.第一类为良性，共条数据，分类正确率为；&
2.第二类为恶性，共条数据，分类正确率为&；&
3.第三类为混合类，共条数据&
因此，对比可以看到，虽然良性的正确率增加了，但是检测出的数据减少了。第三类混合的数量也增多了，说明提出了特种属性较小的属性，可以更加容易区分极端的病例数据，对极端数据的检测更加准确。&
4.主要的Matlab源代码&
1.ReliefF特征提取算法主程序&
1 　　%主函数
3 　　load('matlab.mat')
4 　　D=data(:,2:size(data,2));%
5 　　m =80 ;%抽样次数
6 　　k = 8;
7 　　N=20;%运行次数
8 　　for i =1:N
W(i,:) = ReliefF (D,m,k) ;
10 　　end
11 　　for i = 1:N
%将每次计算的权重进行绘图,绘图N次，看整体效果
plot(1:size(W,2),W(i,:));
14 　　end
15 　　for i = 1:size(W,2)
%计算N次中，每个属性的平均值
result(1,i) = sum(W(:,i))/size(W,1) ;
17 　　end
18 　　xlabel('属性编号');
19 　　ylabel('特征权重');
20 　　title('ReliefF算法计算乳腺癌数据的特征权重');
21 　　axis([1 10 0 0.3])
22 　　%------- 绘制每一种的属性变化趋势
23 　　xlabel('计算次数');
24 　　ylabel('特征权重');
25 　　name =char('块厚度','细胞大小均匀性','细胞形态均匀性','边缘粘附力','单上皮细胞尺寸','裸核','Bland染色质','正常核仁','核分裂');
26 　　name=cellstr(name);
28 　　for i = 1:size(W,2)
plot(1:size(W,1),W(:,i));
xlabel('计算次数') ;
ylabel('特征权重') ;
title([char(name(i))
'(属性' num2Str(i) ')的特征权重变化']);
34 　　end
2.ReliefF函数程序&
1 　　%Relief函数实现
2 　　%D为输入的训练集合,输入集合去掉身份信息项目;k为最近邻样本个数
3 　　function W = ReliefF (D,m,k)
4 　　Rows = size(D,1) ;%样本个数
5 　　Cols = size(D,2) ;%特征熟练,不包括分类列
6 　　type2 = sum((D(:,Cols)==2))/R
7 　　type4 = sum((D(:,Cols)==4))/R
8 　　%先将数据集分为2类，可以加快计算速度
9 　　D1 = zeros(0,Cols) ;%第一类
10 　　D2 = zeros(0,Cols) ;%第二类
11 　　for i = 1:Rows
if D(i,Cols)==2
D1(size(D1,1)+1,:) = D(i,:) ;
elseif D(i,Cols)==4
D2(size(D2,1)+1,:) = D(i,:) ;
17 　　end
18 　　W =zeros(1,Cols-1) ;%初始化特征权重，置0
19 　　for i = 1 : m
%进行m次循环选择操作
%从D中随机选择一个样本R
[R,Dh,Dm] = GetRandSamples(D,D1,D2,k) ;
%更新特征权重值
for j = 1:length(W) %每个特征累计一次，循环
W(1,j)=W(1,j)-sum(Dh(:,j))/(k*m)+sum(Dm(:,j))/(k*m) ;%按照公式更新权重
26 　　end
ReliefF辅助函数寻找最近的样本数
1 %获取随机R 以及找出邻近样本
2 %D：训练集;D1：类别1数据集;D2：类别2数据集;
3 %Dh：与R同类相邻的样本距离;Dm：与R不同类的相邻样本距离
4 function [R,Dh,Dm] = GetRandSamples(D,D1,D2,k)
5 %先产生一个随机数，确定选定的样本R
6 r = ceil(1 + (size(D,1)-1)*rand) ;
7 R=D(r,:); %将第r行选中，赋值给R
8 d1 = zeros(1,0) ;%先置0,d1是与R的距离，是不是同类在下面判断
9 d2 = zeros(1,0) ;%先置0,d2是与R的距离
10 %D1,D2是先传入的参数，在ReliefF函数中已经分类好了
11 for i =1:size(D1,1)
%计算R与D1的距离
d1(1,i) = Distance(R,D1(i,:)) ;
14 for j = 1:size(D2,1)%计算R与D2的距离
d2(1,j) = Distance(R,D2(j,:)) ;
17 [v1,L1] = sort(d1) ;%d1排序，
18 [v2,L2] = sort(d2) ;%d2排序
19 if R(1,size(R,2))==2
%如果R样本=2，是良性
H = D1(L1(1,2:k+1),:) ; %L1中是与R最近的距离的编号，赋值给H。
M = D2(L2(1,1:k),:) ; %v2(1,1:k) ;
H = D1(L1(1,1:k),:);
M = D2(L2(1,2:k+1),:) ;
26 %循环计算每2个样本特征之间的特征距离：(特征1-特征2)/(max-min)
27 for i = 1:size(H,1)
for j =1 :size(H,2)
Dh(i,j) = abs(H(i,j)-R(1,j))/9 ; % 本文数据范围都是1-10，所以max-min=9为固定
Dm(i,j) = abs(M(i,j)-R(1,j))/9 ;
3.K-means算法主程序&
2 　　load('matlab.mat')%加载测试数据
3 　　N0 =1 ;
%从多少列开始的数据进行预测分类
4 　　N1 = size(data,1);%所有数据的行数
5 　　data=data(N0:N1,:);%只选取需要测试的数据
6 　　data1=data(:,[2,3,4,5,6,7,8,9]);% [2,4,7,9]
2:size(data,2)-1
7 　　opts = statset('Display','final');%控制选项
8 　　[idx,ctrs,result,D] = kmeans(data1,2,... %data1为要分类的数据,2为分类的类别数,本文只有2类
'Distance','city',... %选择的距离的计算方式
'Options',opts);
% 控制选项,参考matlab帮助
11 　　t=[data(:,size(data,2)),idx(:,1)];%把测试数据最后一列，也就是分类属性 和 分类结果取出来：列 + 列
12 　　d2 = data(idx==1,11);%提取原始数据中属于第1类的数据的最后一列
13 　　a = sum(d2==2) ;
14 　　b=a/length(d2) ;
15 　　totalSum = 0 ;%总的正确率
16 　　rate1 = 0 ;%第一类的判断正确率.分类类别中数据的正确性
17 　　rate2 = 0 ;%第二类的判断正确率.
18 　　if(b&0.5) %说明第1类属于良性,则a的值就是良性中判断正确的个数
totalSum = totalSum +
rate1 = a/length(d2) ;
%然后加上恶性中判断正确的比例
totalSum = totalSum + sum(data(idx==2,11)==4) ;
rate2 = sum(data(idx==2,11)==4)/length(data(idx==2,11)) ;
24 　　else
%说明第1类属于恶性
totalSum = totalSum + sum(data(idx==1,11)==4) ;
totalSum = totalSum + sum(data(idx==2,11)==2) ;
sum(data(idx==2,11)==2)/length(data(idx==2,11)) ;
sum(data(idx==1,11)==4)/length(data(idx==1,11)) ;
29 　　end
30 　　 x1 =1;%第x1个属性
31 　　x2 =1 ;%第x2个属性
32 　　plot(1:sum(idx==1),data1(idx==1,x1),'r.','MarkerSize',12);
34 　　plot(sum(idx==1)+1:sum(idx==1)+sum(idx==2),data1(idx==2,x1),'b.','MarkerSize',12);
35 　　xlabel('记录数');
36 　　ylabel('属性值');
37 　　title('属性9的值分布');
38 　　legend('第一类','第二类');
39 　　axis([0 640 0 10])
40 　　rate = totalSum/size(t,1)
%总的判断准确率
* 以上用户言论只代表其个人观点，不代表CSDN网站的观点或立场
访问：60006次
排名：千里之外
转载：53篇
(1)(10)(13)(1)(4)(2)(1)(3)(1)(3)(6)(3)(1)(6)(2)Posts - 56,
Articles - 6,
Comments - 35
success belongs to the one with clear and long-term goals!
23:55 by DarkHorse, ... 阅读,
最近看了很多关于LCS(Longest&common&subsequence&problem,最长公共子序列)的文章，大部分问题都只是求出最大公共子序列的长度，或者打印处其中的任意一个最大子序列即可，但是如何快速的打印出所有的最大长度子序列？这个问题好像看到的不多。本文给出了传统的DP(dynamic&programming，动态规划)算法进行求解的过程，并用c语言实现。另外参考一篇论文实现了其中的一种打印所有最大公共子序列的算法，这个算法比起传统的算法而言，时间复杂度大大降低.
一：LCS解析
首先看下什么是子序列？定义就不写了，直接举例一目了然。如对于字符串：&student&，那么su，sud，sudt等都是它的子序列。它可以是连续的也可以不连续出现，如果是连续的出现，比如stud，一般称为子序列串，这里我们只讨论子序列。
什么是公共子序列？很简单，有两个字符串，如果包含共同的子序列，那么这个子序列就被称为公共子序列了。如&student&和&shade&的公共子序列就有&s&或者&sd&或者&sde&等。而其中最长的子序列就是所谓的最长公共子序列(LCS)。当然，最长公共子序列也许不止一个，比如：&ABCBDAB&和&BDCABA&，它们的LCS为&BCBA&，&BCAB&，&BDAB&。知道了这些概念以后就是如何求LCS的问题了。
通常的算法就是动态规划(DP)。假设现在有两个字符串序列：X={x1,x2,...xi...xm},Y={y1,y2,...yj...yn}。如果我们知道了X={x1,x2,...xi-1}和Y={y1,y2,...yj-1}的最大公共子序列L，那么接下来我们可以按递推的方法进行求解：
1)如果xi==yj,那么{L,xi(或yj)}就是新的LCS了，其长度也是len(L)+1。这个好理解，即序列{Xi,Yj}的最优解是由{Xi-1,Yj-1}求得的。
2)如果xiyj,那么可以转换为求两种情况下的LCS。
A:&X={x1,x2,...xi}与Y={y1,y2,...yj-1}的LCS，假设为L1
B:&X={x1,x2,...xi-1}与Y={y1,y2,...yj}的LCS，假设为L2
那么xiyj时的LCS=max{L1,L2},即取最大值。同样，实际上序列{Xi,Yj-1}和{Xi-1,Yj}都可以由{Xi-1,Yj-1}的最优解求得。
怎么样，是不是觉得这种方法很熟悉？当前问题的最优解总是包含了一个同样具有最优解的子问题，这就是典型的DP求解方法。好了，直接给出上面文字描述解法中求LCS长度的公式：
这里用一个二维数组存储LCS的长度信息，i，j分别表示两个字符串序列的下标值。这是求最大公共子序列长度的方法，如果要打印出最大公共子序列怎么办？我们还需要另外一个二维数组来保存求解过程中的路径信息，方便最后进行路径回溯，找到LCS。如果看着很含糊，我下面给出其实现过程。
二：DP实现
很多博文上面都有，基本上是用两个二维数组c[m][n]和b[m][n]，一个用来存储子字符串的LCS长度，一个用来存储路径方向，然后回溯。
其中二维数组b[i][j]的取值为1或2或3，其中取值为1时说明此时xi=yj，c[i][j]=c[i-1][j-1]+1。如果将二维数组看成一个矩阵，那么此时代表了一个从左上角到右下角的路径。如果取值为2，说明xi&yj,且c[i][j]=c[i-1][j]，代表了一个从上到下的路径，同理取值为3代表一个从左到右的路径。
最后我们可以根据c[m][n]的值知道最大公共子序列的长度。然后根据回溯，可以打印一条。其中的坐标点对应的字符同时在两个序列中出现，所以依次回溯这个二维数组就可以找到了。这里给出实现代码：
1 #include &stdio.h&
2 #include &stdlib.h&
3 #include &string.h&
4 #include "stack.h"
5 #define MAX_LEN 1024
6 typedef int **M
8 void GetLCSLen(char *str1, char *str2, Matrix pc, Matrix pb, int nrow, int ncolumn);
9 Matrix GreateMatrix(int nrow, int ncolumn);
10 void DeleteMatrix(Matrix p, int nrow, int ncolumn);
11 void TraceBack(char *str1, Matrix pb, int nrow, int ncolumn);
13 void GetLCSLen(char *str1, char *str2, Matrix pc, Matrix pb, int nrow, int ncolumn)
/************initial the edge***************/
for(i=0; i& i++)
pc[i][0] = 0;
pb[i][0] = 0;
for(j=0; j& j++)
pc[0][j] = 0;
pb[0][j] = 0;
/************DP*****************************/
for(i=1; i& i++)
for(j=1; j& j++)
if(str1[i-1] == str2[j-1])
pc[i][j] = pc[i-1][j-1] + 1;//由左上节点转移而来
pb[i][j] = 1;//标记为1
else if(pc[i-1][j] &= pc[i][j-1])
pc[i][j] = pc[i-1][j];//由上节点转移而来
pb[i][j] = 2;//标记为2
pc[i][j] = pc[i][j-1];//由左节点转移而来
pb[i][j] = 3;//标记为2
50 void TraceBack(char *str1, Matrix pb, int nrow, int ncolumn)
if(str1 == NULL || pb == NULL)
if(nrow == 0 || ncolumn == 0)
ntemp = pb[nrow-1][ncolumn-1];
switch(ntemp)
printf("locate:(%d,%d),%4c\n", nrow-1, ncolumn-1, str1[nrow-2]);//打印公共字符，这里下标是nrow-2，因为矩阵的坐标值(i,j)比字符串的实际下标大1
TraceBack(str1, pb, nrow-1, ncolumn-1);//向左上角递归
TraceBack(str1, pb, nrow-1, ncolumn);//向上方向递归
TraceBack(str1, pb, nrow, ncolumn-1);//向左方向递归
我们给出一个测试：
char str1[MAX_LEN] = "BADCDCBA";
char str2[MAX_LEN] = "ABCDCDAB";
GetLCSLen(str1, str2, C, B, str1len+1, str2len+1);
TraceBack(str1, B, str1len+1, str2len+1);
详细的代码见文章结束处给出的链接。本测试output:BDCDB
问题：上面的方法中，我们没有单独考虑的情况，所以在回溯的时候打印的字符只是其中一条最大公共子序列，如果存在多条公共子序列的情况下。怎么解决？我们对二维数组的取值添加一种可能，等于，这代表了我们说的这种多支情况，那么回溯的时候我们可以根据这个信息打印更多可能的选择。这个过程就不写代码了，其实很简单，以下面的路径回溯图举例，你从（8,8)点开始按b[i][j]的值指示的方向回溯，把所有的路径遍历一遍，如果是能达到起点（0,0)的路径，就是LCS了，有多少条打印多少条。可是，
又出现问题了：你发现没有，在回溯路径的时候，如果采用一般的全搜索，会进行了很多无用功。即重复了很多，且会遍历了一些无效路径，因为这些路径最终不会到达终点，比如节点。因此加大算法复杂度和时间消耗。那么如何解决？看下面的这个方法，正式进入本文正题。
&路径回溯图：
&加入状态4后的状态图：
三：算法改进
上面提到路径回溯过程中，一般的方法就是遍历所有可能的路径，但是一些不可能构成最大公共子序列的跳跃点我们也会去计算。这里先解释下什么叫跳跃点，就是导致公共子序列长度发生变化的节点，即b[i][j]=1对应的节点(i,j)。Ok，接下来的问题是，如何不去考虑这些无效跳跃点，降低算法复杂度？参考论文里提出了这样一种方法：矩形搜索。
首先构造两个栈数据结构store和print。故名思议，一个用来储存节点，一个用来打印节点。栈的定义为：
1 #define MAX_STACK_SIZE 1024
2 typedef struct _Element
8 typedef struct _Stack
Element data[MAX_STACK_SIZE];
栈使用数组实现，并有一个指向顶点的下标值top。为了初始化需要，先构造了一个虚拟节点virtualnode，指向节点(m,n)的右下角，即(m+1,n+1),这个节点的LCS的长度假设为最大公共子序列长度+1。将虚拟节点压入栈store，然后然后执行下面的算法：
1)栈store为空吗？是的话退出。
2)否则从store栈顶弹出节点。
3)如果这个节点为边界元素(行或列的小标为1),则将此节点压入栈print，打印栈print里面的所有节点(除virtualnode外)。查看此时store栈顶节点的LCS长度，并以这个长度为参考，弹出print栈里面所有LCS长度小于等于这个参考值的节点，跳转到第1步。
4)如果不是边界元素，那么以该节点的左上节点(i-1,j-1)为出发点，沿着该出发点指示的方向，找到第一个跳跃点e1(即b[i][j]==1的点)。途中碰到分支节点(b[i][j]==4的点)时，沿着其上节点继续探索。
5)找到第一个跳跃点以后，重新回到第4步的出发点，沿着该节点指示的方向，找到第二个跳跃点e2。途中碰到分支节点(b[i][j]==4的点)时，沿着其左节点继续探索
6)如果e1和e2节点为同一节点，将该节点压入store栈，回到步骤1)。
7）如果不为同一节点，在e1和e2构成的矩形范围内，搜索出所有的跳跃点，并全部压入store栈，回到步骤1)。
不明白？不要紧，我们结合上面的矩阵图一步步按照算法来，看看到底是如何计算的：
第一步：压入虚拟节点6(9,9)到栈store，这里6表示这个节点的LCS长度，（9,9)表示坐标值。&
第二步：store栈不为空，则弹出store栈顶，压入print栈，这时候的两个栈的状态如下面的左图。沿出发点(8,8)出发，这是个分支节点，因为b[8][8]==4,所以选择向上走，搜索到e1跳跃点(7,8)，搜索路径为：(8,8)-&(7,8)。然后回到(8,8)找e2点，这时选择向左走，找到e2跳跃点(8,7)。这两个跳跃点不同，所以以e1，e2为对角线方向围成的矩形内搜索所有跳跃点，这里只有e1,和e2本身两个节点，然后将它们压入栈store。此时两个栈的状态见下面的有图。蓝色底的节点表示有store栈弹出然后压到print栈，绿色底表示新压入到store栈的跳跃点，下面所有的图都这样表示。
第三步：弹出5(7,8)到print栈，搜索到新的两跳跃节点。
第七步：关键步骤来了，因为此时从store栈弹出的节点是边界元素1(1,2)，所以我们打印print栈的所有元素(红色字体节点)，而这些元素恰好构成了一个最长公共子序列(自己揣摩一下)。打印完了以后，我们要对print栈进行清理，去除不需要的节点，按照步骤2，此时store栈顶节点的LCS为1，所以print栈中弹出的节点只有1(1,2)。弹完以后，print栈的状态如图所示。虚线框节点表示已弹出，下同。
第八步：继续弹出store栈顶，发现又是边界元素，继续压入print栈，打印print栈。清理print栈。
第九步：清理完后，继续步骤2.
好了，下面的过程就是重复进行上面的这些步骤了，自己动手画一下就一目了然了。
四：代码实现
用c语言实现关键代码：
1 void TraceBack2(char *str1, Matrix pc, Matrix pb, int nrow, int ncolumn);
2 void PrintStack(Stack *ps, char *str1, int len1);
3 void SearchE(Matrix pb, int curposx, int curposy, int *eposx, int *eposy, int ntype);
5 void PrintStack(Stack *ps, char *str1, int len1)
if(ps == NULL || str1 == NULL)
int ntemp = ps-&
int index = -1;
while((index = ps-&data[ntemp].nrow) &= len1)
printf("%2c",str1[index-1]);
printf("\n");
19 void TraceBack2(char *str1, Matrix pc, Matrix pb, int nrow, int ncolumn)
if(str1 == NULL || pc == NULL || pb == NULL)
Stack store,//构造两个栈store，print
E//store栈的栈顶节点
E//临时变量
E//虚拟节点
int//保存store栈顶节点的LCS长度
int ex1,ey1,ex2,ey2;//矩形搜索的两个节点的坐标
InitialStack(&store);//初始化
InitialStack(&print);
virtualnode = CreateElement(pc[nrow-1][ncolumn-1]+1, nrow, ncolumn);
Push(&store, &virtualnode);//压入虚拟节点到store
while(!IsEmpty(&store))
Pop(&store, &storetop);//从栈顶取出一个节点
if(storetop.nrow == 1 || storetop.ncolumn == 1)//如果是边界节点
Push(&print, &storetop);
PrintStack(&print, str1, nrow-1);//打印print栈里面除虚拟节点之外的所有节点
GetTop(&store, &element);
ntoplen = element.//当前store的栈顶节点的LCS长度
/**********弹出print栈中所有LCS长度小于等于ntoplen的节点**************/
while(GetTop(&print, &element) && element.nlcslen&=ntoplen)
Pop(&print, &element);
Push(&print, &storetop);
SearchE(pb, storetop.nrow-1, storetop.ncolumn-1, &ex1, &ey1, 0);
SearchE(pb, storetop.nrow-1, storetop.ncolumn-1, &ex2, &ey2, 1/*also other value is ok*/);
if(ex1 == ex2 && ey1 ==ey2)
element = CreateElement(pc[ex1][ey1], ex1, ey1);
Push(&store,&element);//压入store栈，回到步骤2
for(i=ey2; i&=ey1; i++)
for(j=ex1; j&=ex2; j++)
if(pb[i][j] == 1)
element = CreateElement(pc[i][j], i, j);
Push(&store, &element);
80 void SearchE(Matrix pb, int curposx, int curposy, int *eposx, int *eposy, int ntype)
switch(pb[curposx][curposy])
SearchE(pb, curposx-1, curposy, eposx, eposy, ntype);
SearchE(pb, curposx, curposy-1, eposx, eposy, ntype);
if(ntype == 0)
SearchE(pb, curposx-1, curposy, eposx, eposy, ntype);//搜索e1点，如过碰到分叉点，向上继续搜索
SearchE(pb, curposx, curposy-1, eposx, eposy, ntype);//搜索e2点，如过碰到分叉点，向左继续搜索
同样的测试，这种算法能打印出全部的最长公共子序列。
该算法能将传统算法的指数级复杂度降低到max{O(mn),O(ck)},k为最大公共子序列的个数。详细证明见论文。因为用两个栈存储了所有有效跳跃点，使得许多重复比较被忽略。栈的有序性也能巧妙的得到任意一条最大公共子序列。
本算法的全部实现代码下载路径：
欢迎讨论。
参考论文：
《利用矩阵搜索求所有最长公共子序列的算法》，宫洁卿，安徽工程科技学院学报，vol23,No.4,Dec.,2008}