科技的发展以及互联网的普及加快了人们生活的节奏,同时也产生了海量的数据和信息.我们正处在一个“数据丰富”而“知识贫乏”的时代,如何从浩瀚的数据中获取有价值的知识成为当务之急.大规模聚类是数据挖掘的主要工具之一,它能高效地组织大量数据,便于用户访问.聚类可以应用在各种场景中,例如网页搜索、图像检索、基因表达分析、推荐系统,以及基于交易数据的市场调研等^[1].

很多文献给出的大规模聚类技术是基于欧氏距离的,并且隐含假设所有数据点位于欧氏几何结构中.基于核的聚类方法突破了这一限制,它将数据点嵌入到一个高维非线性流形中,并且使用非线性的核距离函数度量这些点的相似性^[2].谱聚类就是一种典型的基于核的聚类方法,它把数据点当作无向加权图的节点,这样,聚类问题就变成了如何寻找最佳的图划分,使子图内部的连接权值最大,而子图之间的连接权值最小^[3].但是,图聚类的很多目标函数,如Ratio cut^[4],Normalized cut^[5],都无法在多项式时间内找到最优解,属于NP难问题.传统的求解方法是,将它们写成Rayleigh熵的形式,因为Rayleigh熵的最小值,第二小值,…,最大值分别对应Laplacian矩阵的最小特征值,第二小特征值,…,最大特征值,且极值在相应的特征向量处取得.于是.可以通过求解Laplacian矩阵的特征值和特征向量,将图划分的组合优化问题转化为数值优化问题,使其可以在多项式时间内解决.谱聚类的基本思想是:构造数据点的相似矩阵和对应的Laplacian矩阵,然后对Laplacian矩阵特征分解,再利用得到的特征向量进行聚类^[6].由于用到了矩阵的特征值和特征向量,所以这种算法称为谱聚类.

尽管谱聚类在规模较小的数据集上表现很好,可是如果数据点的个数n很大,谱聚类在计算n个点的成对相似性和存储大的相似矩阵时,会遇到O(n²)级资源占用的瓶颈.此外,求解Laplacian矩阵的前k个特征向量也需要相当大的时间和内存,这些问题都限制了谱聚类算法在大数据中的应用.

解决计算和内存难题最常用的方法是:使相似矩阵的一些元素归零,将矩阵稀疏化;从获得的稀疏相似矩阵,找到对应的Laplacian矩阵,然后调用稀疏的特征求解方法计算特征向量^[7].稀疏表示可以有效解决内存瓶颈,但是一些稀疏化方法仍然需要计算相似矩阵的全部元素.另一个加速谱聚类的重要方法是:利用Nyström近似技术进行特征分解,仅需相似矩阵的一部分,避免使用整个相似矩阵^[8].Kumar等人^[9]分析了Nyström方法的误差界,并指出,集成的Nyström方法比标准的Nyström方法具有更快的收敛速度.这种方法牺牲了准确的相似性值,但换来了更短的计算时间,提高了算法的效率.在不同尺度参数下,平移不变核会从低秩结构变化到块对角结构.基于此,Si等人^[10]提出MEKA算法来逼近平移不变核矩阵.该算法同时考虑了核矩阵的低秩和块结构,在速度、逼近误差和内存使用方面都有良好的表现.此外,Chen等人^[11]研究了并行的谱聚类方法,通过将n个数据点分配到p个分布的机器节点上,实现大规模数据的聚类.

因为对Laplacian矩阵特征分解的时间复杂度也很高,一般为O(n³),Dhillon等人^[12]提出了一种不使用特征向量的聚类方法.该方法源于核k-means与谱聚类的内在联系.核k-means算法是标准k-means算法的一个泛化.与标准k-means相比,核k-means的优势在于:通过将数据隐式映射到一个高维空间,可以发现输入空间中非线性分布的簇^[13].研究结果表明,加权形式的核k-means目标函数与Normalized Cut的目标函数都能写成矩阵迹的形式,可以认为两者在数学上是等价的^[14].这种等价性意味着:我们可以使用加权核k-means算法来优化Normalized Cut的目标函数;反过来,也可以将谱方法应用到加权核k-means中.在计算特征向量很困难的情况下,例如要对一个非常大的矩阵特征分解,加权核k-means算法可能比谱方法更合适.

虽然核距离函数有助于捕捉数据中的非线性结构,但是它要求计算并在内存中存储一个n×n的核矩阵,所以加权核k-means算法的空间复杂度也是O(n²),仍然不适合处理大数据问题.本文中,我们尝试解决这个由大的核矩阵构成的挑战.注意到,加权核k-means之所以需要全部核矩阵,是因为依据Representer定理,类中心是由所有数据点的线性组合表示的^[15].换句话说,类中心位于所有待聚类数据点的生成子空间中.我们可以通过把类中心限制在一个较小的子空间中,避免计算全部核矩阵.基于此,本文提出了求解大规模谱聚类的近似加权核k-means算法,通过随机选取一部分数据点,然后使用由这些数据点生成的子空间中的向量来逼近类中心.该近似仅需要计算和存储整个核矩阵的一部分.理论分析和实验对比均表明,近似加权核k-means与使用整个核矩阵的加权核k-means具有相似的聚类表现.

