牛頓-軟閾值迭代魯棒主成分分析算法

王海鵬，降愛蓮，李鵬翔

（太原理工大學信息與計算機學院，山西晉中 030600）

（?通信作者電子郵箱ailianjiang@126.com）

0 引言

主成分分析（Principal Component Analysis，PCA）作為一種常用的數據降維方法，在數據分析、目標檢測、圖像識別等領域得到了廣泛應用。PCA算法能夠提取目標數據的主要特征達到降維的效果［1］。如果數據含有稀疏噪聲，PCA 提取的特征子空間會向噪聲數據方向偏移。魯棒主成分分析（Robust Principal Component Analysis，RPCA）是一種能夠對具有結構信息但是被稀疏噪聲破壞的數據矩陣進行分解的算法，找到隱藏在原始數據中的低秩成分，實際上就是找到了數據的本質低維空間。RPCA 不僅在特征提取中有著抵抗噪聲的特性，而且在視頻的背景前景分離［2-4］、魯棒人臉識別［5］、移動目標檢測［6］等領域也具有廣闊的應用前景。

2011 年，Candès 等［7］首次定義了RPCA 問題，即將一個含有稀疏噪聲的數據矩陣分解為低秩矩陣與稀疏矩陣兩部分。為了保證稀疏矩陣足夠稀疏，使用l0-范數進行約束，同時要使低秩矩陣的秩盡可能小。然而Gillis等［8］證明二進制矩陣（0-1矩陣）的分解是NP-Hard 問題，進而證明l0-范數約束最優(yōu)化是一個NP-Hard 問題。由于l0-范數無法使用凸優(yōu)化解決，文獻［7］將RPCA 問題的低秩與稀疏約束松弛為核范數與l1-范數約束的目標函數最優(yōu)化問題，并用拉格朗日方向交替法［9］求解。由于該算法每次迭代需要對數據進行奇異值分解求核范數，因此該算法的時間復雜度非常高［4］。Shen 等［10］提出了基于增廣拉格朗日方向交替法的低秩矩陣擬合（Low-Rank Matrix Fitting，LMaFit）算法，該算法用兩個低秩矩陣的乘積表示一個低秩矩陣的方法代替核范數，降低了因使用奇異值分解而產生的時間開銷。這種方法是后續(xù)研究中降低時間開銷的重要方法。Yi 等［11］提出基于梯度下降（Gradient Descent，GD）的RPCA 算法，該算法使用稀疏矩陣估計算子估計出原始數據的稀疏矩陣的近似矩陣。GD 算法在優(yōu)化算法上選擇了梯度下降，提高了算法的運行效率。Lu 等［12］推導出了張量核范數、張量軟閾值算子，并以此構建數學模型，采用方向交替拉格朗日乘子法實現了張量魯棒主成分分析（Tensor Robust Principal Component Analysis，TRPCA）算法，該算法提高了RPCA 算法對彩色圖像處理、視頻背景建模的能力。TRPCA 算法因為使用了張量奇異值分解，因此時間復雜度較高。Chereau 等［13］提出了基于最大相關熵模型的RPCA 算法，他們采用了冪迭代的方式來求解最大相關熵模型，進而實現RPCA 算法。該算法的優(yōu)點是不需要預先設置低秩矩陣的主成分的數量，提高了算法的易用性。由于上述RPCA 算法的計算時間隨著數據規(guī)模提高而顯著提升，處理大規(guī)模數據會消耗大量的時間，因此本文將針對RPCA 問題，通過改進RPCA 模型，降低RPCA 算法的時間復雜度，進而降低數據規(guī)模對RPCA算法的影響，提高RPCA算法的運行速度。

1 RPCA問題及模型

由于主成分分析算法是由歐氏距離最小化推導出的算法，因此稀疏噪聲會因為歐氏距離的二次方特性對主成分分析算法求得結果產生較大的影響［14］。

圖1（a）展示了PCA 算法在數據較均勻的條件下的效果，從圖中可以看出，PCA 能夠準確地對這些分布均勻的數據降維。如果數據含有稀疏噪聲，如圖1（b），PCA 算法的結果會產生較大的偏差。RPCA 算法提取了含有稀疏噪聲數據的主要特征，并且對該含噪數據進行了降維。

圖1 PCA算法與RPCA算法對比Fig.1 Comparison of PCA algorithm and RPCA algorithm

為了將含噪數據分解為低秩矩陣和稀疏矩陣，同時使低秩矩陣的秩與稀疏矩陣的非0 元素的數量盡可能少，文獻［7］中將RPCA模型表示為：

