Hankel矩陣填充的保結(jié)構(gòu)閾值修正算法

張春曉，宋儒瑛

(太原師范學(xué)院 數(shù)學(xué)系,山西 晉中 030619)

0 引言

矩陣填充是近些年以來非常熱的一個研究課題,就是如何在不完備的數(shù)據(jù)下把缺少的數(shù)據(jù)補充完整.它的應(yīng)用相當(dāng)廣泛,比如有圖像修復(fù)、協(xié)同過濾等等.在這里主要研究圖像修復(fù)問題.圖像修復(fù)簡單來說就是通過矩陣填充模型將“打碼”的圖片修復(fù)成原來的圖片.通常使用的是由Candès和Recht[1]首先提出的凸優(yōu)化問題模型

(1)

簡單了解一下這個模型.‖X‖*是一個核范數(shù)，是所求矩陣X∈Rn1×n2的奇異值之和,是rank(X)的最優(yōu)凸近似.{Mij:(i,j)∈Ω} 是秩為r的采樣方陣M∈Rn1×n2里隨機已知m個元素的集合,Ω是已知元素的下標(biāo)集合.

現(xiàn)在用更通俗易懂的語言來描述一下研究內(nèi)容:在現(xiàn)實生活中的大規(guī)模數(shù)據(jù)常常會有部分?jǐn)?shù)據(jù)缺失、數(shù)據(jù)誤差、損壞等問題,這將進一步加大數(shù)據(jù)處理和分析難度.這在實際生活中很常見.例如在人臉識別中,人的臉部在識別時會受到來自外界光照或是別的不可控因素的影響,導(dǎo)致識別到的人臉會有陰影、反光、扭曲等;在運動恢復(fù)結(jié)構(gòu)問題中,提取和匹配特征值的時候,經(jīng)常會存在較大的誤差,這便會使得一些常規(guī)的分析方法和處理手段失效,所以需要提供一種更有效和更實用的算法能夠為人臉識別技術(shù)提供強有力的理論支撐.并且,如果能夠使得一些損毀、殘缺的數(shù)據(jù)得到有效恢復(fù),如果能夠以正確的方法使數(shù)據(jù)變得完整,這些將會對大數(shù)據(jù)的建立、對數(shù)據(jù)的綜合分析和處理產(chǎn)生更大效用.再比如一個“打碼”的圖片(就是失真的圖片)用數(shù)學(xué)語言表示就是一個包含很多0元素的矩陣,每個0元素都是攜帶信息的，我們怎么樣把0元素包含的信息挖掘出來呢?這就需要用到低秩了,這叫低秩矩陣重構(gòu)問題.用一個簡單的模型表述一下:已知矩陣是一個給定的M×N的矩陣A,其中一些元素因為某些原因丟失了,如果沒有其他參考條件,我們確定這些數(shù)據(jù)很困難,但是我們已知A的秩是rank(A)?M且rank(A)?N,那么我們就可以通過矩陣各行(列)之間的線性相關(guān)將丟失的元素求出.

目前已經(jīng)出現(xiàn)了大量的求解此問題的快速算法.比較著名的奇異值閾值算法(簡稱為SVT[2])、増廣Lagrange 乘子算法(簡稱為ALM[3])、不需要奇異值分解的快速奇異值閾值法(簡稱為FastSVT[4])和加速奇異值閾值法(簡稱為ASVT[5])、正交秩1矩陣逼近法[6]、梯度投影算法[7]、交替最速下降算法(ASD[8])等等.

首先給出一個n×n的Hankel矩陣H∈Rn×n:

這個矩陣是由H=VDVT生成的,秩為r.這里V是一個n×r的范德蒙德矩陣,D是一個r×r的主對角矩陣.決定Hankel 矩陣的共有2n-1個元素,就是它的第1列和末1行.填充 Hankel 矩陣這個問題的研究價值極大,是一代代數(shù)學(xué)人薪火相傳不斷攻克的難關(guān)之一.之前Qiao等人演算出針對Hankel 矩陣的快速奇異值分解算法復(fù)雜度較小僅為O(n2logn),而一般矩陣的奇異值分解算法復(fù)雜度就較大為O(n3),兩者之間的差距一目了然.Qiao等人這一算法的提出使整個數(shù)學(xué)界掀起了研究Hankel 矩陣的熱潮.因為通過研究這一算法,大大減少了奇異值分解的時間,實現(xiàn)了很大的突破.眾所周知,奇異值分解消耗的時間在矩陣的整體填充中擁有很大的占比.

填充方法的指向性不強是之前所有Hankel矩陣填充算法的關(guān)鍵所在,因此在該篇文章中主要解決這個問題并且提出了解決此問題的修正算法.這一修正算法吸收前人研究中的精華,并進行了嚴(yán)格的數(shù)值實驗對比,最終證明了修正算法是卓有成效的.

本文結(jié)構(gòu)安排:第2節(jié)簡單回顧前人提出的基本算法,通過這些算法,Hankel矩陣問題將得到有效解決;第3節(jié)細(xì)致地介紹我們提出的l步修正的關(guān)于Hankel 矩陣的填充算法,并對新算法進行實驗驗證;第4節(jié)是對全文的總結(jié).

下面是一些相關(guān)的定義:

定義1設(shè)任意矩陣X=(xij)∈Rn×n,E(X) 定義如下:

定義2(奇異值分解(SVD)) 矩陣X∈Rn1×n2,秩為r,一定存在某個正交矩陣U∈Rn1×r和V∈Rn2×r,使得:

X=U∑rVT,∑r=diag(σ1,…,σr),

其中:σ1≥σ2≥…≥σr>0.

定義3(奇異值閾值算子) 對于任意參數(shù)τ≥0，秩為r的矩陣X∈Rn1×n2，存在奇異值分解X=U∑rVT，奇異值閾值算子Dτ定義為:

Dτ(X):=UDτ(∑)VT,Dτ(∑)=diag({σi-τ}+)

表示一個n×n矩陣,其中l(wèi)=-n+1,…,n-1.PΩ是集合Ω上的正交投影滿足:

1 算法回顧

出于文章的完整性考慮,先做一些簡單回顧.

下面的優(yōu)化模型為低秩狀態(tài)下Hankel矩陣的填充問題:

(2)

X、M均為Hankel矩陣,并且M是低秩，Ω?{-n+1,…,n-1}

Hankel 矩陣是擁有特殊結(jié)構(gòu)的,我們在對其進行填充時,使得每次迭代得到的填充矩陣都保持了Hankel 結(jié)構(gòu),而利用好快速奇異值分解就是利用Lanczos方法和快速傅里葉變換技術(shù).

算法1矩陣填充的奇異值閾值(SVT)算法

奇異值閾值法是如下的優(yōu)化模型:

(3)

其算法如下:

第一步: 給定Ω為下標(biāo)集合,PΩ(M)為樣本矩陣,參數(shù)τ,步長δ，誤差ε,以及初始矩陣Y0=k0δPΩ(M),k:=0.

第二步: 計算矩陣YK比閾值τ大的SVD

[Uk,∑k,Vk]=lansvd(Yk)

令Xk+1=UkDτk(∑k)VkT.

第三步: 若‖PΩ(Xk+1-M)‖F(xiàn)/‖PΩ(M)‖F(xiàn)≤ε,停止; 否則，轉(zhuǎn)第四步.

第四步:Yk+1=PΩ(Yk)+δPΩ(M-Xk+1);

算法2矩陣填充的增廣拉格朗日乘子(ALM)算法

增廣拉格朗日乘子法是如下的優(yōu)化模型:

(4)

其拉格朗日函數(shù)為：

其中Y∈Rn1×n2,其算法如下:

第一步: 給定參數(shù)μ0>0,p>1誤差ε1和ε2,給定樣本矩陣PΩ(M),初始矩陣Y0=0,E0=0,k:=0.

第二步: 計算矩陣YK比閾值μk-1大的SVD

