0 引言

近年来，非负矩阵分解(nonnegative matrix factorization，NMF)算法在数据挖掘和机器学习领域受到越来越多的重视，在很多应用场景中发挥着作用，包括财务数据分析^[1]、基因表达分析^[2]、社交网络分析^[3]、视觉追踪^[4]、图像去噪^[5]、手背静脉识别^[6]、文本聚类^[7]、大规模异构数据联合聚类^[8].

给定非负矩阵X∈R^p×n和r＜＜min(p，n)，NMF算法将X分解为两个非负矩阵F∈R^p×r和G∈R^r×n的乘积，其中F被称为基矩阵，G被称为系数矩阵，这样高维空间(p维)的点x_i被投影到低维空间(r维)的点g_i.由于非负的约束，NMF算法得到的矩阵F和G具有稀疏性，从而提供了基于局部的特征表示.

标准NMF算法^[9]利用F范数(Frobenius norm)衡量矩阵分解的质量，然而数据中具有较大误差的噪音和异常值，能显著地影响目标函数值，算法的聚类性能因此下降.基于L₂₁范数的鲁棒NMF算法^[10](robust NMF using L₂₁ norm，L₂₁NMF算法)和基于CIM的rCIM-NMF算法^[11](row-based Correntropy Induced Metric NMF，rCIM-NMF算法)克服了标准NMF算法对噪音和异常值敏感的缺点，提升了算法的鲁棒性.基于Logdet散度的稀疏NMF算法^[12](Logdet divergence based sparse NMF，LDS-NMF算法)利用L₂₁范数作为矩阵分解的误差函数，同时约束基矩阵F的秩为r，并利用L₁范数正则项进一步增强系数矩阵G的稀疏性.基于线性投影结构的NMF算法^[13](Linear Projection-Based NMF，LPB-NMF算法)和线性投影NMF算法^[14](linear projective NMF，LP-NMF算法)将X∈R^p×n分解为基矩阵F∈R^p×r、变换矩阵Q∈R^r×p和X∈R^p×n三者的乘积.

流形是局部具有欧氏空间性质的空间，流形学习通过发现嵌入到高维空间中的低维流形，达到降维目的.流形学习算法包括等度量映射^[15](Isometric mapping，Isomap)，局部线性嵌入^[16](locally linear embedding，LLE)，拉普拉斯特征嵌入^[17](Laplacian Eigenmap，LE)等，它们利用了局部不变性，即高维空间中距离近的样本点在低维流形中的距离也近.谱聚类算法^[18](spectral clustering，SC算法)计算数据集的亲和矩阵，进一步得到嵌入的低维流形，再利用k均值(k-means)算法进行聚类.图正则的非负矩阵分解算法^[19](graph regularized NMF，GNMF算法)通过亲和矩阵来编码数据的几何结构，构建了能反映该几何结构的NMF算法.几何结构保持的非负矩阵分解算法^[20](structure preseving NMF，SPNMF算法)将近邻样本点间的相似性关系的保持和非近邻样本点间的互斥性关系的保持引入到NMF框架中.结构投影NMF算法^[21](structured projective NMF，SP-NMF算法)将图正则的方法引入到投影NMF模型中.

本文提出了一种基于投影梯度方法的鲁棒流形非负矩阵分解算法(robust nonnegative matrix factorization on manifold via projected garadient method，RMPGM-NMF)，该算法：

1) 使用L₂₁范数衡量矩阵分解的质量，对数据中噪音和异常值不敏感，克服了GNMF算法和标准NMF算法中使用F范数受数据中噪音和异常值影响大的缺点；

2) 利用数据的几何结构并考虑局部不变性，将流形学习与鲁棒NMF相结合，克服了L₂₁NMF算法，rCIM-NMF算法和LDS-NMF算法中未考虑数据内在的几何结构的缺点；

3) 利用投影梯度方法求解得到矩阵F和G的更新规则，避免了上述其它算法使用乘性更新规则时分母可能为0的情况.

1 相关NMF算法和SC算法

在这个部分，简单回顾下标准NMF算法、L₂₁NMF算法、rCIM-NMF算法、LDS-NMF算法、LP-NMF算法、SC算法和GNMF算法.输入矩阵定义为X=(x₁，x₂，…，x_n)，每个p维的列向量x_i表示一个非负的数据样本.

1.1 标准NMF算法

标准NMF算法利用F范数衡量X和FG的差值，定义为

(1)

其中是矩阵X的F范数.该问题通常由如下的乘性更新规则求解：

(2)

(3)

1.2 L₂₁NMF算法

L₂₁NMF算法利用L₂₁范数衡量矩阵分解的质量，定义为

