基于隨機(jī)方差調(diào)整梯度的非負(fù)矩陣分解

史加榮, 白姍姍

(西安建筑科技大學(xué) 理學(xué)院, 西安710055)

非負(fù)矩陣分解(non-negative matrix factorization, NMF)是尋找非負(fù)的基矩陣和系數(shù)矩陣, 即將非負(fù)數(shù)據(jù)表示為非負(fù)基下的非負(fù)線性組合, 以實(shí)現(xiàn)維數(shù)約簡(jiǎn)和特征提取[1]. 實(shí)際應(yīng)用中很多物理特性都是非負(fù)的, 非負(fù)約束使模型具有較強(qiáng)的可解釋性[2]. 近年來, NMF受到廣泛關(guān)注, 已被成功應(yīng)用到數(shù)據(jù)挖掘[3]、 計(jì)算機(jī)視覺[4-5]和模式識(shí)別[6]等領(lǐng)域. 乘性更新(multiplicative update, MU)[7]是求解NMF的一種有效方法, 可保證其解收斂于一個(gè)平穩(wěn)點(diǎn). 但MU計(jì)算復(fù)雜度較高, 故不適用于求解大規(guī)模矩陣分解問題. 隨著數(shù)據(jù)規(guī)模的不斷擴(kuò)大, 隨機(jī)梯度下降(stochastic gradient descent, SGD)算法[8]已成為機(jī)器學(xué)習(xí), 特別是深度學(xué)習(xí)研究的熱點(diǎn)[9-10], 其具有參數(shù)更新過程簡(jiǎn)單、 收斂速度快且計(jì)算復(fù)雜度低等特點(diǎn). 因此, 將SGD算法與MU規(guī)則相結(jié)合是求解NMF的一種有效方法.

Kasai[11]將SGD算法應(yīng)用到NMF的求解中, 提出了隨機(jī)乘性更新(stochastic multiplicative update, SMU)和隨機(jī)方差縮減乘性更新(stochastic variance reduced multiplicative update, SVRMU)規(guī)則, 分別將SGD算法和隨機(jī)方差縮減梯度(stochastic variance reduction gradient, SVRG)[12]算法與乘性更新規(guī)則相結(jié)合, 對(duì)梯度估計(jì)量進(jìn)行更新; 之后，Kasai[13]又提出了隨機(jī)平均梯度乘性更新(stochastic average gradient multiplicative update, SAGMU)規(guī)則. 雖然這些算法取得了較好的效果, 但都未考慮偏差與方差之間的不平衡性差異, 即未對(duì)梯度估計(jì)量進(jìn)行校正, 導(dǎo)致迭代效率較低. 文獻(xiàn)[14]提出了一種性能優(yōu)于SVRG的算法： 隨機(jī)方差調(diào)整梯度(stochastic variance adjusted gradient, SVAG)算法. 該算法在迭代過程中引入?yún)?shù)對(duì)梯度估計(jì)量進(jìn)行不斷調(diào)整, 使偏差與方差達(dá)到平衡. 基于此, 本文提出一種將SVAG算法和MU規(guī)則相結(jié)合的隨機(jī)乘性更新算法： 隨機(jī)方差參數(shù)調(diào)整梯度乘性更新(stochastic variance parameter adjusted gradient multiplicative update, SVPAGMU)算法. 該算法結(jié)合方差縮減策略和乘性更新規(guī)則, 在梯度向量更新的同時(shí), 引入一個(gè)參數(shù)調(diào)整梯度估計(jì)量, 不斷校正梯度下降方向使其偏差與方差達(dá)到平衡, 最終快速逼近局部最優(yōu)解.

1 預(yù)備知識(shí)

1.1 非負(fù)矩陣分解及乘性更新

(1)

其中‖·‖F(xiàn)為矩陣的Frobenius范數(shù),W為基向量組成的矩陣,H為系數(shù)矩陣. 由于優(yōu)化模型(1)是非凸的, 因此求解該模型的全局最優(yōu)解非常困難, 塊坐標(biāo)下降算法[15]是解決該問題的一種主要方法.

針對(duì)最優(yōu)化模型(1), 文獻(xiàn)[7]提出了一種簡(jiǎn)單而有效的乘性更新規(guī)則：

H←H⊙((WTV)./(WTWH)),

(2)

W←W⊙((VHT)./(WHHT)),

(3)

其中⊙表示Hadamard積, 即對(duì)應(yīng)元素之積, ./表示對(duì)應(yīng)元素之商. MU規(guī)則通過交替更新W和H, 使模型(1)中目標(biāo)函數(shù)在每步迭代中都減小, 最終達(dá)到平穩(wěn)點(diǎn).

