非負矩陣分解與非負張量分解：算法與應用

宋 珊，馮 巖，徐常青

（蘇州科技大學 數(shù)學科學學院，江蘇 蘇州 215009）

科學研究和實際問題中的大量觀測數(shù)據(jù)，如距離、體積、壓強、降雨量、溶液濃度和網(wǎng)絡可靠性等度量數(shù)據(jù)都具有非負性，這些數(shù)據(jù)集合經(jīng)適當處理后可表示為非負陣列（包括標量和向量）形式，如信息論中的熵（entropy），索引圖像或視頻流，物理中的矩、速度場和量子場，概率統(tǒng)計學中的概率，網(wǎng)絡中的連通度等[1-2]。這些帶非負特征的數(shù)、向量或矩陣甚至高階張量的每個元素因其非負性而具有可解釋性（物理意義）。但是，人們所熟知的矩陣分解，如矩陣滿秩分解、奇異值分解（SVD）和QR分解等，得到的因子矩陣一般不具備非負性，這一方面限制了其物理應用（分解得到的因子矩陣不可解釋），另一方面也會增大計算誤差和分解難度。

非負矩陣分解（Nonnegative Matrix Factorization，簡稱NMF）產(chǎn)生于1994年，最早由Paatero和Tapper[3]在研究元素非負的正定矩陣分解時提出。1999年，Lee和Seung[4-5]將NMF方法運用于人臉識別技術(shù)，并于2000年提出乘性迭代算法。隨后，在乘性迭代算法的基礎(chǔ)上產(chǎn)生了很多推廣的NMF高效算法，如Cichocki等人[6-7]在2006年提出NMF的推廣智能化算法，并在2007年提出層次交替最小二乘算法，而Lin[8]在2007年提出NMF的投影梯度算法。

非負張量分解（Nonnegative Tensor Factorization，簡稱NTF）最早在2005年由Shashua和Hazan[9]提出，他們是在研究如何提取圖片局部特征時建立了三階非負張量分解模型，并用非負梯度下降法求解該模型。2009年，Zafeiriou[10]提出基于KL散度和廣義Bregman距離的NTF算法，并給出收斂性證明。2011年，劉昶等人[11]提出張量的正交非負CP分解算法并用于圖像表示和識別。2016年，梁秋霞等人[12]提出基于NTF的人臉識別算法。實際情況下，對規(guī)模較大、階數(shù)高（例如10階以上）的非負張量分解，目前并沒有可行的分解算法。2014年，祁力群、徐常青和徐毅通過定義張量的分層對角占優(yōu)給出了一類對稱非負張量的非負分解，并給出了較好的算法[13]。

1 NMF模型及算法

給定觀測值矩陣Y∈R+n×m，其中R+n×m表示由所有元素非負的n×m的矩陣構(gòu)成的集合。Y對應的非負矩陣分解（NMF）Y≈AX等價于非負線性約束下的線性模型

1.1 代價函數(shù)

為了求解因子矩陣A和X，首先需要衡量AX對Y的逼近程度，即誤差矩陣E的大小。Lee和Seung在文獻[5]中選取歐氏距離和KL散度構(gòu)造了兩種代價函數(shù)。

基于歐式距離的代價函數(shù)為

此時因子矩陣求解問題轉(zhuǎn)化為最優(yōu)化問題

基于KL散度的代價函數(shù)為

此時模型問題（1）轉(zhuǎn)化為最優(yōu)化問題

1.2 迭代算法

針對最優(yōu)化問題（3）和（5），常用的迭代算法有：乘性迭代算法（MU）和交替最小二乘算法（ALS）。

1.2.1 乘性迭代算法

Lee和Seung提出的乘性迭代算法是在梯度下降法基礎(chǔ)上賦予了特定步長得到的。該迭代算法具有收斂速度較快和可操作性強等優(yōu)點，且只要保證A和X的初值非負，那么每一次迭代得到的結(jié)果都是非負的。文中以最優(yōu)化問題（3）為例，具體說明基于歐氏距離的乘性迭代算法的求解過程。

設(shè)A ik和X kj的梯度下降算法迭代規(guī)則為

其中

將（7）式代入（6）式，得

（9）式即為基于歐氏距離的乘性迭代規(guī)則（簡稱MUEuc）。在下面的算法中，運用矩陣形式來表述該算法的迭代規(guī)則。先給出算法中相關(guān)符號的說明：*表示兩個大小相同的矩陣按元素對應相乘，/表示兩個大小相同的矩陣按元素對應相除。

算法1 NMF-MUEuc算法

輸入：Y，r，ε，maxiter；

輸出：A，X；

步驟1：隨機初始化矩陣A，X，確保非負，初始化計數(shù)器k=0；

步驟2：計數(shù)器自增k=k+1，并求解（2）式的值DEuc（Y||AX），若DEuc（Y||AX）＜ε或k＞maxiter，則進入步驟4，否

則進入步驟3；

步驟3：A和X按（9）式的矩陣形式進行迭代

迭代后進入步驟2；

步驟4：迭代結(jié)束，返回A，X。