(4)

其乘性更新规则为

(5)

(6)

其中，

(7)

1.3 rCIM-NMF算法

rCIM-NMF算法的目标函数为

(8)

其中f(e，σ)为高斯核函数：

(9)

其乘性更新规则为

(10)

(11)

其中，

(12)

(13)

1.4 LDS-NMF算法

LDS-NMF算法约束矩阵F的秩为r，同时引入L₁正则项增强矩阵G的稀疏性，其目标函数为

(14)

其中，λ和γ为可调参数，I_r为r维的单位矩阵，D_ld(A，A₀)的定义为

(15)

其中m为方阵A的维度.

LDS-NMF算法的乘性更新规则为

(16)

(17)

其中，

(18)

和分别表示矩阵的非负部分和负数部分.

1.5 LP-NMF算法

LP-NMF算法的目标函数为

(19)

其乘性更新规则为

(20)

(21)

1.6 SC算法

SC算法首先计算亲和矩阵A∈R^n×n，

(22)

再计算归一化的亲和矩阵L，

(23)

其中

找到矩阵L的k个最大特征向量y₁，y₂，…，y_k，构成矩阵Y=[y₁，y₂，…，y_k]∈R^n×k，令, 使得矩阵Y每一行向量的L₂范数为1，把矩阵Y的每一行向量视作k维空间的一个样本点(原始数据集的样本点x_i是p维的)，最后利用k-means算法在矩阵Y上聚类.

1.7 GNMF算法

GNMF算法引入拉普拉斯矩阵L，其定义为

(24)

其中L=V-W，W为亲和矩阵，相应的更新规则为

(25)

(26)

2 RMPGM-NMF算法模型

这部分将给出RMPGM-NMF算法的定义，为方便求解出更新规则，本文推导出了等价的目标函数.

2.1 RMPGM-NMF算法的定义

在L₂₁NMF算法中，每个数据样本的目标函数值为

(27)

而在标准NMF算法中，每个数据样本的目标函数值为

(28)

数据集中存在的具有较大误差的噪音和异常值会极大地影响标准NMF算法的目标函数值(28)，从而降低算法的聚类性能，而L₂₁NMF算法的目标函数值(27)受噪音和异常值的影响较小.

结合L₂₁NMF和流行学习，将拉普拉斯矩阵引入到L₂₁NMF算法中，提出了RMPGM-NMF算法，该算法保证鲁棒性的同时，利用数据的几何结构并考虑局部不变性，使得高维空间X中样本的几何结构在进行矩阵分解后的低维空间G中得以保持.其定义为

(29)

其中，L=V-W，L为拉普拉斯矩阵，W为亲和矩阵，，亲和矩阵W定义为

(30)

其中N_m(x)表示距离数据样本x最近的m个数据样本的集合，任意两个数据样本之间的距离定义为

(31)

2.2 等价的目标函数

直接由式(29)求解F和G的更新规则比较困难，本文将给出与式(29)等价的目标函数.

定理1  求解式(29)等价于求解式(32)：

(32)

其中D为对角矩阵，其对角元素为

(33)

证明  若求解式(32)得到更新规则：

(34)

(35)

使得

(36)

(37)

文[10]中的不等式(38)和不等式(39)为

(38)

(39)

由不等式(36)可得：

(40)

将其代入到不等式(38)中可得：

(41)

则：

(42)

由不等式(37)可得：

(43)

将其代入到不等式(39)中可得：

(44)

则：

(45)

不等式(42)和不等式(45)表明式(32)的更新规则(34)和(35)同样也使式(29)中目标函数值单调不递增，从而求解式(29)等价于求解式(32).

若X的某一行或者某一列全为0，乘性更新规则的分母可能为0，而使用投影梯度方法则不会遇到这样的问题.我们将在第3小节中利用投影梯度方法对式(32)进行求解.

3 RMPGM-NMF算法的更新规则

在这部分，先利用投影梯度方法求解一般约束优化问题，进一步利用投影梯度方法求解RMPGM-NMF算法中F和G的更新规则.

3.1 一般约束优化问题的投影梯度方法求解

一般约束优化问题的标准形式为

(46)

其中x是一个p维的列向量.

投影梯度方法使用如下规则更新当前解x^t到x^t+1：

(47)

其中，

(48)

不同的投影梯度方法对步长α的选择不同，算法1展示了一个简单有效的求解式(46)的投影梯度方法，每轮迭代中α的初始值设置为1(第1行)，不等式(49)可以保证每轮迭代中目标函数值f(x)有足够的下降，判断不等式(49)是否成立：若成立(第4行)，则不断增大α的值直至不等式(49)不成立，将α_pre代入式(47)中计算得到x^t+1并返回(第5行~第11行)，若不成立(第13行)，则不断减小α的值直至不等式(49)成立，将α_new代入式(47)中计算得到x^t+1并返回(第14行~第19行).

