基于改进Trace Lasso范数和PALM算法的子空间聚类方法

Subspace Clustering Based on Improved Trace Lasso Norm and PALM Algorithm

摘要

子空间聚类是一类广泛应用的聚类方法，其中最关键的技术是表示矩阵的获取，为使表示矩阵更好地满足块对角结构，提出一种基于改进迹Lasso范数和PALM算法的子空间聚类方法。首先，将原始数据减去噪声所得干净数据作为数据自表示的字典，能够促使表示矩阵更接近块对角结构；其次，提出一种改进的迹Lasso范数，利用非凸FCP范数约束矩阵的奇异值向量，能更好促使矩阵满足低秩性；最后，由于提出模型的非凸非光滑性及约束条件的非线性，利用近端交替线性极小化算法求解模型，具有收敛性保证。在CFP人脸数据集和动物面部图像数据集上进行聚类的数值实验表明，提出的子空间聚类方法相比于普遍应用的K⁃means聚类、谱聚类及稀疏子空间聚类有更好的聚类性能。

Abstract

Subspace clustering is a widely used clustering method， in which the most critical technology is representation matrix acquisition. To make the representation matrix better fit the block-diagonal structure， this study proposes a subspace clustering method based on an improved trace Lasso norm and the proximal alternating linearized minimization （PALM） algorithm. Firstly， the clean data obtained by subtracting the noise from the original data is used as the dictionary of data self-representation， which makes the representation matrix closer to the block-diagonal structure. Secondly， an improved trace Lasso norm is proposed. It utilizes a non-convex FCP norm to constrain the singular value vector of the matrix， so that the matrix can be better promoted to satisfy the low rank. Finally， due to the non-convexity and non-smoothness of the proposed model and the nonlinearity of the constraint conditions， the PALM algorithm is used to solve the model， which ensures convergence. Numerical experiments of clustering on the CFP face dataset and an animal face image dataset show that the proposed subspace clustering method outperforms the commonly used K-means clustering， spectral clustering， and sparse subspace clustering （SSC）.

Graphical abstract

关键词

子空间聚类 / 迹Lasso范数 / 非凸FCP范数 / 近端交替线性最小化

Key words

subspace clustering / trace Lasso norm / non-convex FCP norm / proximal alternating linearized minimization

聚类分析在机器学习^［1］、计算机视觉^［2］、图像处理^［3⁃4］等领域有广泛的应用，其作为一种无监督学习方法，能够在没有预先给定标签的情况下，根据数据点之间的相似性将其归入不同的类别，目前常用的聚类算法有K-means（K均值）聚类^［5］、谱聚类^［6］、层次聚类^［7］、子空间聚类^［8］等。如今数据的高维性使得聚类问题面临着极大挑战，子空间聚类方法针对高维数据提出。现有的聚类方法大致可分为四种：基于统计学的方法、矩阵分解方法、基于代数学方法及基于稀疏化方法，如稀疏子空间聚类^［9］（sparse subspace clustering，SSC）和低秩表示^［10］（low rank representation，LRR），本文主要研究子空间聚类问题。子空间聚类方法分为两个步骤^［11］：获取邻接矩阵和谱聚类，获取邻接矩阵可以转化成求解一个最小化问题；基于邻接矩阵，利用谱聚类算法得到最终的聚类结果。假设数据点{\mathit{x}}_i\in S_i(i\in \mathrm{1,2},\dots ,c)可以表示成其他所有数据的线性组合，对应矩阵形式为

\mathit{X}=\mathit{X}\mathit{Z}， s.t.{\mathit{Z}}_{jj}=\mathrm{0}(j=\mathrm{1,2},\cdots ,n)，（1）

上式中\mathit{X}为数据字典，\mathit{Z}为表示系数矩阵，结构为

\mathit{Z}=\left[\begin{array}{cccc}{\mathit{Z}}^{\mathrm{1}}& \mathrm{0}& \cdots & \mathrm{0}\\ \mathrm{0}& {\mathit{Z}}^{\mathrm{2}}& \cdots & \mathrm{0}\\ ⋮& ⋮& \ddots & ⋮\\ \mathrm{0}& \mathrm{0}& \cdots & {\mathit{Z}}^c\end{array}\right]，（2）

其中{\mathit{Z}}^i\in R^{n_i\times n_i}表示第i个子空间相应的表示系数。

SSC^［12］与LRR^［13］是得到广泛研究与应用的两类子空间聚类方法，SSC通过{\mathrm{𝓁}}_{\mathrm{1}}范数约束系数正则项。LRR方法^{［14⁃15］}通常由核范数代替秩函数约束系数正则项。文献［16］提出潜在低秩表示方法，在低秩约束的前提下利用潜在数据构造字典，通过求解核范数最小化问题近似恢复潜在数据的影响。Zhang等^［17］提出相关结构低秩子空间聚类方法，在模型中加入结构范数，将表示系数矩阵学习与谱聚类结合为一步。2016年Wu等^［18］利用Trace Lasso范数{‖\mathit{X}Diag({\mathit{z}}_j)‖}_{\mathrm{*}}约束系数矩阵\mathit{Z}促使\mathit{X}Diag({\mathit{z}}_j)低秩，给出以下模型：

