聚类旨在根据数据点之间的相似性将其划分到不同的簇，使簇内相似度最大，簇间相似度最小^[1]。传统聚类方法如K-means算法和模糊聚类算法，缺乏处理复杂数据结构的能力，当样本空间非凸时，易陷入局部最优。近年来，谱聚类算法因其可在任意形状空间内进行聚类，并收敛到全局最优，在非凸数据集表现出良好聚类性能，在人脸识别、社区检测、图像分割等领域有着广泛的应用^[2]。

谱聚类将聚类问题转化为图分割问题^[3]，通过寻找最优子图对数据点进行划分。但现存的多数谱聚类算法在处理大规模数据集时，涉及相似度矩阵构造和对应的拉普拉斯矩阵分解，往往需要高昂的时空开销，高计算成本成为制约其在大规模场景应用中的瓶颈^[4]。为实现对大规模数据的高效聚类分析，提升谱聚类算法的扩展性，大量研究方法被提出。Nyström扩展作为一种近似低秩矩阵的有效方法，利用抽样技术，减少特征分解复杂度^[5]。针对抽样策略的选择，Zhan等^[6]设计了一种自适应的Nyström采样方法，通过多次遍历并更新抽样概率，选取更有意义的样本点，以较小的抽样集获得理想的聚类效果。考虑到高斯核参数调整问题，Yang等^[7]提出一种基于层次二部图的谱聚类算法，通过金字塔结构式的多层锚点构造层次二部图降低计算复杂度，此外采用无参数邻接分配策略构造相似度矩阵，避免参数调优过程。蔡登等^[8-9]提出基于地标表示的研究，择取具有代表性的数据点作为地标点，通过地标点的稀疏线性组合近似构造图相似度矩阵，从而有效降低谱嵌入的计算复杂度，该算法一定程度上解决了矩阵存储开销问题，但其随机采样确定地标点的方式会造成聚类结果的不稳定，故而地标点的选取为算法的关键。叶茂等^[10]利用由近似奇异向量的长度计算的抽样概率来进行抽样，以保证聚类结果的稳定性和精度。Zhang等^[11]提出一种基于增量视角的抽样框架，每一个要抽样的点是由先前选定的地标点决定，在地标点之间建立显式关系。但这些改进的方法在择取地标点时都忽略了节点的拓扑属性，其可有效描述亲和图整体结构特性，对于捕捉高维数据空间的拓扑信息有着重要意义。故邓思雨等^[12]将PageRank分值作为样本信息量的度量指标，Rafailid等^[13]提出一种基于地标选择的快速谱聚类算法，根据加权PageRank算法选择亲和图中最重要的节点作为地标点。Liu等^[14]将数据点视为web页面，数据点之间的距离类似为链路间的权重，采用PageRank算法评估数据点的重要性，选择代表点。该方法能够区分球状和非球状的簇，且减少噪声点的负面影响。以上这些改进的谱聚类算法虽然降低了计算复杂度，仍然需要对矩阵进行特征分解。Jia等^[15]设计了一种无需特征分解的大规模近似Ncut谱聚类算法，一方面通过对数据点的采样推断数据集的全局特征，减少归一化切割的空间需求，另一方面利用近似加权核k-均值优化Ncut的目标函数，避免拉普拉斯矩阵的直接特征分解。Tian等^[16]证明了自编码器和谱聚类之间的相似性，利用栈式自编码器代替特征分解，有效降低计算复杂度。Banijamali等^[17]提出将深层结构与地标表示相结合，实现快速且准确的聚类。但上述提及的谱聚类算法均使用传统基于欧氏距离度量方案，并不能完全反映数据复杂的空间分布特性。光俊叶等^[18]提出融合欧氏距离和Kendall Tau距离的谱聚类方法，充分利用不同相似性度量从不同角度获取数据信息的优势，全面反映底层结构信息。

本文在地标点采样^[8-9]的理论基础上，采用基于网页质量的PageRank算法计算节点的重要性，作为抽样策略依据代替随机抽样，选取代表点作为地标点，通过稀疏表示近似获得图相似度矩阵，以降低存储开销。同时为全面反映数据底层结构信息，将欧氏距离与Kendall Tau距离两种度量方案融合，用局部标准差代替特定的缩放尺度，构造融合多度量方式的自适应相似度矩阵，以提高聚类精度。此外，用栈式自编码器代替特征分解，在联合学习框架中进行嵌入表示学习和聚类，避免微调环节覆盖之前所得最优参数，从而进一步提高聚类精度。

1 相关算法理论 1.1 基于QPR的地标选择

PageRank算法采用平均分配思想，将节点权重划分给入链节点，忽略节点质量存在差异性，改进后的加权PageRank算法通过传入链路权重 ${{w}}_{\left( {{{a,b}}} \right)}^{\rm{in} }$ 及传出链路权重 ${{w}}_{\left( {{{a,b}}} \right)}^{\rm{out} }$ 2个参数^[13]，依据网络结构有区分地进行权值分配，但该参数在迭代过程中是固定存在的，存在局限性。基于网页质量的PageRank算法^[19](page quality based PageRank algorithm, QPR)定义了相对质量 ${{Q} } (a)$ ，用于迭代过程协助分配PR值，动态评估每个节点的质量，其计算公式为