其中：Y表示包含稀疏噪聲的數據，L表示數據中的低秩部分，S 表示數據中的稀疏部分，rank(L)表示L 的秩，模型（1）滿足Y=L+S。模型（1）理論上可以達到分離低秩矩陣和稀疏矩陣的目的。然而l0-范數是一種無法使用凸優(yōu)化求解的范數。為了優(yōu)化模型（1），因此使用核范數與l1-范數將模型（1）松弛化為：

2 NSTI模型及其優(yōu)化算法分析

本文改進原有RPCA 模型，提出了一個新的RPCA 模型，使得該模型在能夠表示RPCA 問題的基礎上又能夠采用簡單高效的方法進行優(yōu)化，該優(yōu)化方法稱為牛頓-軟閾值迭代（Newton-Soft Threshold Iteration，NSTI）算法。

2.1 NSTI算法模型

為了降低因計算奇異值分解而導致的巨大的算法運行時間的開銷，LMaFit和GD算法使用了如式（3）所示的方法：

其中U、V 是兩個低秩矩陣，滿足L=UVT。鑒于這種利用兩個低秩矩陣的乘積表示一個低秩矩陣的方法，可避免使用核范數，本文改進模型（2），并給出NSTI算法的數學模型：

式（4）左邊的范數中，低秩矩陣L由兩個低秩矩陣U、V 的乘積表示，因為rank(L)=min(rank(U)，rank(V))，所以能夠保證L 是一個低秩矩陣。式（4）左邊是一個能夠對U、V 進行凸優(yōu)化求解的范數，本文用牛頓法來優(yōu)化左邊的范數。采用牛頓法求解的原因有兩個：第一，牛頓法比梯度下降收斂速度更快；第二，不會因為使用了牛頓法而增加過多的時間開銷。

式（4）右邊稀疏矩陣S 通過l1-范數來表示其稀疏性質。l1-范數是l0-范數的凸近似。最小化l1-范數的目的是為了讓S稀疏化。本文采用軟閾值迭代算法求解稀疏矩陣S。該方案相對于文獻［11］中的GD 算法的稀疏矩陣估計算子，具有計算速度快的優(yōu)點。軟閾值迭代算法的缺點是軟閾值是靜態(tài)的，并且需要在算法運行之前根據數據調整軟閾值的大小。本文提出了軟閾值估計算子，用來彌補軟閾值迭代算法中的不足。

2.2 NSTI算法

NSTI算法分為兩部分，其中一部分是求解稀疏矩陣S。S對于模型（4）的梯度為：

由軟閾值迭代算法［15］，可求出稀疏矩陣S 的近似解。式（6）是軟閾值算子也是S近似解的計算式。

NSTI 算法的另一部分是求解低秩矩陣的兩個因子U、V。U、V對于模型（4）的一階梯度為：

Hessian 矩陣是牛頓法求解的關鍵，通過計算出Hessian矩陣，就能利用牛頓法求出低秩矩陣U、V。U、V 矩陣對應的Hessian矩陣為：

因此，可得到如下兩式：

NSTI算法詳細過程如下。

輸入 數據Y，精確度α，收斂閾值t，秩r；

輸出 U、V、S。

1）首先對Y 進行奇異值分解，選取前r個奇異值與前r個左奇異向量的乘積初始化U，選取前r 個奇異值與前r 個右奇異向量的乘積初始化V。稀疏矩陣S為0矩陣。

2）計算Y'=Y -UVT，并計算軟閾值λ=f(Y'，α)。