(49)

算法1  式(41)的更新规则

输入：x^t

输出：x^t+1

BEGIN

1) α_pre=1，β=0.1，σ=0.01

2) x^t+1=P[x^t-α_pre▽f(x^t)]

3) condition=(f(x^t+1)-f(x^t)≤σ(▽f(x^t))^T(x^t+1 -x^t))

4) IF (condition)

5) WHILE (condition)

6) α_new=α_pre/β

7) x^t+1=P(x^t-α_new▽f(x^t))

8) condition=(f(x^t+1)-f(x^t)≤σ(▽f(x^t))^T(x^t+1 -x^t))

9) END WHILE

10) x^t+1=P[x^t-α_pre▽f(x^t)]

11) return x^t+1

12) END IF

13) ELSE

14) WHILE (! condition)

15) α_new=α_pre×β

16) x^t+1=P(x^t-α_new▽f(x^t))

17) condition=(f(x^t+1)-f(x^t)≤σ(▽f(x^t))^T·(x^t+1-x^t))

18) END WHILE

19) return x^t+1

20) END ELSE

END

3.2 RMPGM-NMF算法的投影梯度方法求解

由第2部分的分析，求解式(29)等价于求解式(32)，记：

(50)

则：

(51)

(52)

于是RMPGM-NMF算法中G和F的更新规则分别如算法2和算法3所示，其中算法2中的〈▽_G^tf(F，G^t)，(G^t+1-G^t)〉表示两个矩阵▽_G^tf(F，G^t)和G^t+1-G^t的内积，算法3中的〈▽_F^tf(F^t，G)，(F^t+1-F^t)〉表示两个矩阵▽_F^tf(F^t，G)和F^t+1-F^t的内积，算法2和算法3中的函数P[·]定义为

(53)

算法2和算法3中最耗时的操作是式(54)的计算：

(54)

在每一轮WHILE循环中(第8~12行或第17~21行)，为了计算condition的值(第11行或第20行)，都需要计算式(54).式(54)中等号右边第1项的浮点乘法操作数为kn²+nk²，浮点加法操作数为kn²+nk²+k，等号右边第2项的浮点乘法操作数为2pkn+pn²+np²，浮点加法操作数为2pkn +pn² +np² +2pn +p，则式(54)的总操作数为2pn² +2np²+2kn²+2nk²+4pkn+2pn+p+k，因为一般k＜＜min(p，n)，则式(54)的时间复杂度为O(pn(p+n))，记WHILE循环次数为#iterations，则算法2和算法3的时间复杂度均为#iterations×O(pn(p+n)).

算法2  update_G

输入：X，F，G^t，D，L

输出：G^t+1

BEGIN

1) α_pre=1，β=0.1，σ=0.01

2) f(F，G^t)=Tr(G^tL(G^t)^T)+Tr(X-FG^t)D(X-FG^t)^T

3) f(F，G^t+1)=Tr(G^t+1L(G^t+1)^T)+Tr(X-FG^t+1)·D(X-FG^t+1)^T

4) ▽_G^tf(F，G^t)=2G^tL-2F^TXD+2F^TFG^tD

5) G^t+1=P[G^t-α_pre▽_G^tf(F，G^t)]

6) condition=(f(F，G^t+1)-f(F，G^t)≤σ〈▽_G^tf(F，G^t)，(G^t+1-G^t)〉)

7) IF (condition)

8) WHILE (condition)

9) α_new=α_pre/β

10) G^t+1=P[G^t-α_new▽_G^tf(F，G^t)]

11) condition=(f(F，G^t+1)-f(F，G^t)≤σ〈▽_G^tf(F，G^t)，(G^t+1-G^t)〉)

12) END WHILE

13) G^t+1=P[G^t-α_pre▽_G^tf(F，G^t)]

14) return G^t+1

15) END IF

16) ELSE

17) WHILE (! condition)

18) α_new=α_pre×β

19) G^t+1=P[G^t-α_new▽_G^tf(F，G^t)]

20) condition=(f(F，G^t+1)-f(F，G^t)≤σ〈▽_G^tf(F，G^t)，(G^t+1-G^t)〉)

21) END WHILE

22) return G^t+1

23) END ELSE

END

算法3  update_F

输入：X，F^t，G，D，L

输出：F^t+1

BEGIN

1) α_pre=1，β=0.1，σ=0.01