$\begin{array}{l} {{Q} } (a) = \dfrac{{{\rm{PR}} (a)}}{{\Bigg(\displaystyle \sum\limits_{b \in {{B}} (a)} {{\rm{PR}} (b)\Bigg)\Bigg/{\rm{MPR}} (a)} }} \\ {\rm{MPR}} (a) = \mathop {\max }\limits_{b \in {{B}} (a)} {\rm{PR}} (b) \end{array} $

(1)

式中： ${{B}} (a)$ 表示节点 $a$ 的入链节点集合； ${\rm{MPR}} (a)$ 表示 $a$ 的入链节点集合中的最大PR值。相对质量 ${{Q} } (a)$ 将迭代过程中 $a$ 的入链集合与其本身的PR值信息结合，对于具有相同PR值的节点，拥有较大入链PR值的节点将得到加权，从而在下一轮迭代中具有更高PR值。随着多轮迭代，节点各自的PR值出现差距，高质量的节点能较快积累PR值，获得较高排名。该算法计算PageRank的公式为

$\begin{array}{c} {\rm{PR}} (a) = d\displaystyle \sum\limits_{b \in {\rm{B}} (a)} {{\rm{PR}} (b)\Pr (b \to a) + \frac{{(1 - d)}}{n}} \\ \Pr (b \to a) = {\rm{Q} } (a)/\displaystyle \sum\limits_{c \in {\rm{F}} (b)} {{\rm{Q} } (c)} \end{array} $

(2)

式中： $\Pr (b \to a)$ 表示节点 $b$ 分配给节点 $a$ 的PR值比重， ${\rm{F}} (b)$ 表示节点 $b$ 的出链节点集合，即 ${\rm{F}} (b) = \{ a|(b,a) \in E\} $ ， $d$ 是一个阻尼因子，通常设置为0.85， $n$ 为节点个数。当所有节点满足式(3)时，迭代过程停止：

$\dfrac{{{\rm{PR}} {{(a)}_{{\rm{iter}} - 1}} - {\rm{PR}} {{(a)}_{{\rm{iter}}}}}}{{{\rm{PR}} {{(a)}_{{\rm{iter}}}}}} \leqslant \beta $

(3)

该式表示上一次迭代与当前迭代之间的归一化差值，其中， $\beta $ 是预先设置的停止阈值。当满足该式时，停止迭代，选择PR值最高的 $p$ 节点作为地标点。

1.2 自适应谱聚类算法

相似度矩阵的构造是谱聚类算法的一个重要步骤，而高斯核函数是构造相似度矩阵最常用的度量方法之一。给定数据集 ${{X}} = [{x_1},{x_2} ,\cdots ,{x_n}] \in {R^{n \times d}}, i = 1,2, \cdots ,n$ ，拟将其划分为 $K$ 簇。谱聚类算法将聚类问题视为无向图 ${{G}} = ({{V}},{{E}})$ 的多路划分问题，其中顶点集 ${{V}} = \{ {x_1},{x_2} ,\cdots ,{x_n}\}$ 为所有样本集合， ${{E}} = \{ {A_{ij}}\} $ 表示顶点间边的权重集合， ${{A}} = {({A_{ij}})} \in {R^{n \times n}}$ ， ${i,j = 1,2, \cdots ,n} $ 即为所需构造的相似度矩阵，其元素被定义为

