薛文娟, 劉 瀟
(1.上海電力大學 數(shù)理學院, 上海 200090; 2.上海理工大學 理學院, 上海 200439)
近年來,利用半正定協(xié)方差矩陣的稀疏逼近在聲納數(shù)據處理[1]、細胞信號恢復[2]、基因聚類[3]、語音信號分類[4-5]等領域都有著廣泛的應用。正如文獻[6]提到的,協(xié)方差矩陣不是尺度不變的。稀疏半正定相關系數(shù)矩陣在基因聚類的實際應用中也發(fā)揮著重要作用[6-7]。因此,尋求一種新的、高效求解半正定相關系數(shù)矩陣稀疏逼近問題的算法是非常必要的。
CUI Y等人[6]首次提出采用對偶方法,通過加速近端梯度來求解稀疏逼近問題的某種對偶問題。該方法是對已被廣泛研究的協(xié)方差矩陣稀疏逼近[1-5]的推廣。本文針對稀疏逼近優(yōu)化問題推導出一個新的對偶問題,并利用鄰近梯度法求解該新對偶問題。
本文首先將引入稀疏逼近優(yōu)化問題的一個新對偶,并給出一階最優(yōu)性條件;然后提出一種新的基于近鄰梯度的對偶方法來求解所得到的新對偶問題,并證明其具有全局收斂性;最后通過數(shù)值對比驗證了該算法的有效性。
本文主要考慮如下形式的稀疏優(yōu)化問題
s.t.〈Ai,X〉=bi,
〈Bj,X〉≤uj,
X-
(1)
若對于矩陣P中某些元素,Pij=0,那么對Xij不要求具有稀疏性。若要求矩陣X的對角線元素為1,即定義Ai的第s行第t列中的元素為
(2)
其中:s=1,2,3,…,n;t=1,2,3,…,n。b的第i個元素bi=1,i=1,2,3,…,n。則可以得到式(1)模型的一個特例,即具有以下形式的相關系數(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)
經過計算可得問題式(5)的對偶問題如下
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ù)的,且當Slater約束規(guī)范條件成立時,問題式(5)和對偶問題式(6)之間的對偶間隙為零。從而,只需求解對偶問題式(6)即可。
由于不等式約束-Pij≤Λij≤Pij,i,j∈[n],λ≥o的特殊結構,Mangasarian-Fromovitz約束規(guī)范條件在對偶問題式(6)的任意可行點都成立[12],于是可以得到對偶問題式(6)的一階最優(yōu)性條件。
(9)
(10)
(11)
(12)
此外,容易驗證一階最優(yōu)性條件式(9)等價于
r(Λk,γk,λk)=0
(13)
因此,KKT誤差rk=r(Λk,γk,λk)可以用來度量最優(yōu)性。
本節(jié)將應用鄰近梯度法求解對偶問題式(6)。
鄰近梯度法[13]一般用于求解可分解的、非光滑的無約束優(yōu)化問題,即
minf1(x)+f2(x)
(14)
其中:f1是凸的且可微的;f2是凸的鄰近映射(但不一定是可微的)。在當前迭代點,鄰近梯度法的迭代公式為
xk+1=proxαkf2(xk-αk?f1(xk))
(15)
其中:αk>0是步長。
min〈?ΛΨ(Λk,γk,λk),Λk+DΛ〉+
(16)
(17)
其中:δΩ(Λ,λ)是Ω的指示函數(shù),即
(18)
注意到式(16)中目標函數(shù)的前4項是連續(xù)可微的,最后1項是凸函數(shù),其鄰近映射可由式(19)計算得到
proxδΩ(Λk+DΛ,λk+dλ)=
ΠΩ(Λk+DΛ,λk+dλ)
(19)
其中:ΠΩ(·)由式(11)計算可得。容易推出式(16)和式(17)閉合形式的解分別如下
λk-?λΨ(Λk,γk,λk))-
(Λk,λk)
(20)
(21)
下面將給出完整的對偶鄰近梯度算法(Dual Proximal Gradient Algorithm,DPGA)。令(Λk,γk,λk)為當前迭代點,通過式(16)和式(17)(或式(20)和式(21))可得鄰近梯度步。定義Ψ(Λ,γ,λ)的預計線性下降量Δk為
(22)
利用線搜索回溯技術,在{1,τ,τ2,…}中尋找最大的αk,使得Ψ(Λ,γ,λ)滿足充分下降條件
Ψ(Λk,γk,λk)+σαkΔk
(23)
其中:τ∈(0,1);σ∈(0,1)是給定的常數(shù)。從而,可得DPGA算法的迭代公式為
(24)
對于算法的終止條件,將類似于文獻[6],用ηgap來表示原始-對偶的間隙,可以衡量原始目標與對偶目標函數(shù)值之間的差異,即
(25)
其中:jprob和jdob分別表示原始目標函數(shù)值和對偶目標函數(shù)值。
當式(26)成立時,DPGA算法停止。
(26)
定理1由DPGA算法生成的序列{(Λk,γk,λk)}的任何聚點都是對偶問題式(6)的最優(yōu)解。
證明:應用鄰近梯度方法收斂性結論立即可得到所提算法的全局收斂性。
為了驗證DPGA算法的有效性,本文將其與文獻[6]中的相關系數(shù)矩陣稀疏估計(Sparse Estimation of the Correlation matrix,SEC)算法進行了比較。為了公平比較,本文還將SEC算法代碼中的EIG(·)替換為Mexeig(·)。
本文實驗中的測試問題形式如下
s.t.Xii=1,i∈[n],
X-
(27)
取ε=10-6。當DPGA算法滿足終止條件式(26),或達到最大迭代次數(shù)(5 000次)時停止。DPGA算法中的其他參數(shù)設置如下:Λ0=O,γ0=o,σ=10-4,τ=0.5。此外,SEC算法采用默認參數(shù)。
本文生成矩陣C有以下兩種方法。
(1) 測試問題E1式(27)中的矩陣C由MATLAB命令實現(xiàn)生成。具體為:C=rand(n,n);C=(C+C′)/2;C(1:n+1:end)=1。
(2) 測試問題E2式(27)中的矩陣C是一個樣本相關系數(shù)矩陣,是由一組N(0,Σ)樣本計算生成的,均值為零。其中:Σ是一個的n×n的協(xié)方差矩陣。特別地,定義
i,j∈[n]
(28)
接下來,將從N(0,Σ)中抽取p個獨立同分布的樣本,用D∈Rp×n來表示。最后,應用zscore方法對D進行歸一化,使得每個特征的均值為零,單位方差為1。在歸一化數(shù)據上,樣本協(xié)方差矩陣與樣本相關系數(shù)矩陣相同。在本文的數(shù)值實驗中,對所有給定的矩陣C,ρ=0.01,將通過求解式(27)來估計相關系數(shù)矩陣。
本文從式(27)對偶問題的殘差、CPU時間、總迭代次數(shù)、原始目標值,以及估計的相關系數(shù)矩陣的稀疏性等5個方面與算法SEC進行比較。實驗結果表明,SEC算法和DPGA算法都能在預設誤差小于10-6和最大迭代次數(shù)內成功地解決所有測試問題。
圖1顯示了SEC和DPGA算法在E1上的數(shù)值結果。其中,φ和φmin分別為目標函數(shù)值和其最小值。由圖1(a)可以看出,在n∈[1 000,2 800]之間所有10個測試問題上,盡管差距很小,但DPGA算法得到的原始目標值都優(yōu)于SEC算法。對于殘差,圖1(b)顯示DPGA算法在所有情況下都比SEC算法更精確。此外,由圖1(c)和圖1(d)可以看出,DPGA算法在E1上的迭代次數(shù)不僅小于SEC算法的1/4,而且CPU時間也少了一半,這說明本文所提出的DPGA算法更具優(yōu)越性。
圖1 SEC和DPGA算法在E1上的數(shù)值結果
此外,通過設置n∈[1 000,2 000]和樣本數(shù)p∈[1 000,2 000],對E2進行了實驗研究。圖2和表1給出了E2的相關數(shù)值結果。由圖2可以看出,與SEC算法相比,DPGA算法在E2上幾乎都能以更少的迭代次數(shù)獲得更高精度的解。 從表1可以觀察到,SEC和DPGA算法在E2上都可以達到幾乎相同的原始目標值及相同的稀疏程度。
圖2 SEC和DPGA算法在E2上的數(shù)值結果
表1 E2上的數(shù)值結果
綜上所述,針對本文所提出的新對偶問題式(6),DPGA算法可以較低的計算量獲得相關系數(shù)矩陣稀疏逼近的最優(yōu)解。
本文提出了一種Frobenius范數(shù)下稀疏矩陣逼近問題的對偶鄰近梯度算法,用于求解由原始問題導出的一種新的對偶問題。該算法利用鄰近梯度產生搜索方向,并使用線搜索技術來尋找合適的步長以確保算法具有全局收斂性。此外,本文還對相關模型隨機生成數(shù)據進行了初步的數(shù)值實驗,并與文獻[6]中的SEC算法進行了比較。通過對數(shù)值結果比較,結果顯示本文提出的求解對偶問題的DPGA算法優(yōu)于SEC算法。