本文第1节介绍Normalized Cut的基本原理,并将其归结为矩阵迹的最大化问题.第2节分析加权核k-means算法,然后给出加权核k-means与Normalized Cut目标函数之间统一的数学关系.第3节提出近似加权核k-means算法,用于解决大数据的谱聚类问题.第4节从理论层面讨论近似加权核k-means算法的计算复杂度和误差范围.第5节在不同的数据集上验证所提出的算法有效性,并与其他算法的聚类表现进行对比.最后归纳本文的主要贡献,以及下一步要做的工作.

1 归一化图划分

一个无向加权图可以表示为G=(V,E,A),其中,V是所有节点的集合;E是连接节点的边的集合,每条边都有一个权值,衡量两点属于同一类的可能性的大小;这些权值构成了亲和矩阵A,通常A是非负的和对称的.

给定数据集X={x₁,x₂,…,x_n},X中包含n个数据点.把这些数据点看作图G的节点,令节点集合V=[n]表示所有待聚类的元素.要把n个点聚成k类,就是将V划分成k个不相交的子集,即,$V = \bigcup\nolimits_{i = 1}^k {{V_i}} {\rm{,}}$,且 $ {V_i} \cap {V_j} = \emptyset ,i \ne j. $ .

假设V₁,V₂⊂V,定义links(V₁,V₂)为V₁和V₂之间总的连接权值:

$ links({V_1},{V_2}) = \sum\limits_{i \in {V_1},j \in {V_2}} {{A_{ij}}} $

(1)

一个子集V₁的度(degree)就是V₁中的点与所有数据点的连接权值之和:

$ degree\left( {{V_1}} \right) = links\left( {{V_1},V} \right) $

(2)

使用度作为归一项,定义linkratio(V₁,V₂)表示V₁与V₂的连接在V₁与全集的连接中所占的比例:

$ linkratio({V_1},{V_2}) = \frac{{links({V_1},{V_2})}}{{degree({V_1})}} $

(3)

有两种特殊的linkratio:一种是linkratio(V₁,V₁),用来衡量V₁内部的连接;另一种是它的补集linkratio(V₁, V\V₁),用来衡量V₁外部的连接.一个好的聚类划分希望类内部的连接比较紧密,同时,类之间的连接比较松散.由这两个目标得到两种归一化的k-way划分准则:Normalized Association和Normalized Cut.它们的目标函数如下:

$ NAssoc({V_1},...,{V_k}) = \sum\limits_{i = 1}^k {linkratio({V_i},{V_i})} $

(4)

$ NCut({V_1},...,{V_k}) = \sum\limits_{i = 1}^k {linkratio({V_i},V\backslash {V_i})} $

(5)

根据links(V₂,V\V₁)=degree(V₁)-links(V₁,V₁),可以得到NAssoc(V₁,…,V_k)+NCut(V₁,…,V_k)=k.所以在最大化Normalized Association的同时,也能最小化Normalized Cut.于是,图的k-way划分可以表示成下面的优化问题:

$ maxNAssoc\left( {{V_1}, \ldots ,{V_k}} \right). $

为了便于计算,定义一个隶属矩阵U∈$\mathbb{R}$^k×n来描述划分结果{V₁,…,V_k}.设U=(u₁,…,u_k)^T,其中, $ u_{i}^{T}\in {{\mathbb{R}}^{1\times n}} $ 是U的第i行,表示第i类V_i中包含哪些元素:对于数据集中的第j个点x_j,若x_j∈V_i,则U_ij=1;若x_j∉V_i,则U_ij=0.由于每个节点只能归到V₁到V_k中的一个类,U的每列只包含一个1,故U^T1_k=1_n,其中,1_n∈$\mathbb{R}$^n×1表示元素全为1的向量.根据对称的亲和矩阵A,定义度矩阵D:

$ D=diag\left( A{{\mathbf{1}}_{n}} \right) $

(6)

其中,diag(·)表示由它的矢量参数构成的对角矩阵.然后,可以重写links和degree如下:

$links({{V}_{i}},{{V}_{i}})=u_{i}^{T}A{{u}_{i}} $

(7)

$ degree({{V}_{i}})=u_{i}^{T}D{{u}_{i}} $

(8)

定义矩阵Z∈R^k×n,令Z=(z₁,…,z_k)^T=(UDU^T)^-1/2U,其中, ${{z}_{i}}={{u}_{i}}{{(u_{i}^{T}D{{u}_{i}})}^{-1/2}} $ ,用tr(·)表示矩阵的迹,根据公式(7)、

公式(8),则有:

$ NAssoc({{V}_{1}},...,{{V}_{k}})=\sum\limits_{i=1}^{k}{\frac{u_{i}^{T}A{{u}_{i}}}{u_{i}^{T}D{{u}_{i}}}}=\sum\limits_{i=1}^{k}{\frac{u_{i}^{T}}{{{(u_{i}^{T}D{{u}_{i}})}^{1/2}}}A\frac{{{u}_{i}}}{{{(u_{i}^{T}D{{u}_{i}})}^{1/2}}}}=\sum\limits_{i=1}^{k}{z_{i}^{T}A{{z}_{i}}}=tr(ZA{{Z}^{T}}). $

注意到ZDZ^T=(UDU^T)^-1/2U·D·U^T(UDU^T)^-1/2=I_k,其中,I_k∈$\mathbb{R}$^k×k是单位矩阵.

令 $\tilde{Z}=Z{{D}^{1/2}}\text{,}$ 则 ${{\tilde{Z}}^{T}}={{D}^{1/2}}{{Z}^{T}},\tilde{Z}{{\tilde{Z}}^{T}}=Z{{D}^{1/2}}\cdot {{D}^{1/2}}{{Z}^{T}}=ZD{{Z}^{T}}={{I}_{k}}\text{,} $ 那么,