\underset{z_j}{\mathrm{m}\mathrm{i}\mathrm{n}}{\displaystyle \sum _{j=\mathrm{1}}^n}{‖x_j-\mathit{X}{z}_j‖}_{\mathrm{1}}+\lambda {\displaystyle \sum _{j=\mathrm{1}}^n}{‖\mathit{X}Diag(z_j)‖}_{\mathrm{*}},s.t.x_j=\mathit{X}{z}_j,j=\mathrm{1,2},\cdots ,n，（3）

迹Lasso范数用核范数约束正则项。有文献表明^{［19⁃21］}，极小化非凸稀疏范数，譬如范数，非凸FCP范数（folded concave penalty，FCP）等，比极小化{\mathrm{𝓁}}_{\mathrm{1}}范数能更好地使向量满足稀疏性。文献［22］在数据无噪声且子空间独立的情况下，正则项满足强制块对角（enforced block diagonal，EBD）条件，最优表示矩阵满足块对角结构。{\mathrm{𝓁}}_{\mathrm{1}}范数约束系数矩阵缺乏考虑数据的全局结构。核范数极小化的目的是通过使矩阵的奇异值向量的{\mathrm{𝓁}}_{\mathrm{1}}范数达到极小从而使奇异值向量稀疏，促使表示矩阵满足低秩性。然而，{\mathrm{𝓁}}_{\mathrm{1}}范数不能很好地使向量稀疏，因此也不能很好地使矩阵满足低秩性；且实际问题中的数据不可避免地含有噪声，将原始噪声数据作为字典，不能有效促使表示矩阵满足块对角结构。

本文提出一种基于改进Trace Lasso范数和PALM算法的子空间聚类方法，提出的模型约束条件满足非线性，文献［23］中给出了一种求解一类非线性矩阵等式约束的非凸非光滑最优化问题的近端交替线性极小化算法（proximal alternating linearized minimization，PALM），适用本文模型。

1 提出的模型及求解算法

1.1 提出的模型

文献［22］指出在理想条件下，即数据是无噪声的，且子空间独立，正则项满足EBD条件，所求的最优表示矩阵满足块对角结构，这一特性确定了对数据的聚类（每个块对应每一类）。设x_j=\mathit{X}{z}_j且z_j=\left[z_{\mathrm{1}j}^{\mathrm{1}},\cdots ,z_{n_{\mathrm{1}}j}^{\mathrm{1}},\cdots ,z_{\mathrm{1}j}^c,\cdots ,z_{n_{k}j}^c\right]，假设x_j\in S_i，x_l^i表示子空间S_i中第l个数据，其中{\mathit{z}}_{\mathit{l}\mathit{j}}^{\mathit{i}}为表达式{\mathit{z}}_j中数据点x_l^i前的表示系数，则有

\mathit{X}Diag(z_j)=\left[x_{\mathrm{1}}^{\mathrm{1}},\cdots ,x_{n_{\mathrm{1}}}^{\mathrm{1}},\cdots ,x_{\mathrm{1}}^c,\cdots ,x_{n_k}^c\right]\left[\begin{array}{ccccccc}{z}_{\mathrm{1}j}^{\mathrm{1}}& & & & & & \\ & \ddots & & & & & \\ & & z_{n_{\mathrm{1}}j}^{\mathrm{1}}& & & & \\ & & & \ddots & & & \\ & & & & z_{\mathrm{1}j}^c& & \\ & & & & & \ddots & \\ & & & & & & z_{n_{k}j}^c\end{array}\right]=[\mathrm{0},\cdots ,\mathrm{0};\cdots ;z_{\mathrm{1}j}^i{x}_{\mathrm{1}}^i,\cdots ,z_{n_{i}j}^i{x}_{n_i}^i;\mathrm{0},\cdots ,\mathrm{0}]，

进一步，由于第i个子空间的维数（基底的个数）远远小于其所含的数据点个数，即\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}(S_i)\ll n_i，故表示系数{\left\{z_{lj}^i\right\}}_{l=\mathrm{1}}^{n_i}中非基底数据点前的系数也均为0，因此\mathit{X}Diag({\mathit{z}}_j)具有很高的低秩性。本文利用FCP范数极小化使矩阵\mathit{X}Diag({\mathit{z}}_j)的奇异值向量满足稀疏性，从而达到使\mathit{X}Diag({\mathit{z}}_j)低秩的目的，秩函数用FCP范数约束的正则项为：{‖\mathit{X}Diag({\mathit{Z}}_j)‖}_{P_{\lambda ,a}}，其中，矩阵的FCP范数定义为

{‖I‖}_{P_{\lambda ,a}}={\displaystyle \sum _{k=\mathrm{1}}^c}{P}_{\lambda ,a}\left({\sigma }_k(I)\right)，（4）

{\sigma }_k(I)是矩阵I\in R^{c\times n}(c\le n)的第k个奇异值，用P_{\lambda ,a}^{\mathrm{M}\mathrm{C}\mathrm{P}}(·)或P_{\lambda ,a}^{\mathrm{S}\mathrm{C}\mathrm{A}\mathrm{D}}(·)作为函数P_{\lambda ,a}(·)，MCP的形式为

（5）

