流形学习算法中邻域大小参数的递增式选取

邵 超，万春红，赵静玉

(河南财经政法大学计算机与信息工程学院，郑州 450002)

1 概述

近年来，人们陆续发现实际中的很多高维数据都分布在某个低维非线性流形上［1-2］，为此提出了多种能有效学习这种结构的流形学习算法如等距映射［3-4］(Isometric Mapping，ISOMAP)和局部线性嵌入［5-6］(Locally Linear Embedding，LLE)。现有的大多数流形学习算法能否成功应用还依赖于其邻域大小参数的选取是否合适，然而，目前该参数在实际中通常还难以高效选取，另外，数据中的噪音对邻域大小参数的合适性也会产生一定的影响［7］，从而限制了它们的实际应用［7-10］。

针对流形学习算法敏感于邻域大小参数的这一问题，目前大多数方法都基于残差来选取合适的邻域大小参数［3，11-12］，但这需要多次运行整个流形学习算法以计算相应的残差(尽管文献［11 -12］都只对其中具有极小重建误差或极小成本函数的少数几个邻域大小参数分别计算相应的残差);另外，残差只是维数约简过程中产生的误差，难以正确指导邻域大小参数的合适选取［13］。

通过对采用不同邻域大小参数得到的多个运行结果进行集成可以在一定程度上减轻流形学习算法对邻域大小参数的敏感程度［9］，但算法的时间复杂度会显著增加。通过识别和删除邻域图中可能出现的“短路”边也能在一定程度上避免邻域大小参数难以高效选取的问题，例如，文献［14］删除每一个数据点中具有最小重建权重的几个近邻点，文献［15］删除每一个数据点中具有最大图距离的几个近邻点，但它们又引入了一个同样难以选取的额外参数。文献［16 -17］通过组合若干个最小生成树来创建邻域图，没有了邻域大小参数，但最小生成树的数量同样难以有效确定。文献［18 -19］可以自适应地为每一个数据点选取不同的邻域大小参数，但计算比较复杂。

根据流形的局部欧氏性，邻域图上的所有邻域都呈线性或近似线性，是邻域大小参数被认为合适的基础，此时所有邻域的线性度量可聚成一类;而邻域大小参数一旦变得不合适，邻域图上的某些邻域就不再呈线性或近似线性［20］，其线性度量也不再能聚成一类。本文用加权主成分分析［21］(Principal Component Analysis，PCA)得到的重建误差对邻域图上所有邻域的线性程度进行度量，然后用贝叶斯信息准则［2，22］(Bayesian Information Criterion，BIC)对其聚类个数进行探测，从而能递增式地选取合适的邻域大小参数。

2 邻域大小的递增式选取

具有合适邻域大小参数的邻域图应能正确表达数据的邻域结构，按照流形的局部欧氏性，所有邻域都应呈线性或近似线性;而邻域大小参数一旦变得不合适，邻域图中开始出现“短路”边，从而使某些邻域不再呈线性或近似线性。线性度量有很多，考虑到鲁棒性，本文采用的是加权PCA［21］重建误差。当邻域大小k 合适时，邻域图上所有邻域的加权PCA重建误差都相对较小，可聚成一类，如图1 所示;而邻域大小k 一旦变得不合适，某些不再线性的邻域其加权PCA 重建误差就会变得很大，这样就使邻域图上所有邻域的加权PCA 重建误差不再能聚成一类，如图2 所示。因此，根据邻域图上所有邻域的加权PCA 重建误差所聚成的类别个数即可递增式地选取合适的邻域大小k。本文采用BIG［2，22］来探测其聚类个数。

图1 当邻域大小k 合适时的加权PCA 重建误差

图2 当邻域大小k 不合适时的加权PCA 重建误差

2.1 加权PCA 重建误差的计算

PCA 重建误差可用来度量数据的线性程度，一般而言，线性程度越大，PCA 重建误差就越小。然而，PCA 对数据中的异常点(outlier)比较敏感［21］，异常点的加入会极大地改变主成分，从而使PCA 重建误差难以准确度量数据的线性程度。如图3 所示，异常点(图中标注为“* ”)的加入使原本以实线标注的主成分变成了以虚线标注的主成分，从而降低了该数据集合的最大PCA 重建误差。

