CN105930308A - 基于低秩恢复的非负矩阵分解方法 - Google Patents

Abstract

本发明属于信息处理技术领域，具体涉及一种基于低秩恢复的非负矩阵分解方法，包括以下步骤：1】将原始数据库中的每个图像样本均转换为向量，构成m×n的原始数据矩阵X；m为图像样本的维数，n为图像样本的个数；2】对原始数据矩阵X进行低秩稀疏分解；2.1】设置低秩矩阵的秩为r，设置稀疏矩阵的稀疏度为k；2.2】采用双边随机投影算法求解原始数据矩阵X的秩为r的低秩矩阵L和稀疏度为k的稀疏矩阵S；3】对步骤2】中求解得到的低秩矩阵L进行非负矩阵分解，得到基矩阵W和编码矩阵H。本发明通过低秩稀疏分解得到数据低秩成分和稀疏成分，并对去除稀疏噪声部分的低秩成分进行非负分解，从而使得非负分解结果免受噪声的干扰。

Description

基于低秩恢复的非负矩阵分解方法

技术领域

本发明属于信息处理技术领域，具体涉及一种基于低秩恢复的非负矩阵分解方法。

背景技术

随着信息化和互联网的发展，高维数据在社会各领域不断涌现。总体来讲，这些数据或者是半结构的或者是无结构的，这使得构建这些数据的特征向量高达上万维甚至更高。数据维数的增加对大规模数据处理带来了困难。非负矩阵分解时基于非监督模式识别的一个研究领域，旨在得到数据稀疏的、基于部分的低维数据表示。非负矩阵分解被广泛应用于数学、最优化、神经计算、模式识别与机器学习、数据挖掘、图像工程与计算机视觉等领域，因此研究非负矩阵分解技术具有非常重要的意义与应用价值。非负矩阵分解技术建立在对因子矩阵的非负约束的基础上，通过非监督的方法学习到原数据的低维表示以及其相应的基矩阵，从而有利于后续相关应用中的处理。

目前，非负矩阵分解方法主要分为三类：

一是基于稀疏约束或正交约束的方法，这种方法致力于通过对因子矩阵(即基矩阵和编码矩阵)施加稀疏或者正交约束来学习到更稀疏和局部化的数据表示。P.Hoyer在文献“Non-Negative Matrix Factorization with Sparseness Constraints,J.MachineLearning Research,vol.5,no.9,pp.1457-1469,2004”提出对编码矩阵施加L1范数约束的非负矩阵分解方法，通过L1范数的稀疏约束来取得更稀疏的数据表示。S.Li等人在文献“Learning spatially localized,parts-based representation,in Proc.IEEEInternational Conference Computer Vision Pattern Recognition,pp.207-212,2001”提出基于正交约束的非负矩阵分解方法，通过正交约束去除了基向量间的的冗余成分。

二是基于判别信息的非负矩阵分解方法。这种方法的核心思想是利用标记样本学习出更具有判别性的低维数据表示。Y.Wang等人在文献“Fisher Non-Negative MatrixFactorization for Learning Local Features,Proc.Asian Conference ComputerVision,pp.27-30,2004”中提出了基于费舍尔判别准则的非负矩阵分解算法，通过引入费舍尔判别准则来使得类内样本分布更紧致，类间样本分布更远。J.Yang等人在文献“Non-negative graph embedding,Proc.IEEE Int’l Conf.Computer Vision and PatternRecognition,pp.1-8,2008”中提出了基于图嵌入的非负矩阵分解方法。该方法构建本质图和惩罚图两个图结构，通过最小化本质图和惩罚图使类内样本分布紧致，类间间隔更大。

三是基于流形学习的非负矩阵分解方法。这种方法通过流行学习算法保持数据在高维空间中的拓扑结构，从而考虑了数据的结构信息。D.Cai等人在文献“Non-negativematrix factorization on manifold,IEEE Transaction Pattern Analysis MachineIntelligence,vol.33,no.8,pp.1548-1560,2011”中提出了基于图正则的非负矩阵分解算法，通过最小化图正则项来保持数据内部的几何分布结构。Q.Gu等人在文献“NeighborhoodPreserving Nonnegative Matrix Factorization,Proc.20th British Machine VisionConf.,2009”中提出基于局部线性嵌入的非负矩阵分解算法，假设数据的局部拓扑结构符合局部线性嵌入假设。这两种方法的不同就在于对于数据的局部拓扑结构的假设。