$tr(ZA{Z^T}) = tr(Z{D^{1/2}} \cdot {D^{ - 1/2}}A{D^{ - 1/2}} \cdot {D^{1/2}}{Z^T}) = tr(\tilde Z{D^{ - 1/2}}A{D^{ - 1/2}}{\tilde Z^T}). $

所以,最大化NAssoc(V₁,…,V_k)等价于矩阵迹的最大化问题,即

$ \mathop {\max }\limits_{\tilde Z} tr(\tilde Z{D^{ - 1/2}}A{D^{ - 1/2}}{\tilde Z^T}) $

(9)

注意到,该问题需要满足约束条件 $ \tilde Z{\tilde Z^T} = {I_k}. $ 一个可行的求解方法是,根据谱图理论,定义Laplacian矩阵L=D^-1/2AD^-1/2,通过对L特征分解,利用其前k个特征向量构造矩阵 $ \tilde Z{\rm{.}} $ 由于特征向量中包含了数据点的类属信息,因此可以在由特征向量组成的低维特征空间中对数据点进行划分.

2 加权核k-means算法

核k-means是经典k-means算法的一个非线性扩展.它把k-means算法中使用的欧氏距离函数d²(x_a,x_b)= ||x_a-x_b||²替换成了一个非线性核距离,定义为

$ d_\kappa ^2({x_a},{x_b}) = \kappa ({x_a},{x_a}) + \kappa ({x_b},{x_b}) - 2\kappa ({x_a},{x_b}), $

其中,x_a∈$\mathbb{R}$^d和x_b∈$\mathbb{R}$^d是两个数据点,κ(·,·):$\mathbb{R}$^d×$\mathbb{R}$^d→R表示核函数.核函数建立了从原(输入)空间到高维核空间的非线性映射,有助于识别输入空间中非线性分布的簇.

令X={x₁,x₂,…,x_n}表示包含n个点的输入数据集,k是设定的类数,K∈R^n×n是核矩阵,K_ij=κ(x_i,x_j).由核函数κ(·,·)诱导出的线性空间称为再生核希尔伯特空间(reproducing kernel Hilbert space,简称RKHS),用H_κ表示.将数

据点划分成k个不相交的类{V₁,…,V_k},每个类中数据点到类中心的距离的平方和称为聚类误差,核k-means的目标是寻找合理的划分,使总的聚类误差最小.

为了将核k-means算法与图的Normalized Cut联系起来,需要假设每个点x_i都有一个权值w_i,引入权值后,核k-means就变成了加权核k-means,其目标函数可以表示为

$ \min J({V_1},...,{V_k}) = \sum\limits_{i = 1}^k {\sum\limits_{{x_j} \in {V_i}} {{w_j}|\kappa ({x_j}, \cdot ) - {c_i}( \cdot )|_{{H_\kappa }}^2} } = \sum\limits_{i = 1}^k {\sum\limits_{j = 1}^n {{w_j}{U_{ij}}|\kappa ({x_j}, \cdot ) - {c_i}( \cdot )|_{{H_\kappa }}^2} } $

(10)

其中, $|...{|_{{H_\kappa }}} $ 表示H_κ的泛函范数,c_i(·)∈H_κ表示类中心,U∈{0,1}^k×n是数据点的隶属矩阵.

定义加权的隶属矩阵Y∈$\mathbb{R}$^k×n:Y=(y₁,…,y_k)^T=UW,其中,W∈$\mathbb{R}$^n×n是由数据点权值构成的对角矩阵,即,W= diag(w₁,…,w_n).令 ${s_i} = y_i^T{1_n} $ 表示Y的第i行之和,s_i的倒数组成另一个对角矩阵L∈R^k×k:L=[diag(s₁,…,s_k)]^-1.利用L将Y归一化,得到 $\hat{Y}\in {{\mathbb{R}}^{k\times n}} $ 和 $\tilde{Y}\in {{\mathbb{R}}^{k\times n}}\text{:}$

$\begin{matrix} \hat{Y}={{({{{\hat{y}}}_{1}},...,{{{\hat{y}}}_{k}})}^{T}}={{[diag({{s}_{1}},...,{{s}_{k}})]}^{-1}}UW=LUW=LY, \\ \tilde{Y}={{({{{\tilde{y}}}_{1}},...,{{{\tilde{y}}}_{k}})}^{T}}={{\left[ diag\left( \sqrt{{{s}_{1}}},...,\sqrt{{{s}_{k}}} \right) \right]}^{-1}}U{{W}^{1/2}}={{L}^{1/2}}U{{W}^{1/2}}. \\ \end{matrix} $

那么,加权核k-means最佳的类中心可以表示为

$ {{c}_{i}}(\cdot )=\frac{1}{{{s}_{i}}}\sum\limits_{{{x}_{j}}\in {{V}_{i}}}{{{w}_{j}}\kappa ({{x}_{j}},\cdot )}=\sum\limits_{j=1}^{n}{{{{\hat{Y}}}_{ij}}\kappa ({{x}_{j}},\cdot )},i\in [k] $

(11)

定理1. 设核矩阵K=(φ₁,…,φ_n)^T, $\varphi _{i}^{T}\in {{\mathbb{R}}^{1\times n}} $ 表示K的第i行.结合权值矩阵W,加权核k-means的目标函数也可以由矩阵的迹来表示:

