数据的多流形结构分析

参赛密码（由组委会填写）全第十二届“中关村青联杯”全国研究生数学建模竞赛参赛密码（由组委会填写）第十二届“中关村青联杯”全国研究生数学建模竞赛题目数据的多流形结构分析摘要在这个信息爆炸的时代，海量的数据不断产生，迫切需要对这些大数据进行有效的分析，以至数据的分析和处理方法成为了诸多问题成功解决的关键，涌现出了大量的数据分析方法。本文是一个研究数据多流形结构分析的问题，以解决图像数据聚类、运动分割、人脸识别等问题。针对问题一，附件1给出了一组多维数据，它采样于两个独立的子空间，运用K-means算法对其进行分类。利用MATLAB软件将这组属于两个独立子空间的数据分成了两类，并用表格的形式表示出了数据分类，其中“1”表示第一类，“2”表示第二类。针对问题二，采用稀疏子空间聚类的方法对(a)图和(b)图进行聚类，采用多流形聚类的方法对(c)图和(d)图进行聚类。对(a)图数据进行处理，用稀疏矩阵表示，然后建立邻接矩阵，最后进行谱聚类分析，通过调整稀疏矩阵的正则参数来改变稀疏矩阵，比较聚类错误率得到最优聚类效果，最后将其分为两类，其聚类缺失率为5.3%；图(b)同样采用(a)的方法，调整正则参数，最终将其分为三类，两条直线分别表示一类，平面表示一类，其聚类缺失率为8.7%；(c)图为两条不相交的二次曲线，因欧氏距离的相似性度量不能完全反映这种空间分布特性，通过欧式距离求得流形距离，根椐流形距离核的相似性度量构造相似性矩阵，再求得拉普拉斯矩阵，继而得到拉普拉斯矩阵的最大特征值对应的特征向量，对其构建矩阵并单位化，最后采用K-means算法将其分为2类，其聚类缺失率为7.3%；图(d)为两条相交的螺旋线，同样采用(c)的算法，将其分为两类，其聚类缺失率为7.1%。针对问题三，采用低秩稀疏子空间聚类方法，对图(a)进行聚类，将数据映射到一个低维的隐空间中，同时在这个低维空间中寻求数据的低秩稀疏系数，优化线性变换矩阵得到最优解，最后利用谱聚类算法，将其分为了“横”和“竖”两类。对于图(b)运动分割，基于特征点的轨迹，使用拉普拉斯特征映射降低维数的方法，将场景中不同运动对应的不同特征点轨迹分割出来，使得同一运动的特征点轨迹在同一个线性流形上，之后采用(a)中的算法，将此特征点轨迹分成三类，并用表格表示分类结果。对于(c)问题，也采用(b)的方法，降低维数到100维，将这20幅图像分成两类。针对问题四，采取基于拓扑距离度量的流形聚类的方法，对图(a)进行聚类，使数据集中的每个数据点产生符合正态分布函数的信号，将所有数据点产生的信号叠加起来，形成一条信号曲面，借助流形曲线来计算两两数据点间的流形距离，即点与点的流形距离为曲线的的弧长和平均曲率之和，基于拓扑距离构造数据集合的相似度矩阵，再使用规范化分割聚类算法(Neut)进行流形聚类，将圆台分成了三类，分别是圆台的顶、底、侧面。(b)是机器工件外部边缘轮廓的图像，采用(a)中方法，将轮廓线分别聚类成两类、三类与四类，因未说明分成的类数，即为无监督的聚类，通过计算每种分类的类间偏差平方和和类内偏差平方和之间的差的绝对值，用此来评价分类结果的优劣，绝对值越小分类效果越好，比较发现，将(b)分成两类，分类效果最佳。最后，本文通过对分类结果进行类间偏差平方和和类内偏差平方和之间的差的绝对值的评价，可以说明分类的结果是可靠的。关键词：K-means谱聚类低秩稀疏子空间聚类拉普拉斯特征映射流形距离拓扑距离目录7420目录 330113一、问题的重述 4269561.1问题的背景 4154041.2问题的提出 414921.3本文所需解决的问题 415921二、符号说明 5141三、问题分析 5135173.1问题一 5115193.2问题二 5124853.3问题三 578853.4问题四 527993四、模型的建立与求解 6270664.1问题一模型的建立与求解 6262964.1.1K-means算法 6232444.1.2模型的求解 6146224.2问题二模型的建立与求解 8152114.2.1稀疏子空间聚类与低秩子空间聚类 8137994.2.2问题二模型的建立 8303964.2.3运行结果 9179374.2.4结果分析 11313644.3问题三模型的建立与求解 12117864.3.1问题三模型的建立 1214411低秩稀疏子空间聚类方法 129407LaplacianEigenmaps降维法 1242654.3.2运行结果 13684.4问题四模型的建立与求解 154774.4.1问题四模型的建立 153207规范化分割聚类算法 1531765基于数据拓扑结构的拓扑距离 161272基于拓扑距离度量的流形聚类 1743774.4.2运行结果 18122954.4.3结果分析 208762五、模型的评价 20214715.1模型的优点 2017045.2模型的缺点 2013780六、参考文献 2113735七、附录 22一、问题的重述1.1问题的背景我们已经进入了一个信息爆炸的时代，海量的数据不断产生，迫切需要对这些大数据进行有效的分析，以至数据的分析和处理方法成为了诸多问题成功解决的关键，涌现出了大量的数据分析方法。几何结构分析是进行数据处理的重要基础,已经被广泛应用在人脸识别、手写体数字识别、图像分类、等模式识别和数据分类问题，以及图象分割、运动分割等计算机视觉问题（人脸识别、图像分类、运动分割等实例见下文）中。更一般地，对于高维数据的相关性分析、聚类分析等基本问题，结构分析也格外重要。1.2问题的提出有些实际问题的数据并不符合混合子空间结构的假设，所以假设数据的结构为混合多流形更具有一般性。由于混合流形不全是子空间的情况，数据往往具有更复杂的结构，分析这种数据具有更大的挑战性。基于谱聚类的方法仍然是处理该类问题的流行方法。1.3本文所需解决的问题(1)当子空间独立时，子空间聚类问题相对容易。附件一中1.mat中有一组高维数据（.mat所存矩阵的每列为一个数据点，以下各题均如此），它采样于两个独立的子空间。我们需要将该组数据分成两类。(2)附件二中是四个低维空间中的子空间聚类问题和多流形聚类问题，如图1所示。图1(a)为两条交点不在原点且互相垂直的两条直线，需要将其分为两类；图1(b)为一个平面和两条直线，这是一个不满足独立子空间的关系的例子，要求我们将其分为三类。图1(c)为两条不相交的二次曲线，需将其分为两类。图1(d)为两条相交的螺旋线，也是需要将其分为两类。(3)问题三需要我们解决以下三个实际应用中的子空间聚类问题，数据见附件三(a)受实际条件的制约，在工业测量中往往需要非接触测量的方式，视觉重建是一类重要的非接触测量方法。特征提取是视觉重建的一个关键环节，如图2(a)所示，其中十字便是特征提取环节中处理得到的，十字上的点的位置信息已经提取出来，为了确定十字的中心位置，一个可行的方法是先将十字中的点按照“横”和“竖”分两类。要求我们使用适当的方法将图2(a)中十字上的点分成两类。(b)运动分割是将视频中有着不同运动的物体分开，是动态场景的理解和重构中是不可缺少的一步。基于特征点轨迹的方法是重要的一类运动分割方法，该方法首先利用标准的追踪方法提取视频中不同运动物体的特征点轨迹，之后把场景中不同运动对应的不同特征点轨迹分割出来。已经有文献指出同一运动的特征点轨迹在同一个线性流形上。图2(b)显示了视频中的一帧，有三个不同运动的特征点轨迹被提取出来保存在了3b.mat文件中，我们所需要做的就是，使用适当方法将这些特征点轨迹分成三类。(c)3c.mat中的数据为两个人在不同光照下的人脸图像共20幅（X变量的每一列为拉成向量的一幅人脸图像），需要将这20幅图像分成两类。(4)请作答如下两个实际应用中的多流形聚类问题图3(a)分别显示了圆台的点云，要求我们将这些点按照其所在的面分开（即圆台按照圆台的顶、底、侧面分成三类）。图3(b)是机器工件外部边缘轮廓的图像，我们需要将轮廓线中不同的直线和圆弧分类。二、符号说明符号符号说明d两点之间欧氏距离表示第一个点的第i维坐标表示第二个点的第i维坐标伸缩因子节点和之间的最短路径距离连接顶点顶点i和j之间的相关性Minkowski度量来定义两个数据向量的距离三、问题分析3.1问题一独立子空间的聚类问题相对容易，采用k-means算法对附件一中的数据进行聚类，这组多维数据是属于两个独立的子空间，我们利用MATLAB软件将这组高维数据分成了两类，并用表格的形式表示出第一类用1表示，第二类用2表示。3.2问题二问题二的数据为四个二维数据，涉及到的问题是子空间聚类问题和多流形子空间聚类问题，采用稀疏子空间聚类的方法对(a)图和(b)图进行聚类，用稀疏矩阵表示，然后建立邻接矩阵，最后进行谱聚类分析，然后通过改变稀疏矩阵的正则参数，来调整分类的效果，利用聚类错误率来辨别聚类效果的优劣，最终将其分为两类，每条直线分别表示一类。对于图(c)，(d)，因为欧氏距离的相似性度量不能完全反映数据聚类复杂的空间分布特性的问题。(c)图为两条不相交的二次曲线，因欧氏距离的相似性度量不能完全反映这种空间分布特性，通过欧式距离求得流形距离，根椐流形距离核的相似性度量构造相似性矩阵，再求得拉普拉斯矩阵，继而得到拉普拉斯矩阵的最大特征值对应的特征向量，对其构建矩阵并单位化，最后采用K-means算法将其分为2类。图(d)为两条相交的螺旋线，同样采用(c)的算法，将其分为两类。3.3问题三解决以下三个实际应用中的子空间聚类问题，对于图(a)采用低秩稀疏子空间聚类方法，将数据映射到一个低维的隐空间中，同时在这个低维空间中寻求数据的低秩稀疏系数，优化线性变换矩阵得到最优解，最后利用谱聚类算法，将其分为了“横”和“竖”两类。对于图(b)与(c)，数据的维度较高，考虑到计算的复杂度与时间，以及数据不丢失的情况下，对数据进行降低维数。对于图(b)运动分割，基于特征点的轨迹，使用拉普拉斯特征映射降低维数的方法，将场景中不同运动对应的不同特征点轨迹分割出来，使得同一运动的特征点轨迹在同一个线性流形上，之后采用(a)中的算法，将此特征点轨迹分成三类，并用表格表示分类结果。对于(c)问题，也采用(b)的方法，降低维数到100维，将这20幅图像分成两类。3.4问题四针对问题四，采取基于拓扑距离度量的流形聚类的方法，对图(a)进行聚类，使数据集中的每个数据点产生符合正态分布函数的信号，将所有数据点产生的信号叠加起来，形成一条信号曲面，借助流形曲线来计算两两数据点间的流形距离，即点与点的流形距离为曲线的的弧长和平均曲率之和，基于拓扑距离构造数据集合的相似度矩阵，再使用规范化分割聚类算法(Neut)进行流形聚类，根据图像数据构造一个带权重的无向图，计算最小的几个特征值及其对应的特征向量。利用已经得到的谱的信息将图一分为二。判断分割后的图像是否还需要再次分割，如需要分割，再对分割后的两部分图像分别重复上述过程进行分割。最后将圆台分成了三类，分别是圆台的顶、底、侧面。(b)是机器工件外部边缘轮廓的图像，采用(a)中方法，将轮廓线分别聚类成两类、三类与四类，并计算每种分类的类间偏差平方和和类内偏差平方和之间的差的绝对值，用来评价分类结果的优劣，绝对值越小分类效果越好。四、模型的建立与求解4.1问题一模型的建立与求解4.1.1K-means算法K-means算法是最简单的一种聚类算法。算法的目的是使各个样本与所在类均值的误差平方和达到最小（这也是评价K-means算法最后聚类效果的评价标准）。K-means算法的一般步骤：1)初始化。先随机选取K个对象作为初始的聚类中心。然后计算每个对象与各个种子聚类中心之间的距离，把每个对象分配给距离它最近的聚类中心。2)进行迭代。聚类中心以及分配给它们的对象就代表一个聚类。一旦全部对象都被分配了，每个聚类的聚类中心会根据聚类中现有的对象被重新计算。3)更新聚类中心。这个过程将不断重复直到满足某个终止条件。终止条件可以是以下任何一个：i)没有（或最小数目）对象被重新分配给不同的聚类。ii)没有（或最小数目）聚类中心再发生变化。iii)误差平方和局部最小。在这里我们用的是欧式距离来描述这两个点之间的远近，通过距离来反映他们的相似程度，距离越近，相似程度越高[1]。(1)4.1.2模型的求解我们使用了K-means聚类算法对附件一中的高维数据进行分类，得出如下的结果：表1独立子空间聚类结果编号1234567891011121314151617181920类别11111111111111111111编号2122232425262728293031323334353637383940类别11111111111111111111编号4142434445464748495051525354555657585960类别22222222222222222222编号6162636465666768697071727374757677787980类别22222222222222222222编号81828384858687888990919293949596979899100类别22222222222222222222编号101102103104105106107108109110111112113114115116117118119120类别22222222222222222222编号121122123124125126127128129130131132133134135136137138139140类别22222222222222222222编号141142143144145146147148149150151152153154155156157158159160类别11111111111111111111编号161162163164165166167168169170171172173174175176177178179180类别11111111111111111111编号181182183184185186187188189190191192193194195196197198199200类别11111111111111111111注：图中的类别1表示第一类，2表示第二类。4.2问题二模型的建立与求解4.2.1稀疏子空间聚类与低秩子空间聚类过去的几十年人们见证了数据的爆炸式增长，这对于数据的处理工作提出了巨大的挑战，特别是这些数据集通常都是高维数据．数据的高维特性不仅增加了计算时间，而且由于噪声和环境空间降低了算法的性能．实际上，这些数据的内在尺度往往比实际空间中小得多，这就促使人们运用一些技术发现高维数据的低维表示，比如低秩近似和稀疏表示等．实际上，在许多问题中，高维空间中的数据往往可以用低维子空间进行表示。计算过程是，首先通过对角稀疏矩阵，求得其相似矩阵，然后对这个相似矩阵进行归一化，然后通过谱聚类的三种聚类方法Unnormalizedmethod、Randomwalk、Normalizedsymmetric进行聚类，最后再用K-means的方法对其进行最终的聚类，聚类结束后再对分类结果的好坏进行判断，用聚类错误率作为评判聚类好坏的标准。稀疏子空间聚类(SSC)[2]，每一个数据点可以表示为其他数据点的稀疏线性组合，通过这些稀疏系数构造亲和矩阵进行子空间聚类。也就是说，给定一个数据集X，希望找到一个系数矩阵C，满足X=XC并且diag(c)=0。可以通过求解(2)式得到解。,(2)当数据集被噪声G污染时，SSC算法假设每个数据点可以表示为X=XC+G，可以通过求解凸优化问题(3)得到解。,(3)低秩表示(LRR)算法[3]和稀疏子空间聚类算法非常类似，区别在于LRR算法的目标是寻找数据的低秩表示，而SSC算法在于寻找数据的稀疏表示。LRR通过求解凸优化问题(4)得到解。(4)当数据集被噪声G污染时，LRR通过求解凸优化问题(5)得到解．(5)最后，通过得到的稀疏矩阵（利用SSC或者LRR），构造亲和矩阵C+CT，在这个亲和矩阵上利用谱聚类算法，就可以得到最终的聚类结果[4]。4.2.2问题二模型的建立开始MATLAB程序流程开始读取附件中的数据读取附件中的数据构造稀疏矩阵构造稀疏矩阵构建点的邻接距离向量矩阵构建点的邻接距离向量矩阵谱聚类谱聚类图1(a)为两条交点不在原点且互相垂直的两条直线，将其分为两类，通过上面的方法，用MATLAB进行编程得到很好的分类结果，将两条线各自分为一类。图1(b)为一个平面和两条直线，这是一个不满足独立子空间的关系的例子，将其分为三类。在MATLAB程序中，(a)、(b)图都是通过改变正则化参数，改变稀疏矩阵，从而来得到更好的分类效果，同类样本的相似度构造的较大，不同类的较小，这类方法一般能得到正确分类结果。针对问题二图(c)和(d)，因为图(c)是从图(a)的两条直线变为了两条不相交的二次曲线，图(d)是两条螺旋线，基于欧氏距离的相似性度量不能完全反映数据聚类复杂的空间分布特性的问题，提出了一种基于流形距离核的谱聚类算法．它能充分挖掘数据集中的内在结构信息，较好地反映局部和全局一致性，并且可以很好地防止噪声点的影响。它可以增大位于同一流形结构的数据对间的相似度，并且减小位于不同流形结构的数据对间的相似度，从而解决了复杂数据的空间分布特性。流形上的线段长度用下面的公式(6)来计算：(6)其中，d(x,y)为数据点x和y间的欧氏距离；称为伸缩因子，可以用来描述聚类的全局一致性，通过调节伸缩因子来放大或缩短两点间线段长度。流形距离测度：将数据点看作是一个加权无向图中的顶点V，边集合E表示的是在每一对数据点间定义的连接权值。令表示图上长度为的连接点与之间的路径，其中边，。令表示连接数据点对的所有路径的集合，与之间的流形距离为：(7)(8)其中，是图G上节点和之间的最短路径距离。是图G上节点到最短路径上任意相邻两点的欧氏距离。不难看出，该距离的定义满足距离测度定义的4个约束条件：自反性：，当且仅当；对称性：；非负性：；三角不等式：．该距离测度可以度量流形上的最短路径，能够很好地反映数据集内在的流形结构，使得位于同一流形上的2点可以用许多短边相连，而位于不同流形上的2点则需要用穿过低密度区域的长边相连，最终达到放大位于不同流形上的数据点间距离，缩短位于同一流形上的数据点间距离的目的．从上面的定义可以看出，该距离测度可以反映数据的局部密度特征。4.2.3运行结果图1(a)分类结果如下图所示：图1(a)的分类结果图注：其中红色的线为一类，蓝色的线为一类。图1(b)分类结果如下图所示：图2(b)的分类结果图注：其中红色的线为一类，蓝色的线为一类，绿色的平面为一类。图1(c)为两条不相交的二次曲线，它是属于多流形聚类的问题，通过MATLAB编程，我们将其分为两类，分类结果如图3所示：图3(c)的分类结果图注：其中红色的曲线为一类，蓝色的曲线为一类。图1(d)为两条相交的螺旋线，这两条螺旋线也是属于多流形聚类的问题，将其分为两类，结果如图4所示：图3(d)的分类结果图注：其中红色的螺旋线为一类，蓝色的螺旋线为一类。4.2.4结果分析我们采用聚类错误率来评价聚类算法的性能：对于图1(a)的聚类结果，通过计算，得到数据的聚类缺失率为0.053，错误率比较低，这说明用稀疏子空间聚类的方法对数据进行分类是有效的。对于图1(b)的聚类结果，用同样的计算方法，得到该组数据的聚类缺失率为0.087，这说明用稀疏子空间聚类的方法对数据进行分类的效果比较好。对于图1(c)的聚类结果，聚类缺失率为0.073，这说明用基于多流形的聚类方法对这组数据进行分类达到的效果比较好。对于图1(d)的聚类结果，聚类缺失率为0.071，这说明用稀疏子空间聚类的方法对数据进行分类的效果比较好。4.3问题三模型的建立与求解4.3.1问题三模型的建立低秩稀疏子空间聚类方法针对问题三图2(a)提出了低秩稀疏子空间聚类算法。本算法是在稀疏子空间聚类和低秩表示的基础上，在聚类的过程中可以对高维数据进行降维，同时在低维空间中利用稀疏表示和低秩表示对数据进行聚类．在运动分割和人脸聚类上，该算法都具有很好的聚类性能。稀疏子空间聚类算法和低秩表示，文中将数据映射到一个低维的空间中，同时在这个低维空间中寻求数据的低秩稀疏系数。令为一个线性变换矩阵，它将数据从原始空间映射到一个维数为t的空间中．通过目标函数的最小化，可以同时得到变换矩阵和数据集的低秩稀疏系数：(9)其中，(10)(10)式的第一项为求取数据集的低秩系数；第二项为求取数据集的稀疏系数；第三项的主要目的是去除噪声影响；最后一项是似于PCA的正则项，主要目的是保证映射变换不能过多丢失一些原始空间的信息；和为非负常数。另外，要求P正交并且归一化,这样就避免了解的退化，并且保证优化方法的计算效率。以注意到，(10)式是能够进行扩展的，这样就可以对位于仿射子空间中的数据进行处理．可以对优化问题(9)增加一个约束条件得到：(11)由于是实际应用中的子空间的聚类，受实际条件的制约，图(a)是视觉重建，所以关键环节是特征的提取，我们通过低秩稀疏子空间聚类的算法可以将十字分为两类，分别为“横”和“竖”两类。针对问题三图2(b)，我们使用拉普拉斯特征映射(LaplacianEigenmaps)降维的方法，LE假设每一点只与它距离最近的一些点相似，再远一些的数据相似程度为0，降维后相近的点尽可能保持相近。LaplacianEigenmaps降维法LaplacianEigenmaps降维的方法，首先要构建图的权值矩阵。构建方法为：步骤11)通过设置一个阈值，相似度在阈值以下的都直接置为零，这相当于在一个领域内考虑局部性；2)对每个点选取k个最接近的点作为邻居，与其他的点的相似性则置为零。3)全连通。步骤2LE把降维考虑为一种高维到低维的映射，则一个好的映射应该使连通的点靠的更近。设映射后的点为，则LE的目标就是最小化以下函数：(12)对目标函数进行整理(13)其中，L=D-W，L就叫做LaplacianMatrix。(14)于是，最小化问题等效为：(15)(16)因此，对于y，我们将特征值从小到大排序后，选择第二小的特征值到第m+1小的特征值对应的特征向量（共m个），组成一个的矩阵，矩阵的每一行就是原来的数据降维后得到的m维特征。原本我们只知道数据与数据之间的相似程度，结果把降维后的特征求出来这里的特征值。如果仅有一个特征值为0，那么这个图一定是全通的。通过降维的方法把(b)图的297维降到100维，这个维度是在不失去很多数据，又可以简化算法复杂度的基础上，得到的最终维度，降维后得到稀疏矩阵，然后接下来的处理的方法和问题二的一样。针对问题三图2(c)，同样使用拉普拉斯特征映射降维的方法将3c.mat中的数据20幅两个人在不同光照下的人脸图像分成两类（X变量的每一列为拉成向量的一幅人脸图像），使用相同的降维方法，将2000维数据降到100维，再通过谱聚类的算法将这20幅人脸分成两类，通过降维的方法我们能够得到很好的结果。4.3.2运行结果通过运行程序，图(a)的数据分成了两类，分类结果如下图所示：图4(a)的分类结果图注：其中红色的“横”为一类，蓝色的“竖”为一类。通过MATLAB软件，运行程序，图(b)的数据分成了三类，分类结果如下表所示：表2问题三图（b)的分类结果编号1234567891011121314151617181920类别22222222222222222222编号2122232425262728293031323334353637383940类别22222222222222222222编号4142434445464748495051525354555657585960类别22222222222222222222编号6162636465666768697071727374757677787980类别22222222222222222222编号81828384858687888990919293949596979899100类别22222222222222222222编号101102103104105106107108109110111112113114115116117118119120类别22222222222222222233编号121122123124125126127128129130131132133134135136137138139140类别33333333333333333333编号141142143144145146147148149150151152153154155156157158159160类别33333333333333133333编号161162163164165166167168169170171172173174175176177178179180类别33333333333333333333编号181182183184185186187188189190191192193194195196197198199200类别33333333333333333333编号201202203204205206207208209210211212213214215216217218219220类别33333333333333111111编号221222223224225226227228229230231232233234235236237238239240类别11111111131311113131编号241242243244245246247248249250251252253254255256257258259260类别11111111113311111111编号261262263264265266267268269270271272273274275276277278279280类别11111113131311111111编号281282283284285286287288289290291292293294295296297类别11111111111111111注：表2中共分为三类，其中1表示第一类，2表示第二类，3表示第三类。图(c)的数据分成了两类，分类结果如下表所示：表3问题三图(c)的分类结果编号1234567891011121314151617181920类别11111222221111122212注：表3中数据共分为2类，其中1表示第一类，2表示第二类。4.4问题四模型的建立与求解4.4.1问题四模型的建立规范化分割聚类算法许多机器学习和模式识别方法都依赖于距离度量。用于聚类的K-means方法等的性能都需要一个合适的距离度量。距离度量的选取必须要与数据特征的差异和不同尺度相对应，不同的距离函数可能使得距离结果差距很大。通常使用Minkowski度量来定义两个数据向量的距离为：(17)特别的，p=1意味着乌距离或者马氏距离，p=2意味着几距离或者欧式距离。理论上，距离函数必须满足一下条件：(18)距离度量还需要满足三角不等式：(19)的MinkoWSki不等式就是最大度量：以上距离度量是传递不变的，但是对于尺度变化不能保持不变。达到尺度不变的一般策略是对数据集进行预处理归一化。尺度不变的距离度量有：Canberra度量：(20)余弦相似性：(21)另一个修正欧式距离使其保持尺度不变的方式是Mahalanobis距离：(22)其中是整个数据集的协方差矩阵：(23)或者某个数据类别的协方差矩阵：(24)马氏距离度量通常应用于基于模型的聚类算法中(比如EM聚类算法)，高斯混合模型的最大似然估计使用马氏距离近似描述后验概率[6]。规范化分割聚类算法(NormalizedCut，Neut)是由Shi和Malik提出的一种图像分割方法，该方法的基本思想是把图像中的每个像素点看作图中的一个顶点，以图像象素点的颜色信息、空间信息和纹理信息来定义其相似度作为顶点间的权重，来构造一个带权重的无向图，在图上建立一个稀疏的关系矩阵，并利用该矩阵的谱信息对图或图像进行划分。具体步骤如下：步骤l：对含有n个像素的图像，构造一个带权重的无向图。V是图中顶点的集合，E是图中连接顶点i和j的边的集合；W是一个n阶对称矩阵，矩阵中的元素表示边上的权值，它表示顶点i和j之间的相关性。步骤2：解广义特征值方程：(25)这里只需要计算最小的几个特征值及其对应的特征向量；D是对角矩阵，对角线上的元素等于权重矩阵第i行上所有元素的和。步骤3：利用已经得到的谱的信息将图一分为二。步骤4：判断分割后的图像是否还需要再次分割，如需要分割，再对分割后的两部分图像分别重复上述过程进行分割。n阶权重矩阵W中元素，表示顶点i和j之间的相关性。权越大，顶点i和j的相关性越强，分在同一个子集中的可能性越大。任一顶点与自己的相关性最大。权的取值中应包含图像的重要特征信息，可以用高斯核定义权：；其中Q(i)是顶点i的空间坐标或其它特征信息。由于图像一般只是局部相关的，为计算方便，当顶点i和j间距离大于r(r<<n)时，取值为0。因此，最终得到的权重矩阵W是一个稀疏的对称矩阵[7]。规范化分割问题最终可以转化为谱问题：这里，X表示特征向量，该向量中隐含了图像分割的信息。第二小特征值对应的特征向量也称作Fielder向量，经常用来表示像素的类型隶属度信息。这种利用关系矩阵的谱信息进行图像分割的方法，通常也叫做谱分割。规范化分割方法的主要特点有：图像在大多数情况下仅仅是局部相关的，所以最后得到的矩阵牙都是稀疏矩阵。通常只需要几个最小的特征向量就可以对图像进行很好分割，不必求出矩阵的全部特征向量。分割对特征向量的精确度要求不高。大多数情况下，只需要其分量的正、负号就可以分割图像。通过挖掘流形数据集的局部的和全局的信息，构造高维流形信号曲面，提出了一种基于数据拓扑结构的拓扑距离度量，该距离度量更好的描述了流形的结构，使得同类别数据更加紧凑，不同类别的数据更加疏远。最后基于这种距离度量给出了一种新的流形聚类的算法框架。基于数据拓扑结构的拓扑距离(A)一维数据的情况对于两类一维流形数其基于数据拓扑结构的拓扑距离计算如下：数据集中的每个数据点依据其周围的邻居信息，产生不同大小的符合正态分布函数的信号：(26)其中，是由点的k近邻的平均距离决定的，平均距离越大，就越大，反之亦然。(2)将所有数据点产生的信号叠加起来，形成一条信号曲线。该信号曲线可以反映数据集合中的拓扑结构信息，纵坐标高的位置对应数据密集的区域，纵坐标低的位置对应数据稀疏的区域，通常稀疏区域也就是数据不同类别的分界区域。(3)、两两数据点的流形距离表示为两点间信号曲线的几何能量。(B)二维数据的情况对于二维数据的情况，以两个同心圆数据集合的聚类为例，来给出我们的拓扑距离算法：(1)数据集中的每个数据点根据其周围的邻域信息，产生了不同大小的符合正态分布函数的信号，信号函数为：(27)其中，是由点的k近邻的平均距离决定的，平均距离越大，就越大，反之亦然。系数是为了避免局部密度较大而导致信号曲面局部过高的现象。(2)将所有数据点产生的信号叠加起来，形成一条信号曲面。该信号曲面可以更好的反映数据间的流形特征，纵坐标高的位置对应数据密集的区域，纵坐标低的位置对应数据稀疏的区域，通常稀疏区域也就是数据不同类别的分界区域。(3)得到数据点间的信号曲线。(4)两两数据点的流形距离表示为两点间信号曲线的几何能量。两点间的流形距离定义为两点间对应曲线的弧长和平均曲率之和：(28)同类别数据点B和C间的距离基本没变，而不同类别数据点A和B的距离被表示为一条曲率较大的曲线段，其距离明显变大。(C)高维数据情况对于高维数据集合，可按照下面的步骤定义点与点之间的流形距离：(1)首先定义点与点之间对应的信号曲线是以从点到点的射线为横轴，以周围点对曲线横轴的信号影响为纵轴。比如，对于点到的射线上的点，其纵轴值为和周围的点产生的符合正态分布的信号在处的叠加值。(2)其次定义点与点的流形距离为曲线的的弧长和平均曲率之和为：(29)常用的欧式距离只涉及到两个点的空间位置关系，马氏距离在计算两点间距离时使用了数据集的协方差矩阵，考虑了数据集的整体分布。本章节提出的拓扑距离则更多的考虑了数据局部的拓扑结构，在计算两点间距离时加入了局部的拓扑结构信息。基于拓扑距离度量的流形聚类流形聚类方法主要是寻找数据间的全局或局部的相关性，寻找潜在的流形结构，再根据数据相对流形结构的相似性进行聚类。本文提出的拓扑距离度量很好的描述了流形的结构，使得同类别数据更加紧凑，不同类别的数据更加疏远。基于拓扑距离，使用来定义点A和B的相似度。使用拓扑距离定义的相似度矩阵和欧式距离定义的相似度矩阵。欧式距离定义的右图不同类别数据间的相似度很大，而拓扑距离定义的相似度使得不同类别数据间的相似度变得很小，从中可以看出本章节的距离度量更好的描述了流形的结构，使得同类别数据的距离紧凑，不同类别的数据的距离疏远[8]。本章节基于前一小节提出的拓扑距离，给出一个基于拓扑结构的流形聚类算法框架：步骤1：数据集中的每个数据点产生符合正态分布函数的信号，将所有数据点产生的信号叠加起来，形成一条信号曲面。步骤2：借助流形曲线来计算两两数据点间的流形距离：(30)步骤3：基于拓扑距离构造数据集合的相似度矩阵，再使用规范化分割聚类算法(Neut)进行流形聚类。4.4.2运行结果图中显示了圆台的点云，图中所存在的点构成曲面与平面，平面与曲面之间的点距离近，按问题三中低秩稀疏子空间聚类的方法来分类，分类的结果聚类错误率较高，分类不明显且效果不好，采取基于拓扑距离度量的流形聚类这一方法，能很好的描述了流形的结构，使得同类别数据更加紧凑，不同类别的数据更加疏远，图5为分类结果，点按照其所在的面分开，即圆台按照圆台的顶、底、侧面分成三类，红色部分代表一类，蓝色部分代表一类，绿色部分代表一类。图5(a)的分类结果图对于图(b)，问题给出是机器工件外部边缘轮廓的图像，需根据轮廓线中不同的直线和圆弧分类。根据上节算法，对所给图像分别分成两类、三类、四类，结果如下：图6(b)的分为两类结果图图6(b)的分为三类结果图图6(b)的分为四类结果图注：不同颜色代表不同的类。4.4.3结果分析对于图(a)，我们运用了不同的算法进行处理，通过计算不同算法的聚类错误率，我们得到如下表格：表4不同算法的聚类错误率比较算法LRRSSCLRSSC本章算法聚类错误率34.3%27.8%15.7%4.2%从上面的表格可以看出，我们所选用的算法的聚类错误率是最低的，因此对这组数据的分类是有效的。对于图(b)，我们选用类内类间距离来评判把数据分为几类是最优的。在点集中，是X的中心，为m个类的数据点集，，为类的中心，，其中为中心点构成，根据类间偏差平方和和类内偏差平方和之间的差的绝对值越小，分类效果越好。类间偏差平方和为，类内偏差平方和为。分别计算将题中图像分为两类、三类、四类的类间偏差平方和和类内偏差平方和之间的差的绝对值，通过计算，发现分成两类的的值最小，所以，分成两类的效果最优。五、模型的评价本文主要对数据进行聚类，通过很多聚类的算法对各种数据类型进行聚类。5.1模型的优点拉普拉斯特征映射降维的算法可以将高维数据进行降低维数，降低了算法的复杂度。谱聚类的算法不用对数据的全局结构作假设，并且具有识别非凸分布聚类的能力，非常适合于许多实际问题。流形距离能更好地反映这种空间的复杂分布特性。基于拓扑距离度量的流形聚类考虑了数据集的整体分布。分类之后的类间偏差平方和和类内偏差平方和之间的差的绝对值能更好地评价分类结果。5.2模型的缺点对许多真实世界中的数据集，由于采样密度稀疏，会使得所寻找的局部邻域结构产生偏差。(2)多流形上存在噪声时，数据点并不是精确地位于某一个流形上，将会影响数据的划分。六、参考文献[1]韩晓红，胡彧.K-means聚类算法的研究[J].太原理工大学学报，2009，40(3).[2]G.Liu,Z.Lin,S.Yan,J.Sun,Y.Yu,andY.Ma.Robustrecoveryofsubspacestructuresbylow-rankrepresentation.IEEETransactionsonPatternAnalysisandMachineIntelligence,35(1):171–184,2013.[3]B.Cheng,G.Liu,J.Wang,Z.Huang,andS.Yan,Multi-tasklowrankaffinitypursuitforimagesegmentation,ICCV,2011.[4]刘建华.基于隐空间的低秩稀疏子空间聚类[J].西北师范大学学报，2015，51(3)[5]XiangyangLiu,HongtaoLu,Multi—ManifoldModelingforHeadPoseEstimation,inProc.IEEEInternationalConferenceonImageProcessing(ICIP),I:3277-3280,HongKong,Sep.2010.[6]郭启勇.流形聚类的算法研究及其应用[D].复旦大学，2009.[7]黄启宏.流形学习方法理论研究及图像中应用[D].电子科技大学，2007.[8]GuillaumeObozinski,BenTaskar,MichaelI.Jordan.Jointcovariateselectionandjointsubspaceselectionformultipleclassificationproblems[J].StatisticsandComputing.2010(2).七、附录问题1主要程序：load('C:\Users\jsw\Downloads\2015研究生试题\2015赛题\B\0915版本\数据\第1题数据\1.mat')>>data=data';>>[ure]=KMeans(data,2);问题2(a)图主要程序：load('C:\Users\jsw\Downloads\2015研究生试题\2015赛题\B\0915版本\数据\第2题数据\2a.mat')>>X=data;CMat=admmLasso_mat_func(X,true,100);C=CMat;CKSym=BuildAdjacency(thrC(C,1));grps=SpectralClustering(CKSym,2);data=data’;fori=1:340if(grps(i,1)==1)plot(data(i,1),data(i,2),'*');holdonelseplot(data(i,1),data(i,2),'r+');holdonendend问题2(b)图主要程序：CMat=admmLasso_mat_func(X,false,20);C=CMat;CKSym=BuildAdjacency(thrC(C,1));grps=SpectralClustering(CKSym,3);data=data’;fori=1:300if(grps(i,1)==1)plot3(data(i,1),data(i,2),data(i,3),'*');holdonelseif(grps(i,1)==2)plot3(data(i,1),data(i,2),data(i,3),'r+');holdonelseplot3(data(i,1),data(i,2),data(i,3),'g.');holdonendendendend问题2(c)图主要程序：>>load('C:\Users\jsw\Downloads\2015研究生试题\2015赛题\B\0915版本\数据\第2题数据\2c.mat')>>CMat=admmLasso_mat_func(X,false,20);C=CMat;CKSym=BuildAdjacency(thrC(C,1));grps=SpectralClustering(CKSym,3);data=data’;d=pdist(data);>>W=squareform(d);>>N=tril(W,0);>>s=sum(N);>>D=diag(s);>>L=D-N;>>[Q,A]=eigs(L,2);>>[ure]=KMeans(Q,2);fori=1:400if(grps(i,1)==1)plot(data(i,1),data(i,2),'*');holdonelseplot(data(i,1),data(i,2),'r+');holdonendend问题2(d)图主要程序：load('C:\Users\jsw\Downloads\2015研究生试题\2015赛题\B\0915版本\数据\第2题数据\2d.mat')>>CMat=admmLasso_mat_func(X,false,20);C=CMat;CKSym=BuildAdjacency(thrC(C,1));grps=SpectralClustering(CKSym,3);data=data’;d=pdist(data);>>W=squareform(d);>>N=tril(W,0);>>s=sum(N);>>D=diag(s);>>L=D-N;>>[Q,A]=eigs(L,2);>>[ure]=KMeans(Q,2);fori=1:400if(grps(i,1)==1)plot(data(i,1),data(i,2),'*');holdonelseplot(data(i,1),data(i,2),'r+');holdonendend问题3(a)图主要程序：load('C:\Users\jsw\Downloads\2015研究生试题\2015赛题\B\0915版本\数据\第3题数据\3a.mat')>>data=data';>>d=pdist(data);>>W=squareform(d);>>N=tril(W,0);>>s=sum(N);>>D=diag(s);>>L=D-N;>>[Q,A]=eigs(L,2);>>[ure]=KMeans(Q,2);问题3(b)图主要程序：>>load('C:\Users\jsw\Downloads\2015研究生试题\2015赛题\B\0915版本\数据\第3题数据\3b.mat')>>>>X=data;CMat=admmLasso_mat_func(X,true,100);C=CMat;CKSym=BuildAdjacency(thrC(C,1));grps=SpectralClustering(CKSym,2);问题3(c)图主要程序：>>load('C:\Users\jsw\Downloads\2015研究生试题\2015赛题\B\0915版本\数据\第3题数据\3c.mat')X=data;CMat=admmLasso_mat_func(X,false,20);C=CMat;CKSym=BuildAdjacency(thrC(C,1));grps=SpectralClustering(CKSym,2);问题4(a)图主要程序：>>load('C:\Users\jsw\Downloads\2015研究生试题\2015赛题\B\0915版本\数据\第4题数据\4a.mat')if(nargin<2)%defaultsubspacesarelinearaffine=false;endif(nargin<3)%defaultregularizarionparametersalpha=800;endif(nargin<4)%defaultcoefficienterrorthresholdtostopADMM%defaultlinearsystemerrorthresholdtostopADMMthr=2*10^-4;endif(nargin<5)%defaultmaximumnumberofiterationsofADMMmaxIter=200;endif(length(alpha)==1)alpha1=alpha(1);alpha2=alpha(1);elseif(length(alpha)==2)alpha1=alpha(1);alpha2=alpha(2);endif(length(thr)==1)thr1=thr(1);thr2=thr(1);elseif(length(thr)==2)thr1=thr(1);thr2=thr(2);endN=size(Y,2);%settingpenaltyparametersfortheADMMmu1=alpha1*1/computeLambda_mat(Y);mu2=alpha2*1;if(~affine)%initializationA=inv(mu1*(Y'*Y)+mu2*eye(N));C1=zeros(N,N);Lambda2=zeros(N,N);err1=10*thr1;err2=10*thr2;i=1;%ADMMiterationswhile(err1(i)>thr1&&i<maxIter)%updatingZZ=A*(mu1*(Y'*Y)+mu2*(C1-Lambda2/mu2));Z=Z-diag(diag(Z));%updatingCC2=max(0,(abs(Z+Lambda2/mu2)-1/mu2*ones(N))).*sign(Z+Lambda2/mu2);C2=C2-diag(diag(C2));%updatingLagrangemultipliersLambda2=Lambda2+mu2*(Z-C2);%computingerrorserr1(i+1)=errorCoef(Z,C2);err2(i+1)=errorLinSys(Y,Z);%C1=C2;i=i+1;endfprintf('err1:%2.4f,err2:%2.4f,iter:%3.0f\n',err1(end),err2(end),i);else%initializationA=inv(mu1*(Y'*Y)+mu2*eye(N)+mu2*ones(N,N));C1=zeros(N,N);Lambda2=zeros(N,N);lambda3=zeros(1,N);err1=10*thr1;err2=10*thr2;err3=10*thr1;i=1;%ADMMiterationswhile((err1(i)>thr1||err3(i)>thr1)&&i<maxIter)%updatingZZ=A*(mu1*(Y'*Y)+mu2*(C1-Lambda2/mu2)+mu2*ones(N,1)*(ones(1,N)-lambda3/mu2));Z=Z-diag(diag(Z));%updatingCC2=max(0,(abs(Z+Lambda2/mu2)-1/mu2*ones(N))).*sign(Z+Lambda2/mu2);C2=C2-diag(diag(C2));%updatingLagrangemultipliersLambda2=Lambda2+mu2*(Z-C2);lambda3=lambda3+mu2*(ones(1,N)*Z-ones(1,N));%computingerrorserr1(i+1)=errorCoef(Z,C2);err2(i+1)=errorLinSys(Y,Z);err3(i+1)=errorCoef(ones(1,N)*Z,ones(1,N));%C1=C2;i=i+1;endendwarningoff;N=size(CKSym,1);MAXiter=1000;%MaximumnumberofiterationsforKMeansREPlic=20;%NumberofreplicationsforKMeans%NormalizedspectralclusteringaccordingtoNg&Jordan&Weiss%usingNormalizedSymmetricLaplacianL=I-D^{-1/2}WD^{-1/2}DN=diag(1./sqrt(sum(CKSym)+eps));LapN=speye(N)-DN*CKSym*DN;[uN,sN,vN]=svd(LapN);kerN=vN(:,N-n+1:N);fori=1:NkerNS(i,:)=kerN(i,:)./norm(kerN(i,:)+eps);endgroups=kmeans(kerNS,n,'maxiter',MAXiter,'replicates',REPlic,'EmptyAction','singleton');问题4(b)图主要程序：>>load('C:\Users\jsw\Downloads\2015研究生试题\2015赛题\B\0915版本\数据\第4题数据\4b.mat')if(nargin<2)%defaultsubspacesarelinearaffine=false;endif(nargin<3)%defaultregularizarionparametersalpha=800;endif(nargin<4)%defaultcoefficienterrorthresholdtostopADMM%defaultlinearsystemerrorthresholdtostopADMMthr=2*10^-4;endif(nargi