一種新的相關(guān)系數(shù)矩陣稀疏估計(jì)的對(duì)偶方法

薛文娟, 劉 瀟

(1.上海電力大學(xué) 數(shù)理學(xué)院, 上海 200090; 2.上海理工大學(xué) 理學(xué)院, 上海 200439)

近年來,利用半正定協(xié)方差矩陣的稀疏逼近在聲納數(shù)據(jù)處理[1]、細(xì)胞信號(hào)恢復(fù)[2]、基因聚類[3]、語音信號(hào)分類[4-5]等領(lǐng)域都有著廣泛的應(yīng)用。正如文獻(xiàn)[6]提到的,協(xié)方差矩陣不是尺度不變的。稀疏半正定相關(guān)系數(shù)矩陣在基因聚類的實(shí)際應(yīng)用中也發(fā)揮著重要作用[6-7]。因此,尋求一種新的、高效求解半正定相關(guān)系數(shù)矩陣稀疏逼近問題的算法是非常必要的。

CUI Y等人[6]首次提出采用對(duì)偶方法,通過加速近端梯度來求解稀疏逼近問題的某種對(duì)偶問題。該方法是對(duì)已被廣泛研究的協(xié)方差矩陣稀疏逼近[1-5]的推廣。本文針對(duì)稀疏逼近優(yōu)化問題推導(dǎo)出一個(gè)新的對(duì)偶問題,并利用鄰近梯度法求解該新對(duì)偶問題。

本文首先將引入稀疏逼近優(yōu)化問題的一個(gè)新對(duì)偶,并給出一階最優(yōu)性條件;然后提出一種新的基于近鄰梯度的對(duì)偶方法來求解所得到的新對(duì)偶問題,并證明其具有全局收斂性;最后通過數(shù)值對(duì)比驗(yàn)證了該算法的有效性。

1 對(duì)偶問題及其最優(yōu)性條件

1.1 稀疏逼近優(yōu)化問題

本文主要考慮如下形式的稀疏優(yōu)化問題

s.t.〈Ai,X〉=bi,

〈Bj,X〉≤uj,

X-

(1)

若對(duì)于矩陣P中某些元素,Pij=0,那么對(duì)Xij不要求具有稀疏性。若要求矩陣X的對(duì)角線元素為1,即定義Ai的第s行第t列中的元素為

(2)

其中:s=1,2,3,…,n;t=1,2,3,…,n。b的第i個(gè)元素bi=1,i=1,2,3,…,n。則可以得到式(1)模型的一個(gè)特例,即具有以下形式的相關(guān)系數(shù)矩陣的稀疏逼近:

s.t.Xii=1,i∈[n],

X-

(3)

其中:[n]表示集合{1,2,3,…,n}。

(4)

s.t.A(X)=b,B(X)≤u,

X-

(5)

1.2 一個(gè)新的對(duì)偶問題

經(jīng)過計(jì)算可得問題式(5)的對(duì)偶問題如下

minΨ(Λ,γ,λ)=

s.t.(Λ,λ)∈Ω

(6)

其中:Ω是可行集,即

Λ≤P,λ≥o}

(7)

?ΛΨ(Λ,γ,λ)=(C+A*(γ)-B*(λ)+

Λ-I)++I

?γΨ(Λ,γ,λ)=A(C+A*(γ)-B*(λ)+

Λ-I)++A(I)-b

?λΨ(Λ,γ,λ)=B(C+A*(γ)-B*(λ)+

Λ-I)++B(I)+u

(8)

式(8)中的梯度都是Lipschitz連續(xù)的,且當(dāng)Slater約束規(guī)范條件成立時(shí),問題式(5)和對(duì)偶問題式(6)之間的對(duì)偶間隙為零。從而,只需求解對(duì)偶問題式(6)即可。

1.3 最優(yōu)性條件

由于不等式約束-Pij≤Λij≤Pij,i,j∈[n],λ≥o的特殊結(jié)構(gòu),Mangasarian-Fromovitz約束規(guī)范條件在對(duì)偶問題式(6)的任意可行點(diǎn)都成立[12],于是可以得到對(duì)偶問題式(6)的一階最優(yōu)性條件。

(9)

(10)

