Random Projections for k-means Clustering（ICML 2010）

Random Projections for k-means Clustering

第四十次写博客，本人数学基础不是太好，如果有幸能得到读者指正，感激不尽，希望能借此机会向大家学习。由于K-Means算法本身的时间复杂度很高，特别是在处理大数据集时，因此这篇文章主要介绍一种使用随机映射（Random Projection）降维方法的近似K-Means算法，其他有关于原型聚类算法的讨论可以移步到该类算法的导航页《原型聚类算法综述（原型聚类算法开篇）》。

文章主要内容？

  这篇文章主要讨论K-Means聚类问题中的维数约减技术，并给出如下结论：对于任意$n$n个$d$d维向量构成的矩阵${\bf{A}}\in{\Bbb{R}^{n\times{d}}}$A∈Rn×d，可以以$O\left(nd\lceil{\varepsilon^{-2}k/\log{\left(d\right)}}\rceil\right)$O(nd⌈ε−2k/log(d)⌉)的时间复杂度，将其维数约减至$t=\Omega\left(k/\varepsilon^2\right)$t=Ω(k/ε2)，通过这个映射可以以某一恒定不变的概率得到近似精度为$2+\varepsilon$2+ε的包含$k$k个簇的聚类结果，其中$\varepsilon\in{\left(0,1/3\right)}$ε∈(0,1/3)。上述映射可以通过在矩阵$\bf{A}$A的右边乘上一个随机矩阵${\bf{R}}\in{\Bbb{R}^{d\times{t}}}$R∈Rd×t得到，其中矩阵$\bf{R}$R中的每个元素的值被等概率的设置为$+1/\sqrt{t}$+1/t​或$-1/\sqrt{t}$−1/t​。
  具体来说，本文关注随机映射（Random Projection，或者称为 Johnson-Lindenstrauss嵌入）技术在K-Means聚类中的应用，以输入是$n$n个$d$d维向量的集合为例，该算法首先将这些向量随机映射到$\tilde{d}$d~（$\tilde{d}\leq{d}$d~≤d）维空间中，然后在这些映射后得到的向量上运行K-Means聚类算法。该算法可以以一定的时间复杂度来计算上述嵌入（随机映射）过程，这个时间复杂度与输入数据的规模（即$n$n）呈线性关系，并且得到的最优聚类与在原始数据集上得到的最优聚类之间的近似度为$2+\varepsilon$2+ε。

之前的研究？

  维度约减技术主要可以分为两大方面：a.特征选择（Feature Selection），即通过选择那些可以真实反映数据集特征的特征，来将数据集嵌入到一个低维空间；b.特征抽取（Feature Extraction），即通过人为构造新的特征，来将数据集嵌入到一个低维空间，例如采用原始特征的线性组合。如果将数据集$\tilde{\bf{A}}$A~的最优聚类结果重新带入数据集$\bf{A}$A中，与直接在$\bf{A}$A上得到的最优聚类结果相比差了一个因子$\phi$ϕ（$\phi\geq{1}$ϕ≥1），其中$\tilde{\bf{A}}$A~是$\bf{A}$A通过嵌入$f$f得到的数据集，那么我们就称该嵌入以因子$\phi$ϕ保存了$\bf{A}$A的聚类结构，其中，嵌入$f:\Bbb{R}^d\rightarrow{\Bbb{R}^{\tilde{d}}}$f:Rd→Rd~（$\tilde{d}&lt;d$d~<d）使得$\bf{A}$A中的任意$d$d维向量满足$f\left({\bf{A}}_{\left(i\right)}\right)=\tilde{\bf{A}}_{\left(i\right)}$f(A(i)​)=A~(i)​。
  之前对于K-Means聚类中降维技术的研究包括：a.奇异值分解（SVD）技术将数据集映射为$\tilde{\bf{A}}={\bf{U}}_k{\Sigma}_k\in{\Bbb{R}^{n\times{k}}}$A~=Uk​Σk​∈Rn×k，这种方法以因子2保存了原始数据集中的聚类结构；b.随机映射（Random Reprojection，简称RP）技术将输入的$d$d维向量映射成为$t=\Omega\left(\log{\left(n\right)/\varepsilon^{2}}\right)$t=Ω(log(n)/ε2)维向量，这种方法以因子$1+\varepsilon$1+ε保存了原始数据集中的聚类结构；c.基于奇异值分解的特征选择技术使用SVD找出了$c=\Omega\left(k\log{\left(k/\varepsilon\right)/\varepsilon^{2}}\right)$c=Ω(klog(k/ε)/ε2)个可以真实反映数据集特征的特征，这种映射方式可以某一恒定不变的概率以因子$2+\varepsilon$2+ε保存原始数据集中的聚类结构，下图所示是这三种降维技术与本文提出的降维技术之间的对比，