3）利用軟閾值迭代求出稀疏矩陣S=soft(Y'，λ)。

4）計算矩陣E=Y'-S，以及U、V 的Hessian 矩陣的逆

NSTI 算法思想借鑒了方向交替法。將模型（4）中要求的3 個變量分為兩部分，通過交替迭代求解模型中的3 個變量。2.1 節(jié)提到牛頓法比梯度下降收斂速度更快，軟閾值迭能夠快速求解稀疏部分，因此NSTI采用牛頓法與軟閾值迭代交替計算的方案來優(yōu)化模型（4）。該方案能夠求解模型（4）中的3個變量并且算法的收斂速度與計算速度都得到了提升。

2.3 軟閾值估計

軟閾值λ 決定著數據Y 中的數據的劃分：λ 值越低，每次迭代，能提取更多的元素成為稀疏矩陣，然而這樣會使得原本屬于低秩矩陣的成分被劃分到稀疏矩陣中，優(yōu)點是能夠提高算法的精確度；反之，λ 值越高，將提高低秩矩陣元素的比例，優(yōu)點是能夠降低算法迭代次數。如果λ 值過高，那么最終收斂的結果會將全部元素劃分到低秩矩陣L中，稀疏矩陣S將會是一個0矩陣；反之，λ值過低，那么全部元素將劃分到稀疏矩陣S中，L矩陣將變成一個0矩陣。軟閾值計算方式為：

對任意屬于Y'(Y'=Y -UVT)的元素a1、a2，通過和閾值λ 比較，小于閾值λ 的則被歸到低秩成分L 中，大于λ 的被歸到稀疏成分S中。如果|a1|＜|a2|，而且一定存在低秩成分，那么a1比a2更有可能是低秩成分，低秩空間UVT的元素被認為是其對應原始數據Y 的元素在低秩空間的投影，原始數據減去其在低秩空間大的投影的絕對值是原始數據元素到低秩子空間的距離。將數據擴展到Y 全部元素，對于任意ai(i=1，2，…，mn -1，mn)，一定存在這樣的一個序列，其中ai與bj是一一對應關系。通過式（13）估計一個λ，將數據中距離平面較近的一部分劃分為一個低秩矩陣，將距離平面較遠的數據部分劃分為一個稀疏矩陣。通過調節(jié)α 的大小控制λ 的大小，進而控制算法的精確度。

2.4 NSTI時間復雜度分析

假設數據Y 是一個m × n 的數據矩陣，Y=L+S，L 是一個m × n 的矩陣，L 矩陣的秩為r，r 遠小于m 和n。假設矩陣U、V 的規(guī)模分別為m × r，n × r，那么L 的秩最大為r。對于稀疏矩陣S，S的規(guī)模為m × n。

基于上述假設，奇異值分解的時間復雜度為O(n2m)，初始化U、V 矩陣的過程為O(mr)，因此初始化過程時間復雜度為O(n2m)。

在算法迭代過程中：第2）、3）步時間復雜度為O(nm)；第4）步求解Hessian 矩陣時，需要計算的時間復雜度為O(mnr)，在求Hessian 矩陣的逆矩陣中，因為Hessian 矩陣的規(guī)模時r ×r，所以Hessian 矩陣的逆運算開銷可以忽略；同理第6）步中的時間復雜度為O(mnr)。因此每次迭代復雜度為O(mnr)，軟閾值估計算子的時間復雜度為O(nm)。假設算法迭代x次收斂，則算法完整的時間復雜度為O(xnmr)，因為r遠小于m 和n，因此算法時間復雜度近似為O(xnm)。

3 實驗與分析

3.1 實驗環(huán)境

實驗選取NSTI 算法、LMaFit 算法和GD 算法三個算法進行對比。LMaFit、GD 是文獻［16-17］評價認可的效率較高的兩個算法。通過與LMaFit算法和GD算法對比來證明NSTI算法時間效率和算法精確度的提高。實驗環(huán)境為Matlab 2019a，實驗機器配置為Intel Core i7 8700，6 核心12 線程，內存為8 GB，操作系統(tǒng)為Windows 10。在實驗中的奇異值分解函數為PROPACK 數據分析工具包中的lansvd 函數?；趌anczos分解的lansvd 奇異值分解函數具有比傳統(tǒng)奇異值分解函數更快的計算速度。

3.2 算法有效性分析

為了驗證算法是有效的，將從算法的整體損失、低秩矩陣損失和稀疏矩陣損失三個角度來說明。算法整體損失為l，l越小說明算法越精確，提取的稀疏矩陣越稀疏。其中l(wèi)為：

決定數據Y 規(guī)模的參數有m、n 和低秩矩陣的秩r 以及稀疏矩陣的非0元素的比例a，其中數據Y=L+S。低秩矩陣L中的元素生成過程如下：首先生成一組基向量，選擇r 個基向量，隨機組合生成一個m × n 的秩為r 的低秩矩陣L。然后利用a 隨機生成一個稀疏矩陣S，該實驗使用Matlab 內置函數sprand。

實驗數據的規(guī)模設為1 000×1 000，低秩矩陣的秩為20，稀疏矩陣的稀疏度為0.2。圖2（a）顯示隨著迭代次數增加，算法整體與各部分逐漸收斂。NSTI 算法借鑒了方向交替法的思想，每次迭代求出一組參數，作為下次迭代的輸入。算法每次迭代，都能夠降低算法的整體損失與各參數的損失，最終達到算法的整體收斂。圖2（b）中驗證了算法能夠正確地求出一個符合輸入數據的稀疏度的稀疏矩陣。算法在迭代初期，計算出的稀疏矩陣包含大量的非0 元素，但隨著算法迭代，稀疏矩陣收斂于原始的稀疏矩陣。

