基于稀疏性约束与指数图正则化的鲁棒非负矩阵分解

Robust NMF method based on sparse constraints and exponential graph regularization

摘要

在数据挖掘与机器学习领域，非负矩阵分解（Non-negative Matrix Factorization，NMF）作为一种高效的数据降维与特征表示方法受到广泛关注。但是，标准NMF只能处理非负数据，并且对噪声或异常值敏感，处理特征与样本失衡的小样本（Small Samples Size，SSS）问题时易出现模型过拟合、泛化性减弱及鲁棒性降低。为克服这些局限、获得更好的聚类性能，本文提出了一种基于稀疏性约束与指数图正则化的鲁棒非负矩阵分解算法。算法借助Semi-NMF处理混合符号数据，采用L_{\mathrm{2,1}}范数缓解噪声和异常值的影响，并引入指数图正则化以保留数据的分布特征与几何结构信息，结合稀疏性和正交性约束减少冗余信息、增强特征间独立性并防止过拟合。在9个公共数据集上的比较实验结果显示，算法的聚类准确性和鲁棒性优于其他8种经典聚类算法，从而验证了算法在小样本聚类任务中的优越性。

Abstract

In the field of data mining and machine learning， non-negative matrix factorization （NMF） has received extensive attention as an efficient method for data dimensionality reduction and feature representation.However， standard NMF is sensitive to noise or outliers and can handle only non-negative data， thus is prone to overfitting and reduced generalization and robustness when dealing with small sample size （SSS） problems where features and samples are imbalanced.To overcome these limitations and achieve better clustering performance， a robust NMF algorithm based on sparsity constraints and exponential graph regularization is proposed.The Semi-NMF is used to handle mixed-sign data， the L_{\mathrm{2,1}} norm is employed to mitigate the influence of noise and outliers， and the exponential graph regularization is introduced to preserve the distribution characteristics and geometric structure information of the data.Meanwhile， the combination of sparsity and orthogonality constraints is used to reduce the redundant information， enhance the independence of features and prevent overfitting.Comparitive experimental results on nine public datasets show that the clustering accuracy and robustness of the proposed algorithm are superior to those of the other 8 classic clustering algorithms， which verifys its effectiveness in small sample clustering tasks.

Graphical abstract

关键词

Semi-NMF / 鲁棒性 / 指数图正则化 / SSS问题 / 聚类

Key words

semi-NMF / robustness / exponential graph regularization / SSS problem / clustering

作为一种重要的数据分析方法，聚类被广泛应用于机器学习、图像分析、信号处理和生物信息学等领域^［1-3］。聚类分析的核心思想是：通过对样本特征间的相似性进行建模，将数据划分为若干簇，以便揭示其潜在结构与模式。常见聚类方法主要有k-means^［4］、密度峰值聚类^［5］、谱聚类^［6］以及基于神经网络的聚类^［7］等。

虽然聚类方法在不同类型数据分析任务中均能取得良好的效果，但随着数据规模和复杂度不断增加，其在高维或非球状数据中常常难以有效区分簇结构，并对初始条件和簇数设定较为敏感。同时，参数设定与分辨率选择也成为制约聚类方法性能的关键因素。因此，如何在高维、小样本和噪声干扰等复杂应用场景下保证聚类方法的鲁棒性与准确性已成为当前聚类研究的重要挑战。对此问题，研究者通常采用低维嵌入技术，通过将原始数据映射到低维空间来提升聚类算法的泛化能力与抗干扰性能。

非负矩阵分解（Non-negative Matrix Factorization， NMF）^［8-10］是低维表示的一种常用方法，主要通过引入非负性约束来挖掘原始数据的潜在信息，使分解结果具有可解释性，目前已被广泛应用于聚类分析中。传统的NMF采用Frobenius范数来最小化重构误差，易受噪声和异常值的干扰。Kong等^［11］提出了鲁棒非负矩阵分解（Robust Non-negative Matrix Factorization， RNMF），并采用L_{\mathrm{2,1}}范数来衡量重构误差、增强算法的鲁棒性。然而，该方法得到的低维矩阵往往仍然较为稠密。

研究表明，在分解过程中对低维矩阵施加稀疏性约束能够引导算法学习更具判别力的特征表示。Lee和Seung^［12］通过添加L_{\mathrm{1}}或L_{\mathrm{2}}范数来促进分解矩阵稀疏化。为优化系数矩阵稀疏性较弱的问题，Hoyer^［13］在NMF中引入显式稀疏正则项，利用L_{\mathrm{1}}与L_{\mathrm{2}}范数联合约束稀疏度，获得了更具判别性的低维表示。Gao等^［14］提出了基于L_{\mathrm{2}}范数的稀疏约束算法，实现对分解结果的可控稀疏性。杨亮东等^［15］针对新增数据增大引起的运算效率降低现象，提出了一种稀疏约束的增量式非负矩阵分解改进算法。杨国亮等^［16］提出一种基于L_{\mathrm{1}/\mathrm{2}}稀疏性和峰度平滑约束非负矩阵分解算法，成功解决了传统高光谱图像解混方法中存在的解混效率低等问题。值得注意的是，以上算法虽然在高维特征学习方面取得了一定成效，但因实际数据之间往往存在复杂关联，其性能依然存在局限性。

在低维表示过程中，将数据结构关系引入目标函数并进行优化能够更加准确地反映数据内在的分布规律与拓扑信息。Cai等^［17］提出了图正则化非负矩阵分解方法，通过引入图嵌入以保持局部流形结构，使NMF在处理具有复杂图结构的数据（如图像或社交网络）时依然能够兼顾几何一致性与聚类判别能力。Li等^［18］提出的图正则化稀疏非负矩阵分解方法在保留流形结构的同时强化了稀疏特征表达。Chen等^［19］提出的稀疏约束正交图正则化非负矩阵分解算法结合正交约束、图正则化与稀疏约束，在保持局部结构的同时提升了分解的判别性。高海燕等^［20］提出的双图正则化鲁棒非负矩阵分解算法引入了样本图与特征图正则化，结合L_{\mathrm{2,1}}范数增强抗噪性，并利用稀疏约束减少异常值干扰，使低维表示既保留了样本间的局部邻域关系，又保持特征间的内在相关性，显著提升了聚类精度与鲁棒性。这些新算法的提出进一步推动了NMF在复杂场景下的泛化性与鲁棒性。

尽管图正则化在一定程度上提升了聚类表现，但在特征维度高、样本量有限的情况下仍不可避免地面临小样本（Small Sample Size， SSS）问题^［21］。为解决此问题，研究者提出了基于矩阵指数的降维与聚类方法（如指数局部保持投影^［22］、指数弹性保持投影^［23］、指数局部判别嵌入^［24］及指数判别分析^［25］等），这些方法不仅在控制理论与马尔可夫过程分析中有重要应用，在机器学习领域也能有效缓解奇异性问题。

本文基于Semi-NMF框架，结合指数图正则化、稀疏性约束及正交性约束条件提出了一种基于稀疏约束与指数图正则化的鲁棒非负矩阵分解聚类算法（Robust NMF via Sparse Constraints and Exponential Graph Regularization， RNMFSCEG）。算法通过引入L_{\mathrm{2,1}}范数重构误差项显著提升了噪声及异常值的干扰，且所设计的指数图正则化项能够刻画数据流形的局部几何结构与全局分布特征，还通过施加稀疏性约束增强了特征子空间的可解释性。在Semi-NMF下，算法进一步拓展了对混合符号型数据的处理能力，同时在小样本SSS场景下也具有更好的适应性。

1 相关工作

1.1 NMF与Semi-NMF

NMF^［8］将数据矩阵\mathit{X}\in {\mathbb{R}}_{+}^{m\times n}分解为两个低秩非负矩阵\mathit{U}和\mathit{V}的乘积，即\mathit{X}\approx \mathit{U}\times {\mathit{V}}^{\mathrm{T}}。采用Frobenius范数度量重构误差，其目标函数为：

J_{\mathrm{N}\mathrm{M}\mathrm{F}}=\underset{\mathit{U},\mathit{V}\ge \mathrm{0}}{\mathrm{m}\mathrm{i}\mathrm{n}}{‖\mathit{X}-\mathit{U}{\mathit{V}}^{\mathrm{T}}‖}_{\mathrm{F}}^{\mathrm{2}}（1）

其中， {‖\cdot ‖}_{\mathrm{F}}为Frobenius范数，基矩阵\mathit{U}\in {\mathbb{R}}_{+}^{m\times k}表示特征空间的潜在结构，系数矩阵\mathit{V}\in {\mathbb{R}}_{+}^{n\times k}表示样本空间的潜在结构，k\ll \left\{m,n\right\}。根据乘法迭代更新规则，其更新公式为：