$ J({{V}_{1}},...,{{V}_{k}})=tr({{W}^{1/2}}K{{W}^{1/2}})-tr(\tilde{Y}{{W}^{1/2}}K{{W}^{1/2}}{{\tilde{Y}}^{T}}). $

注意到,核矩阵K、权值矩阵W都是给定的,故tr(W^1/2KW^1/2)是一个常量,所以最小化加权核k-means的目标函数等价于:

$\underset{{\tilde{Y}}}{\mathop{\max }}\,tr(\tilde{Y}{{W}^{1/2}}K{{W}^{1/2}}{{\tilde{Y}}^{T}}) $

(12)

与 $NAssoc({{V}_{1}},...,{{V}_{k}})=tr(\tilde{Z}{{D}^{-1/2}}A{{D}^{-1/2}}{{\tilde{Z}}^{T}}) $ 对比可以发现:若令加权核k-means的权值矩阵W=D,核矩阵K=D^-1AD^-1,则

$ \tilde{Y}={{L}^{1/2}}U{{W}^{1/2}}={{(UW{{U}^{T}})}^{-1/2}}U{{W}^{1/2}}={{(UD{{U}^{T}})}^{-1/2}}U{{D}^{1/2}}=\tilde{Z}. $

于是, $tr(\tilde{Y}{{W}^{1/2}}K{{W}^{1/2}}{{\tilde{Y}}^{T}})=tr(\tilde{Z}{{D}^{1/2}}\cdot {{D}^{-1}}A{{D}^{-1}}\cdot {{D}^{1/2}}{{\tilde{Z}}^{T}})=tr(\tilde{Z}{{D}^{-1/2}}A{{D}^{-1/2}}{{\tilde{Z}}^{T}}). $

可见,Normalized Cut问题中 $tr(\tilde{Z}{{D}^{-1/2}}A{{D}^{-1/2}}{{\tilde{Z}}^{T}}) $ 的最大化与加权核k-means问题中 $tr(\tilde{Y}{{W}^{1/2}}K{{W}^{1/2}}{{\tilde{Y}}^{T}}) $ 的最大化是等价的,因此可以使用加权核k-means算法来求解Normalized Cut的目标函数,进而得到最佳的划分子图.为了保证加权核k-means算法最后收敛,一般要求核矩阵K是正定的.加权核k-means算法通过多次迭代,不断更新聚类中心和隶属矩阵,使聚类误差逐渐减小.与传统的谱聚类方法相比,利用迭代算法解决图划分问题,避免了对拉普拉斯矩阵特征分解,降低了计算复杂度,同时便于使用局部搜索等技术来提高聚类质量.

由定理1可知:加权核k-means的实现需要计算和存储全部n×n的核矩阵K,使它不适合处理大规模的数据集.Zhang和Rudnicky^[16]通过把核矩阵划分成小块,每次只用其中的一块来降低空间需求.尽管该技术考虑了空间复杂度,它仍然需要计算整个核矩阵.为了同时降低加权核k-means的时间和空间复杂度,本文考虑利用数据集中的部分抽样点计算近似的核矩阵,并设计了近似加权核k-means算法,用于处理大数据的谱聚类问题.

3 近似加权核k-means算法

降低加权核k-means的复杂度的一个简单而直观的方法是:从待聚类的数据集中随机选取m个点,在这些采样点中,使用标准加权k-means算法确定k个类中心;然后对于剩余的每个点,分别计算k个类中心与该点的距离,将该点与最近的类中心归为一类.采用这种方法尽管可以减少运行时间和空间需求,但是由于采样的任意性和局限性,它的表现无法和使用全部核矩阵的原算法相比,除非抽样点的数量足够大.

仔细观察公式(10)可以发现,由于类中心{c_i(·),i∈[k]}是全部待聚类数据点的线性组合,所以加权核k-means需要计算全部核矩阵K.换句话说,类中心位于所有数据点生成的子空间中,即,c_i(·)∈H_κ=span(κ(x₁,·),…,κ(x_n,·)), i∈[k].倘若把类中心的解限制在一个较小的子空间 ${{\hat{H}}_{\kappa }}\subset {{H}_{\kappa }} $ 里,就可以避免计算整个核矩阵.构造的 ${{\hat{H}}_{\kappa }}$ 应该具有如下性质:(1) ${{\hat{H}}_{\kappa }}$应该足够小,以便高效计算;(2) ${{\hat{H}}_{\kappa }}$的覆盖面应该足够大,以产生类似于使用H_κ的聚类结果.基于这个简单而重要的发现,本文提出了一个构造${{\hat{H}}_{\kappa }}$的有效方法,可以显著降低加权核k-means的复杂度.首先,随机选取m个数据点(m<<n),表示为 $\hat{X}=\{{{\hat{x}}_{1}},...,{{\hat{x}}_{m}}\}$ ,用来构造子空间 ${{\hat{H}}_{\kappa }}=span(\kappa ({{\hat{x}}_{1}},\cdot ),...,\kappa ({{\hat{x}}_{m}},\cdot ))$ .给定子空间${{\hat{H}}_{\kappa }}$,将加权核k-means的目标函数写成:

$ \min J({{V}_{1}},...,{{V}_{k}})=\sum\limits_{i=1}^{k}{\sum\limits_{j=1}^{n}{{{w}_{j}}{{U}_{ij}}|\kappa ({{x}_{j}},\cdot )-{{c}_{i}}(\cdot )|_{{{H}_{\kappa }}}^{2}}},{{c}_{i}}(\cdot )\in {{\hat{H}}_{\kappa }} $

(13)

定义两种核矩阵 $\hat{K}\in {{\mathbb{R}}^{m\times m}}$ 和 $\tilde{K}\in {{\mathbb{R}}^{n\times m}},\hat{K} $ 由$\hat{X}$中样本点之间的核相似性构成, $\tilde{K}$ 由X中的数据点与$\hat{X}$中的样本点之间的核相似性构成.下面的引理1给出了限制在子空间${{\hat{H}}_{\kappa }}$中的类中心c_i(·)的最优解:

引理1. 给定隶属矩阵U和权值矩阵W,根据采样点得到的类中心设为 ${{c}_{i}}(\cdot )=\sum\limits_{j=1}^{m}{{{\alpha }_{ij}}\kappa ({{{\hat{x}}}_{j}},\cdot )}\text{,}$ 其中,α∈$\mathbb{R}$^k×m是未知参数.为使目标函数J(V₁,…,V_k)最小,最佳的α应满足 $\alpha =\hat{Y}\tilde{K}{{\hat{K}}^{-1}}\text{,}$ 其中,矩阵$\hat{Y}=LUW.$

定理2. 若加权核k-means的类中心${{c}_{i}}(\cdot )\in {{\hat{H}}_{\kappa }}\text{,}$i∈[k],其目标函数可以归结为矩阵迹的最优化问题:

$J({{V}_{1}},...,{{V}_{k}})=tr({{W}^{1/2}}K{{W}^{1/2}})-tr(\tilde{Y}{{W}^{1/2}}\tilde{K}{{\hat{K}}^{-1}}{{\tilde{K}}^{T}}{{W}^{1/2}}{{\tilde{Y}}^{T}}).$

定理2表明,数据点的隶属矩阵U可以通过矩阵$\tilde{K}$和$\hat{K}$来确定.由于$\hat{K}$是$\tilde{K}$的一部分,事实上只需要计算$\tilde{K}$.当m<<n时,计算$\tilde{K}$的开销远小于计算整个核矩阵K.观察定理1和定理2,定理2中的方法也可以视为用$\tilde{K}{{\hat{K}}^{-1}}{{\tilde{K}}^{T}}$来逼近定理1中的核矩阵K,这与矩阵低秩逼近的Nyström方法是类似的^[17].由定理2得到的算法称为近似加权核k-means,利用该算法来求解大规模谱聚类问题,可以显著降低算法的时间和空间复杂度,提高聚类的效率.详细过程见算法1.

算法1. 近似加权核k-means算法.

输入:

• X={x₁,…,x_n}:待聚类的n个数据点的集合;

• m:随机采样的数据点个数(m<<n);

• k:聚类个数;

• MAXITER:最大迭代次数.

输出:包含最终聚类信息的隶属矩阵U.

Step 1. 从X中随机选取m个点,表示为$\hat{X}=({{\hat{x}}_{1}},...,{{\hat{x}}_{m}})$.

Step 2. 根据核函数K=D^-1AD^-1,计算$\tilde{K}={{[\kappa ({{x}_{i}},{{\hat{x}}_{j}})]}_{n\times m}}$和$\hat{K}={{[\kappa ({{\hat{x}}_{i}},{{\hat{x}}_{j}})]}_{m\times m}}.$

Step 3. 计算$T=\tilde{K}{{\hat{K}}^{-1}}.$

Step 4. 随机初始化隶属矩阵U.

Step 5. 令迭代次数t=0.

Step 6. repeat

Step 7. 令t=t+1.

Step 8. 计算加权隶属矩阵Y=UW=UD和对角矩阵L=[diag(Y1_n)]^-1,并将Y归一化:$\hat{Y}=LY.$

Step 9. 计算$\alpha =\hat{Y}T=\hat{Y}\tilde{K}{{\hat{K}}^{-1}}.$

Step 10. for j=1,…,n do

Step 11. 为每个点x_j找到最近的类中心i^*:

${{i}^{*}}=\underset{i\in [k]}{\mathop{\arg \min }}\,|\kappa ({{x}_{j}},\cdot )-{{c}_{i}}(\cdot )|_{{{H}_{\kappa }}}^{2}=\underset{i\in [k]}{\mathop{\arg \min }}\,(\alpha _{i}^{T}\hat{K}{{\alpha }_{i}}-2\tilde{\varphi }_{j}^{T}{{\alpha }_{i}}),$

其中,$\alpha _{i}^{T}$是α的第i行,$\tilde{\varphi }_{j}^{T}$ 是$\tilde{K}$的第j行.

Step 12. 更新U的第j列,令其第i=i行元素U_ij=1,而其余元素为0.

Step 13. end for

Step 14. until 隶属矩阵U不再变化或t>MAXITER.

4 算法分析

本节首先分析所提出的近似加权核k-means算法的计算复杂度和收敛性,然后与标准加权核k-means算法进行对比,从理论上研究了抽样近似方法对算法1聚类误差的影响.

4.1 计算复杂度