圖2 NSTI算法運行效果Fig.2 Running effect of NSTI algorithm

3.3 算法性能對比

為了驗證NSTI 算法的時間效率的提升。設計實驗對比三種算法在同一數據下經過迭代達到相同的損失l 所需要的運行時間。

圖3 是三種算法在相同規(guī)模的數據下運行時間的對比，說明了NSTI 算法比LMaFit 算法、GD 算法，運行速度要快，運行時間穩(wěn)定。

圖3 三種算法在不同數據規(guī)模下的對比Fig.3 Comparison of three algorithms under different data sizes

LMaFit 的單次迭代時間復雜度為O(n3)，GD 算法的單次迭代時間復雜度為O(n3)，NSTI 算法單次迭代時間復雜度為O(n2)，所以NSTI的計算時間相對LMaFit、GD增長緩慢。在算法機理上，NSTI 選用的牛頓法，在理論上比梯度下降收斂速度快，又因為RPCA 問題的低秩的特性，牛頓法的計算速度由O(n3)變?yōu)镺(n2)，所以牛頓法比梯度下降法計算速度快。NSTI選擇的軟閾值迭代法理論上時間復雜度為O(n2)，在保證算法準確度的同時，保證了算法的時間效率。

圖4 是當非0 元素比例a=0.1，數據規(guī)模為5 000×5 000，通過調整低秩矩陣的秩r，對數據進行矩陣分解，對比三個算法在達到相同的損失l所需要的時間：當r ＜8時，NSTI算法與GD 算法運行時間比較穩(wěn)定；當r ≥8時，三種算法的變化規(guī)律逐漸接近。

由于算法機理的不同，NSTI、GD 與LMaFit 單次迭代的時間隨著r提升而小幅提升，然而LMaFit的迭代次數隨著r的提升而降低，進而降低了計算時間。如圖4，NSTI 與GD 算法略有提升，LMaFit算法有明顯下降。

圖4 三種算法在不同低秩r下運行時間對比Fig.4 Comparison of running time of three algorithms under different low rank r conditions

圖5 展示了NSTI 對比其他算法的時間效率提升百分比，圖5的數據來源是圖4。通過圖5可以看出，NSTI算法在秩為1 時相比較GD、LMaFit 算法，能夠分別提升54.0%與93.5%的時間效率；當秩為20 時，能夠分別提升24.6%與45.5%的時間效率。

圖5 與對比算法相比NSTI的性能提升Fig.5 Performance improvement of NSTI compared with comparison algorithm

3.4 視頻背景/前景提取實驗

為了驗證NSTI算法在視頻前景背景分離的能力，設計實驗對比算法在視頻背景/前景提取中的效果。將視頻看做數據Y，視頻的背景看作是低秩矩陣L，視頻的前景看作是稀疏矩陣S［15］。對視頻的數據矩陣進行RPCA，得到低秩矩陣和稀疏矩陣。低秩矩陣輸出為背景，稀疏矩陣輸出為前景。圖6為該實驗要處理視頻的第40、70、100幀。

圖6 原始視頻Fig.6 Original video

圖7 展示不同算法提取的背景。對比提取的背景，在第70幀和第100幀中原畫汽車的位置，GD算法與LMaFit算法出現黑影，在NSTI 算法對應的位置沒有黑影，色彩與原始圖像的色彩基本一致，這說明NSTI 算法提高了去除噪聲的能力。在右側欄桿位置，LMaFit 算法將欄桿全部劃分到前景，NSTI算法，GD算法將欄桿算作背景。

圖8 中的GD 和LMaFit 算法出現一個NSTI 算法中沒有出現的圖塊。對比9 張?zhí)崛〉那熬皥D片，NSTI 算法相對GD、LMaFit 算法，提取的視頻前景最完整。表1 說明NSTI 算法比其他算法的運行時間效率至少提高了78.7%。

RPCA 算法能夠將數據分解為一個低秩矩陣和一個稀疏矩陣，進而將視頻分解為背景與前景。又因為該視頻數據的背景的秩為1，因此能夠保證NSTI 算法的牛頓法部分的時間效率最高。NSTI 采用了軟閾值迭代法，同時用軟閾值估計算子確定軟閾值，這種改進了的軟閾值迭代算法能夠求出更加精確的稀疏矩陣，因此可以從數據中看出NSTI算法求出的前景更加精確。