{\mathit{U}}_{ij}\leftarrow {\mathit{U}}_{ij}\frac{(\mathit{X}{\mathit{V}}^{\mathrm{T}}{)}_{ij}}{(\mathit{U}\mathit{V}{\mathit{V}}^{\mathrm{T}}{)}_{ij}},\begin{array}{c}\end{array}{\mathit{V}}_{\mathit{i}\mathit{j}}\leftarrow {\mathit{V}}_{\mathit{i}\mathit{j}}\frac{({\mathit{U}}^{\mathrm{T}}{\mathit{X})}_{ij}}{({\mathit{U}}^{\mathrm{T}}{\mathit{U}\mathit{V})}_{ij}}（2）

在实际应用中，常常会出现包含负值的混合数据，如心电图（Electrocardiogram，ECG）信号可能包含负向波峰（如Q波和S波），股票收益率矩阵可能包含负值（如做空收益）等。此时，标准NMF因强制所有矩阵元素非负而失效，需采用Semi-NMF^［26］或其变体。Semi-NMF放松了NMF的约束条件，仅要求系数矩阵\mathit{V}非负，基矩阵\mathit{U}可以包含负值，其目标函数为：

J_{\mathrm{S}\mathrm{e}\mathrm{m}\mathrm{i}‐\mathrm{N}\mathrm{M}\mathrm{F}}=\underset{\mathit{V}\ge \mathrm{0}}{\mathrm{m}\mathrm{i}\mathrm{n}}{‖\mathit{X}-\mathit{U}{\mathit{V}}^{\mathrm{T}}‖}_{\mathrm{F}}^{\mathrm{2}}（3）

由于系数矩阵\mathit{V}\in {\mathbb{R}}_{+}^{n\times k}需要保持非负性，利用

\mathit{Z}={\mathit{Z}}^{+}-{\mathit{Z}}^{-}，{\mathit{Z}}^{+}=\left(\left|\mathit{Z}\right|+\mathit{Z}\right)/\mathrm{2}，{\mathit{Z}}^{-}=\left(\left|\mathit{Z}\right|-\mathit{Z}\right)/\mathrm{2}，采用乘法更新规则优化\mathit{V}，有

{\mathit{V}}_{ik}\leftarrow {\mathit{V}}_{ik}\sqrt[]{\frac{{\left({\mathit{X}}^{\mathrm{T}}\mathit{U}\right)}_{ik}^{+}+{\left[\mathit{V}{\left({\mathit{U}}^{\mathrm{T}}\mathit{U}\right)}^{-}\right]}_{ik}}{{\left({\mathit{X}}^{\mathrm{T}}\mathit{U}\right)}_{ik}^{-}+{\left[\mathit{V}{\left({\mathit{U}}^{\mathrm{T}}\mathit{U}\right)}^{+}\right]}_{ik}}}（4）

进一步，通过最小二乘法可求解得\mathit{U}的更新：

\mathit{U}=\mathit{X}\mathit{V}{\left({\mathit{V}}^{\mathrm{T}}\mathit{V}\right)}^{-\mathrm{1}}（5）

1.2 RNMF

RNMF^［11］采用L_{\mathrm{2,1}}范数重构误差。通过将鲁棒损失函数替代为Frobenius范数，RNMF在含噪声或含异常样本的数据中仍能保持特征提取的准确性与聚类的可靠性。RNMF的目标函数为：

\underset{\mathit{U}\ge \mathrm{0},\mathit{V}\ge \mathrm{0}}{\mathrm{m}\mathrm{i}\mathrm{n}}{‖\mathit{X}-\mathit{U}{\mathit{V}}^{\mathrm{T}}‖}_{\mathrm{2,1}}={\displaystyle \sum _{j=\mathrm{1}}^n}{‖{\left(\mathit{X}-\mathit{U}{\mathit{V}}^{\mathrm{T}}\right)}_j‖}^{\mathrm{2}}（6）

更新规则为：

{\mathit{U}}_{ik}\leftarrow {\mathit{U}}_{ik}\frac{\mathit{X}\mathit{J}{\mathit{V}}_{ik}}{{\left(\mathit{U}{\mathit{V}}^{\mathrm{T}}\mathit{J}\mathit{V}\right)}_{ik}},{\mathit{V}}_{ij}\leftarrow {\mathit{V}}_{ik}\frac{{\left({\mathit{X}}^{\mathrm{T}}\mathit{U}\mathit{J}\right)}_{jk}}{{\left[\mathit{V}{\mathit{U}}^{\mathrm{T}}\mathit{U}\mathit{J}\right]}_{jk}}（7）

其中，\mathit{J}为对角矩阵，{\mathit{J}}_{ii}=\mathrm{1}/{‖{\left(\mathit{X}-\mathit{U}{\mathit{V}}^{\mathrm{T}}\right)}_i‖}_{\mathrm{2}}为其对角元素。

2 算法

传统NMF仅适用于处理非负数据，存在对噪声和异常值敏感的缺陷。同时，在面临特征与样本数量失衡的小样本场景时，容易出现模型过拟合问题，削弱模型的泛化能力和鲁棒性。为解决上述局限、提升聚类性能，本文提出的RNMFSCEG算法通过矩阵指数自适应地调节邻域权重，结合稀疏与正交约束增强低维表示的判别性与稳健性，以便显著提升高维小样本数据的聚类精度与鲁棒性。算法的目标函数为：

\underset{\mathit{V}\ge \mathrm{0}}{\mathrm{m}\mathrm{i}\mathrm{n}}{‖\mathit{X}-\mathit{U}{\mathit{V}}^{\mathrm{T}}‖}_{\mathrm{2,1}}+\frac{\lambda }{\mathrm{2}}\mathrm{T}\mathrm{r}\left({\mathit{V}}^{\mathrm{T}}\mathrm{e}\mathrm{x}\mathrm{p}\left(\mathit{L}\right)\mathit{V}\right)+\alpha {‖\mathit{V}‖}_{\mathrm{1}}+\frac{\beta }{\mathrm{2}}{‖{\mathit{V}}^{\mathrm{T}}\mathit{V}-\mathit{I}‖}_{\mathrm{F}}^{\mathrm{2}}（8）

其中，\mathrm{T}\mathrm{r}\left(\cdot \right)表示矩阵的迹，\alpha为稀疏误差项权重系数，\beta是控制稀疏性与重构精度之间的正则化系数，\lambda表示指数图正则化参数，\mathit{I}为单位矩阵，\mathit{L}=\mathit{D}-\mathit{W}为图拉普拉斯矩阵，\mathit{W}为邻接矩阵，\mathit{D}为度矩阵，{\mathit{D}}_{jj}=\sum _l{\mathit{W}}_{jl}。

鉴于图拉普拉斯矩阵\mathit{L}仅能捕捉样本之间的一阶邻接关系，而在高维小样本条件下有限的样本数量常常导致局部结构信息不足、难以完整表达数据的全局流形特征，本文采用指数图拉普拉斯来增强样本间的全局结构建模能力，定义为：

\mathrm{e}\mathrm{x}\mathrm{p}\left(\mathit{L}\right)={\displaystyle \sum _{k=\mathrm{0}}^{\mathrm{\infty }}}\frac{{\mathit{L}}^k}{k!}（9）

在实际数值实验中，矩阵指数通过有限阶近似进行计算。矩阵指数可看作是将样本之间的高阶连接关系引入建模过程，不仅直接考虑相邻的样本对，而且间接相连的样本对（如两点之间存在多跳连接）也能通过高阶项得到适度关联，以便捕捉到更平滑、更具全局一致性的流形结构。基于指数图拉普拉斯的平滑度测度表示为：

R_{\mathrm{e}\mathrm{x}\mathrm{p}}=\mathrm{T}\mathrm{r}\left[{\mathit{V}}^T\mathrm{e}\mathrm{x}\mathrm{p}\left(\mathit{L}\right)\mathit{V}\right]（10）

与传统的\mathrm{T}\mathrm{r}\left({\mathit{V}}^{\mathrm{T}}\mathit{L}\mathit{V}\right)相比，该形式在谱空间中对不同特征分量进行指数加权，可自适应地强化低频结构（全局相似性）并抑制高频噪声成分，在样本稀缺时仍能保持嵌入结果的平滑性与全局一致性。