SCAD的形式为

（6）

式（5）和（6）中的\lambda和a是参数，t是函数变量。需要指出的是，虽然上述对\mathit{X}Diag({\mathit{z}}_j)的低秩性分析是在对数据矩阵排序的前提下给出的，而在实际问题中，数据矩阵未排序，但排序操作可以通过矩阵乘积实现，设排序矩阵为\mathit{P}，由于\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}\left(\left(\mathit{P}\mathit{X}\right)Diag({\mathit{z}}_j)\right)=\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}\left(\mathit{X}Diag({\mathit{z}}_j)\right)，即对矩阵排序操作不影响矩阵的秩，因此\mathit{X}Diag({\mathit{z}}_j)的低秩性仍成立。

另一方面，文献［12］表明，只有用干净数据作为自表示字典时，才能保证表示矩阵\mathit{Z}更好地满足块对角结构。而实际中，数据矩阵\mathit{X}不可避免含有噪声，用含有噪声的数据矩阵\mathit{X}作为字典，不利于\mathit{Z}块对角结构的形成。设\mathit{E}为数据矩阵中所含的噪声，本文提出利用干净数据\mathit{X}-\mathit{E}作为自表示字典，给出约束条件\mathit{X}-\mathit{E}=(\mathit{X}-\mathit{E})\mathit{Z}。针对不同的噪声情况，对\mathit{E}所使用的范数也不同，如表1。可以看出，对噪声用稀疏范数，可以同时去除奇异样本或野点及随机噪声。因此，由于非凸稀疏范数能更好地使向量满足稀疏性，对噪声也使用非凸FCP范数{‖\mathit{E}‖}_{P_{\lambda ,a}}正则化。

综上，提出的子空间聚类模型为

\underset{Z_j,E}{\mathrm{m}\mathrm{i}\mathrm{n}}{\displaystyle \sum _{j=\mathrm{1}}^n}{‖\mathit{X}Diag({\mathit{Z}}_j)‖}_{P_{\lambda ,a}}+\lambda {‖\mathit{E}‖}_{P_{\lambda ,a}}s.t.\mathit{X}=(\mathit{X}-\mathit{E})\mathit{Z}+\mathit{E}\mathrm{。}（7）

为说明模型中表示矩阵\mathit{Z}的块对角结构，首先，给出以下定理。

定理 考虑k个独立子空间{\left\{S_j\right\}}_{j=\mathrm{1}}^k中的数据集，令\mathit{X}=[{\mathit{X}}_{\mathrm{1}},{\mathit{X}}_{\mathrm{2}},\cdots ,{\mathit{X}}_k]，\mathit{E}为噪声且\mathit{E}=[{\mathit{E}}_{\mathrm{1}},{\mathit{E}}_{\mathrm{2}},\cdots ,{\mathit{E}}_k]，数据集{\mathit{X}}_j-{\mathit{E}}_j\in S_j。对于满足\mathit{X}=(\mathit{X}-\mathit{E})\mathit{Z}+\mathit{E}的可行解{\mathit{Z}}^{\mathrm{*}}，将{\mathit{Z}}^{\mathrm{*}}分解为两部分，即{\mathit{Z}}^{\mathrm{*}}={\mathit{Z}}^B+{\mathit{Z}}^C，其中

{\mathit{Z}}^B=\left[\begin{array}{cccc}{\mathit{Z}}_{\mathrm{1}}^B& \mathrm{0}& \cdots & \mathrm{0}\\ \mathrm{0}& {\mathit{Z}}_{\mathrm{2}}^B& \cdots & \mathrm{0}\\ ⋮& ⋮& \ddots & ⋮\\ \mathrm{0}& \mathrm{0}& \cdots & {\mathit{Z}}_k^B\end{array}\right],{\mathit{Z}}^C=\left[\begin{array}{cccc}\mathrm{0}& \mathrm{*}& \cdots & \mathrm{*}\\ \mathrm{*}& \mathrm{0}& \cdots & \mathrm{*}\\ ⋮& ⋮& \ddots & ⋮\\ \mathrm{*}& \mathrm{*}& \cdots & \mathrm{0}\end{array}\right],

{\mathit{Z}}_j^B对应于子空间S_j，则\mathit{X}=(\mathit{X}-\mathit{E}){\mathit{Z}}^B+\mathit{E}成立或者{\mathit{X}}_j=(\mathit{X}-\mathit{E}){\mathit{Z}}_j^B+\mathit{E}，j=\mathrm{1,2},\cdots ,k，且(\mathit{X}-\mathit{E}){\mathit{Z}}^C=\mathrm{0}。