1.2 隨機(jī)梯度下降法的代表性算法

許多機(jī)器學(xué)習(xí)任務(wù)都可轉(zhuǎn)化為經(jīng)驗(yàn)風(fēng)險(xiǎn)最小化模型：

(4)

其中x∈d為模型參數(shù),fi(x):d→是第i個(gè)樣本關(guān)于參數(shù)x的損失函數(shù). SGD算法是求解大規(guī)模學(xué)習(xí)問題的一種有效方法, 下面簡(jiǎn)單介紹幾種具有代表性的隨機(jī)梯度下降算法, 更多算法可參見文獻(xiàn)[16].

SGD算法以其基本的梯度下降形式更新迭代, 在每輪更新參數(shù)時(shí), 僅隨機(jī)選取一個(gè)或幾個(gè)樣本計(jì)算其梯度, 并以此梯度作為全局梯度的估計(jì)值, 極大降低了計(jì)算復(fù)雜度. SGD算法的參數(shù)更新公式為

xt+1=xt-ηfi(xt),

(5)

其中參數(shù)η表示迭代步長(或?qū)W習(xí)率),i∈[N]([N]={1,2,…,N})表示第t輪迭代中隨機(jī)選取的樣本序號(hào),fi(xt)表示當(dāng)前迭代梯度. 由于梯度方差隨迭代次數(shù)的增加而不斷累加, 因此導(dǎo)致SGD算法收斂速度較慢.

為平衡有偏低方差的SAG算法和無偏高方差的SAGA算法, 文獻(xiàn)[14]引入了SVAG算法, 該算法引入一個(gè)參數(shù)權(quán)值θ對(duì)梯度估計(jì)量進(jìn)行校正, 即需判斷隨機(jī)選取的單個(gè)采樣梯度fit(xt)與該樣例處的最新梯度fit(xtit)之間的相關(guān)程度, 不斷調(diào)整梯度使其偏差與方差達(dá)到平衡, 參數(shù)更新公式為

(6)

其中ti表示第t輪迭代中隨機(jī)選取的樣本序號(hào)i. 特別地,θ=1和θ=N分別對(duì)應(yīng)于SAG算法和SAGA算法. 數(shù)值算例驗(yàn)證了SVAG算法優(yōu)于SVRG算法及其他隨機(jī)梯度下降算法.

2 求解NMF的隨機(jī)方差調(diào)整梯度算法

設(shè)vn和hn分別表示矩陣V和H的第n列, 且n∈[N]是選取的樣本序號(hào). 最優(yōu)化模型(1)等價(jià)于

(7)

圖1 SVPAGMU算法的基本流程Fig.1 Basic flow chart of SVPAGMU algorithm

(8)

(9)

在更新W前, 為使其偏差與方差之間達(dá)到平衡, 校正梯度下降方向, 需判斷隨機(jī)采樣梯度g1與當(dāng)前快照點(diǎn)處梯度g2的相關(guān)程度, 記為

(10)

(11)

若‖g1-g2‖2<ε, 則選擇較大的參數(shù)θ； 否則, 選擇較小的參數(shù)θ. 通過θ取值的不同調(diào)整梯度估計(jì)量, 使其沿著最佳的梯度下降方向更新W. 根據(jù)SVAG算法可知： 參數(shù)θ的選取與樣本數(shù)目N有關(guān), 根據(jù)不同的數(shù)據(jù), 判斷梯度g1與g2的相關(guān)關(guān)系, 選取合適的參數(shù)θ對(duì)梯度估計(jì)量進(jìn)行不斷調(diào)整.

(13)

(14)

(15)

隨著迭代次數(shù)的增加,αs不斷縮減, 即

αs+1=ραs,

(16)

其中ρ∈(0,1). 綜上可得：

算法1隨機(jī)方差調(diào)整梯度乘性更新(SVPAGMU)算法.

步驟2) Fors=0,1,…, do

步驟4) Fort=1,2,…,T

步驟6) 根據(jù)式(9)更新hm；

步驟7) 根據(jù)式(10)計(jì)算g1；

步驟8) 根據(jù)式(11)計(jì)算g2；

步驟9) if ‖g1-g2‖2<ε；θ=θmax； elseθ=θmin； end

步驟14) End For

步驟15) 根據(jù)式(16)更新αs；

步驟17) End For.