(11)

(12)

此外,容易驗(yàn)證一階最優(yōu)性條件式(9)等價(jià)于

r(Λk,γk,λk)=0

(13)

因此,KKT誤差rk=r(Λk,γk,λk)可以用來度量最優(yōu)性。

2 一種新的對(duì)偶方法及其收斂性分析

本節(jié)將應(yīng)用鄰近梯度法求解對(duì)偶問題式(6)。

2.1 求解對(duì)偶問題的鄰近梯度法

鄰近梯度法[13]一般用于求解可分解的、非光滑的無約束優(yōu)化問題,即

minf1(x)+f2(x)

(14)

其中:f1是凸的且可微的;f2是凸的鄰近映射(但不一定是可微的)。在當(dāng)前迭代點(diǎn),鄰近梯度法的迭代公式為

xk+1=proxαkf2(xk-αk?f1(xk))

(15)

其中:αk>0是步長。

min〈?ΛΨ(Λk,γk,λk),Λk+DΛ〉+

(16)

(17)

其中:δΩ(Λ,λ)是Ω的指示函數(shù),即

(18)

注意到式(16)中目標(biāo)函數(shù)的前4項(xiàng)是連續(xù)可微的,最后1項(xiàng)是凸函數(shù),其鄰近映射可由式(19)計(jì)算得到

proxδΩ(Λk+DΛ,λk+dλ)=

ΠΩ(Λk+DΛ,λk+dλ)

(19)

其中:ΠΩ(·)由式(11)計(jì)算可得。容易推出式(16)和式(17)閉合形式的解分別如下

λk-?λΨ(Λk,γk,λk))-

(Λk,λk)

(20)

(21)

2.2 主要算法

下面將給出完整的對(duì)偶鄰近梯度算法(Dual Proximal Gradient Algorithm,DPGA)。令(Λk,γk,λk)為當(dāng)前迭代點(diǎn),通過式(16)和式(17)(或式(20)和式(21))可得鄰近梯度步。定義Ψ(Λ,γ,λ)的預(yù)計(jì)線性下降量Δk為

(22)

利用線搜索回溯技術(shù),在{1,τ,τ2,…}中尋找最大的αk,使得Ψ(Λ,γ,λ)滿足充分下降條件

Ψ(Λk,γk,λk)+σαkΔk

(23)

其中:τ∈(0,1);σ∈(0,1)是給定的常數(shù)。從而,可得DPGA算法的迭代公式為

(24)

對(duì)于算法的終止條件,將類似于文獻(xiàn)[6],用ηgap來表示原始-對(duì)偶的間隙,可以衡量原始目標(biāo)與對(duì)偶目標(biāo)函數(shù)值之間的差異,即

(25)

其中:jprob和jdob分別表示原始目標(biāo)函數(shù)值和對(duì)偶目標(biāo)函數(shù)值。

當(dāng)式(26)成立時(shí),DPGA算法停止。

(26)

定理1由DPGA算法生成的序列{(Λk,γk,λk)}的任何聚點(diǎn)都是對(duì)偶問題式(6)的最優(yōu)解。

證明:應(yīng)用鄰近梯度方法收斂性結(jié)論立即可得到所提算法的全局收斂性。

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

為了驗(yàn)證DPGA算法的有效性,本文將其與文獻(xiàn)[6]中的相關(guān)系數(shù)矩陣稀疏估計(jì)(Sparse Estimation of the Correlation matrix,SEC)算法進(jìn)行了比較。為了公平比較,本文還將SEC算法代碼中的EIG(·)替換為Mexeig(·)。

3.1 設(shè)置問題

本文實(shí)驗(yàn)中的測(cè)試問題形式如下

s.t.Xii=1,i∈[n],

X-

(27)

取ε=10-6。當(dāng)DPGA算法滿足終止條件式(26),或達(dá)到最大迭代次數(shù)(5 000次)時(shí)停止。DPGA算法中的其他參數(shù)設(shè)置如下:Λ0=O,γ0=o,σ=10-4,τ=0.5。此外,SEC算法采用默認(rèn)參數(shù)。

3.2 合成數(shù)據(jù)集

本文生成矩陣C有以下兩種方法。