证明 对于可行解{\mathit{Z}}^{\mathrm{*}}，满足\mathit{X}=(\mathit{X}-\mathit{E}){\mathit{Z}}^{\mathrm{*}}+\mathit{E}，子空间S_j中的数据\mathit{X}-\mathit{E}=[{\mathit{X}}_{\mathrm{1}},{\mathit{X}}_{\mathrm{2}},\cdots ,{\mathit{X}}_k]=(\mathit{X}-\mathit{E})\left(\begin{array}{ccc}{Z}_{\mathrm{1}}^{\mathrm{*}}& & \\ & \ddots & \\ & & Z_k^{\mathrm{*}}\end{array}\right)，那么就有{\mathit{X}}_j-{\mathit{E}}_j=(\mathit{X}-\mathit{E}){\mathit{Z}}_j^{\mathrm{*}}\in S_j。因为{\mathit{Z}}^{\mathrm{*}}={\mathit{Z}}^B+{\mathit{Z}}^C，则{\mathit{X}}_j-{\mathit{E}}_j=(\mathit{X}-\mathit{E})({\mathit{Z}}_j^B+{\mathit{Z}}_j^C)。只有与{\mathit{X}}_{\mathit{j}}-{\mathit{E}}_{\mathit{j}}在同一子空间中数据前的表示系数不为0，即(\mathit{X}-\mathit{E}){\mathit{Z}}_j^B\in S_j，则(\mathit{X}-\mathit{E}){\mathit{Z}}_j^C中与{\mathit{X}}_{\mathit{j}}-{\mathit{E}}_{\mathit{j}}不在同一子空间中数据前的表示系数不为0，因此(\mathit{X}-\mathit{E}){{\mathit{Z}}_j}^C\in \oplus i\ne j{S}_j（\oplus是直和符号，这里表示i\ne j的所有子空间）。另一方面{\mathit{Z}}^C满足：(\mathit{X}-\mathit{E}){\mathit{Z}}_j^C=({\mathit{X}}_j-{\mathit{E}}_j)-(\mathit{X}-\mathit{E}){\mathit{Z}}_j^B\in S_j，则(\mathit{X}-\mathit{E}){{\mathit{Z}}_j}^C\in S_j\bigcap \oplus i\ne j{S}_j，这里假设子空间是独立的，那么S_j\bigcap \oplus i\ne j{S}_j=\{\mathrm{0}\}。考虑上面的过程，有(\mathit{X}-\mathit{E}){{\mathit{Z}}_j}^C=\{\mathrm{0}\}，当j取\mathrm{1,2},\cdots ,k时，满足(\mathit{X}-\mathit{E}){\mathit{Z}}^C=\mathrm{0}，因此(\mathit{X}-\mathit{E}){\mathit{Z}}^B=(\mathit{X}-\mathit{E})-(\mathit{X}-\mathit{E}){\mathit{Z}}^C=\mathit{X}-\mathit{E}成立，即\mathit{X}=(\mathit{X}-\mathit{E}){\mathit{Z}}^B+\mathit{E}。

定理给出了在独立子空间假设下表示矩阵{\mathit{Z}}^{\mathrm{*}}的性质，结果表明只有来自同一子空间S_j中的数据点{\mathit{X}}_j有实际贡献，即\mathit{X}=(\mathit{X}-\mathit{E}){\mathit{Z}}^B+\mathit{E}。表示矩阵的分解{\mathit{Z}}^{\mathrm{*}}={\mathit{Z}}^B+{\mathit{Z}}^C在{\mathit{Z}}^C\ne \mathrm{0}时块对角结构是未知的，同一子空间中数据点的贡献并没有被表示矩阵{\mathit{Z}}^{\mathrm{*}}明确地反映出来。这种情况下，问题（7）约束条件中的{\mathit{Z}}^{\mathrm{*}}不一定服从块对角性质。为了解决这个问题，文献［18］指出，对\mathit{X}Diag({\mathit{Z}}_j)利用低秩约束，从而有效促使{\mathit{Z}}^C=\mathrm{0}，使{\mathit{Z}}^{\mathrm{*}}满足块对角结构。本文利用非凸FCP范数对奇异值向量进行稀疏正则化约束，比利用Trace Lasso范数更利于促使\mathit{X}Diag({\mathit{Z}}_j)低秩。

求解上述模型（7）的困难在于目标函数的非凸非光滑性以及约束条件的非线性，本文提出的模型符合［23］中给出的一般模型的形式，且提出模型的目标函数满足文献［23］中收敛性定理的条件，从而具有收敛性保证的属性。

1.2 求解算法

令\mathit{X}Diag({\mathit{Z}}_j)={\mathit{T}}_j，两边同时乘以({\mathit{X}}^{\mathrm{T}}{\mathit{X})}^{-\mathrm{1}}{\mathit{X}}^{\mathrm{T}}，有({\mathit{X}}^{\mathrm{T}}{\mathit{X})}^{-\mathrm{1}}{\mathit{X}}^{\mathrm{T}}\mathit{X}Diag({\mathit{Z}}_j)=({\mathit{X}}^{\mathrm{T}}{\mathit{X})}^{-\mathrm{1}}{\mathit{X}}^{\mathrm{T}}{\mathit{T}}_j，从而Diag({\mathit{Z}}_j)=({\mathit{X}}^{\mathrm{T}}{\mathit{X})}^{-\mathrm{1}}{\mathit{X}}^{\mathrm{T}}{\mathit{T}}_j，进一步有{\mathit{Z}}_j=({\mathit{X}}^{\mathrm{T}}{\mathit{X})}^{-\mathrm{1}}{\mathit{X}}^{\mathrm{T}}{\mathit{T}}_j\mathit{e}，则问题（7）的形式等价为