以上三类方法虽然从不同的角度对非负矩阵分解进行了改进，但是，都没有考虑原数据本身可能包含的噪声对非负分解带来的不良影响。

发明内容

为了解决非负分解容易受到噪声干扰的技术问题，本发明提供一种基于低秩恢复的非负矩阵分解方法。

本发明的技术解决方案是：一种基于低秩恢复的非负矩阵分解方法，其特殊之处在于：包括以下步骤：

1】将原始数据库中的每个图像样本均转换为向量，构成m×n的原始数据矩阵X；m为图像样本的维数，n为图像样本的个数；

2】对原始数据矩阵X进行低秩稀疏分解；

2.1】设置低秩矩阵的秩为r，设置稀疏矩阵的稀疏度为k；

2.2】采用双边随机投影算法求解原始数据矩阵X的秩为r的低秩矩阵L和稀疏度为k的稀疏矩阵S；

3】对步骤2】中求解得到的低秩矩阵L进行非负矩阵分解，得到基矩阵W和编码矩阵H。

上述步骤2.2】包括以下步骤：

2.2.1】初始化低秩矩阵为L₀＝X，初始化稀疏矩阵为S₀＝0，初始化迭代次数为t＝0；设置重构相对误差阈值ε；

2.2.2】计算其中，q为0或者正整数；

2.2.3】计算的r双边随机投影和其中，A₁是n×r的随机矩阵，A₂是m×r的随机矩阵；

2.2.4】进行迭代更新：t＝t+1，A₂＝Y₁，

2.2.5】计算Y₁和Y₂的QR分解：Y₂＝Q₂R₂，Y₁＝Q₁R₁；

2.2.6】更新低秩矩阵和稀疏矩阵：S_t＝P_Ω(X-L_t)；其中，P_Ω表示将一个矩阵投影到元素集合Ω；

2.2.7】计算重构误差判断重构误差是否小于重构相对误差阈值ε，若重构误差小于ε则执行步骤2.2.8】，若重构误差大于或者等于ε则执行步骤2.2.4】；

2.2.8】得到低秩矩阵L＝L_t和稀疏矩阵S＝S_t。

上述步骤3】包括以下步骤：

3.1】将步骤2】中求得的低秩矩阵L中的非负元素赋值为0；

3.2】初始化基矩阵W₀为m×l的随机矩阵，初始化编码矩阵H₀为l×n的随机矩阵，初始化迭代次数为t＝0；其中，l为样本类个数；设定迭代误差限ε′；

3.3】利用K近邻算法构建近邻图，计算图拉普拉斯矩阵L_ap＝D_ap-S_ap；其中，S_ap是对称的权重矩阵，D_ap是对角矩阵，对角元素是S_ap的列和；

3.4】迭代求解基矩阵W和编码矩阵H；

3.4.1】计算其中α是图正则参数；

3.4.2】计算其中，β是Tiknohov正则参数；

3.4.3】计算非负矩阵分解重构误差；若非负矩阵分解重构误差大于或者等于迭代误差限ε′则执行步骤3.4.1】；若非负矩阵分解重构误差小于迭代误差限ε则执行步骤3.4.4】；

3.4.4】得到基矩阵W＝W_t′+1和编码矩阵H＝H_t′+1。

上述基于低秩恢复的非负矩阵分解方法还包括以下步骤：

4】用k-means聚类算法对编码矩阵H进行聚类；

5】计算聚类结果评判指标聚类精度AC和归一化互信息NMI：

$AC=\frac{{\mathrm{\Σ}}_{i=1}^n\mathrm{\δ}({\mathrm{gnd}}_i,map(z_i))}{n}$

$MI(C,C^{\′})=\underset{c_i\∈C,c_j^{\′}\∈C^{\′}}{\Σ}p(c_i,c_j^{\′})\log \frac{p(c_i,c_j^{\′})}{p\left(c_i\right)p\left(c_j^{\′}\right)}$

$NMI(C,C^{\′})=\frac{MI(C,C^{\′})}{max\left(H\right(C),H\left(C^{\′}\right))}$

其中，n是样本集中的样本个数，gnd_i是已知的样本标签；