下面討論SVPAGMU算法在一次外部迭代過程的計(jì)算復(fù)雜度. 若內(nèi)層循環(huán)最大迭代次數(shù)為T, 則更新系數(shù)矩陣H的計(jì)算復(fù)雜度為O(TMK2), 計(jì)算梯度g1和g2的復(fù)雜度均為O(TMK), 計(jì)算Q和P的復(fù)雜度均為O(TMNK), 計(jì)算S的復(fù)雜度為O(TM2K2), 最終得到更新基矩陣W的計(jì)算復(fù)雜度為O(T(M2K2+MNK)). 因此, 一次迭代總的計(jì)算復(fù)雜度為O(T(M2K2+MNK)).

3 數(shù)值實(shí)驗(yàn)

下面在人臉圖像集和氣象數(shù)據(jù)集上進(jìn)行實(shí)驗(yàn), 并將SVPAGMU算法與SMU,SVRMU和SAGMU三種算法進(jìn)行比較. 使用MATLAB R2019a進(jìn)行編程, 涉及的全部算法程序均在IntelCoreTMi7-8565U CPU @1.80 GHz 1.99 GHz處理器上運(yùn)行, 操作系統(tǒng)為64位Windows 10系統(tǒng).

3.1 實(shí)驗(yàn)設(shè)置

(17)

相對(duì)誤差越小, 逼近性能越好.

3.2 人臉圖像集

通過對(duì)人臉圖像集進(jìn)行測(cè)試, 比較上述4種算法的實(shí)際性能. 實(shí)驗(yàn)選擇ORL數(shù)據(jù)集, 該數(shù)據(jù)集來自于UCI機(jī)器學(xué)習(xí)庫, 是由40個(gè)不同年齡、 不同性別的人在不同時(shí)間獲取的共400幅灰度圖像組成, 即每個(gè)人的人臉圖像10幅, 每個(gè)像素的灰度大小為0～255, 每幅圖像的分辨率為64×64.

圖2 不同算法對(duì)人臉圖像集的迭代效率(A)和時(shí)間效率(B)對(duì)比Fig.2 Comparison of iterative efficiency (A) and time efficiency (B) of different algorithms on face image set

下面考慮K取不同值時(shí)SVPAGMU算法迭代效率和時(shí)間效率的比較. 選取K∈{20,30,40,50}, 實(shí)驗(yàn)結(jié)果如圖3所示. 由圖3(A)可見, 相對(duì)誤差隨迭代次數(shù)的增加而逐漸改善, 并且在同一迭代次數(shù)中K值越大, 相對(duì)誤差越小; 由圖3(B)可見, 隨著K值的增加, 其運(yùn)行時(shí)間逐漸延長. 因此, 若將SVPAGMU算法應(yīng)用于較高維數(shù)據(jù), 需選取合適的K值, 使其在迭代效率和時(shí)間效率上均達(dá)到相對(duì)最優(yōu).

圖3 不同K值下SVPAGMU算法對(duì)人臉圖像集的迭代效率(A)和時(shí)間效率(B)對(duì)比Fig.3 Comparison of iterative efficiency (A) and time efficiency (B) of SVPAGMU algorithm on face image set under different K values

為更直觀地說明SVPAGMU算法引入?yún)?shù)θ的有效性, 圖4展示了目標(biāo)函數(shù)的震蕩軌跡, 其中橫坐標(biāo)表示迭代次數(shù), 縱坐標(biāo)表示當(dāng)前目標(biāo)函數(shù)值與最優(yōu)目標(biāo)函數(shù)值之差, 圖中實(shí)心圓表示快照點(diǎn). 由圖4可見, SVPAGMU算法目標(biāo)函數(shù)的震蕩范圍隨著迭代次數(shù)的增加而不斷減小, 并且震蕩幅度較小, 同時(shí)快照點(diǎn)的變化軌跡一直保持平緩下降, 算法優(yōu)化效率較高. 這是因?yàn)镾VPAGMU算法在梯度更新時(shí), 引入了調(diào)整梯度估計(jì)量的參數(shù)θ, 使其方差與偏差達(dá)到平衡, 不斷地校正梯度下降方向, 沿著最優(yōu)軌跡逼近局部最優(yōu)解, 從而提升了算法性能.

圖4 SVPAGMU算法目標(biāo)函數(shù)的震蕩軌跡Fig.4 Oscillation trajectory of objective function for SVPAGMU algorithm

3.3 氣象數(shù)據(jù)集