2) f(F^t，G)=Tr(GLG^T)+Tr(X-F^tG)D(X-F^tG)^T

3) f(F^t+1，G)=Tr(GLG^T)+Tr(X-F^t+1G)D(X-F^t+1G)^T

4) ▽_F^tf(F^t，G)=-2XDG^T+2F^tGDG^T

5) F^t+1=P[F^t-α_pre▽_F^tf(F^t，G)]

6) condition=(f(F^t+1，G)-f(F^t，G)≤σ〈▽_F^tf(F^t，G)，(F^t+1-F^t)〉)

7) IF (condition)

8) WHILE (condition)

9) α_new=α_pre/β

10) F^t+1=P[F^t-α_new▽_F^tf(F^t，G)]

11) condition=(f(F^t+1，G)-f(F^t，G)≤σ〈▽_F^tf(F^t，G)，(F^t+1-F^t)〉)

12) END WHILE

13) F^t+1=P[F^t-α_pre▽_F^tf(F^t，G)]

14) return F^t+1

15) END IF

16) ELSE

17) WHILE (! condition)

18) α_new=α_pre×β

19) F^t+1=P[F^t-α_new▽_F^tf(F^t，G)]

20) condition=(f(F^t+1，G)-f(F^t，G)≤σ〈▽_F^tf(F^t，G)，(F^t+1-F^t)〉)

21) END WHILE

22) return F^t+1

23) END ELSE

END

4 实验

本文将RMPGM-NMF算法同k-means算法、PCA算法、标准NMF算法^[9]、L₂₁NMF算法^[10]、rCIM-NMF算法^[11]、LDS-NMF算法^[12]和LP-NMF算法^[14]、SC算法^[19]、GNMF算法^[20]进行比较，给出聚类结果，并分析RMPGM-NMF算法的收敛性和算法的性能与参数的关系.

4.1 数据集

使用了5个图片数据集(来源：http://www.vision.caltech.edu/Image_Datasets/Caltech101和http://www.cad.zju.edu.cn/home/dengcai/Data/data.html)，表 1概括了这5个数据集的基本特性.对于每个数据集，本文随机选取K₁个类别的样本点构成第1个子集，再从第1个子集中随机选取K₂个类别的样本点构成第2个子集，各子集的类别数K₁和K₂见表 2.

ORL：包含40个人的400张人脸图片，它们拍摄于不同的时间、不同的光线条件且同一个人的多张图片中面部表情也不同.每张图片被裁剪为32×32像素.

Umist：包含20个人的360张人脸图片，同一个人的18张图片差别很大，因此在人脸识别时具有挑战性.每张图片的像素为28×23.

COIL20：包含1 440张图片，分别属于20个类别，每个类别有72张图片，每张图片中的物体处于不同的角度.每张图片的分辨率为32×32像素.

Caltech：Caltech101中的图片可以划分为101个大类，其中"人脸"类中包含了435张人脸图片，我们从中选取了19个人的共342张图片，平均每个人18张图片，每张图片被裁剪为39×26像素.

MNIST：该数据集包含手写的数字"0"到"9"的共70 000张图片，从每个类别中随机选取了100张图片，因此实验中使用的MNIST数据集样本数为1 000.每张图片裁剪为28×28像素.

4.2 性能指标

本文使用的性能指标包括聚类准确度ACC(clustering accuracy)，归一化的互信息NMI(normalized mutual information)，纯度PUR(purity).

ACC定义为

(55)

l_i和c_i分别是样本点x_i的真实类标签和聚类得到的类标签，最佳匹配函数map(·)通过匈牙利算法(Hungarian algorithm)求解，函数δ(x，y)定义为

(56)

ACC的值越大表示聚类性能越好.

NMI的定义为

(57)

(58)

其中，C表示样本真实类标签的集合，C'表示聚类得到的类标签的集合，MI(C，C')为C和C'的互信息，其中p(c_i)表示任取一个样本点属于类别c_i的概率，p(c_j')表示任取一个样本点属于类别c_j'的概率，p(c_i，c_j')表示任取一个样本点同时属于类别c_i和c_j'的概率. NMI(C，C')表示C和C'归一化的互信息，其中H(C)和H(C')分别表示C和C'的信息熵.

NMI的值越大表示聚类性能越好.

PUR的定义如下：

(59)

其中，S_i表示聚类后类别标签为i的样本点集合，n_i表示该集合包含的样本点数目，n_i^j表示真实类别标签为j的样本点聚类后标签为i的样本点数目.

PUR的值越大表示聚类性能越好.

4.3 实验结果

