張春曉，宋儒瑛

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

0 引言

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

(1)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1 算法回顧

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

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

(2)

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

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

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

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

(3)

其算法如下:

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

第二步: 計(jì)算矩陣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.

第二步: 計(jì)算矩陣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，簡(jiǎn)稱(chēng) F-NSPTA)

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

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

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

令

第三步，計(jì)算

第四步: 給定參數(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ù)值實(shí)驗(yàn)

2.1 新算法

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

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

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

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

令

2) 計(jì)算

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) 計(jì)算矩陣YK，l比閾值τk，l大的奇異值

令

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

2) 計(jì)算光滑算子

第四步:給定參數(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時(shí)就等同于F-NSPTA算法.

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

注：采樣矩陣M產(chǎn)生的Hankel矩陣都是隨機(jī)的,未知元素的位置都是不同的,這就會(huì)導(dǎo)致個(gè)別實(shí)驗(yàn)時(shí)間的波動(dòng).

表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)過(guò)全新修正之后的奇異值閾值算法,是更好的填充算法.經(jīng)過(guò)兩種算法的比較之后,可以發(fā)現(xiàn)我們提出的步修正算法具有更好的收斂性及更高的精準(zhǔn)度.總的來(lái)說(shuō),采樣率與計(jì)算時(shí)間是成反比的.