δ(a,b)是单位冲激函数，当a＝b时δ(a,b)＝1，当a≠b时δ(a,b)＝0；

map(z_i)是重标记函数，z_i表示样本i经过k-means聚类得到的聚类标签；

C是数据库提供的类标签，C′是经过k-means聚类得到的类标签；

p(c_i)表示随机选取的样本x属于C中的c_i类的概率，p(c'_j)表示样本x属于C′中的c_j′类的概率，p(c_i,c'_j)表示样本x同时属于c_i类和c_j′类的概率；

H(C)表示C的熵，H(C′)表示C′的熵。

本发明的有益效果在于：本发明通过原数据矩阵的低秩稀疏分解得到数据低秩成分和稀疏成分，并对去除稀疏噪声部分的低秩成分进行非负分解，从而使得非负分解结果免受噪声的干扰。同时，在矩阵非负分解的同时考虑了数据内的几何分布结构，以取得更好的低维数据表示。最后，对降维的数据表示用k-means算法聚类，通过聚类结果的好坏判断得到的低维数据表示的优劣。

附图说明

图1为本发明较佳实施例的基于低秩恢复的非负矩阵分解流程图。

具体实施方式

参见图1，本发明较佳实施例的实现步骤如下：

步骤1，对原始数据矩阵进行低秩稀疏分解。

(1a)将图像样本集中的每幅图像拉成一个向量，共同构成m×n的原始数据矩阵X，m为每个样本的维数，n为样本个数；

(1b)初始化低秩矩阵为L₀＝X，稀疏矩阵为S₀＝0，低秩矩阵的秩设置为r，稀疏矩阵的稀疏度设置为k，重构相对误差阈值ε，迭代次数t＝0。

(1c)为了防止X的奇异值逐渐退化，从而导致基于双边投影的原数据矩阵X的秩为r的低秩近似矩阵的逼近效果不好，采用幂修正方案，即计算的双边随机投影而不是X的双边随机投影，这样的低秩近似矩阵为其中，是的r双边随机投影，A₁和A₂分别是n×r和m×r的随机矩阵，q可调，q＝0则不采用该方案，这里设置q＝1。迭代更新：

t＝t+1，A₂＝Y₁，

(1d)为了得到原数据矩阵X的秩为r的低秩近似，需要计算Y₁和Y₂的QR分解即：

Y₂＝Q₂R₂，Y₁＝Q₁R₁

然后如下更新低秩矩阵和稀疏矩阵：

S_t＝P_Ω(X-L_t)

其中，P_Ω表示将一个矩阵投影到元素集合Ω，即只保留X-L_t前k个元素。

(1e)计算重构误差：

$\left|\right|X-L_t-S_t|{|}_F^2/|\left|X\right|{|}_F^2$

判断是否小于阈值ε，如果小于ε则迭代停止；否则回到(1c)继续迭代求解。

步骤2，求解低秩恢复矩阵的非负分解。

(2a)对步骤1求得的低秩矩阵L进行非负投影，即将低秩矩阵L中的非负元素赋值为0；

(2b)初始化m×l的基矩阵W₀、l×n的编码矩阵H₀为随机矩阵，l为样本类个数，迭代次数t'＝0；

(2c)利用K近邻算法构建近邻图，并计算图拉普拉斯矩阵L_ap，其中，L_ap＝D_ap-S_ap，S_ap表示对称的权重矩阵，D_ap是对角矩阵，对角元素是S_ap的列和(或者行和，因为S_ap是对称矩阵)；

(2d)为了保障基矩阵W的平滑性，对W施加Tiknohov正则约束：

$\underset{i,j}{\Σ}{W}_{ij}^2=tr\left({\mathrm{WW}}^T\right);$

(2e)迭代求解基矩阵W和编码矩阵H，固定一个求另一个，具体操作如下：

固定W更新H：

固定H更新W：

其中α为图正则参数,β是Tiknohov正则参数，平衡图正则项和Tiknohov正则项间的权重。

(2f)计算非负矩阵分解重构误差，如果小于误差限则停止迭代；否则回到(2e)继续迭代。

步骤3，聚类测试。

将编码矩阵H作为原数据矩阵X的降维表示，用k-means聚类算法对新的样本聚类(编码矩阵H每一列是原数据矩阵X每一列的一个降维表示)；