图1 K-Means聚类中的各种降维技术

其中，RP技术采用随机信号矩阵（Random Sign Matrix）和邮差算法（Mailman Algorithm）来人为构造特征，这些手段会在下文中进行简要介绍。

聚类的最优度和近似度？

  首先引入K-Means聚类问题的另一种表示，给定数据集${\bf{A}}\in{\Bbb{R}^{n\times{d}}}$A∈Rn×d、簇数量$k$k，以及“集群指标矩阵”（Cluster Indicator Matrix），K-Means聚类问题的最优解还可以表示为集群指标矩阵的形式，即
${\bf{X}}_{opt}={\rm{argmin}}_{X\in{}{\chi}}||{\bf{A}}-{\bf{XX}}^{T}{\bf{A}}||^{2}_{F} \tag{1}$Xopt​=argminX∈χ​∣∣A−XXTA∣∣F2​(1)

其中，$\chi$χ是集群指标矩阵的集合，${\bf{X}}_{opt}\in{\Bbb{R}^{n\times{k}}}$Xopt​∈Rn×k是问题得到的最优集群指标矩阵，这样K-Means聚类问题还可以表示为
$F\left({\bf{A}},{\bf{X}}\right)=||{\bf{A}}-{\bf{XX}}^{T}{\bf{A}}||^{2}_{F}$F(A,X)=∣∣A−XXTA∣∣F2​

其中，集群指标矩阵${\bf{X}}\in{\Bbb{R}^{n\times{k}}}$X∈Rn×k的每行只有一个非零元素，该元素用来指示样本点所属的簇，即${\bf{X}}_{ij}=1/\sqrt{z_j}$Xij​=1/zj​​（$i=1,2,\dots,n$i=1,2,…,n、$j=1,2,\dots,k$j=1,2,…,k）表示第$i$i个样本点隶属于第$j$j个簇，$z_j$zj​则是簇中所含样本点总数，明显${\bf{XX}}^{T}={\bf{I}}_{k}\in{diag\left(1,1,\dots,1\right)}$XXT=Ik​∈diag(1,1,…,1)是一个$k$k维单位矩阵。
  上述式（1）即聚类的最优度度量，如果对于数据集$\bf{A}$A和簇数量$k$k，某个算法得到的指示矩阵${\bf{X}}_{\gamma}$Xγ​以至少$1-\delta{\gamma}$1−δγ的概率满足下式，
$||{\bf{A}}-{\bf{X_{\gamma}}}{\bf{X}}_{\gamma}^{T}{\bf{A}}||_{F}^{2}\leq{\gamma\min_{{\bf{X}}\in{\chi}}{||{\bf{A}}-{\bf{XX}}^{T}{\bf{A}}||_{F}^{2}}} \\ =\gamma{||{\bf{A}}-{\bf{X}}_{opt}{\bf{X}}_{opt}^{T}{\bf{A}}||_{F}^{2}} \tag{2}$∣∣A−Xγ​XγT​A∣∣F2​≤γX∈χmin​∣∣A−XXTA∣∣F2​=γ∣∣A−Xopt​XoptT​A∣∣F2​(2)

那么称该算法是K-Means问题的“$\gamma-$γ−近似”（$\gamma-$γ−approximation），其中，$\gamma\geq{1}$γ≥1，$\delta_{\gamma}\in{[0,1)}$δγ​∈[0,1)。这里直接引用《A simplelineartime (1+$\varepsilon$ε)-approximation algorithm fork-means clustering in any dimensions.》 一文中的下述定理，
  对于任意$\varepsilon&#x27;\in{(0,1]}$ε′∈(0,1]，令$\gamma=1+\varepsilon&#x27;$γ=1+ε′，那么K-Means问题的“$\gamma-$γ−近似”算法的时间复杂度是$O\left(2^{\left(k/\varepsilon&#x27;\right)^{O\left(1\right)}}dn\right)$O(2(k/ε′)O(1)dn)。