\underset{T_j,E}{\mathrm{m}\mathrm{i}\mathrm{n}}{\displaystyle \sum _{j=\mathrm{1}}^n}{‖{\mathit{T}}_j‖}_{P_{\lambda ,a}}+\lambda {‖\mathit{E}‖}_{P_{\lambda ,a}}s.t.{\mathit{X}}_j=(\mathit{X}-\mathit{E})({\mathit{X}}^{\mathrm{T}}{\mathit{X})}^{-\mathrm{1}}{\mathit{X}}^{\mathrm{T}}{\mathit{T}}_j\mathit{e}+{\mathit{E}}_j,j=\mathrm{1,2},\cdots ,n\mathrm{。}（8）

其中{\mathit{X}}_j是第j个数据，即输入的第j个图像数据，\mathit{e}是元素全为1的向量{[\mathrm{1,1},\cdots ,\mathrm{1}]}^{\mathrm{T}}，{\mathit{E}}_j是第j个噪声图像数据。首先，用增广拉格朗日算法将问题转化为无约束问题

\underset{T_j,E}{\mathrm{m}\mathrm{i}\mathrm{n}}L({\mathit{T}}_j,\mathit{E})={\displaystyle \sum _{j=\mathrm{1}}^n}{‖{\mathit{T}}_j‖}_{P_{\lambda ,a}}+\lambda {‖\mathit{E}‖}_{P_{\lambda ,a}}+tr\left({\mathit{A}}^{\mathrm{T}}\left[(\mathit{X}-\mathit{E})({\mathit{X}}^{\mathrm{T}}{\mathit{X})}^{-\mathrm{1}}{\mathit{X}}^{\mathrm{T}}{\mathit{T}}_j\mathit{e}+{\mathit{E}}_j-{\mathit{X}}_j\right]\right)+

\frac{\mu }{\mathrm{2}}{\displaystyle \sum _{j=\mathrm{1}}^n}‖(\mathit{X}-\mathit{E})({\mathit{X}}^{\mathrm{T}}{\mathit{X})}^{-\mathrm{1}}{\mathit{X}}^{\mathrm{T}}{\mathit{T}}_j\mathit{e}+{\mathit{E}}_j-{\mathit{X}}_j‖,（9）

其中\mathit{A}\in R^{m\times \mathrm{1}}是拉格朗日乘数，\mu是大于0的参数。优化问题（8）和拉格朗日函数（9）是等价问题。令

H({\mathit{T}}_j,\mathit{E})={\displaystyle \sum _{j=\mathrm{1}}^n}tr\left({\mathit{A}}^{\mathrm{T}}\left[(\mathit{X}-\mathit{E})({\mathit{X}}^{\mathrm{T}}{\mathit{X})}^{-\mathrm{1}}{\mathit{X}}^{\mathrm{T}}{\mathit{T}}_j\mathit{e}+{\mathit{E}}_j-{\mathit{X}}_j\right]\right)+\frac{\mu }{\mathrm{2}}{{\displaystyle \sum _{j=\mathrm{1}}^n}‖\left(\mathit{X}-\mathit{E}\right){\left({\mathit{X}}^{\mathrm{T}}\mathit{X}\right)}^{-\mathrm{1}}{\mathit{X}}^{\mathrm{T}}{\mathit{T}}_j\mathit{e}+{\mathit{E}}_j-{\mathit{X}}_j‖}_{\mathrm{2}}^{\mathrm{2}}\mathrm{。}（10）

利用PALM算法求解（9），k=\mathrm{0,1},\mathrm{2},\cdots，n。

（1）固定\mathit{E}={\mathit{E}}^k，H\left({\mathit{T}}_{\mathit{j}},{\mathit{E}}^k\right)可线性化为

H\left({\mathit{T}}_j,{\mathit{E}}^k\right)\approx H\left({{\mathit{T}}_j}^k,{\mathit{E}}^k\right)+tr\left[{\nabla }_{T_j}H{\left({{\mathit{T}}_j}^k,{\mathit{E}}^k\right)}^{\mathrm{T}}\left({\mathit{T}}_j-{{\mathit{T}}_j}^k\right)\right]+\frac{c_k}{\mathrm{2}}{‖{\mathit{T}}_j-{{\mathit{T}}_j}^{\mathit{k}}‖}_F^{\mathrm{2}}=\frac{c_k}{\mathrm{2}}{‖{\mathit{T}}_j-\left({{\mathit{T}}_j}^k-\frac{{\nabla }_{T_j}H\left({{\mathit{T}}_j}^k,{\mathit{E}}^k\right)}{c_k}\right)‖}_F^{\mathrm{2}},

其中{\nabla }_{T_j}H({{\mathit{T}}_j}^k,{\mathit{E}}^k)=-\mu {\mathit{P}}^{\mathrm{T}}({\mathit{X}}^k-{\mathit{E}}^k{)}^{\mathrm{T}}[{{\mathit{X}}_j}^k-({\mathit{X}}^k-{\mathit{E}}^k)\mathit{P}{{\mathit{T}}_j}^k\mathit{e}-{{\mathit{E}}_j}^k+\frac{{\mathit{A}}^k}{\mu }]{\mathit{e}}^{\mathrm{T}}。

求解以下极小化问题