步骤4，计算聚类结果评判指标聚类精度AC和归一化互信息NMI。

$AC=\frac{{\mathrm{\Σ}}_{i=1}^n\mathrm{\δ}({\mathrm{gnd}}_i,map(z_i))}{n}$

其中，n是样本集中样本个数，gnd_i是事先已知的样本标签，单位冲激函数δ(a,b)＝1，当a＝b；当a≠b其值为0。map(z_i)是重标记函数，可以将算法获得的聚类标签映射到样本集提供的标签上。

$MI(C,C^{\′})=\underset{c_i\∈C,c_j^{\′}\∈C^{\′}}{\Σ}p(c_i,c_j^{\′})\log \frac{p(c_i,c_j^{\′})}{p\left(c_i\right)p\left(c_j^{\′}\right)}$

其中p(c_i)和p(c'_j)分别表示随机选取的样本x属于类C和C'的概率，p(c_i,c'_j)表示样本x同时属于类C和C'的概率。

$NMI(C,C^{\′})=\frac{MI(C,C^{\′})}{max\left(H\right(C),H\left(C^{\′}\right))}$

其中，H(C)和H(C')分别表示类C和C'的熵。NMI度量这两个类别的相似性。为了避免每次聚类结果的随机性，因此统计100次聚类结果中精度排序前10的结果的平均。

本发明的效果可以通过以下实验做进一步的说明。

1.仿真条件

本发明是在中央处理器为Intel(R)Core i3-2130 3.4GHZ、内存16G、WINDOWS 8操作***上，运用MATLAB软件进行的仿真。

实验中使用的图像数据库为美国邮政***手写数字字符数据库USPS和剑桥大学人脸图像数据库ORL。USPS数字字符数据库包括“0”到“9”10个数字的9298张灰度图像，每张图像大小为16×16，实验中每类选择100个样本，共1000个数据样本；ORL人脸数据库包括40个人的面部灰度图像，每个人有10个图像样本，共400张图像，每张图像是大小是32×32，每个人的图像样本具有不同的光照变化、表情变化以及面部细节。

2.仿真内容

首先，在USPS数据库和ORL数据库上，完成本发明算法(基于低秩恢复的非负矩阵分解)的实验。为了证明算法的有效性，我们选取了4个对比方法NMF、GNMF、k-means、PCA进行比较。其中NMF是在文献“D.D.Lee and H.S.Seung,Learning the parts of objects bynonnegative matrix factorization,Nature,vol.401,no.6755,pp.788–791,1999”中提出的。GNMF在“D.Cai,X.He,J.Han,and T.S.Huang,Graph regularized nonnegativematrix factorization for data representation,IEEE Trans.Pattern Analysis andMachine Intelligence,vol.33,no.8,pp.1548–1560,2011”中有详细介绍。PCA是被广泛应用的数据降维算法。K-means是广泛应用的聚类算法。

在USPS数据库中随机选取8类、9类、10类数字字符的图像分别实验；在ORL数据库中随机选取30类、35类、40类人脸图像分别实验。用k-means算法对NMF、PCA、GNMF以及我们提出的算法学***均值。

在本发明实施例中对于USPS数据库和ORL数据库分别设置低秩矩阵的秩r为原数据矩阵满秩的9％和10％，即min(m,n)×9％和min(m,n)×10％，稀疏矩阵的稀疏度k为原数据矩阵规模的1％，即m×n×1％，对于USPS数据库和ORL数据库l分别取值8、9或10和30、35或40，图正则项的正则参数α＝100，Tiknohov正则参数β＝0.0001，构建图的近邻个数为5。

实验测试结果如表1和表2所示。

表1 USPS数据库上的聚类结果

表2 ORL数据库上的聚类结果

从表1，2可见，本发明的聚类效果比PCA、NMF和GNMF三种数据降维表示方法都要好，同时比不经降维表示直接用k-means聚类的效果好。

Claims (4)

1.一种基于低秩恢复的非负矩阵分解方法，其特征在于：包括以下步骤：

1】将原始数据库中的每个图像样本均转换为向量，构成m×n的原始数据矩阵X；m为图像样本的维数，n为图像样本的个数；

2】对原始数据矩阵X进行低秩稀疏分解；

2.1】设置低秩矩阵的秩为r，设置稀疏矩阵的稀疏度为k；

2.2】采用双边随机投影算法求解原始数据矩阵X的秩为r的低秩矩阵L和稀疏度为k的稀疏矩阵S；

3】对步骤2】中求解得到的低秩矩阵L进行非负矩阵分解，得到基矩阵W和编码矩阵H。

2.根据权利要求1所述的基于低秩恢复的非负矩阵分解方法，其特征在于：所述步骤2.2】包括以下步骤：

2.2.1】初始化低秩矩阵为L₀＝X，初始化稀疏矩阵为S₀＝0，初始化迭代次数为t＝0；设置重构相对误差阈值ε；

2.2.2】计算其中，q为0或者正整数；

2.2.3】计算的r双边随机投影和其中，A₁是n×r的随机矩阵，A₂是m×r的随机矩阵；

2.2.4】进行迭代更新：t＝t+1，A₂＝Y₁，

2.2.5】计算Y₁和Y₂的QR分解：Y₂＝Q₂R₂，Y₁＝Q₁R₁；

2.2.6】更新低秩矩阵和稀疏矩阵：S_t＝P_Ω(X-L_t)；其中，P_Ω表示将一个矩阵投影到元素集合Ω；

2.2.7】计算重构误差判断重构误差是否小于重构相对误差阈值ε，若重构误差小于ε则执行步骤2.2.8】，若重构误差大于或者等于ε则执行步骤2.2.4】；

2.2.8】得到低秩矩阵L＝L_t和稀疏矩阵S＝S_t。

3.根据权利要求2所述的基于低秩恢复的非负矩阵分解方法，其特征在于：所述步骤3】包括以下步骤：

3.1】将步骤2】中求得的低秩矩阵L中的非负元素赋值为0；

3.2】初始化基矩阵W₀为m×l的随机矩阵，初始化编码矩阵H₀为l×n的随机矩阵，初始化迭代次数为t＝0；其中，l为样本类个数；设定迭代误差限ε′；

3.3】利用K近邻算法构建近邻图，计算图拉普拉斯矩阵L_ap＝D_ap-S_ap；其中，S_ap是对称的权重矩阵，D_ap是对角矩阵，对角元素是S_ap的列和；

3.4】迭代求解基矩阵W和编码矩阵H；

3.4.1】计算其中α是图正则参数；

3.4.2】计算t'＝t'+1；其中，β是Tiknohov正则参数；

3.4.3】计算非负矩阵分解重构误差；若非负矩阵分解重构误差大于或者等于迭代误差限ε′则执行步骤3.4.1】；若非负矩阵分解重构误差小于迭代误差限ε则执行步骤3.4.4】；

3.4.4】得到基矩阵W＝W_t′+1和编码矩阵H＝H_t′+1。

4.根据权利要求1-3中任一所述的基于低秩恢复的非负矩阵分解方法，其特征在于：还包括以下步骤：

4】用k-means聚类算法对编码矩阵H进行聚类；

5】计算聚类结果评判指标聚类精度AC和归一化互信息NMI：

$AC=\frac{{\mathrm{\Σ}}_{i=1}^n\mathrm{\δ}({\mathrm{gnd}}_i,map(z_i\left)\right)}{n}$

$MI(C,C^{\′})=\underset{c_i\∈C,c_j^{\′}\∈C^{\′}}{\mathrm{\Σ}}p(c_i,c_j^{\′})\log \frac{p(c_i,c_j^{\′})}{p\left(c_i\right)p\left(c_j^{\′}\right)}$

$NMI(C,C^{\′})=\frac{MI(C,C^{\′})}{max\left(H\right(C),H(C^{\′}\left)\right)}$

其中，n是样本集中的样本个数，gnd_i是已知的样本标签；

δ(a,b)是单位冲激函数，当a＝b时δ(a,b)＝1，当a≠b时δ(a,b)＝0；

map(z_i)是重标记函数，z_i表示样本i经过k-means聚类得到的聚类标签；

C是数据库提供的类标签，C′是经过k-means聚类得到的类标签；

p(c_i)表示随机选取的样本x属于C中的c_i类的概率，p(c'_j)表示样本x属于C′中的c_j′类的概率，p(c_i,c'_j)表示样本x同时属于c_i类和c_j′类的概率；

H(C)表示C的熵，H(C')表示C′的熵。