图3 异常点对PCA 的影响

邻域大小参数一旦变得不合适，邻域图上就会有某些邻域由线性变成非线性，和之前线性邻域中的那些数据点相比，新加入的导致邻域非线性的数据点就类似于异常点，其PCA 重建误差(通常也是该邻域中的最大PCA 重建误差)用来度量该邻域线性程度随邻域大小的变化趋势是比较合适的。然而，为准确计算这些“异常点”的PCA 重建误差，需要尽量降低异常点对PCA 的影响，因此，本文采用加权PCA 算法［21］，通过迭代给异常点赋以较小的权重，从而得到尽可能接近于原始主成分(图3 中标注为实线)的主成分(图3 中标注为点划线)。

加权主成分分析算法可简要描述如下［21］:

(1)采用标准PCA 算法获得数据集合X 的初始数据中心μ(0)和初始主成分p(0);

(2)t=0;

(3)DO

1)t=t+1;

2)按照式(1)计算每一个数据点Xi在数据中心μ(t-1)和主成分p(t-1)下的PCA 重建误差

3)按照式(2)计算每一个数据点Xi的权重wi:

4)重新计算数据中心μ(t)=

Until μ(t)和p(t)都相对稳定下来。

如上所述，本文用每一个邻域Nk(Xi)(即数据点Xi及其k 个最近邻数据点的集合，也就是加权PCA 算法中的数据集合X，因此有s=k +1)的最大加权PCA 重建误差(一个标量值，维数为1)来度量其线性程度，进而探测该邻域中是否出现了“短路”边。

2.2 贝叶斯信息准则的计算

由于具有不同聚类个数的聚类模型可以看作是不同的模型，因此BIC 作为模型选择准则就可以用来探测数据的聚类个数［2，22］。BIC 使用一个带修正项的对数似然统计量对不同的模型进行打分，BIC分值越大的模型就越适合于对应的数据集合。

模型Mj的BIC 分值BIC(Mj)如式(3)所示［2，22］:

本文采用K-均值聚类算法对邻域图上所有邻域的线性度量(即最大加权PCA 重建误差max(norm进行聚类，因此，数据集合D 就由邻域图上所有邻域的最大加权PCA 重建误差组成，即D=，则有m=n，pj=K(R +1)，其中，K 为聚类个数;R=1 为D 的维数。

不同的聚类个数K 就对应了不同的K-均值聚类模型，因此，分别用BIC(K=1)和BIC(K=2)表示聚类个数为1 和2 时对应K-均值聚类模型的BIC 分值。例如，图4 所示的2 个数据集合，其聚类个数分别为1和2，与之相适应的K-均值聚类模型的BIC 分值都相对较大(图4(a)中BIC(K=1)=-354.373 6，BIC(K=2)=-458.869 9，图4(b)中BIC(K=1)=-558.274 3，BIC(K=2)=-507.237 9)。

图4 BIC 用于探测聚类个数的数据集合

2.3 算法流程

该方法可简要描述如下:

(1)使用广度优先搜索算法获得能使邻域图连通的最小邻域大小参数kmin，作为递增式选取邻域大小参数的起点;

(2)k=kmin-1;

(3)DO

1)k=k+1;

2)对邻域图上的每一个邻域Nk(Xi)执行加权PCA 算法，获得邻域图上每一个邻域的最大加权PCA 重建误差，即式(4)中的数据集合D;

3)计算BIC(K=1)(由于只有一类，所以数据中心即为聚类中心);

4)计算BIC(K=2)，这需要对数据集合D 执行K-均值聚类算法将其聚成2 类;

Until BIC(K=1)＜BIC(K=2)

(4)kselected=k-1。

通过该方法的递增式选取，最终的kselected即为合适的邻域大小参数。和传统的基于残差的参数选取方法相比，该方法只需递增式地执行加权PCA 算法和K-均值聚类算法，并计算相应的BIC(K=1)和BIC(K=2)，而无需运行整个流形学习算法和计算相应的残差。