同理，可以得到基于KL散度的乘性迭代規(guī)則（簡稱MUKL）

并用矩陣形式寫成算法。

算法2 NMF-MUKL算法

輸入：Y，r，ε，maxiter；

輸出：A，X；

步驟1：隨機初始化矩陣A，X，確保非負，初始化計數(shù)器k=0；

步驟2：計數(shù)器自增k=k+1，并求解（4）式的值DKL（Y||AX），若DKL（Y||AX）＜ε或k＞maxiter，則進入步驟4，否則

進入步驟3；

步驟3：A和X按（11）式的矩陣形式進行迭代

這里1表示n×m的全1矩陣。迭代后進入步驟2；

步驟4：迭代結(jié)束，返回A，X。

1.2.2 交替最小二乘算法

交替最小二乘算法（Alternative Least Square，縮寫為ALS）主要用于求解最優(yōu)化問題（3）。它的基本思想是：先固定A并通過最小化目標函數(shù)（2）求非負因子矩陣X，再代入計算得到的X，并通過最小化目標函數(shù)（2）進一步更新非負因子矩陣A，循環(huán)往復，直至代價函數(shù)（2）產(chǎn)生的兩個矩陣序列收斂為止。根據(jù)最小二乘法，易知當A固定時等價于矩陣方程ATAX=ATY在已知A，Y的情形下估計X。同樣固定X時等價于矩陣方程XXTAT=XYT在已知X，Y的情形下估計A。如此便可求得交替最小二乘算法的迭代公式

下面將上述流程寫成具體算法：

算法3 NMF-ALS算法

輸入：Y，r，ε，maxiter；

輸出：A，X；

步驟1：隨機初始化矩陣A，X，確保非負，初始化計數(shù)器k=0；

步驟2：計數(shù)器自增k=k+1，并求解（2）式的值DEuc（Y||AX），若DEuc（Y||AX）＜ε或k＞maxiter，則進入步驟4，否

則進入步驟3；

步驟3：X和A按（13）式進行迭代，迭代后進入步驟2；

步驟4：迭代結(jié)束，返回A，X。

1.3 算法比較

下面通過數(shù)值實驗比較基于歐氏距離的乘性迭代算法（算法1）和交替最小二乘算法（算法3）的優(yōu)劣。隨機生成一個500×100的矩陣，元素介于0～1之間，取特征維數(shù)r=5，用算法1和算法3對生成矩陣進行分解，畫出代價函數(shù)（2）在不同迭代次數(shù)和運行時間下的示意圖（圖1（a）和圖1（b））。從圖中可以看出，最開始的時候（Epoch＜10，Time＜0.01），在相同的迭代次數(shù)和運行時間下，NMF-MUEuc算法得到的代價函數(shù)值比NMF-ALS算法小，隨著迭代次數(shù)和運行時間的增長，NMF-MUEuc算法得到的代價函數(shù)值又比NMF-ALS算法大，最后兩者的代價函數(shù)值趨于一致。這說明NMF-MUEuc算法開始的時候收斂速度比NMF-ALS算法快，之后比NMF-ALS算法慢，但兩者最終都收斂到大致相同的平穩(wěn)點。

圖1 不同迭代次數(shù)和運行時間下代價函數(shù)的變化

2 NTF模型及算法

NMF模型（1）可推廣到非負張量分解（NTF）模型，該文著重于三階非負張量的分解。給定找到非負因子矩陣使得

2.1 代價函數(shù)

與NMF情形相類似，為了求解因子矩陣U，V，W，需要先衡量[U，V，W]對X的逼近程度，即E的大小。Shashua和Hazan在文獻[9]中同樣選取歐氏距離和KL散度構(gòu)造了兩種代價函數(shù)?；跉W式距離的代價函數(shù)為

此時因子矩陣求解問題轉(zhuǎn)化為最優(yōu)化問題

基于KL散度的代價函數(shù)為

此時因子矩陣求解問題轉(zhuǎn)化為最優(yōu)化問題

2.2 迭代算法

針對最優(yōu)化問題（16）和（18），文中介紹最常用的乘性迭代算法，它是其他所有推廣算法的基礎(chǔ)。與NMF的乘性迭代算法類似，NTF的乘性迭代算法同樣是在梯度下降法的基礎(chǔ)上賦予了特定步長得到的。下面以最優(yōu)化問題（16）為例，具體說明基于歐氏距離的乘性迭代算法的求解過程。

設(shè)uij的梯度下降迭代規(guī)則為

所以

其中e i表示I階單位矩陣的第i個列向量。將（21）式代入（19）式得

同理，可得vsj和wtj的乘性迭代公式

（23）、（24）和（25）構(gòu)成NTF的基于歐式距離的乘性迭代規(guī)則。在此規(guī)則下，只要U，V，W的初值非負，那么每次迭代得到的結(jié)果都是非負的。Shashua和Hazan在文獻[9]中給出了算法的收斂性證明。在下面的算法中，將上述迭代規(guī)則運用矩陣的Khatri-Rao積和張量的矩陣化進行改寫。矩陣A＝[a1，a2，…，an]∈Rm×n和B＝[b1，b2，…，bn]∈Rp×n的Khatri-Rao積為