[Uk,∑k,Vk]μk-1=svd(D-Ek+μk-1Yk)

第三步:令

Ak+1=UkDμk-1(∑k)VkT

第四步: 若‖D-Ak+1-Ek+1‖F(xiàn)/‖D‖F(xiàn)<ε1,且μk‖Ek+1-Ek‖F(xiàn)/‖D‖F(xiàn)<ε2停止,否則轉(zhuǎn)第五步;

第五步: 給定參數(shù),令Yk+1=YK+μk(D-Ak+1-Ek+1),如果μk‖Ek+1-Ek‖F(xiàn)/‖D‖F(xiàn)<ε2,

令μk+1=ρμk;否則,應(yīng)該轉(zhuǎn)第二步.

算法3

基于F-模的Hankel 矩陣填充的保結(jié)構(gòu)閾值算法(structure-preserving thresholding algorithm based on F-norm for Hankel matrix completion，簡稱 F-NSPTA)

第一步:給定Ω為下標(biāo)集合,PΩ(M)為樣本矩陣,參數(shù)τ0，0

第二步:計算矩陣YK比閾值τk大的奇異值

[Uk,∑k,Vk]=lansvd(Yk)

令

第三步，計算

第四步: 給定參數(shù),若‖Yk+1-Yk‖F(xiàn)/‖Yk‖F(xiàn)<ε, 就停止 ; 否則改變參數(shù)τ,選擇τk+1使?jié)M足 ‖Yk+1-Xk+1‖F(xiàn)≤‖Yk-Xk‖F(xiàn),轉(zhuǎn)第五步.

第五步: 給定參數(shù)k:=k+1,轉(zhuǎn)第二步.

2 新算法及數(shù)值實驗

2.1 新算法

l步修正的Hankel矩陣填充的保結(jié)構(gòu)閾值算法(structure-preserving threshold algorithm for L-Step modified Hankel matrix completion，簡稱l-NSPTA)

第一步:給定Ω為下標(biāo)集合,PΩ(M)為樣本矩陣，參數(shù)τ0，0

第二步:前l(fā)-1步迭代

1) 計算矩陣YK,q比閾值τk,q大的奇異值

令

2) 計算

Yk+1，q+1=Xk,q+PΩ(M).

3) 給定參數(shù),若‖Yk+1,q+1-Yk,q‖F(xiàn)/‖Yk,q‖F(xiàn)<ε,就停止;否則改變參數(shù)τ，選擇τk+1,q+1使?jié)M足‖Yk+1,q+q-Xk+1,q+1‖F(xiàn)≤‖Yk,q-Xk,q‖F(xiàn).

第三步:第l步迭代

1) 計算矩陣YK，l比閾值τk，l大的奇異值

令

Xk，l=Uk，lDτk，l(∑k，l)Vk，lT.

2) 計算光滑算子

第四步:給定參數(shù)，若‖Yk+l-Yk‖F(xiàn)/‖Yk‖F(xiàn)<ε， 就停止;否則改變參數(shù)τ， 選擇τk+1使?jié)M足‖Yk+l-Xk+l‖F(xiàn)≤‖Yk-Xk‖F(xiàn)，轉(zhuǎn)第五步.

第五步:給定參數(shù)k:=k+1，q:=q+1，轉(zhuǎn)第二步.

很明顯，此算法是F-NSPTA算法的加速算法，如果l=1時就等同于F-NSPTA算法.

2.2 數(shù)值實驗

注：采樣矩陣M產(chǎn)生的Hankel矩陣都是隨機的,未知元素的位置都是不同的,這就會導(dǎo)致個別實驗時間的波動.

表1 p=0.2

表2 p=0.3 續(xù)表2

表3 p=0.4

表4 p=0.6

表5 p=0.6

3 總結(jié)

經(jīng)過全新修正之后的奇異值閾值算法,是更好的填充算法.經(jīng)過兩種算法的比較之后,可以發(fā)現(xiàn)我們提出的步修正算法具有更好的收斂性及更高的精準(zhǔn)度.總的來說,采樣率與計算時間是成反比的.