近似加权核k-means算法最费时的操作是矩阵求逆${{\hat{K}}^{-1}}$和计算$T=\tilde{K}{{\hat{K}}^{-1}}\text{,}$该过程的时间复杂度是O(m³+m²n).计算α和更新隶属矩阵U的计算开销是O(mnkt),其中,t是算法收敛需要迭代的次数.所以,算法1总的计算复杂度为O(m³+m²n+mnkt).为了避免矩阵求逆${{\hat{K}}^{-1}}\text{,}$进一步提高算法的运行效率,可以把计算$\alpha =\hat{Y}\tilde{K}{{\hat{K}}^{-1}}$转化成下面的优化问题:

$\begin{align} & {{\alpha }^{*}}=\underset{\alpha \in {{\mathbb{R}}^{k\times m}}}{\mathop{\arg \min }}\,\left[ \sum\limits_{i=1}^{k}{({{s}_{i}}\alpha _{i}^{T}\hat{K}{{\alpha }_{i}}-2y_{i}^{T}\tilde{K}{{\alpha }_{i}})} \right]=\underset{\alpha \in {{\mathbb{R}}^{k\times m}}}{\mathop{\arg \min }}\,[tr({{L}^{-1}}\alpha \hat{K}{{\alpha }^{T}})-2tr(Y\tilde{K}{{\alpha }^{T}})] \\ & =\underset{\alpha \in {{\mathbb{R}}^{k\times m}}}{\mathop{\arg \min }}\,\left[ \frac{1}{2}tr(\alpha \hat{K}{{\alpha }^{T}})-tr(\hat{Y}\tilde{K}{{\alpha }^{T}}) \right]. \\ \end{align}$

如果$\hat{K}$是正定的,即,其最小的特征值显著大于0,该优化问题可以用梯度下降法来求解,收敛速度是O(log(1/ε)),其中,zε是期望的精度.由于梯度下降法的每一步需要花费O(m²k)的代价,所以整个求解过程的计算开销是O(m²klog(1/ε)).当klog(1/ε)<<m时,O(m²klog(1/ε))<<O(m³).使用这一技巧,可以将算法1的时间复杂度降低到O(m²klog(1/ε)+m²n+mnkt).因为算法1中需要存储的最大矩阵是$\tilde{K}$,所以其空间复杂度为O(mn).与标准加权核k-means算法O(n²)的复杂度相比,该方法大幅度降低了大数据聚类过程中的时间和空间需求.

4.2 算法收敛性

文献[14]中指出:如果核矩阵是正定的,就可以保证加权核k-means算法收敛.利用加权核k-means优化Normalized Cut的目标函数时,需要定义核矩阵K=D^-1AD^-1,权值矩阵W=D.但是,如果A是一个任意的邻接矩阵, D^-1AD^-1不一定是正定的,所以加权核k-means也不一定会收敛.本节通过为核矩阵K增加恰当的对角偏移量来避免这个问题.该解决方法可以看作对Roth等人^[18]工作的推广,他们将对角偏移方法用在了非加权的情况中.

给定亲和矩阵A,定义K'=σD^-1+D^-1AD^-1,其中,σ是一个正的足够大的常数,以保证K'是正定的.

因为D^-1是一个正的对角矩阵,加上σD^-1就为D^-1AD^-1的对角线元素增加了正的偏移量.用K'代替公式(12)中的K,可以得到:

$tr(\tilde{Y}{{W}^{1/2}}{K}'{{W}^{1/2}}{{\tilde{Y}}^{T}})=tr(\tilde{Z}{{D}^{1/2}}\sigma {{D}^{-1}}{{D}^{1/2}}{{\tilde{Z}}^{T}})+tr(\tilde{Z}{{D}^{-1/2}}A{{D}^{-1/2}}{{\tilde{Z}}^{T}})=\sigma k+tr(\tilde{Z}{{D}^{-1/2}}A{{D}^{-1/2}}{{\tilde{Z}}^{T}}).$

因此,使用K'最大化$\tilde{Y}$与公式(9)中的Normalized Association问题是等价的,只是K'已经被构造成一个正定矩阵.运行基于K'的近似加权核k-means,可以单调优化Normalized Association的目标函数,保证算法的收敛性.

4.3 误差范围

与使用全部核矩阵的加权核k-means相比,近似加权核k-means的不同之处仅在于:其类中心限制在一个较小的子空间${{\hat{H}}_{\kappa }}$中,而${{\hat{H}}_{\kappa }}$是基于采样点构成的.下面就从该约束出发,探究使用这种近似的方法后,算法1的期望误差的范围.

命题1. 给定隶属矩阵U和权值矩阵W,设加权隶属矩阵Y=(y₁,…,y_k)^T=UW.定义二值随机变量$\xi $=($\xi $₁,$\xi $₂,…, $\xi $_n)^T∈{0,1}^n×1来描述随机采样过程:如果数据点x_i被选中构造子空间,则×_i=1;否则,×_i=0.于是,近似加权核k-means的聚类误差可以表示为

$L(Y,\xi )=tr({{W}^{1/2}}K{{W}^{1/2}})+\sum\limits_{i=1}^{k}{{{L}_{i}}(Y,\xi )}$

(14)

其中,${{L}_{i}}(Y,\xi )=\underset{{{\alpha }_{i}}\in {{\mathbb{R}}^{n\times 1}}}{\mathop{\min }}\,[{{s}_{i}}{{({{\alpha }_{i}}\circ \xi )}^{T}}K({{\alpha }_{i}}\circ \xi )-2y_{i}^{T}K({{\alpha }_{i}}\circ \xi )].$

注意到,$\xi ={{\mathbf{1}}_{n}}$=1_n时(1_n是元素全为1的向量),说明选择了所有数据点来构造子空间${{\hat{H}}_{\kappa }}$,等价于使用全部核矩阵的加权核k-means算法,因此,L(Y,1_n)是标准加权核k-means算法的聚类误差.下面的引理2给出了L(Y,×)期望的范围.

引理2. 已知加权隶属矩阵Y,L(Y,×)的期望满足下面的不等式:

${{E}_{\xi }}[L(Y,\xi )]\le L(Y,{{1}_{n}})+tr\left( \tilde{Y}{{\left[ {{({{W}^{1/2}}K{{W}^{1/2}})}^{-1}}+\frac{m}{n}{{[diag({{W}^{1/2}}K{{W}^{1/2}})]}^{-1}} \right]}^{-1}}{{{\tilde{Y}}}^{T}} \right),$

其中,$L(Y,{{1}_{n}})=tr({{W}^{1/2}}K{{W}^{1/2}})-tr(\tilde{Y}{{W}^{1/2}}K{{W}^{1/2}}{{\tilde{Y}}^{T}}).$

推论1. 假设对于任意x都有κ(x,x)≤1,权值w(x)≤1,令$\lambda $₁≥$\lambda $₂≥…≥$\lambda $_n≥0为加权核矩阵W^1/2KW^1/2的特征值.给定加权隶属矩阵Y,有下面的不等式成立:

$\frac{{{E}_{\xi }}[L(Y,\xi )]-L(Y,{{1}_{n}})}{L(Y,{{1}_{n}})}\le \frac{k/m}{\sum\limits_{i=k+1}^{n}{{{\lambda }_{i}}/n}}.$

为了说明推论1的结果,考虑一个特殊的加权核矩阵W^1/2KW^1/2,它的前a个特征值等于n/a,剩余的特征值为0,即$\lambda $₁=…=$\lambda $_a=n/a,$\lambda $_a+1=…=$\lambda $_n=0.进一步假设a≥2k,即,W^1/2KW^1/2的非0特征值的个数大于类数的2倍,那么根据推论1,有下式成立:

$\frac{{{E}_{\xi }}[L(Y,\xi )]-L(Y,{{1}_{n}})}{L(Y,{{1}_{n}})}\le \frac{k/m}{\sum\limits_{i=k+1}^{n}{{{\lambda }_{i}}/n}}=\frac{k/m}{{{a}^{-1}}(a-k)}\le \frac{k/m}{{{a}^{-1}}(a-a/2)}=\frac{2k}{m}.$

上式表明:当W^1/2KW^1/2的非0特征值的个数显著的大于类数时,使用逼近策略的近似加权核k-means的聚类误差与采样个数m是息息相关的,其与标准加权核k-means聚类误差的差别会随着样本点的增加以O(1/m)的速率减小.

5 实验与分析

将基于近似加权核k-means的谱聚算法记作AWKK-SC,为了测试AWKK-SC算法的有效性,本节对比了该算法与另外4种典型的聚类算法在基准数据集上的聚类性能.对照算法分别是:标准谱聚类算法(NJW-SC)^[19]、基于标准加权核k-means的谱聚类算法(WKK-SC)^[12]、基于Nyström低秩近似的谱聚类算法(Nyström-SC)^[8]、基于MEKA低秩近似的核聚类算法(MEKA-KC)^[10].所有实验都是在一台高性能惠普工作站上进行的,其配置为:Intel Xeon E5-1620 3.60GHz处理器,18G内存,Windows 7 64位操作系统,开发工具是MATLAB 2012b.表 1给出了实验中使用的基准数据集.

Waveform数据集^[20]包含3类波形,每类各占33%,其40维属性中的后19维都是噪声数据,噪声的均值为0,方差是1.Ringnorm数据集^[20]中的两类样本分别呈现两种不同的正态分布,但是这两种分布也有相互重叠的地方,不易区分.USPS和MNIST都是手写数字数据集^[21],它们各自含有10种不同类型的手写数字图片,每幅图片由一个256维或784维的特征向量表示.Forest Cover Type数据集^[22]来自美国地质调查局(USGS)和美国林务局(USFS),可分成7类数据,每类代表一种森林植被类型.

得到聚类结果后,可以依据簇标签与真实的类标签之间的归一化互信息(normalized mutual information,简称NMI)来衡量算法的聚类准确度^[23].

令${{U}_{c}}={{(u_{1}^{c},...,u_{k}^{c})}^{T}}$表示与聚类得到的簇标签对应的隶属矩阵,${{U}_{t}}={{(u_{1}^{t},...,u_{k}^{t})}^{T}}$表示与真实的类标签对应的隶属矩阵,这两种变量之间的归一化互信息定义为

$NMI({{U}_{c}},{{U}_{t}})=\frac{I({{U}_{c}},{{U}_{t}})}{\sqrt{H({{U}_{c}})\cdot H({{U}_{t}})}},$

其中,I(U_c,U_t)是U_c和U_t之间的互信息,H(U_c)和H(U_t)是信息熵,用于对互信息归一化,使其位于区间[0, 1]内.实践中,常使用下面的公式来估计NMI的值:

$NMI({{U}_{c}},{{U}_{t}})=\frac{\sum\limits_{i=1}^{k}{\sum\limits_{j=1}^{k}{n_{i,j}^{c,t}\log \left( \frac{n\cdot n_{i,j}^{c,t}}{n_{i}^{c}\cdot n_{j}^{t}} \right)}}}{\sqrt{\left( \sum\limits_{i=1}^{k}{n_{i}^{c}\log \frac{n_{i}^{c}}{n}} \right)\left( \sum\limits_{j=1}^{k}{n_{j}^{t}\log \frac{n_{j}^{t}}{n}} \right)}}$

(15)

其中,$n_{i}^{c}={{(u_{i}^{c})}^{T}}{{1}_{n}}$表示簇i中数据点的个数,$n_{j}^{t}={{(u_{j}^{t})}^{T}}{{1}_{n}}$表示类j中数据点的个数,$n_{i,j}^{c,t}={{(u_{i}^{c})}^{T}}u_{j}^{t}$表示属于类j但是被划分到簇i中的数据点的个数.如果聚类结果与真实的类标签完全吻合,则NMI值为1;如果数据被随意划分,则NMI值趋近于0.NMI的值越高,表示聚类的质量越好.为了客观地对比实验结果,各算法统一采用高斯核函数计算数据点之间的相似性,最大迭代次数设为100.由于是随机初始化,算法的每次聚类结果会有小幅波动,因此,将4种算法在每个数据集上都运行20次,并计算其NMI指标的平均值,统计结果见表 2(“-”表示内存不足,实验无法进行).

从表 2中可以看出,NJW-SC算法和WKK-SC算法在Waveform,Ringnorm,USPS这些较小的数据集上可以正常运行;但是当处理MNIST,Forest Cover Type这些大规模的数据集时,会提示内存不足而无法聚类.因为这两种算法都要使用完整的核相似矩阵,空间复杂度都是O(n²),当数据量很大时,需要很大的内存空间来存储相似矩阵.假设每对数据点的相似性都由一个双精度浮点数表示,占8个字节,MNIST数据集有70 000个数据点,这些数据点构成n×n的相似矩阵,大约需要36.5GB的内存;而Forest Cover Type数据集有581 012个数据点,整个相似矩阵占用的内存约为2515.1GB.Nyström-SC算法、MEKA-KC算法和AWKK-SC算法由于采用了近似计算的策略,仅需要使用相似矩阵的一部分,空间复杂度大幅度降低,因而可以在有限的内存里对MNIST,Forest Cover Type进行聚类.而且随着采样点个数的增多,这3种算法的NMI值也逐渐提高.总体来看,Nyström-SC算法在Waveform和USPS数据集上表现较好,MEKA-KC算法对Ringnorm数据集的聚类精度最高,AWKK-SC算法在MNIST和Forest Cover Type数据集上优势明显,这在一定程度上说明了AWKK-SC算法能够较好地处理大规模的数据集.为了进一步对比算法的运行效率,表 3给出了每种算法在各个数据集上20次聚类的平均时间(“-”表示内存不足,实验无法进行).

表 3中,NJW-SC算法的运行时间最长.因为该算法需要构造Laplacian矩阵,并对其特征分解,整个过程的时间复杂度很高.WKK-SC算法利用加权核k-means来优化Normalized Cut的目标函数,不必计算特征向量,所以聚类效率比NJW-SC算法提高不少.但是WKK-SC算法要求使用全部核矩阵进行运算,如果待处理的数据点很多,依然需要花费大量时间.相比之下,仅用部分核矩阵进行近似计算的Nyström-SC,MEKA-KC和AWKK-SC算法在各个数据集上都能很快得到聚类结果,尤其是当数据点的规模达到几十万时,仍然可以流畅运行.Nyström- SC算法与AWKK-SC算法都采用随机抽样策略,虽然提高抽样比例可以改善聚类准确率,但是也会增加程序的计算量,延长聚类时间.仔细观察可以发现:Nyström-SC算法由于需要计算近似的特征向量,在样本点逐渐增多的过程中,其聚类时间的增幅较大;MEKA-KC算法需要花费较长时间来求解最佳的k秩近似核矩阵,当数据集的规模很大时,算法的运行效率较低;而AWKK-SC算法的聚类时间主要与样本数和迭代次数有关,其变化趋势相对平缓,而且大多数情况下能够在更短的时间内完成聚类任务.这也说明AWKK-SC算法的聚类效率更高,适合大数据环境下的数据挖掘工作.

6 结束语

求解谱聚类目标函数的传统方法依赖于计算特征向量,时间和空间复杂性较高,无法处理非常大的数据集.本文从数学上讨论了谱聚类、核k-means与加权核k-means三者之间的联系,证明了Normalized Cut的目标函数与加权核k-means的目标函数是等价的.利用这一等价性,本文设计了一种适用于大数据谱聚类问题的近似加权核k-means算法.该算法一方面采用迭代的方式优化谱聚类的目标函数,避免了对Laplacian矩阵特征分解;另一方面,通过把类中心限制在一个由随机抽样点生成的较小的子空间中,也不用计算整个核矩阵,因此同时降低了谱聚类的计算复杂性和空间需求.理论分析表明:与使用全部核矩阵的加权核k-means相比,近似加权核k-means的期望误差会随着采样点个数的增加而逐渐减小.最后,通过在真实的数据集上进行测试,验证了所提出算法的有效性.虽然采用近似的策略会损失部分精度,但是可以大幅度提高聚类的效率.下一步考虑通过启发式采样或半监督技术来进一步改善所提出算法的性能.

致谢 在此,我们向对本文的工作给予支持和建议的同行表示感谢.