圖7 不同算法提取的背景對比Fig.7 Comparison of background extracted by different algorithms

圖8 不同算法提取的前景對比Fig.8 Comparison of foreground extracted by different algorithms

表1 各算法消耗的時間 單位：sTab.1 Time consumed by different algorithms unit：s

3.5 圖像降噪實驗

為了證明算法在圖像降噪的能力的提升，本實驗對一組a=0.2含噪圖片進行降噪。實驗數據共有10張含噪圖片。圖像降噪實驗中，原始圖像被視為低秩矩陣L，噪聲被視作稀疏矩陣S，通過RPCA 處理達到降噪效果。本實驗將對比算法的運行時間以及算法的殘差e，利用殘差對比算法的精確度高低，e越小精確度越高。e的計算方式為：

其中：L表示圖片的數據矩陣，Lsrc表示原始圖片，Li表示第i幅含噪圖片降噪后的結果。該實驗采用峰值信噪比（Peak Signal-to-Noise Ratio，PSNR）［12］來評價算法的降噪能力。

圖9 的數據說明了NSTI 算法的輸出結果噪點最少，相比GD和LMaFit，NSTI的圖像更加清晰。表2的數據說明了NSTI的殘差值比其他算法至少降低了45.3%。NSTI算法的時間效率比其他兩個算法至少高出64.3%。NSTI算法計算得到的圖片的PSNR 值比LMaFit 算法高出18.8%，比GD 算法高出24.8%。因此可以說明NSTI算法處理的圖像更加清晰。

表2 圖像降噪實驗Tab.2 Image denoising experiment

圖9 圖像降噪實驗Fig.9 Image denoising experiment

因為NSTI算法提高了算法效率，所以在處理圖像數據更加高效，NSTI 算法比LMaFit 算法提高的時間效率為78.29%，與圖5的數據相符合。NSTI算法采用了改進的軟閾值迭代算法，不需要預估稀疏度，相比較稀疏矩陣估計算子，提高了圖像降噪的水平，PSNR高達31.337 5 dB。

3.6 圖像恢復實驗

為了驗證NSTI 算法在圖像恢復方面的能力設計人臉恢復實驗，選取的對比算法有GD、LMaFit、TRPCA 算法，數據集為耶魯大學的人臉數據集。耶魯大學的人臉數據包含面部的不同表情、光影、方向及佩飾的圖像。

圖10 4種算法的圖像恢復效果對比Fig.10 Image recovery effects comparison among four algorithms

表3展示了各個算法分別處理得到的圖像的PSNR 的值，PSNR 值越大，說明算法的恢復圖像能力越強。圖10 驗證了NSTI算法能夠有效地恢復圖像。與LMAFit相比，NSTI的時間效率與PSNR 值比較高，這一點在表3 中得到驗證。NSTI，GD與TRPCA算法的PSNR值在表3中無明顯區(qū)別，這說明對于耶魯大學數據集，NSTI、GD 與TRPCA 三種算法的圖像恢復能力基本一致。因為NSTI算法使用了牛頓法與軟閾值迭代法，所以在時間效率上，NSTI比另外三種算法要高，能夠快速地計算出人臉的低秩特征。表3證明了NSTI有著最高的時間效率。

NSTI 算法比其他算法的時間效率高的原因有兩點：第一，GD 算法每次迭代需要求解稀疏矩陣估計算子，該算子的時間復雜度為O(n3)；LMaFit算法每次迭代都需要計算奇異值分解，奇異值分解的時間復雜度為O(n3)；TRPCA 算法存在張量核范數的運算，張量核范數的時間復雜度O(n3)。第二，NSTI采用了牛頓法與軟閾值迭代法提高計算速度，NSTI每次迭代時間復雜度為O(n2)。

表3 不同算法圖像恢復性能對比Tab.3 Comparison of image restoration performance of different algorithms

4 結語

本文研究了RPCA 問題，利用牛頓法和軟閾值迭代法改進了RPCA模型，并提出了NSTI算法。經實驗證明，本文提出的NSTI算法，與其他RPCA算法相比，在時間效率與算法精確度上有明顯提升。本文提出的軟閾值估計算子能夠計算出一個動態(tài)且精確的軟閾值，能夠提高NSTI算法的收斂速度和精確度。NSTI算法也存在著不足，像GD 算法一樣，參數中需要輸入低秩r，來決定低秩矩陣的秩r，而沒辦法在算法中直接估計出r的大小。