\begin{array}{l}{{\mathit{T}}_j}^{k+\mathrm{1}}=\underset{{\mathit{T}}_j}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}\left\{{‖{\mathit{T}}_j‖}_{P_{\lambda ,a}}+\frac{c_k}{\mathrm{2}}{‖{\mathit{T}}_j-\left({{\mathit{T}}_j}^k-\frac{{\nabla }_{T_j}H({{\mathit{T}}_{\mathit{j}}}^k,{\mathit{E}}^k)}{c_k}\right)‖}_F^{\mathrm{2}}\right\}=\\ pro{x}_{\frac{\mathrm{1}}{c_k}}^{{‖.‖}_{P_{\lambda ,a}}}\left({{\mathit{T}}_j}^k-\frac{{\nabla }_{{\mathit{T}}_j}H({{\mathit{T}}_j}^k,{\mathit{E}}^k)}{c_k}\right),j=\mathrm{1,2},\cdots ,n\mathrm{。}\end{array}（11）

（2）固定{\mathit{T}}_j={{\mathit{T}}_j}^{k+\mathrm{1}}，H({{\mathit{T}}_j}^{k+\mathrm{1}},\mathit{E})可线性化为

H({{\mathit{T}}_j}^{k+\mathrm{1}},\mathit{E})\approx H({{\mathit{T}}_j}^k,{\mathit{E}}^k)+tr\left[{\nabla }_{\mathit{E}}H({{\mathit{T}}_j}^k,{\mathit{E}}^k{)}^{\mathrm{T}}(\mathit{E}-{\mathit{E}}^k)\right]+\frac{d_k}{\mathrm{2}}{‖\mathit{E}-{\mathit{E}}^k‖}_F^{\mathrm{2}}=\frac{d_k}{\mathrm{2}}‖\mathit{E}-\left({\mathit{E}}^k-{\left.\frac{{\nabla }_{\mathit{E}}H({{\mathit{T}}_j}^k,{\mathit{E}}^k)}{d_k}\right)‖}_F^{\mathrm{2}}\right.，

其中{\nabla }_{E}H({{\mathit{T}}_j}^{k+\mathrm{1}},{\mathit{E}}^k)=\mu \mathit{I}-\mu {\mathit{P}}^T({{\mathit{T}}_j}^{k+\mathrm{1}}{)}^{\mathrm{T}}({\mathit{X}}^k-{\mathit{E}}^k{)}^{\mathrm{T}}\left[{{\mathit{X}}_j}^k-({\mathit{X}}^k-{\mathit{E}}^k)\mathit{P}{{\mathit{T}}_j}^{k+\mathrm{1}}\mathit{e}-\right.\left.{{\mathit{E}}_j}^k+\frac{{\mathit{A}}^k}{\mu }\right]{\mathit{e}}^{\mathrm{T}}。

求解以下极小化问题

\begin{array}{l}{\mathit{E}}^{k+\mathrm{1}}=\underset{\mathit{E}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}\left\{\lambda {‖\mathit{E}‖}_{P_{\lambda ,a}}+\frac{d_k}{\mathrm{2}}{‖\mathit{E}-\left({\mathit{E}}^k-\frac{{\nabla }_{\mathit{E}}H({{\mathit{T}}_j}^{k+\mathrm{1}},{\mathit{E}}^k)}{d_k}\right)‖}_F^{\mathrm{2}}\right\}=\\ pro{x}_{\frac{\mathrm{1}}{d_k}}^{\lambda {‖\mathit{E}‖}_{P_{\lambda ,a}}}\left({\mathit{E}}^k-\frac{{\nabla }_{\mathit{E}}H({{\mathit{T}}_j}^{k+\mathrm{1}},{\mathit{E}}^k)}{d_k}\right),j=\mathrm{1,2},\cdots ,n\mathrm{。}\end{array}（12）

利用文献［21］给出的算法可求解（11）和（12），其中\mathit{P}=({\mathit{X}}^{\mathrm{T}}{\mathit{X})}^{-\mathrm{1}}{\mathit{X}}^{\mathrm{T}}。综上，求解模型（8）的算法步骤如下：

（1）输入：观测数据\mathit{X}\in R^{m\times n}，参数\lambda ,\mu ,{\gamma }_{\mathrm{1}}>\mathrm{1},{\gamma }_{\mathrm{2}}>\mathrm{1}，初始化{{\mathit{T}}_j}^{\mathrm{0}}\in R^{m\times n}(j=\mathrm{1,2},\cdots ,n),{\mathit{E}}^{\mathrm{0}}\in R^{m\times n},{\mathit{A}}^{\mathrm{0}}\in R^{m\times \mathrm{1}}；