2.4 时间复杂度分析

该方法主要涉及以下3 个部分的计算:

(1)为确定kmin，需要从k=1 开始反复执行广度优先搜索算法(共执行kmin次)。广度优先搜索算法的时间复杂度为O(n +e)=O(kn)(e 为邻域图的边数)，因此，这一部分的时间复杂度为通常情况下，kmin都比较小。

(2)对于邻域大小k 的每一次迭代(共迭代(kselected-kmin+2)次)，在邻域图的每一个邻域(包括k+1 个d 维数据点，d 即为数据的维数)上分别执行加权主成分分析算法(共执行n 次)。加权主成分分析算法的时间复杂度为O(((k +1)d2+d3)c)(c为加权主成分分析算法的迭代次数，通常都比较小)，因此，这一部分的时间复杂度为O(((k +1)d2+d3)cn(kselected-kmin+2))。

(3)对于邻域大小k 的每一次迭代(共迭代(kselected-kmin+2)次)，对数据集合D(包括n 个1 维数据点)执行K-均值聚类算法将其聚成两类。K-均值聚类算法的时间复杂度为O(nKpl)=O(2nl)(n和p=1 分别为D 中数据点的个数及其维数，K=2为聚类个数，l 为K-均值聚类算法的迭代次数，通常都比较小)，因此，这一部分的时间复杂度为O(2nl(kselected-kmin+2))。

3 实验结果与分析

为验证该方法的有效性，本文采用具有2 000 个随机样本点的Swiss roll［3］和S-curve［5］数据集合(为了降低计算量，本文采用K-均值算法分别从中选取了n=500 个代表点，如图5 所示)，该方法的运行结果如图6 所示。

从图6 可以看出，对于Swiss roll 数据集合，该方法选取的邻域大小参数为k=8;对于S-curve 数据集合，该方法选取的邻域大小参数为k=16。它们的合适性可通过ISOMAP 算法的运行结果得以证实。

众所周知，数据中的噪音对邻域大小参数的合适性也会产生一定的影响［7］。为了验证该方法的鲁棒性，本文在以上2 个数据集合上分别加入一定的噪音［7］(如图7 所示)，该方法的运行结果如图8 所示。

图5 实验中的2 个数据集合

图6 在不同数据集合上的运行结果

图7 加入了噪音的2 个数据集合

图8 本文方法在不同数据集合上的运行结果

从图8 可以看出，对于加入了噪音的Swiss roll数据集合，该方法选取的邻域大小参数为k=6;对于加入了噪音的S-curve 数据集合，该方法选取的邻域大小参数为k=11。它们的合适性同样可通过ISOMAP 算法的运行结果得以证实。

此外，该方法基于的是流形的局部欧氏性，那么在由于采样不均匀而出现“空洞”的数据集合上是否依然有效?为了验证这一点，本文采用带有“空洞”的Swiss roll 和S-curve 数据集合(如图9 所示)，该方法的运行结果如图10 所示。

图9 带有“空洞”的2 个数据集合

从图10 可以看出，对于带有“空洞”的Swiss roll数据集合，该方法选取的邻域大小参数为k=8;对于带有“空洞”的S-curve 数据集合，该方法选取的邻域大小参数为k=11。它们的合适性同样可通过ISOMAP 算法的运行结果得以证实。

与传统的基于残差的参数选取方法(括号内指定的是邻域大小的选取范围)相比，本文方法不但能选取合适的邻域大小(分别如图6、图8 和图10 所示)，而且还具有较高的运行效率，如表1 所示(计算机配置为Intel i3-2120 CPU，主频为3.30 GHz，内存容量为4 GB)。

图10 本文方法在不同数据集合上的运行结果

表1 2 种方法在运行时间上的对比 s

4 结束语

按照流形的局部欧氏性，本文方法通过对邻域图上每一个邻域的线性程度进行度量，并对其聚类个数进行探测，能够高效地选取合适的邻域大小参数，且无需任何额外参数(在数据维数不是特别高的情况下)。下一步的研究方向是当数据的维数非常高时，如何高效地选取合适的邻域大小。