${A_{ij}} = \left\{ \begin{array}{l} \exp \Bigg(\dfrac{{ - {\rho ^2}({x_i},{x_j})}}{{2{\sigma ^2}}}\Bigg),\;\;\;i \ne j \\ 0,\;\;\;i = j \end{array} \right.$

式中： $\;\rho (x{}_i,{x_j})$ 指样本 ${x_i}$ 与 ${x_j}$ 之间的一种距离度量方法，参数 $\sigma $ 是手动给出的固定缩放参数，并不具有自适应性，难以体现数据的真实分布^[20]。自适应谱聚类算法^[21]针对全局尺度参数不能有效反映数据集真实分布信息，提出了局部尺度参数的概念，根据样本 ${x_i}$ 的第 $k$ 个近邻定义样本 ${x_i}$ 的局部尺度参数 ${\sigma _i} = \rho ({x_i},{x_k})$ ，该算法定义相似度矩阵元素为

${A_{ij}} = \left\{ \begin{array}{l} \exp \Bigg(\dfrac{{ - {\rho ^2}({x_i},{x_j})}}{{2{\sigma _i}{\sigma _j}}}\Bigg),\;\;\;i \ne j \\ 0,\;\;\;i = j \end{array} \right.$

特定的缩放参数 ${\sigma _i}$ 可根据样本 ${x_i}$ 和 ${x_j}$ 的邻域点进行自调。该算法虽然在一定程度上克服了自适应谱聚类算法的缺点，但也有一定的局限性。其聚类结果可能受到异常值影响，因为局部尺度参数 ${\sigma _i}$ 可能扭曲了离群值。为避免该问题，文献[22]提出一种局部标准差谱聚类算法，定义局部标准差 ${\sigma _{{\rm{std}}\_i}} = \sqrt {\dfrac{1}{{k - 1}}\displaystyle \sum\limits_{t = 1}^k {{\rho ^2}({x_i},{x_t})} }$ ，其表示样本 ${x_i}$ 的前 $k$ 个近邻的局部标准差，相似度矩阵的元素则如式(4)所示，尽可能地反映数据集的原始分布：

${A_{ij}} = \left\{ \begin{array}{l} \exp \Bigg(\dfrac{{ - {\rho ^2}({x_i},{x_j})}}{{2{\sigma _{{\rm{std}}\_i}}{\sigma _{{\rm{std}}\_j}}}}\Bigg),\;\;\;i \ne j \\ 0,\;\;\;i = j \end{array} \right.$

(4)

2 本文算法 2.1 结合度量融合与地标表示的相似度矩阵

由1.2节可以看出谱聚类算法中所构造的相似度矩阵维度为 $n \times n$ ，对于大规模数据集无法实现高效的谱嵌入，空间计算复杂度较高。对此文献[8-9]提出选取 $p \ll n$ 个地标点作为特征代表点，将所有样本近似表示为这些地标点的线性组合，有效利用稀疏表示矩阵近似获得图相似度矩阵。本文采用1.1节的基于QPR的地标选择方式，选择PR值最高的 $p$ 个节点作为地标点。对于任意给定样本 ${x_i}$ ，其近似样本 ${\hat x_i}$ 可表示为

${\hat x_i} = \displaystyle \sum\limits_{j = 1}^p {{Z_{ji}}{y_j}} $

式中： ${y_j}$ 是地标点矩阵 ${{Y}} \in {R^{d \times p}}$ 的列向量， ${Z_{ij}}$ 为稀疏表示矩阵 ${{Z}} \in {R^{p \times n}}$ 的第i行、第j列元素，表示所选取的 $p$ 个地标点与所有样本点之间的成对相似性。根据稀疏编码策略， ${x_i}$ 与 ${y_j}$ 越接近， ${Z_{ji}}$ 越大，而若 ${y_j}$ 不在样本 ${x_i}$ 的 $k$ 近邻内，则将 ${Z_{ji}}$ 置为0，通过引入 $k$ 近邻点来度量点对之间的局部亲和力。相应稀疏表示矩阵 ${{Z}}$ 的元素为

${Z_{ji}} = \left\{ \begin{array}{l} \dfrac{{{A_{ji}}}}{{\displaystyle \sum\limits_{m \in {\rm{KNN}}({x_i})} {{A_{mi}}} }},\;\;\;j \in {\rm{KNN}}({x_i})\\ 0,\;\;\;j \notin {\rm{KNN}}({x_i}) \end{array} \right.$

(5)

式中： ${A_{ji}}$ 表示样本 ${x_i}$ 与地标点 ${y_j}$ 之间的相似度，采用式(4)计算方式，尽可能反映数据集的原始分布。在计算点对相似性时往往采用欧氏距离传统度量方式，一些邻域有用信息容易被忽略。Kendall Tau距离是衡量两个等级列表之间两两不一致的数量，即求两个排列之间的逆序数，反映了两个排列的相似程度。对于两排列 ${L_1} = ({L_{11}},{L_{21}}, \cdots , {L_{n1}})$ 及 ${L_2} = ({L_{12}},{L_{22}}, \cdots ,{L_{n2}})$ ，它们之间的Kendall Tau距离定义为

$ {\rm{KT}}({L_1},{L_2}) = \left| {\{ (i,j):i<j,(({L_{i1}}<{L_{j1}}) \wedge ({L_{i2}}<{L_{j2}})){\rm{or}}(({L_{i1}}>{L_{j1}}) \wedge ({L_{i2}}>{L_{j2}}))\} } \right| $

(6)

式中： ${L_{i1}}$ 和 ${L_{i2}}$ 表示 ${L_1}$ 和 ${L_2}$ 中第 $i$ 个元素各自的序号， $\left| \cdot \right|$ 为集合中元素的个数。若两排列相同，则 ${\rm{KT}} ({L_1},{L_2})$ 为0；若两排列互为逆序， ${\rm{KT}} ({L_1},{L_2})$ 为 $n(n - 1)/2$ 。通过除以 $n(n - 1)/2$ 进行规范化处理，使 ${\rm{KT}} ({L_1},{L_2})$ 位于区间 $[0,1]$ 内。

若将欧氏距离及Kendall Tau距离两种度量方式融合，生成增强的相似度矩阵，能在反映点对之间地理相似性的同时考虑点之间的拓扑逻辑相似性。将点 ${x_i}$ 和 ${x_j}$ 之间的欧氏距离表示为 ${\rho _1}({x_i},{x_j})$ ,Kendall Tau距离表示为 ${\rho _2}({x_i},{x_j})$ 。在基于欧氏距离对样本进行排序，然后利用Kendall Tau距离度量两个样本之间的全局关系，全面反映数据的底层结构信息。对于样本 ${x_i}$ ，基于欧氏距离计算与 $p$ 个地标点之间的相似度，可得一有序序列： ${{\rm{List}} _i} = [{\rho _1}({x_1},{x_i}),{\rho _1}({x_2},{x_i}), \cdots , {\rho _1}({x_m},{x_i}), \cdots ,{\rho _1}({x_p},{x_i})]$ ，同样的对于样本 ${x_j}$ 可得有序序列 ${{\rm{List}} _j}$ ，那么结合式(6)，样本 ${x_i}$ 和 ${x_j}$ 之间的Kendall Tau距离则为

${\rho _2}({x_i},{x_j}) = {\rm{KT}} ({{\rm{List}} _i},{{\rm{List}} _j})$

从度量融合的角度来看，大多数主流的思想都是通过线性组合多个亲和图来获得互补的相似信息^[23]，然而这个线性组合对分配给每个度量的权重很敏感，且来自多个亲和图的互补信息并非线性相关的。为探索两种距离度量方案相互融合的思想，综合考虑相邻点的几何分布和拓扑分布，受到联合训练算法的启发，采用基于消息传递理论和交叉扩散过程的融合方法。以完全相似矩阵作为初始状态矩阵，即对式(5)不考虑 $k$ 近邻点而构建相似矩阵，该初始状态矩阵 ${{F}}$ 携带关于每个地标点与其他所有样本点之间相似性的完整信息，对应元素为

${F_{ji}} = \frac{{{A_{ji}}}}{{\displaystyle \sum\limits_{m = 1}^p {{A_{mi}}} }}$

(7)

式(7)分别采用欧氏距离、Kendall Tau距离度量方法计算获得的完全相似矩阵表示为 ${{{F}}^{(1)}}$ 、 ${{{F}}^{(2)}}$ 。稀疏表示矩阵实际上是完全相似矩阵的KNN图，在扩散过程中为提高计算效率，采用稀疏表示矩阵作为核矩阵进行交叉扩散。将稀疏表示矩阵 ${{Z}}$ 采用欧氏距离、Kendall Tau距离分别表示为 ${{{Z}}^{(1)}}$ 和 ${{{Z}}^{(2)}}$ 。首先将 ${{F}}_{h = 0}^{(1)} = {{{F}}^{(1)}}$ 和 ${{F}}_{h = 0}^{(2)} = {{{F}}^{(2)}}$ 视为迭代过程第一步中的初始两个状态矩阵( $h = 0$ )，而度量融合的关键步骤是迭代更新对应于每个度量的稀疏表示矩阵

$\begin{array}{l} {{F}}_{h + 1}^{(1)} = {{{Z}}^{(1)}} \times {{F}}_h^{(2)} \times {({{{Z}}^{(1)}})^{\mathop{\rm T}\nolimits} }\\ {{F}}_{h + 1}^{(2)} = {{{Z}}^{(2)}} \times {{F}}_h^{(1)} \times {({{{Z}}^{(2)}})^{\rm{T}}} \end{array}$

(8)

式中： ${{F}}_h^{(1)}$ 表示基于欧氏距离度量的稀疏表示矩阵在 $h$ 步迭代更新后的状态矩阵， ${{F}}_h^{(2)}$ 表示基于Kendall Tau距离度量的矩阵在 $h$ 步迭代更新后的状态矩阵。此步骤每次更新状态矩阵时，生成两个并行交换扩散过程，在 $h$ 步后，融合两种度量方式的稀疏表示矩阵 ${{Z}}$ ：

${{Z}} = \dfrac{{{{F}}_h^{(1)} + {{F}}_h^{(2)}}}{2}$

(9)

从概率角度分析该交叉扩散过程，对于状态矩阵 ${{F}}_h^{(1)}$ ，依据扩散映射^[24]可定义进行 $h$ 步迭代时的扩散距离 $M_h^{(1)}(i,j) = \left\| {{{F}}_h^{(1)}(i,:) - {{F}}_h^{(1)}(j,:)} \right\|$ ，扩散过程即是将数据空间映射到扩散空间 $\Re _h^{(1)}$ 的过程，本质上每个样本都可由自身与其他数据点间相似度表示。为实现核矩阵之间的融合，对于样本 $x_h^{(1)} \in \Re _h^{(1)}$ ，引入线性算子 ${{{Z}}^{(2)}}$ ，经过 $x_{h + 1}^{(2)} = {{{Z}}^{(2)}}x_h^{(1)} + \varepsilon $ ( $\varepsilon $ 为白噪声)线性运算后，由扩散映射性质可获得 $x_{h{{ + }}1}^{(2)}$ 的边缘分布^[23]，同理可得 $x_{h + 1}^{(1)}$ 。该交叉扩散过程的实质并非在原始数据空间中进行线性投影，而是在扩散空间中进行迭代线性运算的，这样对样本的噪声和规模具有较强鲁棒性，且融合了整个样本集的相似流形的内在结构。

在由式(9)获得稀疏表示矩阵 ${{Z}}$ 后，通过该矩阵来近似获得数据 ${{X}}$ 的规范化图亲和矩阵，即： ${{\hat G}} = {{{\hat Z}}^{\rm{T}} }{{\hat Z}} \in {R^{n \times n}}$ ，其中 ${{\hat Z}} = {{{D}}^{ - 1/2}}{{Z}}$ ， ${{D}}$ 为度矩阵，其元素一般为 ${D_{ii}} = \displaystyle\sum\nolimits_j {{Z_{ij}}} $ ，为进一步降低计算复杂度，采用文献[2]计算度矩阵元素 ${D_{ii}} = {\rm{diag}} ({{\hat Z}}_i^{\rm{T}} {{{\hat Z}}^s})$ ，其中 ${{{\hat Z}}^s} = \displaystyle\sum\limits_{j = 1}^n {{{\hat Z}_j}} $ ，为 $p \times 1$ 的向量，其第 $k$ 个元素为 ${{\hat Z}}$ 的第 $k$ 行元素之和，该计算 ${{D}}$ 的时间复杂度为 $O(np)$ 。结合度量融合与稀疏表示的拉普拉斯矩阵表示为

${{L}} = {{{D}}^{ - \frac{1}{2}}}{{\hat G}}{{{D}}^{ - \frac{1}{2}}} = {{{D}}^{ - \frac{1}{2}}}{{{\hat Z}}^{\rm{T}} }{{\hat Z}}{{{D}}^{ - \frac{1}{2}}}$

(10)

然后将 ${{{D}}^{ - \frac{1}{2}}}{{{\hat Z}}^{\rm{T}}}$ 作为栈式自编码器的输入，通过编码和解码重构学习数据潜在特征表示。

2.2 深度嵌入谱聚类

假定聚类任务为 $N$ 个样本， ${{X}} = [{x_1},{x_2} ,\cdots ,{x_n}], {x_i} \in {R^{{d_x}}},i = 1,2, \cdots ,n$ 划分为 $K$ 簇，每个都由聚类中心 ${\mu _j},j = 1,2, \cdots ,K$ 表示。为避免直接在样本空间 ${{X}}$ 上进行聚类，首先使用非线性映射进行数据转换，嵌入函数 ${\varphi _W}:{{X}} \to {{H}}$ ，将原始数据映射到潜在特征空间 ${{H}} = [{h_1},{h_2} ,\cdots ,{h_n}],{h_i} \in {R^{{d_h}}}$ ，与输入样本相比具有较低维度( ${d_h} \ll {d_x}$ )，再经过非线性映射进行数据重构，解码器映射函数 ${g_{{W'}}}:{{H}} \to {{\hat X}}$ , ${{\hat X}} = [{\hat x_1},{\hat x_2} ,\cdots ,{\hat x_n}]$ ， ${\hat x_i}$ 表示样本 ${x_i}$ 的数据重构。本文旨在寻找良好的嵌入特征表示 $\{ {h_i}\} _{i = 1}^n$ 使其更适合于聚类，所构造的损失目标函数包含自编码器重构损失及聚类损失，利用自编码器以无监督方式学习表示，所习得特征能最大限度地保留数据的固有局部结构，聚类损失用于分散嵌入点。与需要分层预训练以及非联合嵌入聚类的学习模型不同，同时避免微调步骤覆盖预训练获得的参数，采用联合聚类和重构损失函数的目标函数，对所有网络层进行端到端优化训练，迭代优化自编码器参数及聚类中心，损失目标函数 $J$ 定义为

$J = {J_r} + \gamma {J_c}$

(11)

其中 ${J_r}$ 和 ${J_c}$ 分别表示为重构损失和聚类损失， $\gamma (0 < \gamma < 1)$ 用来调节潜在特征空间的失真程度。基于2.1节，将构造的相似度矩阵作为栈式自编码的输入，在进行聚类之前，对嵌入特征表示 $\{ {h_i}\} _{i = 1}^n$ 进行预训练，执行k-means算法，获得初始自编码器参数 $\{ W,{W{'}}\}$ 及聚类中心 $\left\{ {{\mu _j}} \right\}_{j = 1}^K$ 。

对于训练所获得的嵌入特征表示，在文献[23]引入了由t-SNE衡量的软标签分布 ${p_{ik}}$ ，用来衡量潜在空间样本 ${h_i}$ 与聚类中心 ${\mu _j}$ 之间的相似度：

${p_{ik}} = \dfrac{{{{(1 + {{\left\| {{h_i} - {m_j}} \right\|}^2})}^{ - 1}}}}{{\displaystyle \sum\limits_{j = 1}^K {{{(1 + {{\left\| {{h_i} - {m_j}} \right\|}^2})}^{ - 1}}} }},\quad j = 1,2, \cdots ,K$

(12)

式中 ${p_{ik}}$ 表示潜在空间样本 ${h_i}$ 属于 ${\mu _j}$ 所在簇的概率。结合最小化相关熵(KL散度)来降低软标签分布 ${{{P}}}$ 与辅助目标分布 ${{{Q}}}$ 之间的差异，将聚类损失目标函数 ${J_c}$ 定义为

${J_c} = {\rm{KL}} ({{Q}}\parallel {{P}}) = \dfrac{1}{N}\displaystyle \sum\limits_i {\displaystyle \sum\limits_k {{q_{ik}}\lg \frac{{{q_{ik}}}}{{{p_{ik}}}}} } $

(13)

其中，辅助目标分布分量 ${q_{ik}}$ 是由 ${p_{ik}}$ 得到的，通过计算 $\dfrac{{\partial {J_c}}}{{\partial {p_{ik}}}}$ 使其为零，求取相应封闭解获得 ${q_{ik}}$ ：

${q_{ik}} = \dfrac{{p_{ik}^2/\displaystyle \sum\limits_{{i^{'}}} {{p_{{i^{'}}k}}} }}{{\displaystyle \sum\limits_{{k^{'}}} {\Bigg(p_{i{k^{'}}}^2/\displaystyle \sum\limits_{{i^{'}}} {{p_{{i^{'}}{k^{'}}}}\Bigg)} } }}$

(14)

最小化聚类损失目标函数 ${J_c}$ 可视为自训练过程，此外，重构误差函数 ${J_r}$ 采用均值误差测量(mean squared error, MSE)：

${J_r} =\displaystyle \sum\limits_{i = 1}^n {\left\| {{x_i} - {g_{{W'}}}({h_i})} \right\|_2^2} $

(15)

关于自编码器参数 $\{ W,{W'}\} $ 及聚类中心 $\left\{ {{\mu _j}} \right\}_{j = 1}^K$ 的更新，采用mini-batch随机梯度算法反向传播逐层训练网络。 ${J_c}$ 关于嵌入特征表示 ${h_i}$ 和聚类中心 ${\mu _j}$ 的梯度计算为

$\dfrac{{\partial {J_c}}}{{\partial {h_i}}} = 2\displaystyle \sum\limits_{j = 1}^K {{{(1 + {{\left\| {{h_i} - {\mu _j}} \right\|}^2})}^{ - 1}}({q_{ij}} - {p_{ij}})(} {h_i} - {\mu _j})$

$\dfrac{{\partial {J_c}}}{{\partial {\mu _i}}} = 2\displaystyle \sum\limits_{j = 1}^K {{{(1 + {{\left\| {{h_i} - {\mu _j}} \right\|}^2})}^{ - 1}}({p_{ij}} - {q_{ij}})(} {h_i} - {\mu _j})$

将 $\dfrac{{\partial {J_c}}}{{\partial {h_i}}}$ 反向传播给栈式自编码器，利用下层神经元梯度由上层神经元残差导出的规律，自上而下反向逐层计算出 $\dfrac{{\partial {J_c}}}{{\partial W}}$ 及 $\dfrac{{\partial {J_c}}}{{\partial {W'}}}$ 。给定小批量样本数目 $m$ 及学习率 $\lambda $ ，编码器参数 $W$ 、解码器参数 ${W'}$ 及聚类中心 ${\mu _j}$ 更新公式为

$W = W - \dfrac{\lambda }{m}\displaystyle \sum\limits_{i = 1}^m {\Bigg(\dfrac{{\partial {J_r}}}{{\partial W}}} + \gamma \frac{{\partial {J_c}}}{{\partial W}}\Bigg)$

(16)

${W'} = {W'} - \dfrac{\lambda }{m}\displaystyle \sum\limits_{i = 1}^m {\dfrac{{\partial {J_r}}}{{\partial {W'}}}} $

(17)

${\mu _j} = {\mu _j} - \dfrac{\lambda }{m}\displaystyle \sum\limits_{j = 1}^k {\dfrac{{\partial {J_r}}}{{\partial {\mu _j}}}} $

(18)

2.3 算法步骤

算法　结合度量融合和地标表示的自编码谱聚类算法

输入　输入数据X，聚类数目 $K$ ，地标点个数 $p$ ，停止阈值 $\;\beta $ ，近邻点数目 $k$ ，学习率 $\lambda $ ，迭代次数 $h$ ，最大迭代次数 ${\rm{maxiter}}$ ，停止收敛阈值 ${\rm{Th}}$ 。

输出　自编码器参数 $\{ W,{W{'}}\}$ 及聚类中心 $\left\{ {{\mu _j}} \right\}_{j = 1}^K$ 及样本类标签Label。

1)根据式(1)、(2)计算节点PR值，当所有节点满足式(3)时，停止迭代，PR值降序排列，择取前 $p$ 个高PR值节点作为地标点；

2)依据式(4)成对相似度计算，引入 $k$ 近邻点，以欧氏距离、Kendall Tau距离2种度量方式，通过式(5)分别计算相应的稀疏表示矩阵 ${{{Z}}^{(1)}}$ 和 ${{{Z}}^{(2)}}$ 。同时，通过式(7)计算对应于2种度量方式的完全相似矩阵 ${{{F}}^{(1)}}$ 、 ${{{F}}^{(2)}}$ ；

3)将 ${{{F}}^{(1)}}$ 、 ${{{F}}^{(2)}}$ 作为初始状态矩阵， ${{{Z}}^{(1)}}$ 、 ${{{Z}}^{(2)}}$ 作为核矩阵，根据式(8)进行矩阵融合，迭代更新对应于每个度量的稀疏表示矩阵，在 $h$ 步后结合式(9)获得融合2种度量方式的稀疏表示矩阵 ${{Z}}$ ；

4)通过式(10)构造拉普拉斯矩阵，将 ${{{D}}^{ - \frac{1}{2}}}{{{\hat Z}}^{\rm{T}}}$ 作为栈式自编码器的输入，前向训练网络获得初始化数据潜在特征表示 $\{ {h_i}\} _{i = 1}^n$ ；

5)对嵌入特征表示 $\{ {h_i}\} _{i = 1}^n$ 进行预训练，执行k-means算法，获得初始自编码器参数 $\{ W,{W'}\} $ 及聚类中心 $\left\{ {{\mu _j}} \right\}_{j = 1}^K$ ；

6)若迭代次数大于 ${\rm{maxiter}}$ ，转至10)，否则重复执行7) ~ 9)；

7)根据式(12)和式(14)计算 ${p_{ik}}$ 和 ${q_{ik}}$ ；

8)根据式(11)、(13)和(15)计算损失目标函数；

9)采用mini-batch随机梯度算法反向传播逐层训练网络，根据式(16)~(18)更新整个网络参数 $\{ W,{W'}\} $ 及聚类中心 $\left\{ {{\mu _j}} \right\}_{j = 1}^K$ ；

10)根据软标签分布 ${{{P}}}$ 来分配样本类标签： ${\rm{Label}} (i) = \mathop {\arg \max }\limits_{k \in 1,2, \cdots ,K} {p_{ik}}$ 。

2.4 算法复杂度分析

首先，从 $n$ 个数据样本中选择 $p$ 个地标点，该步骤时间复杂度为 $O(tpn)$ ， $t$ 为迭代次数。相似度矩阵构造部分涉及了2种度量方案的融合以及构造稀疏表示矩阵^[17]，时间复杂度为 $O(pn{\log _2}n)$ 。本文算法用栈式自编码器代替拉普拉斯矩阵特征分解部分，时间复杂度为 $O(n{D^2} + nvK)$ ， $K$ 为聚类数目， $D$ 为设置的隐藏层最大单元数， $v$ 为数据潜在特征表示维度，一般 $K \leqslant v \leqslant D$ ，所提算法整体时间复杂度为 $O(tpn + pn{\log _2}n + n{D^2})$ 。基于地标表示的快速谱聚类算法^[9]时间复杂度为 $O(tpn + pn + {p^3} + {p^2}n)$ ，由于 $p \ll n$ ，时间复杂度可写为 $O(tpn + pn + {p^2}n)$ 。 $p$ 和 $D$ 取值相当且远小于 $n$ ，因而本文算法复杂度略高于基于地标表示的快速谱聚类算法时间复杂度，但远远低于基于自编码器的谱聚类算法DEC^[25]的时间复杂度 $O({n^2} + n{D^2})$ 、GraEn ^[16]的时间复杂度 $O({n^2} + nKD)$ 。

3 实验及分析结果 3.1 实验环境

为证明本文算法适用于大规模数据集，选取如下5个数据集：MNIST、USPS、COIL100、CoverType和CIFAR10。这5个数据集包含不同类型的图像，如手写数字、英文字母表、物体图像和森林植被类等，样本数较大，可以有效验证本文算法在大规模数据集上的聚类性能。表1为相关数据集的具体信息。仿真实验硬件环境为Intel(R) Core(TM) i7-6700 CPU @3.40 GHz；16 GB RAM；操作系统：Windows7；编程语言：Python 3.5。

本文算法与文献[8-9]基于地标表示的快速谱聚类算法LSC-R和LSC-K，未考虑地标点的基于自编码的谱聚类算法DEC^[25]及GraEn^[16]、文献[17]中自编码器与地标点结合的快速谱聚类算法SCAL-R及SCAL-K进行对比。LSC-R及SCAL-R是通过随机采样获取 $p$ 个地标点，LSC-K及SCAL-K是采用k-means算法获得 $p$ 个地标点的。

为保证算法之间的公平对比，具体参数如下设置：文献[8-9]、[17]及本文算法涉及的地标点个数 $p$ 统一设置为1 000，稀疏表示矩阵构造涉及的近邻点数目 $k$ 设置为5，所有算法涉及k-means部分迭代次数为500。文献[16]、[17]、[25]及本文算法涉及的自编码器部分，编码器维度为 $p - 500 - 500 - 2\;000 - 10$ ，解码器则为 $10{\rm{ - }}2\;000{\rm{ - }}500{\rm{ - }} 500{\rm{ - }}p$ ，非线性激活函数统一采用线性校正单元(rectified linear units，ReLU)。本文算法中涉及的停止阈值 $\beta = {10^{ - 3}}$ ，距离度量融合步骤迭代次数设置为 $h = 10$ ，参数优化更新部分采用mini-batch随机梯度算法，小批量样本数目 $m$ 为256，学习率 $\lambda = 0.1$ ；在预训练阶段，执行k-means算法20次选取最佳结果，初始化自编码器参数和聚类中心；收敛阈值设置为 ${\rm{Th}} = 0.001$ ，最大迭代次数 ${\rm{maxiter}}$ 为20 000，调节潜在特征空间的失真程度 $\gamma = 0.1$ 。

3.2 实验分析

为了比较这些算法性能，通过聚类结果和样本真实标签进行对比。采用聚类准确率(ACC)和标准化互信息(NMI)来作为度量指标进行评估比较。ACC和NMI这2种度量标准在[0,1]之间取值，值越大聚类性能越好。各种算法在5个大规模数据集的实验结果如表2、3所示。

通过表2、3两种指标下的各种算法性能度量对比，可以看出本文提出的结合度量融合和地标表示的自编码谱聚类算法在MNIST、LetterRec、COIL100、CoverType及CIFAR10这5个大规模数据集上，相对于基于地标表示的快速谱聚类算法LSC-R和LSC-K，未考虑地标点的基于自编码的谱聚类算法DEC及GraEn，自编码器与地标点结合的快速谱聚类算法SCAL-R及SCAL-K，在ACC和NMI两种指标下均表现出了更为理想的聚类性能。同时，还可以得出如下结论：

1)本文算法在稀疏表示矩阵构造部分由于引入网页质量进行地标点择取，并进行度量融合计算点对相似度，在算法复杂度方面，相对于旨在提升谱聚类算法效率的LSC-R和LSC-K处于劣势。从运行时间上看，在MNIST数据集上，LSC-R及LSC-K算法仅需10 s左右，本文算法则需200~300 s的运行时间，但相对于DEC及GraEn，本文算法运行时间降低了1 000~2 000 s，有了很大的提升。且从表2中可以看出，就ACC，本文算法与LSC-K算法相比，在MNIST、LetterRec、COIL100、CoverType及CIFAR10F分别提升了17.42%、16.57%、34.47%、7.75%、10.89%，与随机抽样择取地标点的LSC-R算法相比，优势更大。同样地，从表3看出，本文算法NMI值也有了明显的提升。而与SCAL-R及SCAL-K算法相比，本文融合两种距离度量多角度描述数据集的结构信息，使其在ACC、NMI两种指标上均处于优势。

2)相比未考虑地标点的基于自编码的谱聚类算法DEC及GraEn，本文算法ACC和NMI也有了一定的提升，如与DEC算法对比，本文算法在LetterRec数据集上ACC、NMI分别提升了17.15%、5.84%。说明本文构造的结合度量融合和地标表示的相似度矩阵更全面地反映了数据之间的结构信息，且运行效率大大得到提升。

3)此外，间接可以看出，以LSC-R、LSC-K为代表的引入地标点算法与DEC、GraEn为代表的基于深度自编码算法相比，在MNIST、COIL100及CIFAR10这些维度相对较高的大规模数据集中，聚类性能处于劣势，但在LetterRec及CoverType维度较低的数据集上则处于优势，可见基于自编码器的算法对于高维数据集较为友好，地标点的择取有利于大规模数据集的聚类性能。

本文算法采用欧氏距离与Kendall Tau距离度量融合的方式，挖掘数据的真实分布，从上述结论1)可以看出算法性能会受到距离度量方式的影响。为深入分析，将本文算法稀疏表示矩阵构造部分，分别采用基于欧氏距离(SCE)与基于Kendall Tau距离(SCK)计算点对相似度，采用ACC和NMI两个指标与本文算法融合度量方案进行性能对比。表4显示3种算法在3个数据集上的算法性能。从表4可以看出，2种距离度量方式的融合一定程度上能更全面地反映数据之间的结构信息，从而获得更为理想的聚类性能。欧氏距离与Kendall Tau距离这两种度量方式的优劣，并不能准确地从数据中评判出来，但结合表2、3可以看出，基于本文算法的SCE和SCK算法在3个数据集上与其他几种算法相比，相对处于优势。特别是与SCAL_R及SCAL_K算法对应ACC和NMI对比可以得出，本文通过基于网页质量的PageRank算法择取地标点总体上要优于随机选择与k-means算法的方法，但也有低于这两种方法的情况出现，为深入分析，更改地标点个数 $p$ 从100~1 000，在3个数据集上，进行本文算法与SCAL-R、SCAL-K的ACC值的对比实验。图1、2分别显示了SCAL-R、SCAL-K及本文算法在选取的3个数据集上，地标点个数从100~1 000之间变化，所对应的ACC及NMI值。图中可直观看出，本文算法在3个数据集上的ACC、NMI普遍比LSC-K、SCAL-K要高，除去在MNIST数据集上，地标点个数较少时，本文算法与SCAL-R算法相比略处于劣势，这是由于本文构造相似度矩阵采用的度量融合方案，在地标点个数相对多时，能有效地挖掘数据的结构信息，地标点个数从300开始，本文算法便展现出优势。同时从折线走势可以看出，地标点个数的增加，3种算法的ACC、NMI值也随之提高，且本文算法呈现稳步上升趋势，这表明两个聚类性能指标受地标点个数的选取影响很大。总体来说，在不同地标点的个数选取下，本文算法展现出优于其他算法的聚类性能。

4 结束语

随着数据规模的增大，结构信息的复杂度提高，在聚类过程中往往会耗费大量时间，存在相似度矩阵存储开销大及矩阵分解复杂度高的问题，且聚类精度也会受到影响。为此，本文提出一种结合度量融合和地标表示的自编码谱聚类算法，通过引入节点相对质量概念作为地标点选择的依据，以地标点与其他样本点之间相似度构造图相似度矩阵，以降低存储开销。同时，融合欧氏距离与Kendall Tau距离作为相似度度量方式，充分挖掘数据底层结构信息，且以栈式自编码器代替拉普拉斯矩阵分解步骤，通过联合学习框架进一步提高聚类精度。实验表明在几种大规模数据集上本文算法具有较好的聚类性能，但由于本文算法聚类精度受地标点个数影响较大，且选取方式略微增加了算法复杂度，未来将致力于寻求更为有效的地标点选择方式，在保证聚类精度的同时进一步降低算法复杂度。