算法在Semi-NMF框架基础上引入指数图正则化、稀疏性约束和正交性约束，以增强算法在样本有限情况下的鲁棒性与判别性。在式（8）中，第1项{‖\mathit{X}-\mathit{U}{\mathit{V}}^{\mathrm{T}}‖}_{\mathrm{2,1}}为重构误差项，这里的L_{\mathrm{2,1}}范数能有效抑制噪声和异常值，在样本量较少时保持稳定重构。第2项\mathrm{T}\mathrm{r}\left({\mathit{V}}^{\mathrm{T}}\mathrm{e}\mathrm{x}\mathrm{p}\left(\mathit{L}\right)\mathit{V}\right)为指数图正则化项，旨在通过引入图拉普拉斯矩阵的指数形式刻画样本间的非线性流形结构，弥补小样本条件下的结构信息不足。第3项{‖\mathit{V}‖}_{\mathrm{1}}为稀疏约束项，促使低维表示更加紧凑，有助降低维度冗余并防止过拟合。第4项{‖{\mathit{V}}^{\mathrm{T}}\mathit{V}-\mathit{I}‖}_{\mathrm{F}}^{\mathrm{2}}为正交约束项，用于约束特征表示之间的相关性，使不同特征维度尽可能独立，提升低维空间的判别能力和信息利用效率。

2.1 更新规则

在（8）式中，目标函数的增广拉格朗函数为：

\begin{array}{l}\mathrm{ℒ}=\mathrm{T}\mathrm{r}[(\mathit{X}-\mathit{U}{\mathit{V}}^{\mathrm{T}})\mathit{H}(\mathit{X}-\mathit{U}{\mathit{V}}^{\mathrm{T}}{)}^{\mathrm{T}}]+\\ \frac{\lambda }{\mathrm{2}}\mathrm{T}\mathrm{r}[{\mathit{V}}^{\mathrm{T}}\mathrm{e}\mathrm{x}\mathrm{p}(\mathit{L})\mathit{V}]+\alpha {‖\mathit{V}‖}_{\mathrm{1}}+\\ \frac{\beta }{\mathrm{2}}\mathrm{T}\mathrm{r}({\mathit{V}}^{\mathrm{T}}\mathit{V}{\mathit{V}}^{\mathrm{T}}\mathit{V})-\beta \mathrm{T}\mathrm{r}({\mathit{V}}^{\mathrm{T}}\mathit{V})+\\ \frac{\mathrm{1}}{\mathrm{2}}\beta \mathrm{T}\mathrm{r}(\mathit{I}{\mathit{I}}^{\mathrm{T}})-\mathrm{T}\mathrm{r}(\mathrm{\Psi }{\mathit{V}}^{\mathrm{T}})=\\ \mathrm{T}\mathrm{r}(\mathit{X}\mathit{H}{\mathit{X}}^{\mathrm{T}})-\mathrm{2}\mathrm{T}\mathrm{r}(\mathit{U}{\mathit{V}}^{\mathrm{T}}\mathit{H}{\mathit{X}}^{\mathrm{T}})+\\ \mathrm{T}\mathrm{r}(\mathit{U}{\mathit{V}}^{\mathrm{T}}\mathit{H}\mathit{V}{\mathit{U}}^{\mathrm{T}})+\mathrm{2}\lambda \mathrm{T}\mathrm{r}({\mathit{V}}^{\mathrm{T}}\mathrm{e}\mathrm{x}\mathrm{p}(\mathit{L})\mathit{V})+\\ \alpha {‖\mathit{V}‖}_{\mathrm{1}}+\frac{\beta }{\mathrm{2}}\mathrm{T}\mathrm{r}({\mathit{V}}^{\mathrm{T}}\mathit{V}{\mathit{V}}^{\mathrm{T}}\mathit{V})-\\ \beta \mathrm{T}\mathrm{r}({\mathit{V}}^{\mathrm{T}}\mathit{V})+\frac{\beta }{\mathrm{2}}\mathrm{T}\mathrm{r}(\mathit{I}{\mathit{I}}^{\mathrm{T}})-\mathrm{T}\mathrm{r}(\mathrm{\Psi }{\mathit{V}}^{\mathrm{T}})\end{array}（11）

其中，\mathit{H}为对角矩阵，第i个对角元素{\mathit{H}}_{ii}=\mathrm{1}/{‖\mathit{X}-\mathit{U}{\mathit{V}}^{\mathrm{T}}‖}_{\mathrm{2}}。对式（11）分别关于\mathit{U}和\mathit{V}求导数，有

\begin{array}{l}\frac{\partial \mathrm{ℒ}}{\partial \mathit{V}}=-\mathrm{2}\mathit{H}{\mathit{X}}^{\mathrm{T}}\mathit{U}+\mathrm{2}\mathit{H}\mathit{V}{\mathit{U}}^{\mathrm{T}}\mathit{U}+\lambda \mathrm{e}\mathrm{x}\mathrm{p}(\mathit{L})\mathit{V}+\\ \alpha \mathit{I}+\mathrm{2}\beta (\mathit{V}{\mathit{V}}^{\mathrm{T}}\mathit{V})-\mathrm{2}\beta \mathit{V}-\mathrm{\Psi }-\\ \mathrm{2}\mathit{H}{\mathit{X}}^{\mathrm{T}}({\mathit{U}}^{+}-{\mathit{U}}^{-})+\mathrm{2}\mathit{H}\mathit{V}({\mathit{U}}^{\mathrm{T}}{\mathit{U})}^{+}-\\ \mathrm{2}\mathit{H}\mathit{V}({\mathit{U}}^{\mathrm{T}}{\mathit{U})}^{-}+\lambda \mathrm{e}\mathrm{x}\mathrm{p}(\mathit{L})\mathit{V}+\alpha \mathit{I}+\\ \mathrm{2}\beta (\mathit{V}{\mathit{V}}^{\mathrm{T}}\mathit{V})-\mathrm{2}\beta \mathit{V}-\mathrm{\Psi }\end{array}（12）

\frac{\partial \mathrm{ℒ}}{\partial \mathit{U}}=-\mathrm{2}\mathit{X}\mathit{H}\mathit{V}+\mathrm{2}\mathit{U}{\mathit{V}}^{\mathrm{T}}\mathit{H}\mathit{V}（13）

其中，{\mathit{U}}^{+}=\left(\left|\mathit{U}\right|+\mathit{U}\right)/\mathrm{2}，{\mathit{U}}^{-}=\left(\left|\mathit{U}\right|-\mathit{U}\right)/\mathrm{2}。利用KKT条件有{\mathrm{\Psi }}_{jk}{\mathit{V}}_{jk}=\mathrm{0}。令\frac{\partial \mathrm{ℒ}}{\partial \mathit{V}}=\mathrm{0}，得：

\begin{array}{c}\begin{array}{l}(-\mathrm{2}\mathit{H}{\mathit{X}}^{\mathrm{T}}({\mathit{U}}^{+}-{\mathit{U}}^{-})+\mathrm{2}\mathit{H}\mathit{V}({\mathit{U}}^{\mathrm{T}}{\mathit{U})}^{+}-\\ \mathrm{2}\mathit{H}\mathit{V}({\mathit{U}}^{\mathrm{T}}{\mathit{U})}^{-}{)}_{jk}{\mathit{V}}_{jk}+(\lambda \mathrm{e}\mathrm{x}\mathrm{p}(\mathit{L})\mathit{V}+\\ \alpha \mathit{I}+\mathrm{2}\beta (\mathit{V}{\mathit{V}}^{\mathrm{T}}{\mathit{V})-\mathrm{2}\beta \mathit{V})}_{jk}{\mathit{V}}_{jk}=\mathrm{0}\end{array}\end{array}（14）

本文通过以下方式更新\mathit{V}：

{\mathit{V}}_{jk}\leftarrow {\mathit{V}}_{\mathit{j}k}\frac{[\mathrm{2}\mathit{H}{\mathit{X}}^{\mathrm{T}}{\mathit{U}}^{+}+\mathrm{2}\mathit{H}\mathit{V}({\mathit{U}}^{\mathrm{T}}{\mathit{U})}^{-}{+\mathrm{2}\beta \mathit{V}]}_{\mathit{j}\mathit{k}}}{[\mathrm{2}\mathit{H}{\mathit{X}}^{\mathrm{T}}{\mathit{U}}^{-}+\mathrm{2}\mathit{H}\mathit{V}({\mathit{U}}^{\mathrm{T}}{\mathit{U})}^{+}+\lambda \mathrm{e}\mathrm{x}\mathrm{p}(\mathit{L})\mathit{V}+\mathrm{2}\beta (\mathit{V}{\mathit{V}}^{\mathrm{T}}{\mathit{V})]}_{\mathit{j}\mathit{k}}+\alpha }（15）

为加速算法收敛，修改公式（14）为：