［1］Seung H S，Lee D D.The Manifold Ways of Perception［J］.Science，2000，290(5500):2268-2269.

［2］杨 剑，李伏欣，王 珏.一种改进的局部切空间排列算法［J］.软件学报，2005，16(9):1584-1590.

［3］Tenenbaum J B，de Silva V，Langford J C.A Global Geometric Framework for Nonlinear Dimensionality Reduction［J］.Science，2000，290(5500):2319-2323.

［4］王耀南，张 莹，李春生.基于核矩阵的Isomap 增量学习算法研究［J］.计算机研究与发展，2009，46(9):1515-1522.

［5］Roweis S T，Saul L K.Nonlinear Dimensionality Reduction by Locally Linear Embedding［J］.Science，2000，290(5500):2323-2326.

［6］Zhang S.Enhanced Supervised Locally Linear Embedding［J］.Pattern Recognition Letters，2009，30(13):1208-1218.

［7］Balasubramanian M，Shwartz E L，Tenenbaum J B，et al.The ISOMAP Algorithm and Topological Stability［J］.Science，2002，295(5552):7-17.

［8］Saul L K，Roweis S T.Think Globally，Fit Locally:Unsupervised Learning of Low Dimensional Manifolds［J］.Journal of Machine Learning Research，2003，4(1):119-155.

［9］詹德川，周志华.基于集成的流形学习可视化［J］.计算机研究与发展，2005，42(9):1533-1537.

［10］邵 超，黄厚宽，赵连伟.一种更具拓扑稳定性的ISOMAP 算法［J］.软件学报，2007，18(4):869-877.

［11］Kouropteva O，Okun O，Pietikainen M.Selection of the Optimal Parameter Value for the Locally Linear Embedding Algorithm［C］//Proc.of the 1st International Conference on Fuzzy Systems and Knowledge Discovery，Orchid Country Club.Singapore:IEEE Press，2002:359-363.

［12］Samko O，Marshall A D，Rosin P L.Selection of the Optimal Parameter Value for the Isomap Algorithm［J］.Pattern Recognition Letters，2006，27(1):968-979.

［13］黄启宏，刘 钊.流形学习中非线性维数约简方法概述［J］.计算机应用研究，2007，24(11):19-25.

［14］Saxena A，Gupta A，Mukerjee A.Non-linear Dimensionality Reduction by Locally Linear Isomaps［C］//Proc.of the 11th International Conference on Neural Information Processing.Calcutta，India:Springer，2004:1038-1043.

［15］Wen G，Jiang L，Shadbolt N R.Using Graph Algebra to Optimize Neighborhood for Isometric Mapping［C］//Proc.of the 20th International Joint Conference on Artificial Intelligence.Hyderabad，India:AAAI Press，2007，2398-2403.

［16］Carreira-Perpinan M A，Zemel R S.Proximity Graphs for Clustering and Manifold Learning［C］//Proc.of the 18th Annual Conference on Neural Information Processing Systems.Vancouver，Canada:MIT Press，2004:225-232.

［17］Yang L.Building k Edge-disjoint Spanning Trees of Minimum Total Length for Isometric Data Embedding［J］.IEEE Transactions on Pattern Analysis and Machine Intelligence，2005，27(10):1680-1683.

［18］Zhang Z，Wang J，Zha H.Adaptive Manifold Learning［J］.IEEE Transactions on Pattern Analysis and Machine Intelligence，2012，34(2):1473-1480.

［19］文贵华，江丽君，文 军.邻域参数动态变化的局部线性嵌入［J］.软件学报，2008，19(7):1666-1673.

［20］邵 超，张 斌，万春红.流形学习中邻域大小参数的合适性判定［J］.计算机工程与应用，2010，46(20):172-175.

［21］Chang H，Yeung D.Robust Locally Linear Embedding［J］.Pattern Recognition，2006，39(6):1053-1065.

［22］Pelleg D，Moore A.X-means:Extending K-means With Efficient Estimation of the Number of Clusters［C］//Proc.of the 17th International Conference on Machine Learning.San Francisco，USA:Morgan Kaufmann Publishers，2000:727-734.