其中?為Kronecker乘積，A?B表示A的每個元素乘以B。三階張量沿模-1，模-2，模-3方向的矩陣化結(jié)果記為X（1），X（2），X（3），其中X的元素xi1i2i3與矩陣X（k），k=1，2，3的元素xik j是一一對應的，并且有

算法4 NTF-MUEuc算法

輸入：X，r，ε，maxiter；

輸出：U，V，W；

步驟1：隨機初始化矩陣U，V，W，確保非負，初始化計數(shù)器k=0；

步驟2：計數(shù)器自增k=k+1，計算（15）式的值DEuc（X||[U，V，W]），若DEuc（X||[U，V，W]）＜ε或k＞maxiter，則進入

步驟4，否則進入步驟3；

步驟3：U，V，W按（23）、（24）、（25）式的矩陣形式進行迭代

其中A=W⊙V，B=W⊙U，C=V⊙U。迭代后進入步驟2；

步驟4：迭代結(jié)束，返回U，V，W。

同樣，可以得到NTF的基于KL散度的乘性迭代規(guī)則

該乘性迭代規(guī)則的收斂性由Zafeiriou在文獻[10]中用構(gòu)造輔助函數(shù)的思想證明。下面將迭代規(guī)則用矩陣形

式寫成算法。

算法5 NTF-MUKL算法

輸入：X，r，ε，maxiter；

輸出：U，V，W；

步驟1：隨機初始化矩陣U，V，W，確保非負，初始化計數(shù)器k=0；

步驟2：計數(shù)器自增k=k+1，計算（17）式的值DKL（X||[U，V，W]），若DKL（X||[U，V，W]）＜ε或k＞maxiter，則進入

步驟4，否則進入步驟3；

步驟3：U，V，W按（29）、（30）、（31）式的矩陣形式進行迭代

迭代后進入步驟2；

步驟4：迭代結(jié)束，返回U，V，W。

3 實驗模擬

下面將非負矩陣分解和非負張量分解的基于歐式距離的乘性迭代算法（算法1和算法4）用于人臉識別的特征提取，通過識別準確率比較它們之間的優(yōu)劣。

3.1 實驗條件和參數(shù)設(shè)定

文中實驗環(huán)境為PC機，CPU2.10 GHz，Windows10，Matlab（2014a）。實驗數(shù)據(jù)來自AR數(shù)據(jù)庫和ORL數(shù)據(jù)庫，從AR數(shù)據(jù)庫中下載100個人的像素為165×120的面部灰度圖像，共2 600張，其中每個人有26種不同的姿態(tài)（可看成一類）。從ORL數(shù)據(jù)庫中下載40個人的像素為112×92的面部灰度圖像，共400張，其中每個人有10種不同的姿態(tài)（可看成一類）。在實驗開始前，為了使實驗數(shù)據(jù)滿足文中算法的要求，對面部灰度圖像進行預處理。首先，把面部圖像裁剪成60×60像素的圖像，再把灰度值介于0～255之間的uint8類型的數(shù)據(jù)轉(zhuǎn)換為灰度值介于0～1之間的雙精度（2double）數(shù)據(jù)。實驗過程中所有的實驗參數(shù)都視情況而定。

3.2 實驗步驟

對AR數(shù)據(jù)庫作兩組實驗：實驗一是從每類中隨機選取5張圖像作為訓練樣本，其余為測試樣本；實驗二是從每類中隨機選取10張圖像作為訓練樣本，其余為測試樣本。對ORL數(shù)據(jù)庫也作兩組實驗：實驗一是從每類中隨機選取4張圖像作為訓練樣本，其余為測試樣本；實驗二是從每類中隨機選取6張圖像作為訓練樣本，其余為測試樣本。最后通過最近鄰分類器實現(xiàn)人臉識別。

3.3 結(jié)果分析

表1和表2分別是AR數(shù)據(jù)庫和ORL數(shù)據(jù)庫的識別率比較表。從表中可以看出，不管是AR數(shù)據(jù)庫還是ORL數(shù)據(jù)庫，隨著特征維數(shù)和訓練樣本數(shù)的增加，NMF-MUEuc算法和NTF-MUEuc算法的識別率都在逐步提高，但當特征維數(shù)和訓練樣本數(shù)相同時，NTF-MUEuc算法的識別率總比NMF-MUEuc算法的識別率高，說明NTF-MUEuc算法對人臉特征提取的更好。

表1 AR數(shù)據(jù)庫識別率比較表

表2 ORL數(shù)據(jù)庫識別率比較表

4 結(jié)語

首先，介紹NMF模型在不同代價函數(shù)下的乘性迭代算法和交替最小二乘算法，并通過實驗比較兩種算法的不同；其次，介紹NTF模型在不同代價函數(shù)下的乘性迭代算法；最后，用具體的人臉識別實驗說明NTF算法比NMF算法在人臉特征提取上更有優(yōu)勢。