\begin{array}{l}\mathrm{【}-\mathrm{2}\mathit{H}{\mathit{X}}^{\mathrm{T}}({\mathit{U}}^{+}-{\mathit{U}}^{-})+\mathrm{2}\mathit{H}\mathit{V}({\mathit{U}}^{\mathrm{T}}{\mathit{U})}^{+}-\mathrm{2}\mathit{H}\mathit{V}({\mathit{U}}^{\mathrm{T}}{\mathit{U})}^{-}{]}_{jk}{\mathit{V}}_{{}_{\mathit{j}\mathit{k}}}^{\mathrm{2}}+\\ [\lambda \mathrm{e}\mathrm{x}\mathrm{p}(\mathit{L})\mathit{V}+\alpha \mathit{I}+\mathrm{2}\beta (\mathit{V}{\mathit{V}}^{\mathrm{T}}{\mathit{V})-\mathrm{2}\beta \mathit{V}]}_{jk}{\mathit{V}}_{{}_{jk}}^{\mathrm{2}}=\mathrm{0}\end{array}（16）

这样，通过式（16）更新\mathit{V}的更新公式为：

{\mathit{V}}_{\mathit{j}\mathit{k}}\leftarrow {\mathit{V}}_{\mathit{j}\mathit{k}}\sqrt[]{\frac{[\mathrm{2}\mathit{H}{\mathit{X}}^{\mathrm{T}}{\mathit{U}}^{+}+\mathrm{2}\mathit{H}\mathit{V}({\mathit{U}}^{\mathrm{T}}{\mathit{U})}^{-}{+\mathrm{2}\beta \mathit{V}]}_{\mathit{j}\mathit{k}}}{[\mathrm{2}\mathit{H}{\mathit{X}}^{\mathrm{T}}{\mathit{U}}^{-}+\mathrm{2}\mathit{H}\mathit{V}({\mathit{U}}^{\mathrm{T}}{\mathit{U})}^{+}+\lambda \mathrm{e}\mathrm{x}\mathrm{p}(\mathit{L})\mathit{V}+\mathrm{2}\beta (\mathit{V}{\mathit{V}}^{\mathrm{T}}{\mathit{V})]}_{\mathit{j}\mathit{k}}+\alpha }}（17）

进一步，令\frac{\partial \mathrm{ℒ}}{\partial \mathit{U}}=\mathrm{0}，得\mathit{U}的更新公式：

\mathit{U}=\mathit{X}\mathit{H}\mathit{V}{\left({\mathit{V}}^{\mathrm{T}}\mathit{H}\mathit{V}\right)}^{-\mathrm{1}}(18)

综上，算法的求解步骤如以下算法所示。

算法1 RNMFSCEG算法

输入： 初始矩阵\mathit{X}\in {\mathbb{R}}_{}^{m\times n}，类别数k，正则化参数\lambda，稀疏性参数\alpha，正交约束系数\beta，{\epsilon }_{\mathrm{0}}={\mathrm{10}}^{-\mathrm{7}}以及最大迭代次数t_{\mathrm{m}\mathrm{a}\mathrm{x}}=\mathrm{200}；

输出： 基矩阵\mathit{U}，系数矩阵\mathit{V}。

1） 初始化：随机初始化基矩阵\mathit{U}和系数矩阵\mathit{V}；

2） 固定\mathit{U}，根据式（17）更新系数矩阵\mathit{V}；

3） 固定\mathit{V}，根据式（18）更新基矩阵\mathit{U}；

4） 如果收敛条件{‖\mathit{X}-\mathit{U}{\mathit{V}}^{\mathrm{T}}‖}_{\mathrm{2,1}}\le {\epsilon }_{\mathrm{0}}和迭代次数达到预期；

5） 结束；

6） 得到聚类指示矩阵\mathit{V}，并运用k-means对系数矩阵\mathit{V}进行聚类。

2.2 收敛性分析

本节利用辅助函数法证明目标函数（8）的收敛性。

定理2.1 目标函数（8）在更新规则（17）下不增。

证明 设G\left({\mathit{V}}_{jk},{\mathit{V}}_{jk}^t\right)是\mathrm{ℒ}\left({\mathit{V}}_{jk}\right)的辅助函数，满足条件：

G\left({\mathit{V}}_{jk},{\mathit{V}}_{jk}^t\right)\ge \mathrm{ℒ}\left({\mathit{V}}_{jk}\right),G\left({\mathit{V}}_{jk},{\mathit{V}}_{jk}^{}\right)=\mathrm{ℒ}\left({\mathit{V}}_{jk}\right)（19）

其中，辅助函数G\left({\mathit{V}}_{jk},{\mathit{V}}_{jk}^t\right)定义为：

\begin{array}{l}G\left({\mathit{V}}_{jk},{\mathit{V}}_{jk}^t\right)=\sum _{j,k}{\mathit{X}}_{jk}{\mathit{H}}_{jk}{\mathit{X}}_{jk}-\\ \begin{array}{cc}& \end{array}\mathrm{2}\sum _{j,k}{\left(\mathit{H}{\mathit{X}}^{\mathrm{T}}\mathit{U}\right)}_{jk}^{+}{\mathit{V}}_{jk}^t\left(\mathrm{1}+\mathrm{l}\mathrm{o}\mathrm{g}\frac{{\mathit{V}}_{jk}}{{\mathit{V}}_{jk}^t}\right)+\\ \begin{array}{l}\begin{array}{cc}& \end{array}\sum _{j,k}{\left(\mathit{H}{\mathit{X}}^{\mathrm{T}}\mathit{U}\right)}_{jk}^{-}\frac{{\mathit{V}}_{jk}^{\mathrm{2}}+{\left({\mathit{V}}_{jk}^t\right)}^{\mathrm{2}}}{{\mathit{V}}_{jk}^t}+\\ \begin{array}{cc}& \end{array}\sum _{j,k}{\left(\mathit{H}\mathit{V}{\mathit{U}}^{\mathrm{T}}\mathit{U}\right)}_{jk}^{+}\frac{{\mathit{V}}_{jk}^{\mathrm{2}}}{{\mathit{V}}_{jk}^t}-\end{array}\\ \begin{array}{l}\begin{array}{cc}& \end{array}\sum _{j,k}{\left(\mathit{H}\mathit{V}{\mathit{U}}^{\mathrm{T}}\mathit{U}\right)}_{jk}^{-}{\mathit{V}}_{jk}^t\left(\mathrm{1}+\mathrm{l}\mathrm{o}\mathrm{g}\frac{{\mathit{V}}_{jk}}{{\mathit{V}}_{jk}^t}\right)+\\ \begin{array}{cc}& \end{array}\frac{\lambda }{\mathrm{2}}\sum _{j,k}{\left(\mathrm{e}\mathrm{x}\mathrm{p}\left(\mathit{L}\right){\mathit{V}}^t\right)}_{jk}\frac{{\mathit{V}}_{jk}^{\mathrm{2}}}{{\mathit{V}}_{jk}^t}+\end{array}\\ \begin{array}{cc}& \end{array}\alpha \sum _{j,k}\frac{{\mathit{V}}_{\mathit{j}\mathit{k}}^{\mathrm{2}}+{\left({\mathit{V}}_{jk}^t\right)}^{\mathrm{2}}}{\mathrm{2}{\mathit{V}}_{jk}^t}+\frac{\beta }{\mathrm{2}}\sum _{j,k}{\left(\mathit{V}{\mathit{V}}^{\mathrm{T}}\mathit{V}\right)}_{jk}\frac{{\mathit{V}}_{j,k}^{\mathrm{2}}}{{\mathit{V}}_{jk}^t}-\\ \begin{array}{l}\begin{array}{cc}& \end{array}\beta \sum _{j,k}{\mathit{V}}_{jk}{\mathit{V}}_{jk}^t\left(\mathrm{1}+\mathrm{l}\mathrm{o}\mathrm{g}\frac{{\mathit{V}}_{jk}}{{\mathit{V}}_{jk}^t}\right)+\frac{\beta }{\mathrm{2}}\mathrm{T}\mathrm{r}\left(\mathit{I}{\mathit{I}}^{\mathrm{T}}\right)-\\ \mathrm{T}\mathrm{r}\left(\mathrm{\Psi }{\mathit{V}}_{jk}^{\mathrm{T}}\right)\end{array}\end{array}（20）

关于\mathit{V}的增广拉格朗日函数\mathrm{ℒ}\left({\mathit{V}}_{jk}\right)可写为：

\begin{array}{l}\mathrm{ℒ}\left({\mathit{V}}_{\mathit{j}\mathit{k}}\right)=\mathrm{T}\mathrm{r}(\mathit{X}\mathit{H}{\mathit{X}}^{\mathrm{T}})-\mathrm{2}\mathrm{T}\mathrm{r}(\mathit{U}{\mathit{V}}^{\mathrm{T}}\mathit{H}{\mathit{X}}^{\mathrm{T}}{)}^{+}-\\ \mathrm{T}\mathrm{r}(\mathit{U}{\mathit{V}}^{\mathrm{T}}\mathit{H}\mathit{V}{\mathit{U}}^{\mathrm{T}}{)}^{-}+\mathrm{2}\mathrm{T}\mathrm{r}(\mathit{U}{\mathit{V}}^{\mathrm{T}}\mathit{H}{\mathit{X}}^{\mathrm{T}}{)}^{-}+\end{array}