（2）设置步长c_k={\gamma }_{\mathrm{1}}{L}_{\mathrm{1}}({\mathit{E}}^k)={\gamma }_{\mathrm{1}}{‖\mu ({\mathit{P}}^{\mathrm{T}}{(\mathit{X}-\mathit{E})}^{\mathrm{T}}(\mathit{X}-\mathit{E})\mathit{P}‖}_F，通过公式（11）计算变量{\mathit{T}}_j；设置步长d_k={\gamma }_{\mathrm{2}}{L}_{\mathrm{2}}\left({{\mathit{T}}_j}^{k+\mathrm{1}}\right)={\gamma }_{\mathrm{2}}{‖{\mathit{P}}^{\mathrm{T}}{\mathit{T}}_j\mathit{e}{\mathit{e}}^{\mathrm{T}}‖}_F，通过公式（12）计算变量{\mathit{E}}_j；

（3）更新拉格朗日乘子：{\mathit{A}}^{k+\mathrm{1}}={\mathit{A}}^k+\mu \left[(\mathit{X}-{\mathit{E}}^{k+\mathrm{1}})({\mathit{X}}^{\mathrm{T}}{\mathit{X})}^{-\mathrm{1}}{\mathit{X}}^{\mathrm{T}}{{\mathit{T}}_j}^{k+\mathrm{1}}\mathit{e}+{{\mathit{E}}_j}^{k+\mathrm{1}}-{\mathit{X}}_j\right]；

（4）输出：{{\mathit{Z}}_j}^{\mathrm{*}}=({\mathit{X}}^{\mathrm{T}}{\mathit{X})}^{-\mathrm{1}}{\mathit{X}}^{\mathrm{T}}{{\mathit{T}}_j}^{\mathrm{*}}\mathit{e}(j=\mathrm{1,2},\cdots ,n)。

2 数值实验

利用所求的表示系数矩阵构造相似度矩阵\mathit{W}=\frac{{\mathit{Z}}^{\mathrm{*}}+({\mathit{Z}}^{\mathrm{*}}{)}^{\mathrm{T}}}{\mathrm{2}}，最后利用谱聚类算法得到最终的聚类结果。

2.1 评价指标

将提出的子空间聚类算法和K-means、谱聚类等聚类方法进行对比分析，使用轮廓系数作为评价指标，轮廓系数是衡量类内稠密、类间稀疏程度的指标。轮廓系数的公式为

S(i)=\frac{b(i)-a(i)}{\mathrm{m}\mathrm{a}\mathrm{x}\{a(i),b(i)\}},（13）

其中b(i)是类间距离，a(i)是类内距离。轮廓系数的取值范围为[-\mathrm{1,1}]，轮廓系数越大，聚类效果越好。

2.2 数据集

本文选取可公开访问的两种数据集：CFP人脸数据集和动物面部数据集，机器学习的学习性能和泛化能力在很大程度上取决于数据集的特征。本文对CFP人脸数据集选取10类数据，每类数据包含同一个人的10张正脸和4张侧脸，如图1所示；动物面部数据集选取4类不同表情动物共250张图像进行学习，如图2所示。

2.3 算法收敛性

图3为提出模型的目标函数随迭代次数的变化曲线图，其中横轴表示迭代次数，纵轴表示目标函数值，由图3可知，目标函数值随迭代次数的增加而减少并逐渐趋于稳定。图3（a）是在CFP人脸数据集的实验结果，当迭代次数达到48时，函数值逐渐趋于一个稳定的值。图3（b）是在动物图像数据集的实验结果，当迭代次数达到43时，函数值逐渐趋于一个稳定的值。

2.4 聚类结果

对聚类结果进行可视化处理，为了更直观地展示聚类效果，利用t分布⁃随机近邻嵌入（t-SNE，t-distributed stochastic neighbor embedding）降维方法将输入的高维数据投影到二维平面上，根据其预测标签用不同的颜色标记相应数据点的聚类结果，对聚类结果可视化。然而聚类好坏难以单一地用肉眼直接观察，因此需要对聚类效果进行量化，本文对聚类结果使用轮廓系数进行量化度量，确定最佳的聚类簇数并对聚类效果进行评估，如表2所示。由轮廓系数得到的CFP人脸数据集最佳聚类簇数为10，轮廓系数为0.673；动物面部数据集最佳聚类簇数为4，轮廓系数为0.739。

2.4.1 可视化结果

图4和图5中的K-means和谱聚类是直接对两组数据集使用t-SNE降维后的数据做聚类，SSC和本文模型是根据降维的数据矩阵求出稀疏表示矩阵和低秩表示矩阵后，用所求的表示矩阵做邻接矩阵进一步使用谱聚类算法进行聚类。

2.4.2 聚类结果定量评估

对K-means聚类、谱聚类、SSC以及本文提出的聚类模型进行聚类精度比较结果见表3。根据表中的数值结果可得到以下结论：本文所提出的模型在两个数据集上的实验结果均具有较好的聚类性能；针对CFP人脸数据集轮廓系数为0.673，均高于其他三类聚类算法；针对动物图像数据集，得到的轮廓系数为0.739，均高于其他三种聚类算法。

3 结论

本文提出了基于改进Trace Lasso范数和PALM算法的子空间聚类模型，主要创新点和贡献有：（1）针对约束系数矩阵低秩，本文将Trace Lasso范数和非凸FCP范数结合，对\mathit{X}Diag({\mathit{z}}_{\mathit{j}})的奇异值向量利用非凸稀疏范数约束，能更好地使\mathit{X}Diag({\mathit{Z}}_{\mathit{j}})低秩，进而促使系数矩阵\mathit{Z}低秩；（2）针对噪声数据，对表示字典引入噪声项，能有效促使表示矩阵满足块对角结构。本模型考虑到输入数据包含噪声的随机性和稀疏性及表示矩阵的低秩性，用干净数据作为字典表示并对正则项的低秩约束项作了改进，选用ALM算法和PALM算法求解，PALM针对非凸非光滑组合优化问题，适用于提出的模型并有收敛性保证。在数值实验部分给出了算法收敛性分析并作了聚类效果比较，结果表明，所提出的子空间聚类模型有更好的聚类准确率。然而模型的参数是通过手动调试的，在后续的研究中，将更关注于从数据中自适应地学习最优参数的研究。

参考文献

[1] CHEN Z D， LI D Z， WAN H X， et al. Unsupervised machine learning methods for polymer nanocomposites data via molecular dynamics simulation［J］.Molecular Simulation， 2020， 46 （18）： 1509-1521.

[2] BHATT M S， PATALIA T. Computer vision systems for content-based image classification［J］.Statistics， Optimization & Information Computing， 2019， 7 （4）： 840-853.

[3] PENG Z H， LIU H， JIA Y H， et al.Adaptive attribute and structure subspace clustering network［J］.IEEE Transactions on Image Processing， 2022， 31：3430-3439.

[4] YANG A Y， WRIGHT J， MA Y， et al.Unsupervised segmentation of natural images via lossy data compression［J］.Computer Vision and Image Understanding， 2008， 110（2）：212-225.

[5] 周蓉蓉，陈栋，刘思远.基于K均值聚类算法的生鲜运输路径优化模型［J］.农业大数据学报， 2022， 4（1）： 89-97.

[6] 王丙参，魏艳华，李旭.基于独立与距离相似度的谱聚类比较及数据模拟［J］.统计与决策， 2022，38（15）： 22-28.

[7] 邢涛，雍毅，侯江，等.基于主成分分析与层次聚类分析的水质综合评价［J］.四川环境， 2022，41（4）： 131-139.

[8] LI B H， LU H C， ZHANG Y， et al. Subspace clustering under complex noise［J］. IEEE Transactions on Circuits and Systems for Video Technology， 2019， 29（4）： 930-940.

[9] ELHAMIFAR E， VIDAL R. Sparse subspace clustering： algorithm， theory， and applications［J］. IEEE Transactions on Pattern Analysis and Machine Intelligence， 2013， 35（11）： 2765-2781.

[10] LIU G， LIN Z， YAN S， et al.Robust recovery of subspace structures by low-rank representation［J］.Pattern Analysis and Machine Intelligence， 2010， 35（1）：171-184.

[11] 王卫卫，李小平，冯象初，等.稀疏子空间聚类综述［J］.自动化学报， 2015，41（8）：1373-1384.

[12] ZOU X Y， XUE J Q， LI X T， et al. High-fidelity sEMG signals recorded by an on-skin electrode based on AgNWs for hand gesture classification using machine learning［J］. ACS Applied Materials & Interfaces， 2023， 15（15）： 19374-19383.

[13] LU G F， YU Q R， WANG Y，et al.Hyper-Laplacian regularized multi-view subspace clustering with low-rank tensor constraint［J］.Neural Networks， 2020， 125（6）：214-223.

[14] 于茜茜.稀疏低秩子空间聚类算法研究［D］.大连：大连海事大学，2020.

[15] 李欢，唐科威.基于低秩张量表示的多视图子空间聚类［J］.理论数学， 2023， 13（10）：2877-2887.

[16] ZHANG H Y， LIN Z C， ZHANG C，et al.Robust latent low rank representation for subspace clustering［J］.Neurocomputing， 2014， 145（18）：369-373.

[17] ZHANG X， SUN F， LIU G，et al.Fast low-rank subspace segmentation［J］.IEEE Transactions on Knowledge and Data Engineering，2014（5）：1293-1297.

[18] WU L， WANG Y， PAN S.Exploiting attribute correlations： A novel Trace Lasso-based weakly supervised dictionary learning method［J］.IEEE Transactions on Cybernetics， 2017，47（12）：4497-4508.

[19] XU Z B， CHANG X Y， XU F M，et al.L1/2 regularization： A thresholding representation theory and a fast solver［J］.IEEE Transactions on Neural Networks and Learning Systems， 2012， 23（7）：1013-1027.

[20] CHARTRAND R. Exact reconstruction of sparse signals via nonconvex minimization［J］.IEEE Signal Processing Letters， 2007， 14（10）：707-710.

[21] ZHANG W J， FENG X C， XIAO F， et al. A folded concave penalty regularized subspace clustering method to integrate affinity and clustering［J］. Mathematical Problems in Engineering， 2021（3）： 6641180.

[22] LU C Y， FENG J S， LIN Z C， et al. Subspace clustering by block diagonal representation［J］. IEEE Transactions on Pattern Analysis and Machine Intelligence， 2019， 41（2）： 487-501.

[23] ZHANG W J， FENG X C， XIAO F，et al.A class of bilinear matrix constraint optimization problem and its applications［J］.Knowledge-Based Systems， 2021， 231（8）：107429.

基金资助

国家自然科学基金资助项目“面向视觉注意选择的复杂场景层次化语义理解建模研究”(62171361)