文章符号定义与线性代数基础？

  矩阵${\bf{A}}\in{{\Bbb{R}}^{n\times{d}}}$A∈Rn×d，聚类数量$k&lt;\min{\{n,d\}}$k<min{n,d}，对于矩阵${\bf{A}}$A的奇异值分解${\bf{A}}={\bf{U}}\Sigma{\bf{V}}^{T}$A=UΣVT，令${\bf{U}}_{K}\in{\Bbb{R}^{n\times{d}}}$UK​∈Rn×d表示$\bf{U}$U中与最大的$k$k个奇异值对应的列向量构成的矩阵，${\bf{V}}_{K}\in{\Bbb{R}^{d\times{k}}}$VK​∈Rd×k表示$\bf{V}$V中与最大的$k$k个奇异值对应的列向量构成的矩阵，$\Sigma_{k}$Σk​是由这$k$k个最大的奇异值构成的对角矩阵。如果$rank\left({\bf{A}}\right)=\rho$rank(A)=ρ，那么${\bf{A}}_{\rho-k}={\bf{A}}-{\bf{A}}_k$Aρ−k​=A−Ak​，其中${\bf{A}}_k={\bf{U}}_k\Sigma_k{\bf{V}}_k^{T}$Ak​=Uk​Σk​VkT​。${\bf{A}}_{\left(i\right)}$A(i)​代表$\bf{A}$A的第$i$i行，$i\in[n]$i∈[n]代表$i\in{\{1,2,\dots,n\}}$i∈{1,2,…,n}，$\bf{A}$A的$\rho$ρ个非负奇异值可以表示为$\sigma_{\left(i\right)}\left({\bf{A}}\right)$σ(i)​(A)（$i\in[\rho]$i∈[ρ]），$||{\bf{A}}||_2$∣∣A∣∣2​和$||{\bf{A}}||_F$∣∣A∣∣F​分别表示矩阵$\bf{A}$A的谱范数（2-范数）和弗罗贝尼乌斯单数（F-范数），${\bf{A}}^{\dagger}$A†表示$\bf{A}$A的伪逆矩阵，即唯一满足下述各个不等式的矩阵，

${\bf{A}}{\bf{A}}^{\dagger}{\bf{A}}={\bf{A}}、{\bf{A}}^{\dagger}{\bf{A}}{\bf{A}}^{\dagger}={\bf{A}}^{\dagger}、\left({\bf{A}}{\bf{A}}^{\dagger}\right)^{T}={\bf{A}}{\bf{A}}^{\dagger}、\left({\bf{A}}^{\dagger}{\bf{A}}\right)^{T}={\bf{A}}^{\dagger}{\bf{A}}$AA†A=A、A†AA†=A†、(AA†)T=AA†、(A†A)T=A†A

${\bf{A}}^{\dagger}$A†的谱范数（即${\bf{A}}^{\dagger}$A†的最大奇异值）与${\bf{A}}$A的谱范数互为倒数，还有一条关于矩阵范数的重要性质，即任意满足矩阵乘法的矩阵${\bf{C}}$C与${\bf{T}}$T，满足不等式$||{\bf{C}}{\bf{T}}||_{F}\leq{||{\bf{C}}||_{F}||{\bf{T}}||_{2}}$∣∣CT∣∣F​≤∣∣C∣∣F​∣∣T∣∣2​。我们将满足${\bf{P}}^{2}={\bf{P}}$P2=P的方阵${\bf{P}}$P作为投影矩阵，令$E[Y]$E[Y]和$Var[Y]$Var[Y]作为随机变量$Y$Y的期望和方差，$P\left(e\right)$P(e)是事件$e$e发生的可能性，并且将“独立同分布”（independent identically distributed）缩写为“i.i.d.”，将“以一定概率”（with probability）缩写为“w.p.”，所有对数都以2为底。

随机投影与随机信号矩阵？

  随机投影（Random Projection，或者称为 Johnson-Lindenstrauss嵌入）定理的一个重要结论是：
  对于任意矩阵${\bf{A}}\in{\Bbb{R}^{n\times{d}}}$A∈Rn×d，可以将其中的$n$n个$d$d维向量线性的投影到$t=\Omega\left(\log{\left(n\right)/\varepsilon^{2}}\right)$t=Ω(log(n)/ε2)维空间中，这种投影使用随机标准正交矩阵，并以因子$1+\varepsilon$1+ε保存了原始空间中任意两点之间的距离。
