宋 珊，馮 巖，徐常青

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

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

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

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

1 NMF模型及算法

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

1.1 代價(jià)函數(shù)

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

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

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

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

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

1.2 迭代算法

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

1.2.1 乘性迭代算法

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

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

其中

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

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

算法1 NMF-MUEuc算法

輸入：Y，r，ε，maxiter；

輸出：A，X；

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

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

則進(jìn)入步驟3；

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

迭代后進(jìn)入步驟2；

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

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

并用矩陣形式寫成算法。

算法2 NMF-MUKL算法

輸入：Y，r，ε，maxiter；

輸出：A，X；

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

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

進(jìn)入步驟3；

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

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

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

1.2.2 交替最小二乘算法

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

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

算法3 NMF-ALS算法

輸入：Y，r，ε，maxiter；

輸出：A，X；

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

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

則進(jìn)入步驟3；

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

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

1.3 算法比較

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

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

2 NTF模型及算法

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

2.1 代價(jià)函數(shù)

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

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

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

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

2.2 迭代算法

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

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

所以

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

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

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

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

算法4 NTF-MUEuc算法

輸入：X，r，ε，maxiter；

輸出：U，V，W；

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

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

步驟4，否則進(jìn)入步驟3；

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

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

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

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

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

式寫成算法。

算法5 NTF-MUKL算法

輸入：X，r，ε，maxiter；

輸出：U，V，W；

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

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

步驟4，否則進(jìn)入步驟3；

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

迭代后進(jìn)入步驟2；

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

3 實(shí)驗(yàn)?zāi)M

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

3.1 實(shí)驗(yàn)條件和參數(shù)設(shè)定

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

3.2 實(shí)驗(yàn)步驟

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

3.3 結(jié)果分析

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

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

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

4 結(jié)語

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