\begin{array}{c}\begin{array}{l}\mathrm{T}\mathrm{r}(\mathit{U}{\mathit{V}}^{\mathrm{T}}\mathit{H}\mathit{V}{\mathit{U}}^{\mathrm{T}}{)}^{+}+\frac{\lambda }{\mathrm{2}}\mathrm{T}\mathrm{r}({\mathit{V}}^{\mathrm{T}}\mathrm{e}\mathrm{x}\mathrm{p}(\mathit{L})\mathit{V})+\\ \alpha {‖\mathit{V}‖}_{\mathrm{1}}+\frac{\beta }{\mathrm{2}}\mathrm{T}\mathrm{r}({\mathit{V}}^{\mathrm{T}}\mathit{V}{\mathit{V}}^{\mathrm{T}}\mathit{V})-\beta \mathrm{T}\mathrm{r}({\mathit{V}}^{\mathrm{T}}\mathit{V})+\\ \frac{\beta }{\mathrm{2}}\mathrm{T}\mathrm{r}(\mathit{I}{\mathit{I}}^{\mathrm{T}})-\mathrm{T}\mathrm{r}(\mathrm{\Psi }{\mathit{V}}^{\mathrm{T}})\end{array}\end{array}（21）

利用式（19），有

\begin{array}{l}\mathrm{ℒ}\left({\mathit{V}}_{jk}^{t+\mathrm{1}}\right)\le \mathit{G}\left({\mathit{V}}_{jk}^{t+\mathrm{1}},{\mathit{V}}_{jk}^t\right)\le \\ \mathit{G}\left({\mathit{V}}_{jk}^t,{\mathit{V}}_{jk}^t\right)=\mathrm{ℒ}\left({\mathit{V}}_{jk}^t\right),\end{array}

从而，\mathrm{ℒ}\left({\mathit{V}}_{jk}\right)从t迭代到t+\mathrm{1}的更新式：

{\mathit{V}}_{jk}^{t+\mathrm{1}}=\mathrm{a}\mathrm{r}\mathrm{g}\underset{\mathit{V}}{\mathrm{m}\mathrm{i}\mathrm{n}}G\left({\mathit{V}}_{jk},{\mathit{V}}_{jk}^t\right)（22）

不增。

下面证明G\left({\mathit{V}}_{jk},{\mathit{V}}_{jk}^t\right)是\mathrm{ℒ}\left({\mathit{V}}_{jk}\right)的辅助函数。显然，

G\left({\mathit{V}}_{jk},{\mathit{V}}_{jk}\right)=\mathrm{ℒ}\left({\mathit{V}}_{jk}\right)。

只需证明G\left({\mathit{V}}_{jk},{\mathit{V}}_{jk}^t\right)\ge \mathrm{ℒ}\left({\mathit{V}}_{jk}\right)。当\mathit{A}\in {\mathbb{R}}_{+}^{n\times n}、\mathit{B}\in {\mathbb{R}}_{+}^{p\times p}是对角矩阵时，以下不等式成立：

\mathrm{T}\mathrm{r}\left({\mathit{S}}^{\mathrm{T}}\mathit{A}\mathit{S}\mathit{B}\right)\le {\displaystyle \sum _{i=\mathrm{1}}^n}{\displaystyle \sum _{j=\mathrm{1}}^p}\left(\mathit{A}{\mathit{S}}^{\mathrm{T}}\mathit{B}\right)\frac{{\mathit{S}}_{ij}^{\mathrm{2}}}{{\mathit{S}}_{ij}^t}（23）

利用不等式（23）可得：