選取中國661個(gè)氣象臺(tái)站2004—2013年共10年的相對(duì)濕度數(shù)據(jù)[19]. 為避免氣象數(shù)據(jù)存在異常波動(dòng)的影響, 使用每個(gè)臺(tái)站10年的日平均數(shù)據(jù). 因此, 對(duì)氣象變量中的相對(duì)濕度, 每個(gè)臺(tái)站的觀測(cè)值由一個(gè)d=365維列向量表示, 則所有數(shù)據(jù)可表示為365×661維的非負(fù)矩陣V. 比較SMU,SVRMU,SAGMU和SVPAGMU四種算法的性能, 分別選K=15,K=30進(jìn)行對(duì)比, 實(shí)驗(yàn)結(jié)果如圖5和圖6所示.

圖5 不同算法在不同K值下對(duì)平均相對(duì)濕度的迭代效率對(duì)比Fig.5 Comparison of iterative efficiency of different algorithms on average relative humidity under different K values

由圖5可見, 各算法的實(shí)際性能趨勢(shì)與人臉圖像集下的性能基本保持一致. SVPAGMU算法性能明顯優(yōu)于其他3種對(duì)比算法, 相對(duì)誤差隨迭代次數(shù)的增加而迅速減小. 當(dāng)K=30時(shí), SVPAGMU算法相對(duì)誤差已達(dá)0.112 5, 說明該算法矩陣分解結(jié)果更準(zhǔn)確.

由圖6可見, SAGMU算法所耗費(fèi)的時(shí)間最長, 可能是由于添加了梯度平均, 增加了計(jì)算復(fù)雜度. 隨著K值的增加, 每個(gè)算法的相對(duì)誤差都在減少, 但運(yùn)行時(shí)間也逐漸延長. 以SVPAGMU算法為例, 當(dāng)K=15和K=30時(shí), 相對(duì)誤差分別為0.137 9,0.127 5, 運(yùn)行時(shí)間分別為1.63,2.47 s. 因此, 針對(duì)不同的數(shù)據(jù)需選取合適的K值, 使其在迭代效率和時(shí)間效率上都能夠達(dá)到相對(duì)最優(yōu).

圖6 不同算法在不同K值下對(duì)平均相對(duì)濕度的時(shí)間效率對(duì)比Fig.6 Comparison of time efficiency of different algorithms on average relative humidity under different K values

為檢驗(yàn)矩陣分解的逼近性能, 先對(duì)平均相對(duì)濕度數(shù)據(jù)V應(yīng)用SVPAGMU算法進(jìn)行非負(fù)矩陣分解, 再將其分解得到的矩陣進(jìn)行重構(gòu). 通過實(shí)驗(yàn)對(duì)比恢復(fù)結(jié)果與真實(shí)數(shù)據(jù)的差異, 參數(shù)設(shè)置如下： 最大迭代次數(shù)為50,K=30, 其他參數(shù)與之前設(shè)置保持一致. 部分實(shí)驗(yàn)結(jié)果如圖7所示. 圖7中隨機(jī)選取的4個(gè)氣象臺(tái)站(臺(tái)站號(hào))分別為： 內(nèi)蒙古朱日和鎮(zhèn)(53276)、 西藏類烏齊縣(56128)、 浙江麗水市(58646)和新疆巴音布魯克(51542). 由圖7可見, 重構(gòu)數(shù)據(jù)與原始數(shù)據(jù)非常接近, 且重構(gòu)數(shù)據(jù)振幅小, 這可能是因?yàn)镹MF有效地去除了噪聲. 此外, 4個(gè)氣象臺(tái)站的平均相對(duì)誤差分別為0.084 3,0.077 2,0.089 7,0.082 3, 表明該算法逼近性能較好.

圖7 部分氣象臺(tái)站相對(duì)濕度的原始數(shù)據(jù)與重構(gòu)數(shù)據(jù)對(duì)比Fig.7 Comparison between original data and reconstructed data of relative humidity on some meteorological stations

綜上所述, 針對(duì)非負(fù)矩陣分解問題, 本文提出了一種隨機(jī)方差參數(shù)調(diào)整梯度乘性更新算法. 該算法將SVAG算法與MU規(guī)則相結(jié)合, 通過引入對(duì)梯度估計(jì)量進(jìn)行不斷調(diào)整的參數(shù)θ, 使其偏差與方差達(dá)到平衡, 從而有效提高了矩陣分解效果. 選取已有的3種算法SMU,SVRMU,SAGMU與本文SVPAGMU算法進(jìn)行實(shí)驗(yàn)比較, 同時(shí)選取不同的K值進(jìn)行分析, 結(jié)果驗(yàn)證了本文算法的可行性與高效性.