(1) 測(cè)試問題E1式(27)中的矩陣C由MATLAB命令實(shí)現(xiàn)生成。具體為:C=rand(n,n);C=(C+C′)/2;C(1:n+1:end)=1。

(2) 測(cè)試問題E2式(27)中的矩陣C是一個(gè)樣本相關(guān)系數(shù)矩陣,是由一組N(0,Σ)樣本計(jì)算生成的,均值為零。其中:Σ是一個(gè)的n×n的協(xié)方差矩陣。特別地,定義

i,j∈[n]

(28)

接下來,將從N(0,Σ)中抽取p個(gè)獨(dú)立同分布的樣本,用D∈Rp×n來表示。最后,應(yīng)用zscore方法對(duì)D進(jìn)行歸一化,使得每個(gè)特征的均值為零,單位方差為1。在歸一化數(shù)據(jù)上,樣本協(xié)方差矩陣與樣本相關(guān)系數(shù)矩陣相同。在本文的數(shù)值實(shí)驗(yàn)中,對(duì)所有給定的矩陣C,ρ=0.01,將通過求解式(27)來估計(jì)相關(guān)系數(shù)矩陣。

本文從式(27)對(duì)偶問題的殘差、CPU時(shí)間、總迭代次數(shù)、原始目標(biāo)值,以及估計(jì)的相關(guān)系數(shù)矩陣的稀疏性等5個(gè)方面與算法SEC進(jìn)行比較。實(shí)驗(yàn)結(jié)果表明,SEC算法和DPGA算法都能在預(yù)設(shè)誤差小于10-6和最大迭代次數(shù)內(nèi)成功地解決所有測(cè)試問題。

圖1顯示了SEC和DPGA算法在E1上的數(shù)值結(jié)果。其中,φ和φmin分別為目標(biāo)函數(shù)值和其最小值。由圖1(a)可以看出,在n∈[1 000,2 800]之間所有10個(gè)測(cè)試問題上,盡管差距很小,但DPGA算法得到的原始目標(biāo)值都優(yōu)于SEC算法。對(duì)于殘差,圖1(b)顯示DPGA算法在所有情況下都比SEC算法更精確。此外,由圖1(c)和圖1(d)可以看出,DPGA算法在E1上的迭代次數(shù)不僅小于SEC算法的1/4,而且CPU時(shí)間也少了一半,這說明本文所提出的DPGA算法更具優(yōu)越性。

圖1 SEC和DPGA算法在E1上的數(shù)值結(jié)果

此外,通過設(shè)置n∈[1 000,2 000]和樣本數(shù)p∈[1 000,2 000],對(duì)E2進(jìn)行了實(shí)驗(yàn)研究。圖2和表1給出了E2的相關(guān)數(shù)值結(jié)果。由圖2可以看出,與SEC算法相比,DPGA算法在E2上幾乎都能以更少的迭代次數(shù)獲得更高精度的解。 從表1可以觀察到,SEC和DPGA算法在E2上都可以達(dá)到幾乎相同的原始目標(biāo)值及相同的稀疏程度。

圖2 SEC和DPGA算法在E2上的數(shù)值結(jié)果

表1 E2上的數(shù)值結(jié)果

綜上所述,針對(duì)本文所提出的新對(duì)偶問題式(6),DPGA算法可以較低的計(jì)算量獲得相關(guān)系數(shù)矩陣稀疏逼近的最優(yōu)解。

4 結(jié) 語

本文提出了一種Frobenius范數(shù)下稀疏矩陣逼近問題的對(duì)偶鄰近梯度算法,用于求解由原始問題導(dǎo)出的一種新的對(duì)偶問題。該算法利用鄰近梯度產(chǎn)生搜索方向,并使用線搜索技術(shù)來尋找合適的步長以確保算法具有全局收斂性。此外,本文還對(duì)相關(guān)模型隨機(jī)生成數(shù)據(jù)進(jìn)行了初步的數(shù)值實(shí)驗(yàn),并與文獻(xiàn)[6]中的SEC算法進(jìn)行了比較。通過對(duì)數(shù)值結(jié)果比較,結(jié)果顯示本文提出的求解對(duì)偶問題的DPGA算法優(yōu)于SEC算法。