\left\{\begin{array}{l}\mathrm{T}\mathrm{r}({\mathit{V}}^{\mathrm{T}}\mathit{V}{\mathit{V}}^{\mathrm{T}}\mathit{V})\le \sum _{j,k}{\left(\mathit{V}{\mathit{V}}^{\mathrm{T}}\mathit{V}\right)}_{jk}\frac{{\mathit{V}}_{jk}^{\mathrm{2}}}{{\mathit{V}}_{jk}^t},\begin{array}{cc}& \end{array}\\ \mathrm{T}\mathrm{r}({\mathit{V}}^{\mathrm{T}}\mathrm{e}\mathrm{x}\mathrm{p}(\mathit{L})\mathit{V})\le \sum _{j,k}{\left(\mathrm{e}\mathrm{x}\mathrm{p}\left(\mathit{L}\right){\mathit{V}}^t\right)}_{jk}\frac{{\mathit{V}}_{jk}^{\mathrm{2}}}{{\mathit{V}}_{jk}^t},\\ \mathrm{T}\mathrm{r}{\left({\mathit{V}}^{\mathrm{T}}\mathit{H}\mathit{V}{\mathit{U}}^{\mathrm{T}}\mathit{U}\right)}^{+}\le \sum _{j,k}{\left(\mathit{H}\mathit{V}{\mathit{U}}^{\mathrm{T}}\mathit{U}\right)}_{jk}^{+}\frac{{\mathit{V}}_{jk}^{\mathrm{2}}}{{\mathit{V}}_{jk}^t}\end{array}\right.（24）

由不等式{\mathit{V}}_{\mathit{j}\mathit{k}}\mathrm{1}+\mathrm{l}\mathrm{o}\mathrm{g}{\mathit{V}}_{\mathit{j}\mathit{k}}\left({\mathit{V}}_{\mathit{j}\mathit{k}}\mathrm{0}\right)可得：

\left\{\begin{array}{l}\mathrm{T}\mathrm{r}({\mathit{V}}^{\mathrm{T}}\mathit{V})\ge \beta \sum _{j,k}{\mathit{V}}_{jk}{\mathit{V}}_{jk}^t\left(\mathrm{1}+\mathrm{l}\mathrm{o}\mathrm{g}\frac{{\mathit{V}}_{jk}}{{\mathit{V}}_{jk}^t}\right),\\ \begin{array}{l}\mathrm{T}\mathrm{r}(\mathit{U}{\mathit{V}}^{\mathrm{T}}\mathit{H}{\mathit{X}}^{\mathrm{T}}{)}^{+}\ge \\ \sum _{j,k}{\left(\mathit{H}{\mathit{X}}^{\mathrm{T}}\mathit{U}\right)}_{jk}^{+}{\mathit{V}}_{jk}^t\left(\mathrm{1}+\mathrm{l}\mathrm{o}\mathrm{g}\frac{{\mathit{V}}_{jk}}{{\mathit{V}}_{jk}^t}\right),\end{array}\\ \begin{array}{l}\mathrm{T}\mathrm{r}(\mathit{U}{\mathit{V}}^{\mathrm{T}}\mathit{H}\mathit{V}{\mathit{U}}^{\mathrm{T}}{)}^{-}\ge \\ \sum _{j,k}{\left(\mathit{H}\mathit{V}{\mathit{U}}^{\mathrm{T}}\mathit{U}\right)}_{jk}^{-}{\mathit{V}}_{jk}^t\left(\mathrm{1}+\mathrm{l}\mathrm{o}\mathrm{g}\frac{{\mathit{V}}_{jk}}{{\mathit{V}}_{jk}^t}\right)\end{array}\end{array}\right.（25）

利用不等式a^{\mathrm{2}}+b^{\mathrm{2}}\ge \mathrm{2}ab\left(a,b\ge \mathrm{0}\right)可得：

\mathrm{T}\mathrm{r}(\mathit{U}{\mathit{V}}^{\mathrm{T}}\mathit{H}{\mathit{X}}^{\mathrm{T}}{)}^{-}\le \sum _{j,k}{\left(\mathit{H}{\mathit{X}}^{\mathrm{T}}\mathit{U}\right)}_{jk}^{-}\frac{{\mathit{V}}_{jk}^{\mathrm{2}}+{\left({\mathit{V}}_{jk}^t\right)}^{\mathrm{2}}}{{\mathit{V}}_{jk}^t}（26）

{‖\mathit{V}‖}_{\mathrm{1}}=\sum _{j,k}{\mathit{V}}_{jk}\le \sum _{j,k}\frac{{\mathit{V}}_{jk}^{\mathrm{2}}+{\left({\mathit{V}}_{jk}^t\right)}^{\mathrm{2}}}{\mathrm{2}{\mathit{V}}_{jk}^t}（27）

综合式（24）~式（27）可得G\left({\mathit{V}}_{jk},{\mathit{V}}_{jk}^t\right)\ge \mathrm{ℒ}\left({\mathit{V}}_{jk}\right)，即证得G\left({\mathit{V}}_{jk},{\mathit{V}}_{jk}^t\right)是\mathrm{ℒ}\left({\mathit{V}}_{jk}\right)的辅助函数。

考虑到G\left({\mathit{V}}_{jk},{\mathit{V}}_{jk}^t\right)是一个凸优化，令\partial G\left({\mathit{V}}_{jk},{\mathit{V}}_{jk}^t\right)/\partial {\mathit{V}}_{jk}=\mathrm{0}可得{\mathit{V}}_{jk}的迭代规则为：

{\mathit{V}}_{jk}\leftarrow {\mathit{V}}_{jk}\sqrt[]{\frac{[\mathrm{2}\mathit{H}{\mathit{X}}^{\mathrm{T}}{\mathit{U}}^{+}+\mathrm{2}\mathit{H}\mathit{V}({\mathit{U}}^{\mathrm{T}}{\mathit{U})}^{-}{+\mathrm{2}\beta \mathit{V}]}_{jk}}{[\mathrm{2}\mathit{H}{\mathit{X}}^{\mathrm{T}}{\mathit{U}}^{-}+\mathrm{2}\mathit{H}\mathit{V}({\mathit{U}}^{\mathrm{T}}{\mathit{U})}^{+}+\lambda \mathrm{e}\mathrm{x}\mathrm{p}(\mathit{L})\mathit{V}+\mathrm{2}\beta (\mathit{V}{\mathit{V}}^{\mathrm{T}}{\mathit{V})]}_{jk}+\alpha }}（28）

综上，目标函数式（8）在更新规则式（17）下不增。证毕。

2.3 复杂度分析

本节给出算法的计算复杂度。该算法的复杂度主要取决于各项计算的矩阵操作。设m为特征维度，n为样本数量，k为分解后子矩阵维度，t为最大迭代次数，则算法整体复杂度为​O(t{n}^{\mathrm{3}}+t{n}^{\mathrm{2}}k+tn{k}^{\mathrm{2}}+tmnk)，其中更新\mathit{U}的计算复杂度为O(k^{\mathrm{3}}+n{k}^{\mathrm{2}}+mnk)，更新\mathit{V}的计算复杂度为O(n^{\mathrm{2}}k+n{k}^{\mathrm{2}}+mnk)。表1给出了作为比较对象的其他8种算法的时间复杂度。可以看到，相较于NMF和Semi-NMF等算法，本文算法的计算量较大。

3 实验验证

为了验证算法的聚类性能，本文采用聚类纯度（Purity，PUR）、调整兰德指数（Adjusted Rand Index，ARI）、聚类精度（Clustering Accuracy，ACC）和归一化互信息（Normalized Mutual Information，NMI）等4个评价指标，在9个公开数据集上与NMF^［8］、Semi-NMF^［26］、GNMF^［17］、REGNMF^［27］、ONMF^［28］、RNMF^［11］、OGNMFSC_UU^［19］和NMFOS^［29］等8个聚类算法进行比较。本文使用Python 3.8进行实验，计算环境为Intel Celeron N5095，64位Windows操作系统，运行内存为8 GB。

3.1 数据集

选取9个公共数据集进行实验，所有数据集均来源自UCI机器学习库，数据网址为http：//archive.ics.uci.edu/ml/datasets.html。鉴于本文算法能够处理混合数据，本文选取的数据包括8个非负数据集及1个包含负值数的数据集。9个数据集的基本情况如表2所示。需要说明的是，虽然Sonar、Ecg和Haberman并非严格意义上的高维小样本数据集，但由于样本量有限，在模型训练时同样容易出现过拟合和泛化性能不足的问题。例如，虽然Ecg数据集的特征维度处于中等水平，但相较于有限样本量，其数据仍呈现出典型的小样本特征。为直观展示9个数据集的分布特征，本文绘制了数据集的部分特征散点图，如图1所示。可以看到，数据集中存在噪声和异常值。

3.2 评价指标

为了更好地评价本文算法的性能和鲁棒性，我们采用PUR、ARI、ACC和NMI等4个指标。

PUR是计算聚类结果中每个簇中最多数目的样本所属的类别，用于衡量聚类结果中样本被正确分类到其真实类别中的比例，取值范围为\mathrm{0}\le V_{\mathrm{P}\mathrm{U}\mathrm{R}}\le \mathrm{1}，定义为：

V_{\mathrm{P}\mathrm{U}\mathrm{R}}=\frac{\mathrm{1}}{n}{\displaystyle \sum _{j=\mathrm{1}}^n}\underset{j}{\mathrm{m}\mathrm{a}\mathrm{x}}\left(n_i^j\right)（29）

其中，n_i^j表示聚类i中同样属于原始类j的样本。

ARI被用于衡量两个数据分组之间的一致性，通过ARI以避免随机性偏差，计算公式为：

V_{\mathrm{A}\mathrm{R}\mathrm{I}}=\frac{\mathrm{I}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{x}-\mathrm{E}\mathrm{x}\mathrm{p}\mathrm{e}\mathrm{t}\mathrm{c}\mathrm{t}\mathrm{e}\mathrm{d}\mathrm{I}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{x}}{\mathrm{M}\mathrm{a}\mathrm{x}\mathrm{I}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{x}-\mathrm{E}\mathrm{x}\mathrm{p}\mathrm{e}\mathrm{t}\mathrm{c}\mathrm{t}\mathrm{e}\mathrm{d}\mathrm{I}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{x}}（30）

其中，\mathrm{I}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{x}表示实际聚类和真实标签的相似性计数，\mathrm{E}\mathrm{x}\mathrm{p}\mathrm{e}\mathrm{t}\mathrm{c}\mathrm{t}\mathrm{e}\mathrm{d}\mathrm{I}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{x}表示随机情况下的期望相似性计数，\mathrm{M}\mathrm{a}\mathrm{x}\mathrm{I}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{x}为理论上的最大相似性计数。-\mathrm{1}\le V_{\mathrm{A}\mathrm{R}\mathrm{I}}\le \mathrm{1}，绝对值越大代表聚类效果越好。

给定实际真值标签l_i和学习到的聚类标签r_i，定义：

V_{\mathrm{A}\mathrm{C}\mathrm{C}}=\frac{{\displaystyle \sum _{i=\mathrm{1}}^n}\delta \left(l_i,\mathrm{m}\mathrm{a}\mathrm{p}(r_i)\right)}{n}（31）

其中，\mathrm{0}\le V_{\mathrm{A}\mathrm{C}\mathrm{C}}\le \mathrm{1}。如果x=y，则\delta \left(x,y\right)=\mathrm{1}，否则\delta \left(x,y\right)=\mathrm{0}。映射函数\mathrm{m}\mathrm{a}\mathrm{p}(\cdot )由Hungarian算法求解。

NMI值用于衡量两个数据分布的吻合程度，即两个事件集合之间的相关性，互信息越大，词条和类别的相关程度也越大，取值范围为\mathrm{0}\le V_{\mathrm{N}\mathrm{M}\mathrm{I}}\le \mathrm{1}。给定数据的真实标签C和通过算法聚类得到的标签\overline{C}，NMI的计算公式如下：

V_{\mathrm{N}\mathrm{M}\mathrm{I}}=\frac{\mathrm{M}\mathrm{I}\left(C,\overline{C}\right)}{\mathrm{m}\mathrm{a}\mathrm{x}\left(H(C)+H(\overline{C})\right)}（32）

其中，\mathrm{M}\mathrm{I}\left(C,\overline{C}\right)为簇C和\overline{C}的互信息，H(C)和H(\overline{C})分别为C和\overline{C}的熵。

3.3 基线算法

参与比较实验的主要有8种算法，分述如下。

1） NMF^［8］：NMF通过构建全局的重构误差函数衡量原始矩阵与分解后矩阵乘积之间的差异，在此基础上优化得到最优的因子矩阵解，以尽可能保留原始数据的主要结构和特征信息。

2） Semi-NMF^［26］：Semi-NMF是NMF的一个变形，要求系数矩阵是非负的而基矩阵是任意实数矩阵；该方法在保持非负数据特性的同时可以更好地逼近原始数据。

3） GNMF^［17］：GNMF在标准NMF的基础上引入了图正则项，利用图结构信息使低维表示更加贴近原始数据结构。

4） REGNMF^［27］。在GNMF基础上，REGNMF引入了指数图正则化项，将原拉普拉斯矩阵指数化，提高了算法的鲁棒性和聚类性能。

5） ONMF^［28］：ONMF在标准NMF基础上引入了正交约束，且要求分解子矩阵中的某一部分正交，以增强聚类能力。

6） RNMF^［11］：RNMF用稀疏项L_{\mathrm{2,1}}代替Frobenius范数来衡量重构误差，缓解了噪声和异常值的影响，增强了算法的鲁棒性。

7） OGNMFSC_UU^［19］：该算法在矩阵分解过程中引入了正交性、稀疏性等结构性约束，通过构建多目标联合优化模型，有效增强了因子矩阵的判别性与可解释性，在保留数据内在结构特征的同时显著提升了算法在聚类的鲁棒性。

8） NMFOS^［29］：在NMF中，会有冗余信息的干扰，为了提高数据间的独立性， NMFOS将正交约束引入到了目标函数中，提高了算法的性能。

对于含有图正则化项的算法，其近邻参数k统一设置为5。

3.4 聚类结果可视化

本节将直观展示9种算法的聚类性能。同时，为了验证本文算法在子空间学习上的有效性，本实验将算法和其他8种聚类算法学习到的低维数据投影到PCA降维后的2维空间中。限于篇幅，本文以Wdbc、Breath-cancer、Spiral、Wine等4个数据集为代表展示聚类效果，如图2~图5所示。其中，同一类别的数据节点以相同颜色标识，以增强对比效果。为降低随机因素的影响，实验重复进行30次并记录聚类性能的平均值及标准误差。

3.5 SSS实验

本节验证算法处理小样本问题的有效性。采用经过处理的SSS数据集，在保持原始特征维度的前提下，将K1a reduced样本量缩减至1325个，Wine reduced样本量缩减至12个，Breath-cancer reduced样本量缩减至30个。基于这些数据集，本文将研究算法与其他8种算法在样本量受限条件下的分类或聚类任务的准确性与鲁棒性进行了对比，结果如表3所示。主要评估指标包括ACC、NMI、PUR和ARI值。

可以看到，在Wine reduced数据集上，本文算法的ACC、ARI、NMI和PUR值均最高，显著优于其他算8种算法，表明算法在小样本情况下仍能准确区分类别并保持良好的PUR。在Breath-cancer reduced数据集上，本文算法在各项指标上也保持优先，其中ACC、ARI、NMI和PUR值分别为0.701 50、 0.105 19、 0.157 68和0.701 50，在高维小样本环境下具有稳定性与鲁棒性。对于K1a reduced数据集，本文算法依然在4项指标上优于NMF、Semi-NMF、GNMF、REGNMF和ONMF算法，尤其是在ACC和NMI值上提升明显，说明算法在复杂SSS条件下具有更强的聚类能力。综上，本文算法在ACC、NMI和PUR值指标上表现突出，能够更有效地保留样本的类别结构信息，有效解决SSS聚类问题。

3.6 鲁棒性分析

为进一步验证本文算法的鲁棒性，本文设计了多噪声级别的对比实验。以小样本Wine数据集为例，本文分别在原始数据上添加5%、10%和15%等3种不同强度的高斯噪声，并与鲁棒算法RNMF和REGNMF比较，以评估算法在噪声递增环境下的性能表现。为保证实验结果的可靠性，每种方法均独立运行30次取平均值。在噪声干扰下，3种算法的聚类性能如表4和表5所示。可以看到，本文算法在各噪声水平下的4个评价指标均最高。在5%噪声下，RNMF算法的表现最弱，凸显出算法具有良好的鲁棒性。

3.7 旋转不变性

本文以Wdbc数据集为例对原始数据集施加随机正交变换，然后对比本文算法与RNMF、REGNMF算法在原始数据及旋转数据上的性能。本文通过评价指标ARI的差值来说明本文算法的特征旋转不变性，30次独立实验的平均结果如表6所示，图6则展示了Wdbc数据集旋转前后的可视化对比。

由表6可知，旋转实验揭示了本文算法具有旋转不变性，即在旋转扰动下几乎不受影响。在Wdbc数据集上，ΔARI值为0.038 4，保持率91.46%，比其他两种经典算法更优。

3.8 聚类结果分析

为了验证算法的聚类性能，本文在9个数据集上分别比较9种算法。聚类性能的定量对比结果参见表7~表10，其中最佳结果以加粗形式突出显示，“—”表示无运行结果。

从实验结果可以得到以下结论。

1） 从传统NMF算法以及Semi-NMF算法的结果来看，本文算法在SSS条件下展现出显著的性能优势。在Wdbc数据集上，传统NMF算法依赖于欧氏距离进行数据重构，对噪声样本和样本稀疏性较为敏感，本文算法则较NMF算法在ACC值上分别由0.808 34提升至0.838 85，同时ARI和NMI值等指标也明显提高。对比Semi-NMF算法，本文算法在负数集ECG上的ARI值提升4.6%。本文算法通过引入L_{\mathrm{2,1}}范数增强了鲁棒性，并结合指数图正则项强化了局部结构保持能力，从而使其在样本有限的情况下仍能准确挖掘潜在特征分布，实现更高的ACC。

2） 与结构保持类算法GNMF相比较，本文算法在小样本数据集上的表现更为稳健。GNMF算法在Sonar与Haberman等数据集上表现出一定程度的性能下降，而本文算法的ACC和NMI值分则别提升6.17%和1.82%。这种优势源于本文算法采用的指数图正则化形式使样本间相似度权重可自适应调整，更充分地反映了真实流形结构。相比GNMF算法依赖固定邻接图的方式，本文算法能够在样本量有限、分布稀疏的条件下优化邻域关系，有效缓解SSS导致的结构失真。

3） 在与鲁棒型分解算法的比较中，本文算法同样展现出更强的抗噪能力与泛化性。在K1a和Spiral数据集上，其ACC值分别较RNMF提升7.6%和5.1%，同时ARI与PUR值也同步提高。RNMF算法仅依靠鲁棒范数削弱异常值影响，而本文算法则在此基础上引入了指数映射的图正则化，使模型在高维低样本比环境下能够保持簇内紧凑性与簇间分离性，提升其整体聚类稳定性。这一改进充分体现了本文算法应对SSS问题时的理论优势与适用性。

4） 与集成多约束信息的复杂模型相比，本文算法在综合性能上依旧占据明显优势。以Australian和Wine数据集为例，其ACC值从OGNMFSC_UU算法的0.559 67提高至0.600 29，且NMI与PUR值分别提升4.42%与3.93%。尽管OGNMFSC_UU与NMFOS等算法在特征选择上引入了多重约束机制，但过度复杂的正则化结构在小样本情况下容易造成模型过拟合。而本文算法则通过指数图正则与鲁棒误差项的协同作用，更高效地学习到潜在子空间表示，从而确保模型在有限样本下仍具较强的特征表达与抗噪能力。

5） 在高维大样本数据集K1a上，本文算法相比OGNMFSC_UU算法的ACC值由0.302 34提升至0.403 11，且NMI和PUR值也同步提升，体现出算法在高维复杂数据环境中的适应性和鲁棒性。总体来看，本文算法在SSS数据集中有效缓解了小样本条件下的维数灾难问题，在高维大样本数据上同样取得了显著优势，具有广泛的适用性。

3.9 参数敏感性分析

本节讨论本文算法的聚类结果对超参数的敏感性。这些参数主要包括稀疏性参数\alpha、图正则化参数\lambda以及正交约束参数\beta，取值分别为\left\{{\mathrm{10}}^{-\mathrm{3}},{\mathrm{10}}^{-\mathrm{2}},{\mathrm{10}}^{-\mathrm{1}},{\mathrm{10}}^{\mathrm{0}},{\mathrm{10}}^{\mathrm{1}},{\mathrm{10}}^{\mathrm{2}},{\mathrm{10}}^{\mathrm{3}}\right\}。不同参数取值组合下的ACC、ARI 、NMI及PUR值如图7~图9所示。

鉴于图7~图9中参数的变化趋势是类似的，这里仅以图7为例进行分析。可以看到，对于参数\alpha和\beta，本文算法在Australian数据集上的ACC值大致集中在0.559左右，当\alpha =\mathrm{0.1}或\alpha =\mathrm{1}且\beta =\mathrm{0.001}以及\alpha =\mathrm{1}且\beta =\mathrm{10}时，ACC达到最高值0.562。整体来看，参数变化对算法在该数据集上的聚类性能影响较小，仅在\alpha =\mathrm{0.1}或1和\beta =\mathrm{0.001}时略有提升。在Breast-cancer数据集上，本文算法表现极其稳定，几乎所有参数组合下的ACC值均为0.626，仅有一个组合略低些，为0.624，表明该数据集对超参数\alpha和\beta调整不敏感，表现出高度鲁棒性。在Haberman数据集上，随参数的变化ACC值的波动幅度相较Breast-cancer数据集较小，最优ACC达到0.520。对于Sonar数据集，参数变化对ACC值的影响同样较为有限，取值主要集中在0.562~0.588之间，整体表现稳定。最后，在Spiral数据集上，当\alpha =\mathrm{0.01}且\beta =\mathrm{1000}时，聚类性能最优取值为0.413。

3.10 经验收敛性分析

本文算法的目标函数值呈单调递减趋势，在大部分数据集上迭代50次后，目标函数的值急剧下降并能够收敛到一个稳定水平，然后随着迭代次数的增加缓慢下降，最终达到收敛。本文算法在9个数据集上目标函数收敛速度的变化如图10所示。

4 结论

本文基于Semi-NMF、指数图正则化、稀疏性和正交性约束提出了一种稀疏约束与指数图正则化的RNMFSCEG。9个公共数据集上的算例分析结果表明，采用L_{\mathrm{2,1}}范数重构误差能够增强对噪声和异常值的鲁棒性，设计指数图正则化项能够有效捕捉数据的分布特征与局部几何结构，而引入稀疏约束则能够提升特征可解释性。同时，使用Semi-NMF框架扩展了算法处理混合符号数据的能力。本文算法在聚类准确性和鲁棒性等方面表现出色，并且在SSS场景下也展现出显著优势，在复杂小样本聚类任务中具有明显有效性与适用性。

该算法仍存在一定局限性。指数图正则化的指数矩阵计算增加了时间与存储成本，在大规模数据上可能影响算法效率。同时，模型中包含多个正则化参数可能需要调优，以平衡重构精度与泛化能力。未来研究可从降低指数图运算复杂度、自适应参数选择优化以及结合深度非负矩阵分解等方向进一步增强算法在更复杂SSS场景下的适应性。

参考文献

[1] Dong Y F， Deng Y H， Dong Y， et al.Survey of clustering based on deep learning ［J］.Journal of Compiter Applications，2022， 42（4）： 1021-1028.

[2] 董永峰， 邓亚晗， 董瑶， 等.基于深度学习的聚类综述［J］.计算机应用， 2022， 42（4）： 1021-1028.

[3] Xu M M， Hou X M.Density estimation clustering method based on reverse nearest neighbor ［J］.Computer Engineering and Applications， 2025， 61（1）： 165-173.

[4] 许梅梅， 侯新民.基于反向最近邻的密度估计聚类算法［J］.计算机工程与应用， 2025， 61（1）： 165-173.

[5] Yao X H， Guo M. Deep embedding clustering algorithm based on convolutional residual autoencoder ［J/OL］.Journal of Chongqing Technology and Business University （Natural Science Edition），2025-10-23.

[6] 姚晓红， 郭苗.基于卷积残差自编码器的深度嵌入聚类算法［J/OL］.重庆工商大学学报（自然科学版），2025-10-23.

[7] Saba T， Rehman A， Mujahid M， et al.Smart defense based on explainable stacked machine learning architecture for securing internet of health things with K-means clustering ［J］.Sci Rep， 2025， 15（1）： 37241.

[8] Yan S K， Ding H L， Qiu J Y， et al.Low-voltage distribution network topology identification method based on manifold density peak clustering ［J］.Journal of Electric Power Science and Technology， 2025， 40（4）： 61-71.

[9] 严绍奎，丁海丽，邱嘉怡，等.基于流形密度峰值聚类的低压配电网拓扑辨识方法［J］.电力科学与技术学报， 2025， 40（4）： 61-71.

[10] Luxburg U.A tutorial on spectral clustering ［J］.Stat Comput， 2007， 17（4）： 395-416.

[11] Fu X， Gu Z L， Qi B J， et al. Structural modal parameter automatic identification and order determination method based on clustering and neural networks ［J/OL］.Journal of Vibration Engineering，2025-08-06.

[12] 付兴， 顾政力， 祁宝金， 等.基于聚类和神经网络的结构模态参数自动识别与定阶方法［J/OL］.振动工程学报，2025-08-06.

[13] Lee D D， Seung H S.Learning the parts of objects by non-negative matrix factorization ［J］.Nature， 1999， 401（6755）： 788-791.

[14] Zhou J.Research of SWNMF with new iteration rules for facial feature extraction and recognition ［J］.Symmetry， 2019， 11（3）： 354.

[15] Qian Y， Tan C， Ding D， et al.Fast and secure distributed nonnegative matrix factorization ［J］.IEEE Trans Knowl Data Eng， 2020， 34（2）： 653-666.

[16] Kong D， Ding C， Huang H.Robust nonnegative matrix factorization using L₂₁-norm ［C］//Proceedings of the 20th ACM International Conference on Information and Knowledge Management.New York： ACM， 2011： 673-682.

[17] Lee D D， Seung H S.Algorithms for non-negative matrix factorization ［C］//Advances in Neural Information Processing Systems.Cambridge： MIT Press， 2001： 556-562.

[18] Hoyer P O.Nonnegative matrix factorization with sparseness constraints ［J］.J Mach Learn Res， 2004， 5（9）： 1457-1469.

[19] Gao Y， Church G.Improving molecular cancer class discovery through sparse non-negative matrix factorization ［J］.Bioinformatics， 2005， 21（21）： 3970-3975.

[20] Yang L D， Zhao Y J， Pan Z H.Research on L₂₁ incremental non-negative matrix factorization with sparsity constraints ［J］.Science Technology Information， 2024， 22（12）： 240-244.

[21] 杨亮东， 赵妍杰， 潘正红.稀疏约束的L₂₁增量式非负矩阵分解研究［J］.科技资讯， 2024， 22（12）： 240-244.

[22] Yang G L， Zhang J Q， Sheng Y Y.Based on L_1/2 sparsity and kurtosis smoothing constrained non-negative matrix factorization ［J］.Modern Information Technology， 2025， 9（5）： 45-50.

[23] 杨国亮， 张佳琦， 盛杨杨.基于L_1/2稀疏性和峰度平滑约束非负矩阵分解的高光谱图像解混［J］.现代信息科技， 2025， 9（5）： 45-50.

[24] Cai D， He X， Han J， et al.Graph regularized non-negative matrix factorization for data representation ［J］.IEEE Trans Pattern Anal Mach Intell， 2011， 33（8）： 1548-1560.

[25] Li S， Lu L， Liu Q， et al.Graph-regularized， sparsity-constrained non-negative matrix factorization with earth mover’s distance metric ［J］.Mathematics （Basel）， 2023， 11（8）： 1894.

[26] Chen Y， Qu G， Zhao J.Orthogonal graph regularized non-negative matrix factorization under sparse constraints for clustering ［J］.Expert Syst Appl， 2024， 249（Part C）： 123797.

[27] Gao H Y， Liu M S， Zhou G G， et al.Robust non-negative matrix factorization clustering algorithm with dual-graph regularization ［J］.Journal of Frontiers of Computer Science and Technology， 2026， 20（4）： 1061-1078.

[28] 高海燕， 刘孟淑， 周改改， 等.基于双图正则化的鲁棒非负矩阵分解聚类算法［J］.计算机科学与探索， 2026， 20（4）： 1061-1078.

[29] Kuo B C， Chang K Y.Feature extractions for small sample size classification problem ［J］.IEEE Trans Geosci Remote Sens， 2007， 45： 756-764.

[30] Wang S J， Chen H L， Peng X J， et al.Exponential locality preserving projections for small sample size problem ［J］.Neurocomputing， 2011， 74（17）： 3654-3662.

[31] Yuan S， Mao X.Exponential elastic preserving projections for facial expression recognition ［J］.Neurocomputing， 2017， 275： 711-724.

[32] Dornaika F， Bosaghzadeh A.Exponential local discriminant embedding and its application to face recognition ［J］.IEEE Trans Cybern， 2013， 43（3）： 921-934.

[33] Zhang T， Fang B， Tang Y Y， et al.Generalized discriminant analysis： A matrix exponential approach ［J］.IEEE Trans Systems Man Cybernetics Part B， 2010， 40（1）： 186-197.

[34] Ding C H Q， Li T， et al.Convex and semi-nonnegative matrix factorizations ［J］.IEEE Trans Pattern Anal Mach Intell， 2010， 32（1）： 45-55.

[35] Wan M， Cai M， Yang G.Robust exponential graph regularization non-negative matrix factorization technology for feature extraction ［J］.Mathematics （Basel）， 2023， 11（7）： 1716.

[36] Ding C H Q， Li T， Peng W， et al.Orthogonal nonnegative matrix tri-factorizations for clustering ［C］//Proceedings of the 12th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining.New York： ACM， 2006： 126-135.

[37] Li Z， Wu X， Peng H.Nonnegative matrix factorization on orthogonal subspace ［J］.Pattern Recognition Letters， 2010， 31（9）： 905-911.

基金资助

国家社会科学基金(21BTJ042)

甘肃省自然科学基金(23JRRA1186)

全国统计科学研究重点项目(2025LZ007)

甘肃省高校青年博士支持项目(2025QB-058)