随后的研究对上述结论的证明过程进行了简化，并证明使用任意高斯随机矩阵（例如矩阵中所有元素均是服从均值为零，方差为$1/\sqrt{t}$1/t​的独立的高斯随机分布）也可以完成随机映射，具体地说，对于任意高斯随机矩阵$\bf{R}$R如下不等式均以较高的概率成立，

$\left(1-\varepsilon\right)||{\bf{A}}_{(i)}-{\bf{A}}_{(j)}||_{2}\leq{||{\bf{A}}_{}{\bf{R}}-{\bf{A}}_{}{\bf{R}}||_{2}\leq{\left(1+\varepsilon\right)||{\bf{A}}_{(i)}-{\bf{A}}_{(j)}||_{2}}}$(1−ε)∣∣A(i)​−A(j)​∣∣2​≤∣∣A​R−A​R∣∣2​≤(1+ε)∣∣A(i)​−A(j)​∣∣2​

由于映射$\tilde{\bf{A}}={\bf{AR}}$A~=AR具有上述性质，因此也以因子$1+\varepsilon$1+ε保存了在原矩阵$\bf{A}$A上的最优K-Means目标函数，Achlioptas证明了，即使使用一个（经过重新放缩后的）“随机信号矩阵”（Random Sign Matrix）进行投影也可以保证上述结论成立，本文使用这个结论对投影过程中涉及到的计算进行计简化。

本文提出的算法？

图2 K-Means聚类的随机投影算法

算法时间复杂度分析？

  假设投影矩阵$\bf{R}$R是通过图1中步骤2所示方法构造的，那么，步骤3中使用邮差算法（Mailman Algorithm）计算矩阵乘法的时间复杂度是$O(nd\lceil\varepsilon^{-2}k/\log{(d)}\rceil)$O(nd⌈ε−2k/log(d)⌉)，当$k=O(\log{(d)})$k=O(log(d))时，上述步骤的时间复杂度近乎为线性$O(nd/\varepsilon^{2})$O(nd/ε2)，这对于处理大规模数据是非常有效的。虽然采用不同的方法构造投影矩阵会使得上述步骤的时间复杂度产生变化，例如采用均值为零、方差为$1/\sqrt{t}$1/t​的随机高斯矩阵得到的时间复杂度是$O(knd/\varepsilon^{2})$O(knd/ε2)，但使用图1所示算法最终产生的聚类结果是近似的。另外，本文使用MATLAB的matrix-matrix BLAS工具来执行步骤3中的矩阵相乘。
  在算法步骤4中，假设$\gamma=1+\varepsilon$γ=1+ε，那么对于任意$\varepsilon\in{(0,1/3)}$ε∈(0,1/3)，“$\gamma-$γ−近似”算法将会以$O\left(2^{\left(k/\varepsilon\right)^{O\left(1\right)}}kn/\varepsilon^{2}\right)$O(2(k/ε)O(1)kn/ε2)的时间复杂度产生一个以因子$2+\varepsilon$2+ε保留原始聚类结构的个解，那么整个算法的时间复杂度可以表示为$O(nd\lceil\varepsilon^{-2}k/\log{(d)}\rceil+2^{\left(k/\varepsilon\right)^{O\left(1\right)}}kn/\varepsilon^{2})$O(nd⌈ε−2k/log(d)⌉+2(k/ε)O(1)kn/ε2)。即使不采用“$\gamma-$γ−近似”算法，而直接在$\tilde{\bf{A}}$A~上运行标准K-Means聚类，步骤4的时间复杂度$O(nk^{2}/\varepsilon^{2})$O(nk2/ε2)也远远小于直接在原始数据集$\bf{A}$A上运行标准K-Means聚类的时间复杂度$O(nkd)$O(nkd)。

---

参考资料

【1】 Arthur, D. . “k-means++ : The advantages of careful seeding.” Proceedings of the eighteenth annual ACM-SIAM symposium on Discrete algorithms, 2007 Society for Industrial and Applied Mathematics, 2007.