实验中，NMF算法、L₂₁NMF算法、rCIM-NMF算法、LDS-NMF算法、GNMF算法和RMPGM-NMF算法的矩阵F和G的初始值均为相同的随机生成的矩阵，迭代次数均为1 000，在最终得到的系数矩阵G上运行k-means算法获得聚类结果. LP-NMF算法的迭代次数为1 000次，在最终得到的矩阵QX上运行k-means算法获得聚类结果.

实验中10个算法均运行5次，并计算5次运行得到的ACC、NMI、PUR的平均值.其中每次调用k-means算法，均执行20次，返回的聚类结果使得各个样本点到其类中心的距离之和最小.

SC算法中的σ参数设置为σ²=0.5×10⁶. LDS-NMF算法有2个可调的参数，即参数λ和γ，实验中设定λ=0.01，γ依次变化为{0.001，0.01，0.1，1，10，100，1 000}；再设定γ=1，λ依次变化为{0.000 01，0.000 1，0.001，0.01，0.1，1，10}. RMPGM-NMF算法和GNMF算法中有2个可调的参数，即参数λ和构造亲和矩阵W使用的最近邻居样本点的数目m.实验中先设定λ=100，m依次变化为{3，5，7，9，11，13，15}；再设定m=5，λ依次变化为{0.001，0.01，0.1，1，10，100，1 000}.不同的参数组合导致不同的聚类结果，选取最好的聚类结果衡量该算法的聚类性能.

10个算法在5个数据集及其子集上的聚类结果如表 3~表 7所示，每个算法的聚类数目r设置为各个数据集真实的类别数，即r=K或r=K₁或r=K₂，各性能指标的最优值用粗体显示，可以看出，RMPGM-NMF算法的性能优于其它9种算法.

4.4 收敛性

分析RMPGM-NMF算法的收敛性，以数据集ORL和COIL20为例，结果如图 1所示，横坐标为迭代的次数，纵坐标为目标函数(29)的以2为底的对数值，算法参数设置为m=5，λ=100，F和G的初始值随机生成. 图 1显示RMPGM-NMF算法可以快速收敛.

4.5 参数m和λ的分析

RMPGM-NMF算法有2个可调的参数m和λ，以数据集ORL和COIL20为例，图 2和图 3分别显示了RMPGM-NMF算法的性能随参数m和λ变化的情况.

图 2中λ=100，m依次变化为{3，5，7，9，11，13，15}，可以看出RMPGM-NMF算法的性能随着m的增大而呈下降的趋势，这是因为本文使用m近邻的亲和矩阵刻画数据集的几何结构，即任一样本点只与其最近的m个样本点建立邻接关系，显然随着m的增大，这样建立的亲和矩阵的不合理性增大.当m取值为3~7时，RMPGM-NMF算法的性能较优.

图 3中m=5，λ依次变化为{0.001，0.01，0.1，1，10，100，1 000}，可以看出对于数据集ORL，RMPGM-NMF算法的性能随着λ的变化基本稳定，对于数据集COIL20，当λ取值过小(小于0.1)时会导致RMPGM-NMF算法的性能降低.总体上，当λ取值为0.1~100时RMPGM-NMF算法的性能受λ的影响不大.

4.6 聚类数目r的分析

在4.3节中，聚类数目r设置为各个数据集真实的类别数，在本节中分析RMPGM-NMF算法的性能随聚类数目r变化的情况，以数据集ORL和COIL20为例，它们的真实类别数K分别为40和20，r依次变化为{2，4，6，8，10，15，20，30，40，50，60}，以PUR为衡量算法性能的指标，图 4显示随着r的变化，除了数据集ORL上r=2时，RMPGM-NMF算法的性能均优于其它算法.

5 结束语

本文提出了RMPGM-NMF算法，它是基于L₂₁范数的鲁棒流形NMF算法，利用L₂₁范数衡量矩阵分解的质量，降低了噪音和异常值对算法聚类性能的影响，同时引入了拉普拉斯矩阵，利用数据的几何结构并考虑局部不变性，使得低维空间G能反映高维空间X中样本的几何结构，进一步提升了算法的聚类性能.本文使用了投影梯度方法求解出基矩阵F和系数矩阵G的更新规则，分析了算法的时间复杂度.在5个数据集及其子集上运行了RMPGM-NMF算法，并与k-means算法、PCA算法、SC算法、标准NMF算法及其它5种NMF相关算法进行了比较，表明了RMPGM-NMF算法的性能优于其它算法，同时分析了RMPGM-NMF算法的收敛性和算法的性能随参数m、λ及聚类数目r变化的情况.如何将RMPGM-NMF算法扩展到半监督聚类是下一步研究的重点.